R version 3.6.3 (2020-02-29) -- "Holding the Windsock"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

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.

  Natural language support but running in an English locale

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-2017 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.05/
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.05)  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 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.05)  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.04s 
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.05)  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.05s 
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.05)  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.06s 
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.51    0.59    5.35