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
> #
> # 15  Diagnostics for Multilevel Models . . . . . . . . . . . . . . . . .225
> #
> #     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)
> 
> 
> # 15.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .225
> 
> data(diag1, package = "R2MLwiN")
> summary(diag1)
     school            sex              vrq                 ilea      
 Min.   : 1.000   Min.   :0.0000   Min.   :-27.54000   Min.   : 0.00  
 1st Qu.: 4.000   1st Qu.:0.0000   1st Qu.: -9.54000   1st Qu.: 8.50  
 Median : 8.000   Median :0.0000   Median : -0.54000   Median :21.00  
 Mean   : 8.656   Mean   :0.4741   Mean   :  0.00024   Mean   :21.24  
 3rd Qu.:13.000   3rd Qu.:1.0000   3rd Qu.:  9.46000   3rd Qu.:31.00  
 Max.   :18.000   Max.   :1.0000   Max.   : 42.46000   Max.   :68.00  
            type         pupil            cons       n_ilea        
 Comprehensive:864   Min.   :  1.0   Min.   :1   Min.   :-1.43198  
 Grammar      : 43   1st Qu.:227.5   1st Qu.:1   1st Qu.:-0.67559  
                     Median :454.0   Median :1   Median : 0.01658  
                     Mean   :454.0   Mean   :1   Mean   : 0.01727  
                     3rd Qu.:680.5   3rd Qu.:1   3rd Qu.: 0.65125  
                     Max.   :907.0   Max.   :1   Max.   : 3.06113  
     n_vrq          
 Min.   :-2.329253  
 1st Qu.:-0.659817  
 Median : 0.023493  
 Mean   : 0.001778  
 3rd Qu.: 0.664979  
 Max.   : 3.262956  
> 
> (mymodel1 <- runMLwiN(n_ilea ~ 1 + n_vrq + (1 + n_vrq | school) + (1 | pupil), estoptions = list(resi.store = TRUE, 
+   resioptions = c("standardised", "leverage", "influence", "deletion")), data = diag1))
/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 18  21 50.38889  79         18           21      50.38889           79
Estimation algorithm:  IGLS        Elapsed time : 0.11s 
Number of obs:  907 (from total 907)        The model converged after 4 iterations.
Log likelihood:      -847.1 
Deviance statistic:  1694.3 
--------------------------------------------------------------------------------------------------- 
The model formula:
n_ilea ~ 1 + n_vrq + (1 + n_vrq | school) + (1 | pupil)
Level 2: school     Level 1: pupil      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
              Coef.   Std. Err.       z    Pr(>|z|)         [95% Conf.   Interval] 
Intercept   0.00679     0.03082    0.22      0.8257           -0.05362     0.06720 
n_vrq       0.72543     0.03632   19.97   9.717e-89   ***      0.65424     0.79662 
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.00900     0.00564 
cov_Intercept_n_vrq   0.00991     0.00524 
var_n_vrq             0.01559     0.00789 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the pupil level: 
                  Coef.   Std. Err. 
var_Intercept   0.36740     0.01757 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> xb <- predict(mymodel1)
> 
> u0 <- mymodel1@residual$lev_2_resi_est_Intercept
> u1 <- mymodel1@residual$lev_2_resi_est_n_vrq
> 
> yhat <- xb + u0[mymodel1@data$school] + u1[mymodel1@data$school] * mymodel1@data$n_vrq
> 
> plot(mymodel1@data$n_vrq, yhat, type = "n")
> for (i in 1:18) {
+   lines(mymodel1@data$n_vrq[mymodel1@data$school == i], yhat[mymodel1@data$school == i], col = 1)
+ }
> lines(mymodel1@data$n_vrq[mymodel1@data$school == 17], yhat[mymodel1@data$school == 17], col = 2, lwd = 3)
> 
> yhatold <- yhat
> 
> plot(mymodel1@data$n_vrq[mymodel1@data$school != 17], mymodel1@data$n_ilea[mymodel1@data$school != 17], type = "p")
> points(mymodel1@data$n_vrq[mymodel1@data$school == 17], mymodel1@data$n_ilea[mymodel1@data$school == 17], col = 2, 
+   cex = 2)
> 
> u0se <- mymodel1@residual$lev_2_resi_se_Intercept
> u1se <- mymodel1@residual$lev_2_resi_se_n_vrq
> 
> u0CI95 <- 1.96 * u0se
> u0rank <- rank(u0)
> u0rankhi <- u0 + u0CI95
> u0ranklo <- u0 - u0CI95
> u0rankno <- order(u0rank)
> 
> u1CI95 <- 1.96 * u0se
> u1rank <- rank(u1)
> u1rankhi <- u1 + u1CI95
> u1ranklo <- u1 - u1CI95
> u1rankno <- order(u1rank)
> 
> sch17 <- which(levels(as.factor(mymodel1@data$school)) == 17)
> 
> plot(1:18, u0[u0rankno], ylim = c(-0.3, 0.3), pch = 15, xlab = "Rank", ylab = "u0 residual estimate")
> points(1:18, u0rankhi[u0rankno], pch = 24, bg = "grey")
> points(1:18, u0ranklo[u0rankno], pch = 25, bg = "grey")
> for (i in 1:18) lines(rep(i, 2), c(u0ranklo[u0rankno[i]], u0rankhi[u0rankno[i]]))
> points(x = which(u0rankno == sch17), y = u0[u0rankno[which(u0rankno == sch17)]], pch = 22, bg = 2, cex = 2)
> 
> plot(1:18, u1[u1rankno], ylim = c(-0.3, 0.3), pch = 15, xlab = "Rank", ylab = "u1 residual estimate")
> points(1:18, u1rankhi[u1rankno], pch = 24, bg = "grey")
> points(1:18, u1ranklo[u1rankno], pch = 25, bg = "grey")
> for (i in 1:18) lines(rep(i, 2), c(u1ranklo[u1rankno[i]], u1rankhi[u1rankno[i]]))
> points(x = which(u1rankno == sch17), y = u1[u1rankno[which(u1rankno == sch17)]], pch = 22, bg = 2, cex = 2)
> 
> plot(u0[-17], u1[-17], ylim = c(-0.3, 0.3), xlim = c(-0.1, 0.2), type = "p", xlab = "u0 residual estimate", ylab = "u1 residual estimate")
> points(u0[17], u1[17], bg = 2, col = 2, cex = 2)
> 
> # 15.2 Diagnostics plotting: Deletion residuals, influence and leverage . 231
> 
> hist(u0, breaks = seq(min(u0) - 0.01, max(u0) + 0.01, 0.01))
> 
> u0std <- mymodel1@residual$lev_2_std_resi_est_Intercept
> hist(u0std, breaks = seq(min(u0std) - 0.1, max(u0std) + 0.1, 0.1))
> 
> u0lev <- mymodel1@residual$lev_2_resi_leverage_Intercept
> hist(u0lev, breaks = seq(min(u0lev) - 0.01, max(u0lev) + 0.01, 0.01))
> 
> u0inf <- mymodel1@residual$lev_2_resi_influence_Intercept
> hist(u0inf, breaks = seq(min(u0inf) - 0.025, max(u0inf) + 0.025, 0.025))
> 
> u0del <- mymodel1@residual$lev_2_resi_deletion_Intercept
> hist(u0del, breaks = seq(min(u0del) - 0.1, max(u0del) + 0.1, 0.1))
> 
> plot(u0std[-17], u0lev[-17], ylim = c(0.15, 0.4), xlim = c(-0.2, 2.2), type = "p", xlab = "u0 standardised residual", 
+   ylab = "u0 leverage residual")
> points(u0std[17], u0lev[17], bg = 2, col = 2, cex = 2)
> 
> hist(u1, breaks = seq(min(u1) - 0.01, max(u1) + 0.01, 0.01))
> 
> u1std <- mymodel1@residual$lev_2_std_resi_est_n_vrq
> hist(u1std, breaks = seq(min(u1std) - 0.1, max(u1std) + 0.1, 0.1))
> 
> u1lev <- mymodel1@residual$lev_2_resi_leverage_n_vrq
> hist(u1lev, breaks = seq(min(u1lev) - 0.01, max(u1lev) + 0.01, 0.01))
> 
> u1inf <- mymodel1@residual$lev_2_resi_influence_n_vrq
> hist(u1inf, breaks = seq(min(u1inf) - 0.025, max(u1inf) + 0.025, 0.025))
> 
> u1del <- mymodel1@residual$lev_2_resi_deletion_n_vrq
> hist(u1del, breaks = seq(min(u1del) - 0.1, max(u1del) + 0.1, 0.1))
> 
> plot(u1std[-17], u1lev[-17], ylim = c(0.1, 0.35), xlim = c(-1.4, 2.4), type = "p", xlab = "u1 standardised residual", 
+   ylab = "u1 leverage residual")
> points(u1std[17], u1lev[17], bg = 2, col = 2, cex = 2)
> 
> e0 <- mymodel1@residual$lev_1_resi_est_Intercept
> e0rank <- rank(e0)
> e0std <- (e0 - mean(e0))/sd(e0)
> e0uniform <- e0rank/(length(e0rank) + 1)
> e0nscore <- qnorm(e0uniform)
> 
> plot(e0nscore[mymodel1@data$school != 17], e0std[mymodel1@data$school != 17], ylim = c(-4, 5), xlim = c(-4, 4), type = "p", 
+   xlab = "e0nscore", ylab = "e0std")
> points(e0nscore[mymodel1@data$school == 17], e0std[mymodel1@data$school == 17], bg = 2, col = 2, cex = 2)
> 
> diag1$pupilnumber <- unlist(by(diag1$school, diag1$school, function(x) 1:length(x)))
> 
> diag1[order(e0)[1], c("school", "pupil", "pupilnumber")]
    school pupil pupilnumber
886     17   886          22
> 
> diag1$s17p22 <- as.integer(diag1$school == 17 & diag1$pupilnumber == 22)
> 
> (mymodel2 <- runMLwiN(n_ilea ~ 1 + n_vrq + s17p22 + (1 + n_vrq | school) + (1 | pupil), estoptions = list(resi.store = TRUE, 
+   startval = list(FP.b = mymodel1@FP, FP.v = mymodel1@FP.cov, RP.b = mymodel1@RP, RP.v = mymodel1@RP.cov)), data = diag1))
/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 0.00678814289217072 0.725429729819978 0 '_FP_b'
JOIN 0.000950040114645168 0.000559145844151685 0.00131935737838271 0 0 0 '_FP_v'
JOIN 0.00900486942250067 0.00990715226493238 0.0155909529719536 0.367397093219737 '_RP_b'
JOIN 3.17831572748176e-05 1.95036244206765e-05 2.7436609719712e-05 1.01505178752961e-05 2.65508582965033e-05 6.23168796461516e-05 -5.62507380372694e-06 -6.94461782252319e-07 -5.95683596423323e-06 0.000308600156539668 '_RP_v'
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4

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 18  21 50.38889  79         18           21      50.38889           79
Estimation algorithm:  IGLS        Elapsed time : 0.1s 
Number of obs:  907 (from total 907)        The model converged after 5 iterations.
Log likelihood:      -843.1 
Deviance statistic:  1686.2 
--------------------------------------------------------------------------------------------------- 
The model formula:
n_ilea ~ 1 + n_vrq + s17p22 + (1 + n_vrq | school) + (1 | pupil)
Level 2: school     Level 1: pupil      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
               Coef.   Std. Err.       z    Pr(>|z|)         [95% Conf.   Interval] 
Intercept    0.00781     0.03139    0.25      0.8034           -0.05370     0.06933 
n_vrq        0.72807     0.03748   19.43   4.529e-84   ***      0.65462     0.80153 
s17p22      -1.76586     0.61792   -2.86    0.004266   **      -2.97695    -0.55477 
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.00968     0.00587 
cov_Intercept_n_vrq   0.01085     0.00558 
var_n_vrq             0.01717     0.00842 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the pupil level: 
                  Coef.   Std. Err. 
var_Intercept   0.36351     0.01738 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> u0 <- mymodel2@residual$lev_2_resi_est_Intercept
> u1 <- mymodel2@residual$lev_2_resi_est_n_vrq
> 
> u0se <- sqrt(mymodel2@residual$lev_2_resi_var_Intercept)
> u1se <- sqrt(mymodel2@residual$lev_2_resi_var_n_vrq)
> 
> u0CI95 <- 1.96 * u0se
> u0rank <- rank(u0)
> u0rankhi <- u0 + u0CI95
> u0ranklo <- u0 - u0CI95
> u0rankno <- order(u0rank)
> 
> u1CI95 <- 1.96 * u1se
> u1rank <- rank(u1)
> u1rankhi <- u1 + u1CI95
> u1ranklo <- u1 - u1CI95
> u1rankno <- order(u1rank)
> 
> plot(1:18, u0[u0rankno], ylim = c(-0.3, 0.3), pch = 15, xlab = "Rank", ylab = "u0 residual estimate")
> points(1:18, u0rankhi[u0rankno], pch = 24, bg = "grey")
> points(1:18, u0ranklo[u0rankno], pch = 25, bg = "grey")
> for (i in 1:18) lines(rep(i, 2), c(u0ranklo[u0rankno[i]], u0rankhi[u0rankno[i]]))
> points(x = which(u0rankno == sch17), y = u0[u0rankno[which(u0rankno == sch17)]], pch = 22, bg = 2, cex = 2)
> 
> plot(1:18, u1[u1rankno], ylim = c(-0.3, 0.3), pch = 15, xlab = "Rank", ylab = "u1 residual estimate")
> points(1:18, u1rankhi[u1rankno], pch = 24, bg = "grey")
> points(1:18, u1ranklo[u1rankno], pch = 25, bg = "grey")
> for (i in 1:18) lines(rep(i, 2), c(u1ranklo[u1rankno[i]], u1rankhi[u1rankno[i]]))
> points(x = which(u1rankno == sch17), y = u1[u1rankno[which(u1rankno == sch17)]], pch = 22, bg = 2, cex = 2)
> 
> diag1$s17 <- as.integer(diag1$school == 17)
> diag1$s17Xn_vrq <- diag1$s17 * diag1$n_vrq
> 
> (mymodel3 <- runMLwiN(n_ilea ~ 1 + n_vrq + s17p22 + s17 + s17Xn_vrq + (1 + n_vrq | school) + (1 | pupil), estoptions = list(startval = list(FP.b = mymodel2@FP, 
+   FP.v = mymodel2@FP.cov, RP.b = mymodel2@RP, RP.v = mymodel2@RP.cov)), data = diag1))
/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 0.00781355596708382 0.728074202559495 -1.76585970936219 0 0 '_FP_b'
JOIN 0.000985050537488226 0.000611808950430028 0.00140450397250478 -0.000296089693912399 -0.000546871723951445 0.381819278042617 0 0 0 0 0 0 0 0 0 '_FP_v'
JOIN 0.00968086421081849 0.0108481219117231 0.0171700307410536 0.363513347691856 '_RP_b'
JOIN 3.44107180127172e-05 2.22894176495363e-05 3.1085943038604e-05 1.25771820910112e-05 3.11855684771216e-05 7.0818712210117e-05 -5.5377593090256e-06 -6.93373041165234e-07 -5.85926118146108e-06 0.000302117404489703 '_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
school 18  21 50.38889  79         18           21      50.38889           79
Estimation algorithm:  IGLS        Elapsed time : 0.08s 
Number of obs:  907 (from total 907)        The model converged after 6 iterations.
Log likelihood:      -837.8 
Deviance statistic:  1675.5 
--------------------------------------------------------------------------------------------------- 
The model formula:
n_ilea ~ 1 + n_vrq + s17p22 + s17 + s17Xn_vrq + (1 + n_vrq | 
    school) + (1 | pupil)
Level 2: school     Level 1: pupil      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
               Coef.   Std. Err.       z    Pr(>|z|)         [95% Conf.   Interval] 
Intercept   -0.00687     0.02812   -0.24      0.8069           -0.06200     0.04825 
n_vrq        0.70848     0.03343   21.19   1.116e-99   ***      0.64296     0.77400 
s17p22      -1.76047     0.62183   -2.83    0.004639   **      -2.97923    -0.54170 
s17          1.25280     0.47913    2.61     0.00893   **       0.31372     2.19189 
s17Xn_vrq   -0.27708     0.30028   -0.92      0.3561           -0.86561     0.31145 
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.00584     0.00452 
cov_Intercept_n_vrq   0.00597     0.00398 
var_n_vrq             0.01098     0.00630 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the pupil level: 
                  Coef.   Std. Err. 
var_Intercept   0.36183     0.01730 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> (mymodel4 <- runMLwiN(n_ilea ~ 1 + n_vrq + s17p22 + s17 + (1 + n_vrq | school) + (1 | pupil), estoptions = list(resi.store = TRUE, 
+   startval = list(FP.b = mymodel3@FP, FP.v = mymodel3@FP.cov, RP.b = mymodel3@RP, RP.v = mymodel3@RP.cov)), data = diag1))
/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 -0.00687493689426086 0.708478528516249 -1.76046506874345 1.25280447388533 '_FP_b'
JOIN 0.000790992066755186 0.000377064818875047 0.00111758909440556 1.35665290713598e-19 5.62684607613758e-19 0.386672795136772 -0.000790992066755184 -0.000377064818875042 0.0223397865766697 0.229569985677798 '_FP_v'
JOIN 0.0058384238455603 0.00597469259084182 0.0109832541821904 0.36183120387269 '_RP_b'
JOIN 2.03971349519669e-05 9.66068247406285e-06 1.5859927794676e-05 3.01139677063176e-06 1.27522725775572e-05 3.97199013437828e-05 -5.38848173218383e-06 -5.82295398950229e-07 -5.67986100231665e-06 0.000299234566709773 '_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
school 18  21 50.38889  79         18           21      50.38889           79
Estimation algorithm:  IGLS        Elapsed time : 0.1s 
Number of obs:  907 (from total 907)        The model converged after 6 iterations.
Log likelihood:      -838.2 
Deviance statistic:  1676.4 
--------------------------------------------------------------------------------------------------- 
The model formula:
n_ilea ~ 1 + n_vrq + s17p22 + s17 + (1 + n_vrq | school) + (1 | 
    pupil)
Level 2: school     Level 1: pupil      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
               Coef.   Std. Err.       z    Pr(>|z|)         [95% Conf.   Interval] 
Intercept   -0.00806     0.02813   -0.29      0.7744           -0.06320     0.04707 
n_vrq        0.70508     0.03337   21.13   4.378e-99   ***      0.63967     0.77049 
s17p22      -1.83525     0.61676   -2.98    0.002924   **      -3.04408    -0.62643 
s17          0.88372     0.26347    3.35   0.0007962   ***      0.36732     1.40012 
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.00586     0.00453 
cov_Intercept_n_vrq   0.00604     0.00401 
var_n_vrq             0.01115     0.00636 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the pupil level: 
                  Coef.   Std. Err. 
var_Intercept   0.36211     0.01731 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> xb <- predict(mymodel4)
> 
> u0 <- mymodel4@residual$lev_2_resi_est_Intercept
> u1 <- mymodel4@residual$lev_2_resi_est_n_vrq
> 
> yhat <- xb + u0[mymodel4@data$school] + u1[mymodel4@data$school] * mymodel4@data$n_vrq
> 
> plot(mymodel4@data$n_vrq, yhat, type = "n")
> for (i in 1:18) {
+   lines(mymodel4@data$n_vrq[mymodel4@data$school == i], yhatold[mymodel4@data$school == i], col = 1)
+ }
> lines(mymodel4@data$n_vrq[mymodel4@data$school == 17], yhatold[mymodel4@data$school == 17], col = 2, lwd = 3)
> 
> plot(mymodel4@data$n_vrq, yhat, type = "n")
> for (i in 1:18) {
+   points(mymodel4@data$n_vrq[mymodel4@data$school == i], yhat[mymodel4@data$school == i], col = 1)
+ }
> points(mymodel4@data$n_vrq[mymodel4@data$school == 17], yhat[mymodel4@data$school == 17], col = 2, lwd = 3)
> 
> # 15.3 A general approach to data exploration . . . . . . . . . . . . . .240
> 
> # Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . .240
> 
> ############################################################################
> 
> proc.time()
   user  system elapsed 
   3.87    0.42    5.78