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 > # > # 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-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="") > # 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.05) 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.09s 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.05) 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.1s 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 4.07 0.45 5.25