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
> #
> # 5   Graphical Procedures for Exploring the Model . . . . . . . . . . . .65
> #
> #     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)
> 
> 
> # 5.1 Displaying multiple graphs . . . . . . . . . . . . . . . . . . . . .65
> 
> data(tutorial, package = "R2MLwiN")
> 
> (mymodel1 <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), estoptions = list(resi.store = TRUE), 
+   data = tutorial))
/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

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
school 65   2 62.44615 198         65            2      62.44615          198
Estimation algorithm:  IGLS        Elapsed time : 0.11s 
Number of obs:  4059 (from total 4059)        The model converged after 4 iterations.
Log likelihood:      -4678.6 
Deviance statistic:  9357.2 
--------------------------------------------------------------------------------------------------- 
The model formula:
normexam ~ 1 + standlrt + (1 | school) + (1 | student)
Level 2: school     Level 1: student      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
              Coef.   Std. Err.       z   Pr(>|z|)         [95% Conf.   Interval] 
Intercept   0.00239     0.04002    0.06     0.9524           -0.07605     0.08083 
standlrt    0.56337     0.01247   45.19          0   ***      0.53894     0.58780 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the school level: 
                  Coef.   Std. Err. 
var_Intercept   0.09213     0.01815 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                  Coef.   Std. Err. 
var_Intercept   0.56573     0.01266 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> u0 <- mymodel1@residual$lev_2_resi_est_Intercept
> u0se <- sqrt(mymodel1@residual$lev_2_resi_var_Intercept)
> u0CI95 <- 1.96 * u0se
> 
> u0rank <- rank(u0)
> u0rankhi <- u0 + u0CI95
> u0ranklo <- u0 - u0CI95
> u0rankno <- order(u0rank)
> 
> plot(1:65, u0[u0rankno], ylim = c(-1, 1), pch = 15, xlab = "Rank", ylab = "u0 residual estimate")
> points(1:65, u0rankhi[u0rankno], pch = 24, bg = "grey")
> points(1:65, u0ranklo[u0rankno], pch = 25, bg = "grey")
> for (i in 1:65) {
+   lines(rep(i, 2), c(u0ranklo[u0rankno[i]], u0rankhi[u0rankno[i]]))
+ }
> 
> plot(mymodel1@data$standlrt, mymodel1@data$normexam, asp = 1)
> 
> # 5.2 Highlighting in graphs . . . . . . . . . . . . . . . . . . . . . . .68
> 
> xb <- predict(mymodel1)
> 
> xbu <- xb + u0[mymodel1@data$school]
> 
> pred <- as.data.frame(cbind(mymodel1@data$school, mymodel1@data$standlrt, xb, xbu)[order(mymodel1@data$school, mymodel1@data$standlrt), 
+   ])
> 
> colnames(pred) <- c("school", "standlrt", "xb", "xbu")
> 
> if (!require(lattice)) {
+   warning("package lattice required to run this example")
+ } else {
+   xyplot(xbu ~ standlrt, type = "l", group = school, data = pred)
+ 
+   xyplot(xbu ~ standlrt, panel = function(x, y, subscripts) {
+     panel.xyplot(x, y, type = "l", groups = pred$school, subscripts = subscripts)
+     panel.xyplot(pred$standlrt, pred$xb, type = "l", lwd = 3, color = "black")
+   }, data = pred)
+ }
Loading required package: lattice
> 
> c(unique(mymodel1@data$school)[u0rank == 65], u0rank[u0rank == 65])
[1] 53 65
> 
> sch53 <- which(levels(as.factor(mymodel1@data$school)) == 53)
> 
> plot(1:65, u0[u0rankno], ylim = c(-1, 1), pch = 15, xlab = "Rank", ylab = "u0 residual estimate")
> points(1:65, u0rankhi[u0rankno], pch = 24, bg = "grey")
> points(1:65, u0ranklo[u0rankno], pch = 25, bg = "grey")
> for (i in 1:65) lines(rep(i, 2), c(u0ranklo[u0rankno[i]], u0rankhi[u0rankno[i]]))
> points(x = which(u0rankno == sch53), y = u0[u0rankno[which(u0rankno == sch53)]], pch = 22, bg = 2)
> legend(5, 1, "School 53", pch = 22, pt.bg = 2, col = 2)
> 
> plot(mymodel1@data$standlrt, xbu, type = "n")
> for (i in 1:65) {
+   lines(mymodel1@data$standlrt[mymodel1@data$school == i], xbu[mymodel1@data$school == i], col = "blue")
+ }
> lines(mymodel1@data$standlrt, xb, col = 1, lwd = 3)
> lines(mymodel1@data$standlrt[mymodel1@data$school == 53], xbu[mymodel1@data$school == 53], col = "red")
> legend(-3, 2, "School 53", lty = 1, col = "red")
> 
> plot(mymodel1@data$standlrt, mymodel1@data$normexam, type = "p")
> points(mymodel1@data$standlrt[mymodel1@data$school == 53], mymodel1@data$normexam[mymodel1@data$school == 53], col = "red")
> legend(-3, 3, "School 53", lty = 1, col = "red")
> 
> schid <- as.vector(by(mymodel1@data$school, mymodel1@data$school, function(x) x[1]))
> schsize <- as.vector(by(mymodel1@data$school, mymodel1@data$school, length))
> 
> cbind(schid, schsize, u0rank)[u0rank >= 28 & u0rank <= 32, ]
     schid schsize u0rank
[1,]     8     102     32
[2,]    47      82     30
[3,]    48       2     29
[4,]    51      58     28
[5,]    61      64     31
> 
> sch48 <- which(levels(as.factor(mymodel1@data$school)) == 48)
> 
> plot(1:65, u0[u0rankno], ylim = c(-1, 1), pch = 15, xlab = "Rank", ylab = "u0 residual estimate")
> points(1:65, u0rankhi[u0rankno], pch = 24, bg = "grey")
> points(1:65, u0ranklo[u0rankno], pch = 25, bg = "grey")
> for (i in 1:65) {
+   lines(rep(i, 2), c(u0ranklo[u0rankno[i]], u0rankhi[u0rankno[i]]))
+ }
> points(x = which(u0rankno == sch53), y = u0[u0rankno[which(u0rankno == sch53)]], pch = 22, bg = 2)
> points(x = which(u0rankno == sch48), y = u0[u0rankno[which(u0rankno == sch48)]], pch = 22, bg = 3)
> legend(5, 1, c("School 53", "School 48"), pch = 22, pt.bg = c(2, 3), col = c(2, 3))
> 
> plot(mymodel1@data$standlrt, xbu, type = "n")
> for (i in 1:65) {
+   lines(mymodel1@data$standlrt[mymodel1@data$school == i], xbu[mymodel1@data$school == i], col = "blue")
+ }
> lines(mymodel1@data$standlrt, xb, col = 1, lwd = 3)
> lines(mymodel1@data$standlrt[mymodel1@data$school == 53], xbu[mymodel1@data$school == 53], col = "red")
> lines(mymodel1@data$standlrt[mymodel1@data$school == 48], xbu[mymodel1@data$school == 48], col = "green")
> legend(-3, 2, c("The average school", "School 53", "School 48"), lty = 1, col = c("black", "red", "green"))
> 
> plot(mymodel1@data$standlrt, mymodel1@data$normexam, type = "p")
> points(mymodel1@data$standlrt[mymodel1@data$school == 53], mymodel1@data$normexam[mymodel1@data$school == 53], col = "red")
> points(mymodel1@data$standlrt[mymodel1@data$school == 48], mymodel1@data$normexam[mymodel1@data$school == 48], col = "green")
> legend(-3, 3, c("School 53", "School 48"), lty = 1, col = c("red", "green"))
> 
> cbind(schid, u0rank)[u0rank == 1, ]
 schid u0rank 
    59      1 
> 
> sch59 <- which(levels(as.factor(mymodel1@data$school)) == 59)
> 
> plot(1:65, u0[u0rankno], ylim = c(-1, 1), pch = 15, xlab = "Rank", ylab = "u0 residual estimate")
> points(1:65, u0rankhi[u0rankno], pch = 24, bg = "grey")
> points(1:65, u0ranklo[u0rankno], pch = 25, bg = "grey")
> for (i in 1:65) lines(rep(i, 2), c(u0ranklo[u0rankno[i]], u0rankhi[u0rankno[i]]))
> points(x = which(u0rankno == sch53), y = u0[u0rankno[which(u0rankno == sch53)]], pch = 22, bg = 2)
> points(x = which(u0rankno == sch59), y = u0[u0rankno[which(u0rankno == sch59)]], pch = 22, bg = 3)
> legend(5, 1, c("School 53", "School 59"), pch = 22, pt.bg = c(2, 3), col = c(2, 3))
> 
> plot(mymodel1@data$standlrt, xbu, type = "n")
> for (i in 1:65) {
+   lines(mymodel1@data$standlrt[mymodel1@data$school == i], xbu[mymodel1@data$school == i], col = "blue")
+ }
> lines(mymodel1@data$standlrt, xb, col = 1, lwd = 3)
> lines(mymodel1@data$standlrt[mymodel1@data$school == 53], xbu[mymodel1@data$school == 53], col = "red")
> lines(mymodel1@data$standlrt[mymodel1@data$school == 59], xbu[mymodel1@data$school == 59], col = "green")
> legend(-3, 2, c("The average school", "School 53", "School 59"), lty = 1, col = c("black", "red", "green"))
> 
> plot(mymodel1@data$standlrt, mymodel1@data$normexam, type = "p")
> points(mymodel1@data$standlrt[mymodel1@data$school == 53], mymodel1@data$normexam[mymodel1@data$school == 53], col = "red")
> points(mymodel1@data$standlrt[mymodel1@data$school == 59], mymodel1@data$normexam[mymodel1@data$school == 59], col = "green")
> legend(-3, 3, c("School 53", "School 59"), lty = 1, col = c("red", "green"))
> 
> (mymodel2 <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student), estoptions = list(resi.store = TRUE, 
+   resioptions = c("estimates", "sampling")), data = tutorial))
/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

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
school 65   2 62.44615 198         65            2      62.44615          198
Estimation algorithm:  IGLS        Elapsed time : 0.14s 
Number of obs:  4059 (from total 4059)        The model converged after 4 iterations.
Log likelihood:      -4658.4 
Deviance statistic:  9316.9 
--------------------------------------------------------------------------------------------------- 
The model formula:
normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student)
Level 2: school     Level 1: student      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
               Coef.   Std. Err.       z     Pr(>|z|)         [95% Conf.   Interval] 
Intercept   -0.01151     0.03978   -0.29       0.7724           -0.08948     0.06647 
standlrt     0.55673     0.01994   27.92   1.344e-171   ***      0.51765     0.59581 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the school level: 
                           Coef.   Std. Err. 
var_Intercept            0.09044     0.01792 
cov_Intercept_standlrt   0.01804     0.00672 
var_standlrt             0.01454     0.00441 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                  Coef.   Std. Err. 
var_Intercept   0.55366     0.01248 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> xb <- predict(mymodel2)
> 
> u0 <- mymodel2@residual$lev_2_resi_est_Intercept
> u1 <- mymodel2@residual$lev_2_resi_est_standlrt
> 
> xbu <- xb + u0[mymodel2@data$school] + u1[mymodel2@data$school] * mymodel2@data$standlrt
> 
> plot(u1, u0, xlab = "Slope", ylab = "Intercept")
> points(u1[schid == 53], u0[schid == 53], col = "red")
> points(u1[schid == 59], u0[schid == 59], col = "green")
> legend(0.2, -0.5, c("School 53", "School 59"), pch = 22, pt.bg = c("red", "green"), col = c("red", "green"))
> 
> plot(mymodel2@data$standlrt, xbu, type = "n")
> for (i in 1:65) {
+   lines(mymodel2@data$standlrt[mymodel2@data$school == i], xbu[mymodel2@data$school == i], col = "blue")
+ }
> lines(mymodel2@data$standlrt, xb, col = 1, lwd = 3)
> lines(mymodel2@data$standlrt[mymodel2@data$school == 53], xbu[mymodel2@data$school == 53], col = "red")
> lines(mymodel2@data$standlrt[mymodel2@data$school == 59], xbu[mymodel2@data$school == 59], col = "green")
> legend(-3, 2, c("School 53", "School 59"), lty = 1, col = c("red", "green"))
> 
> plot(mymodel2@data$standlrt, mymodel2@data$normexam, type = "p")
> points(mymodel2@data$standlrt[mymodel2@data$school == 53], mymodel2@data$normexam[mymodel2@data$school == 53], col = "red")
> points(mymodel2@data$standlrt[mymodel2@data$school == 59], mymodel2@data$normexam[mymodel2@data$school == 59], col = "green")
> legend(-3, 3, c("School 53", "School 59"), lty = 1, col = c("red", "green"))
> 
> u0var <- mymodel2@residual$lev_2_resi_var_Intercept
> u0u1cov <- mymodel2@residual$lev_2_resi_cov_standlrt_Intercep
> u1var <- mymodel2@residual$lev_2_resi_var_standlrt
> 
> xbu_lo <- xbu - 1.96 * sqrt(u0var[mymodel2@data$school] + 2 * u0u1cov[mymodel2@data$school] * mymodel2@data$standlrt + 
+   u1var[mymodel2@data$school] * mymodel2@data$standlrt^2)
> xbu_hi <- xbu + 1.96 * sqrt(u0var[mymodel2@data$school] + 2 * u0u1cov[mymodel2@data$school] * mymodel2@data$standlrt + 
+   u1var[mymodel2@data$school] * mymodel2@data$standlrt^2)
> 
> plotdata <- as.data.frame(cbind(mymodel2@data$standlrt, mymodel2@data$school, xb, xbu, xbu_lo, xbu_hi)[order(mymodel2@data$standlrt), 
+   ])
> colnames(plotdata) <- c("standlrt", "school", "xb", "xbu", "xbu_lo", "xbu_hi")
> 
> plot(plotdata$standlrt, plotdata$xb, xlim = c(-4, 3), ylim = c(-2.5, 3), type = "l")
> lines(plotdata$standlrt[plotdata$school == 53], plotdata$xbu[plotdata$school == 53], col = "red")
> lines(plotdata$standlrt[plotdata$school == 53], plotdata$xbu_lo[plotdata$school == 53], lty = 3, col = "red")
> lines(plotdata$standlrt[plotdata$school == 53], plotdata$xbu_hi[plotdata$school == 53], lty = 3, col = "red")
> lines(plotdata$standlrt[plotdata$school == 59], plotdata$xbu[plotdata$school == 59], col = "green")
> lines(plotdata$standlrt[plotdata$school == 59], plotdata$xbu_lo[plotdata$school == 59], lty = 3, col = "green")
> lines(plotdata$standlrt[plotdata$school == 59], plotdata$xbu_hi[plotdata$school == 59], lty = 3, col = "green")
> legend(-4, 3, c("School 53", "School 59"), lty = 1, col = c("red", "green"))
> 
> #     Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . 77
> 
> ############################################################################
> 
> proc.time()
   user  system elapsed 
   5.31    0.42    6.51