R version 4.4.1 (2024-06-14 ucrt) -- "Race for Your Life"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ############################################################################
> #     MLwiN User Manual
> #
> # 13  Fitting Models to Repeated Measures Data . . . . . . . . . . . . . 191
> #
> #     Rasbash, J., Steele, F., Browne, W. J. and Goldstein, H. (2012).
> #     A User's Guide to MLwiN, v2.26. Centre for Multilevel Modelling,
> #     University of Bristol.
> ############################################################################
> #     R script to replicate all analyses using R2MLwiN
> #
> #     Zhang, Z., Charlton, C., Parker, R, Leckie, G., and Browne, W.J.
> #     Centre for Multilevel Modelling, 2012
> #     http://www.bristol.ac.uk/cmm/software/R2MLwiN/
> ############################################################################
> 
> library(R2MLwiN)

R2MLwiN: A package to run models implemented in MLwiN from R
Copyright 2013-2024 Zhengzheng Zhang, Christopher M. J. Charlton, Richard M. A. Parker, William J. Browne and George Leckie

Support provided by the Economic and Social Research Council (ESRC)
(Grants RES-149-25-1084, RES-576-25-0032 and ES/K007246/1)

To cite R2MLwiN in publications use:

  Zhengzheng Zhang, Richard M. A. Parker, Christopher M. J. Charlton,
  George Leckie, William J. Browne (2016). R2MLwiN: A Package to Run
  MLwiN from within R. Journal of Statistical Software, 72(10), 1-43.
  doi:10.18637/jss.v072.i10

A BibTeX entry for LaTeX users is

  @Article{,
    title = {{R2MLwiN}: A Package to Run {MLwiN} from within {R}},
    author = {Zhengzheng Zhang and Richard M. A. Parker and Christopher M. J. Charlton and George Leckie and William J. Browne},
    journal = {Journal of Statistical Software},
    year = {2016},
    volume = {72},
    number = {10},
    pages = {1--43},
    doi = {10.18637/jss.v072.i10},
  }

The MLwiN_path option is currently set to C:/Program Files/MLwiN v3.11/
To change this use: options(MLwiN_path="<path to MLwiN>")
> # MLwiN folder
> mlwin <- getOption("MLwiN_path")
> while (!file.access(mlwin, mode = 1) == 0) {
+   cat("Please specify the root MLwiN folder or the full path to the MLwiN executable:\n")
+   mlwin <- scan(what = character(0), sep = "\n")
+   mlwin <- gsub("\\", "/", mlwin, fixed = TRUE)
+ }
> options(MLwiN_path = mlwin)
> 
> 
> # 13.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .191
> 
> # 13.2 A basic model . . . . . . . . . . . . . . . . . . . . . . . . . . 194
> 
> data(reading1, package = "R2MLwiN")
> summary(reading1)
       id             age1             read1                age2         
 Min.   :  1.0   Min.   :-10.000   Min.   :-10.00000   Min.   :-10.0000  
 1st Qu.:149.0   1st Qu.:-10.000   1st Qu.:-10.00000   1st Qu.: -1.5804  
 Median :363.0   Median : -2.500   Median :  4.26330   Median : -1.4704  
 Mean   :342.6   Mean   : -4.834   Mean   :  0.01828   Mean   : -2.1946  
 3rd Qu.:520.5   3rd Qu.: -2.380   3rd Qu.:  4.72640   3rd Qu.: -1.3604  
 Max.   :751.0   Max.   : -1.700   Max.   :  7.50500   Max.   : -0.6804  
     read2              age3              read3              age4         
 Min.   :-10.000   Min.   :-10.0000   Min.   :-10.000   Min.   :-10.0000  
 1st Qu.:  4.894   1st Qu.: -0.6104   1st Qu.:  5.743   1st Qu.:-10.0000  
 Median :  5.490   Median : -0.5004   Median :  6.531   Median :  0.4196  
 Mean   :  4.302   Mean   : -1.6561   Mean   :  4.581   Mean   : -2.4247  
 3rd Qu.:  6.214   3rd Qu.: -0.4004   3rd Qu.:  6.981   3rd Qu.:  0.5396  
 Max.   :  7.150   Max.   :  0.2896   Max.   :  9.907   Max.   :  1.2496  
     read4              age5             read5              age6        
 Min.   :-10.000   Min.   :-10.000   Min.   :-10.000   Min.   :-10.000  
 1st Qu.:-10.000   1st Qu.:-10.000   1st Qu.:-10.000   1st Qu.:-10.000  
 Median :  7.031   Median :  1.390   Median :  7.603   Median :-10.000  
 Mean   :  2.728   Mean   : -2.607   Mean   :  1.980   Mean   : -3.089  
 3rd Qu.:  7.853   3rd Qu.:  1.520   3rd Qu.:  9.068   3rd Qu.:  4.210  
 Max.   :  9.908   Max.   :  2.250   Max.   : 10.420   Max.   :  4.510  
     read6         
 Min.   :-10.0000  
 1st Qu.:-10.0000  
 Median :-10.0000  
 Mean   :  0.3622  
 3rd Qu.: 11.4330  
 Max.   : 13.8690  
> 
> reading1[reading1 == -10] <- NA
> 
> summary(reading1)
       id             age1           read1            age2        
 Min.   :  1.0   Min.   :-2.71   Min.   :3.893   Min.   :-2.0104  
 1st Qu.:149.0   1st Qu.:-2.51   1st Qu.:4.263   1st Qu.:-1.5404  
 Median :363.0   Median :-2.43   Median :4.541   Median :-1.4504  
 Mean   :342.6   Mean   :-2.41   Mean   :4.720   Mean   :-1.4372  
 3rd Qu.:520.5   3rd Qu.:-2.34   3rd Qu.:5.004   3rd Qu.:-1.3504  
 Max.   :751.0   Max.   :-1.70   Max.   :7.505   Max.   :-0.6804  
                 NA's   :130     NA's   :130     NA's   :36       
     read2            age3             read3            age4        
 Min.   :4.341   Min.   :-1.0604   Min.   :4.955   Min.   :-0.0404  
 1st Qu.:5.107   1st Qu.:-0.5704   1st Qu.:6.080   1st Qu.: 0.3996  
 Median :5.575   Median :-0.4754   Median :6.643   Median : 0.4896  
 Mean   :5.690   Mean   :-0.4608   Mean   :6.670   Mean   : 0.4869  
 3rd Qu.:6.299   3rd Qu.:-0.3804   3rd Qu.:7.093   3rd Qu.: 0.5696  
 Max.   :7.150   Max.   : 0.2896   Max.   :9.907   Max.   : 1.2496  
 NA's   :36      NA's   :51        NA's   :51      NA's   :113      
     read4            age5            read5             age6      
 Min.   :5.490   Min.   :0.9596   Min.   : 5.462   Min.   :3.660  
 1st Qu.:6.928   1st Qu.:1.3996   1st Qu.: 7.744   1st Qu.:4.130  
 Median :7.545   Median :1.4896   Median : 8.617   Median :4.210  
 Mean   :7.620   Mean   :1.4838   Mean   : 8.610   Mean   :4.207  
 3rd Qu.:8.367   3rd Qu.:1.5696   3rd Qu.: 9.519   3rd Qu.:4.290  
 Max.   :9.908   Max.   :2.2496   Max.   :10.420   Max.   :4.510  
 NA's   :113     NA's   :145      NA's   :145      NA's   :209    
     read6       
 Min.   : 6.451  
 1st Qu.:10.575  
 Median :11.544  
 Mean   :11.300  
 3rd Qu.:12.319  
 Max.   :13.869  
 NA's   :209     
> 
> reading <- reshape(reading1, idvar = "student", timevar = "id", varying = c("read1", "age1", "read2", "age2", "read3", 
+   "age3", "read4", "age4", "read5", "age5", "read6", "age6"), sep = "", direction = "long")
> 
> reading <- reading[c("student", "id", "age", "read")]
> reading <- reading[order(reading$student, reading$id), ]
> 
> colnames(reading) <- c("student", "occasion", "age", "reading")
> rownames(reading) <- NULL
> 
> 
> summary(reading)
    student       occasion        age             reading      
 Min.   :  1   Min.   :1.0   Min.   :-2.7104   Min.   : 3.893  
 1st Qu.:102   1st Qu.:2.0   1st Qu.:-1.4704   1st Qu.: 5.497  
 Median :204   Median :3.5   Median :-0.4304   Median : 6.756  
 Mean   :204   Mean   :3.5   Mean   : 0.0000   Mean   : 7.125  
 3rd Qu.:306   3rd Qu.:5.0   3rd Qu.: 1.3096   3rd Qu.: 8.279  
 Max.   :407   Max.   :6.0   Max.   : 4.5096   Max.   :13.869  
                             NA's   :684       NA's   :684     
> 
> head(reading, 5)
  student occasion     age reading
1       1        1 -2.6504  6.3009
2       1        2 -1.6104  6.8948
3       1        3 -0.6204  7.7685
4       1        4  0.3396  8.8804
5       1        5  1.3396 10.0820
> 
> tab <- aggregate(reading ~ occasion, reading, function(x) c(N = length(x), mean = mean(x), sd = sd(x)))
> tab <- rbind(tab, c(NA, NA))
> tab$reading[7, ] <- c(length(na.omit(reading$reading)), mean(na.omit(reading$reading)), sd(na.omit(reading$reading)))
> rownames(tab)[7] <- "Total"
> tab
      occasion    reading.N reading.mean   reading.sd
1            1  277.0000000    4.7200072    0.6099982
2            2  371.0000000    5.6899981    0.7399990
3            3  356.0000000    6.6699997    0.8699961
4            4  294.0000000    7.6200003    0.9900053
5            5  262.0000000    8.6099938    1.1199879
6            6  198.0000000   11.3000247    1.4700124
Total       NA 1758.0000000    7.1254008    2.1544282
> 
> tab <- aggregate(age ~ occasion, reading, function(x) c(N = length(x), mean = mean(x), sd = sd(x)))
> tab <- rbind(tab, c(NA, NA))
> tab$age[7, ] <- c(length(na.omit(reading$age)), mean(na.omit(reading$age)), sd(na.omit(reading$age)))
> rownames(tab)[7] <- "Total"
> tab
      occasion         age.N      age.mean        age.sd
1            1  2.770000e+02 -2.410003e+00  1.449145e-01
2            2  3.710000e+02 -1.437192e+00  1.539733e-01
3            3  3.560000e+02 -4.607933e-01  1.571502e-01
4            4  2.940000e+02  4.868789e-01  1.269735e-01
5            5  2.620000e+02  1.483837e+00  1.307948e-01
6            6  1.980000e+02  4.206671e+00  1.150215e-01
Total       NA  1.758000e+03  9.554091e-06  1.943460e+00
> 
> (mymodel1 <- runMLwiN(reading ~ 1 + (1 | student) + (1 | occasion), data = reading))
/nogui option ignored
ECHO    0


Echoing is ON
BATC 1

Batch mode is ON
MAXI 2
STAR
iteration 0
iteration 1

Convergence not achieved
TOLE 2
MAXI 20
NEXT
iteration 2

Convergence achieved
ECHO 0
Execution completed


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.11)  multilevel model (Normal) 
          N min mean max N_complete min_complete mean_complete max_complete
student 407   6    6   6        407            1       4.31941            6
Estimation algorithm:  IGLS        Elapsed time : 0.11s 
Number of obs:  1758 (from total 2442)        The model converged after 3 iterations.
Log likelihood:      -3842.9 
Deviance statistic:  7685.7 
--------------------------------------------------------------------------------------------------- 
The model formula:
reading ~ 1 + (1 | student) + (1 | occasion)
Level 2: student     Level 1: occasion      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
              Coef.   Std. Err.        z   Pr(>|z|)         [95% Conf.   Interval] 
Intercept   7.11491     0.05303   134.17          0   ***      7.01097     7.21885 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                  Coef.   Std. Err. 
var_Intercept   0.07761     0.08309 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the occasion level: 
                  Coef.   Std. Err. 
var_Intercept   4.56149     0.17245 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> 
> # 13.3 A linear growth curve model . . . . . . . . . . . . . . . . . . . 201
> 
> (mymodel2 <- runMLwiN(reading ~ 1 + age + (1 | student) + (1 | occasion), estoptions = list(startval = list(FP.b = mymodel1@FP, 
+   FP.v = mymodel1@FP.cov, RP.b = mymodel1@RP, RP.v = mymodel1@RP.cov)), data = reading))
/nogui option ignored
ECHO    0


Echoing is ON
BATC 1

Batch mode is ON
MAXI 2
STAR
iteration 0
iteration 1

Convergence not achieved
JOIN 7.11491148444759 0 '_FP_b'
JOIN 0.0028122591557062 0 0 '_FP_v'
JOIN 0.0776119399796296 4.56148954336666 '_RP_b'
JOIN 0.00690456124179504 -0.0060668859565199 0.0297395862616497 '_RP_v'
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
iteration 5

Convergence achieved
ECHO 0
Execution completed


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.11)  multilevel model (Normal) 
          N min mean max N_complete min_complete mean_complete max_complete
student 407   6    6   6        407            1       4.31941            6
Estimation algorithm:  IGLS        Elapsed time : 0.16s 
Number of obs:  1758 (from total 2442)        The model converged after 6 iterations.
Log likelihood:      -1897.8 
Deviance statistic:  3795.6 
--------------------------------------------------------------------------------------------------- 
The model formula:
reading ~ 1 + age + (1 | student) + (1 | occasion)
Level 2: student     Level 1: occasion      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
              Coef.   Std. Err.        z   Pr(>|z|)         [95% Conf.   Interval] 
Intercept   7.11710     0.04128   172.41          0   ***      7.03619     7.19801 
age         0.99676     0.00722   138.07          0   ***      0.98261     1.01091 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                  Coef.   Std. Err. 
var_Intercept   0.60298     0.04844 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the occasion level: 
                  Coef.   Std. Err. 
var_Intercept   0.30724     0.01181 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> (mymodel3 <- runMLwiN(reading ~ 1 + age + (1 + age | student) + (1 | occasion), estoptions = list(resi.store = TRUE, 
+   startval = list(FP.b = mymodel2@FP, FP.v = mymodel2@FP.cov, RP.b = mymodel2@RP, RP.v = mymodel2@RP.cov)), data = reading))
/nogui option ignored
ECHO    0


Echoing is ON
BATC 1

Batch mode is ON
MAXI 2
STAR
iteration 0
iteration 1

Convergence not achieved
JOIN 7.11710098696708 0.996758938786093 '_FP_b'
JOIN 0.00170404929632015 1.11939939725668e-05 5.21184505439158e-05 '_FP_v'
JOIN 0.602981523429859 0 0 0.307237233390386 '_RP_b'
JOIN 0.00234657001927917 0 0 0 0 0 -3.83095220098178e-05 0 0 0.000139481811276294 '_RP_v'
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
iteration 5

Convergence achieved
ECHO 0
Execution completed


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.11)  multilevel model (Normal) 
          N min mean max N_complete min_complete mean_complete max_complete
student 407   6    6   6        407            1       4.31941            6
Estimation algorithm:  IGLS        Elapsed time : 0.23s 
Number of obs:  1758 (from total 2442)        The model converged after 6 iterations.
Log likelihood:      -1604.7 
Deviance statistic:  3209.4 
--------------------------------------------------------------------------------------------------- 
The model formula:
reading ~ 1 + age + (1 + age | student) + (1 | occasion)
Level 2: student     Level 1: occasion      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
              Coef.   Std. Err.        z   Pr(>|z|)         [95% Conf.   Interval] 
Intercept   7.11726     0.04335   164.19          0   ***      7.03230     7.20222 
age         0.99506     0.01219    81.63          0   ***      0.97117     1.01896 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                      Coef.   Std. Err. 
var_Intercept       0.68263     0.05302 
cov_Intercept_age   0.12324     0.01237 
var_age             0.03735     0.00386 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the occasion level: 
                  Coef.   Std. Err. 
var_Intercept   0.16101     0.00705 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> u0 <- mymodel3@residual$lev_2_resi_est_Intercept
> u0std <- (u0 - mean(u0))/sd(u0)
> 
> u1 <- mymodel3@residual$lev_2_resi_est_age
> u1std <- (u1 - mean(u1))/sd(u1)
> 
> plot(u0std, u1std, asp = 1)
> 
> e0 <- na.omit(mymodel3@residual$lev_1_resi_est_Intercept)
> e0std <- (e0 - mean(e0))/sd(e0)
> e0rank <- rank(e0)
> e0uniform <- (e0rank - 0.5)/length(e0rank)
> e0nscore <- qnorm(e0uniform)
> 
> plot(e0std, e0nscore, asp = 1)
> 
> # 13.4 Complex level 1 variation . . . . . . . . . . . . . . . . . . . . 204
> 
> (mymodel4 <- runMLwiN(reading ~ 1 + age + (1 + age | student) + (1 + age | occasion), data = reading))
/nogui option ignored
ECHO    0


Echoing is ON
BATC 1

Batch mode is ON
MAXI 2
STAR
iteration 0
iteration 1

Convergence not achieved
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 7
iteration 8

Convergence achieved
ECHO 0
Execution completed


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.11)  multilevel model (Normal) 
          N min mean max N_complete min_complete mean_complete max_complete
student 407   6    6   6        407            1       4.31941            6
Estimation algorithm:  IGLS        Elapsed time : 0.36s 
Number of obs:  1758 (from total 2442)        The model converged after 9 iterations.
Log likelihood:      -1588.5 
Deviance statistic:  3177 
--------------------------------------------------------------------------------------------------- 
The model formula:
reading ~ 1 + age + (1 + age | student) + (1 + age | occasion)
Level 2: student     Level 1: occasion      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
              Coef.   Std. Err.        z   Pr(>|z|)         [95% Conf.   Interval] 
Intercept   7.11696     0.04390   162.12          0   ***      7.03092     7.20300 
age         0.99505     0.01232    80.75          0   ***      0.97089     1.01920 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                      Coef.   Std. Err. 
var_Intercept       0.70489     0.05437 
cov_Intercept_age   0.12797     0.01278 
var_age             0.03546     0.00404 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the occasion level: 
                      Coef.   Std. Err. 
var_Intercept       0.12268     0.00825 
cov_Intercept_age   0.00067     0.00273 
var_age             0.01451     0.00305 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> # 13.5 Repeated measures modelling of non-linear polynomial growth . . . 205
> 
> (mymodel5 <- runMLwiN(reading ~ 1 + age + I(age^2) + (1 + age + I(age^2) | student) + (1 + age | occasion), estoptions = list(resi.store = TRUE), 
+   data = reading))
/nogui option ignored
ECHO    0


Echoing is ON
BATC 1

Batch mode is ON
MAXI 2
STAR
iteration 0
iteration 1

Convergence not achieved
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 7

Convergence achieved
ECHO 0
Execution completed


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.11)  multilevel model (Normal) 
          N min mean max N_complete min_complete mean_complete max_complete
student 407   6    6   6        407            1       4.31941            6
Estimation algorithm:  IGLS        Elapsed time : 0.49s 
Number of obs:  1758 (from total 2442)        The model converged after 8 iterations.
Log likelihood:      -1566 
Deviance statistic:  3132 
--------------------------------------------------------------------------------------------------- 
The model formula:
reading ~ 1 + age + I(age^2) + (1 + age + I(age^2) | student) + 
    (1 + age | occasion)
Level 2: student     Level 1: occasion      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
              Coef.   Std. Err.        z   Pr(>|z|)         [95% Conf.   Interval] 
Intercept   7.11503     0.04629   153.72          0   ***      7.02431     7.20575 
age         0.99454     0.01259    79.00          0   ***      0.96987     1.01921 
I(age^2)    0.00084     0.00316     0.26     0.7916           -0.00537     0.00704 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                            Coef.   Std. Err. 
var_Intercept             0.76685     0.06009 
cov_Intercept_age         0.14083     0.01390 
var_age                   0.03939     0.00427 
cov_Intercept_I(age^2)   -0.01435     0.00306 
cov_age_I(age^2)         -0.00178     0.00077 
var_I(age^2)              0.00133     0.00031 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the occasion level: 
                       Coef.   Std. Err. 
var_Intercept        0.12908     0.00857 
cov_Intercept_age   -0.00608     0.00267 
var_age              0.00039     0.00392 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> l2varfn <- mymodel5@RP["RP2_var_Intercept"] + 2 * mymodel5@RP["RP2_cov_Intercept_age"] * mymodel5@data$age + mymodel5@RP["RP2_var_age"] * 
+   mymodel5@data$age^2 + 2 * mymodel5@RP["RP2_cov_Intercept_I(age^2)"] * mymodel5@data$age^2 + 2 * mymodel5@RP["RP2_cov_age_I(age^2)"] * 
+   mymodel5@data$age * mymodel5@data[["I(age^2)"]] + mymodel5@RP["RP2_var_I(age^2)"] * mymodel5@data[["I(age^2)"]]^2
> 
> l1varfn <- mymodel5@RP["RP1_var_Intercept"] + 2 * mymodel5@RP["RP1_cov_Intercept_age"] * mymodel5@data$age + mymodel5@RP["RP1_var_age"] * 
+   mymodel5@data$age^2
> 
> totvarfn <- l2varfn + l1varfn
> 
> plot(mymodel5@data$age, totvarfn)
> 
> xb <- predict(mymodel5)
> 
> u0 <- mymodel5@residual$lev_2_resi_est_Intercept
> u1 <- mymodel5@residual$lev_2_resi_est_age
> u2 <- mymodel5@residual[["lev_2_resi_est_I(age^2)"]]
> 
> yhat <- xb + u0[mymodel5@data$student] + u1[mymodel5@data$student] * mymodel5@data$age + u2[mymodel5@data$student] * 
+   mymodel5@data[["I(age^2)"]]
> 
> plot(mymodel5@data$age, yhat, type = "n")
> lines(mymodel5@data$age[mymodel5@data$student == 1], yhat[mymodel5@data$student == 1], col = 1)
> lines(mymodel5@data$age[mymodel5@data$student == 2], yhat[mymodel5@data$student == 2], col = 2)
> lines(mymodel5@data$age[mymodel5@data$student == 3], yhat[mymodel5@data$student == 3], col = 3)
> lines(mymodel5@data$age[mymodel5@data$student == 4], yhat[mymodel5@data$student == 4], col = 4)
> 
> plot(mymodel5@data$age, yhat, type = "n")
> lines(mymodel5@data$age[mymodel5@data$student == 10], yhat[mymodel5@data$student == 10], col = 1)
> lines(mymodel5@data$age[mymodel5@data$student == 11], yhat[mymodel5@data$student == 11], col = 2)
> lines(mymodel5@data$age[mymodel5@data$student == 12], yhat[mymodel5@data$student == 12], col = 3)
> lines(mymodel5@data$age[mymodel5@data$student == 13], yhat[mymodel5@data$student == 13], col = 4)
> lines(mymodel5@data$age[mymodel5@data$student == 14], yhat[mymodel5@data$student == 14], col = 4)
> 
> # Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . .209
> 
> ############################################################################
> 
> proc.time()
   user  system elapsed 
   4.39    0.35    7.20