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 > # > # 7 Modelling the Variance as a Function of Explanatory Variables . . . 89 > # > # 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) > > > # 7.1 A level 1 variance function for two groups . . . . . . . . . . . . .89 > > data(tutorial, package = "R2MLwiN") > > covmatrix <- matrix(, nrow = 3, ncol = 1) > covmatrix[1, 1] <- 1 > covmatrix[2, 1] <- "sexboy" > covmatrix[3, 1] <- "sexgirl" > > contrasts(tutorial$sex, 2) <- diag(2) > > (mymodel1 <- runMLwiN(normexam ~ 0 + sex + (0 + sex | student), estoptions = list(clre = covmatrix), data = tutorial)) /nogui option ignored ECHO 0 Echoing is ON BATC 1 Batch mode is ON MAXI 2 STAR iteration 0 iteration 1 Convergence achieved TOLE 2 MAXI 20 NEXT iteration 2 Convergence achieved ECHO 0 Execution completed -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- MLwiN (version: 3.05) multilevel model (Normal) Estimation algorithm: IGLS Elapsed time : 0.04s Number of obs: 4059 (from total 4059) The model converged after 3 iterations. Log likelihood: -5724.8 Deviance statistic: 11449.5 --------------------------------------------------------------------------------------------------- The model formula: normexam ~ 0 + sex + (0 + sex | student) Level 1: student --------------------------------------------------------------------------------------------------- The fixed part estimates: Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] sexboy -0.14035 0.02545 -5.51 3.504e-08 *** -0.19024 -0.09046 sexgirl 0.09332 0.01964 4.75 2.028e-06 *** 0.05482 0.13182 Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 --------------------------------------------------------------------------------------------------- The random part estimates at the student level: Coef. Std. Err. var_sexboy 1.05144 0.03691 var_sexgirl 0.93997 0.02693 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- > > contrasts(tutorial$sex, 1) <- c(0, 1) > > # 7.2 Variance functions at level 2 . . . . . . . . . . . . . . . . . . . 95 > > (mymodel2 <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student), 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.07s 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 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- > > l2varfn <- mymodel2@RP["RP2_var_Intercept"] + (2 * mymodel2@RP["RP2_cov_Intercept_standlrt"] * mymodel2@data$standlrt) + + (mymodel2@RP["RP2_var_standlrt"] * mymodel2@data$standlrt^2) > > varfndata <- as.data.frame(cbind(mymodel2@data$standlrt, l2varfn)[order(mymodel2@data$standlrt), ]) > colnames(varfndata) <- c("standlrt", "l2varfn") > > plot(varfndata$standlrt, varfndata$l2varfn, type = "l") > > # 7.3 Further elaborating the model for the student-level variance . . . .99 > > (mymodel3 <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 + standlrt | student), 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 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 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 6 iterations. Log likelihood: -4655.8 Deviance statistic: 9311.6 --------------------------------------------------------------------------------------------------- The model formula: normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 + standlrt | student) Level 2: school Level 1: student --------------------------------------------------------------------------------------------------- The fixed part estimates: Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] Intercept -0.01174 0.03985 -0.29 0.7683 -0.08984 0.06636 standlrt 0.55787 0.01982 28.14 3.092e-174 *** 0.51901 0.59672 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.09082 0.01799 cov_Intercept_standlrt 0.01863 0.00671 var_standlrt 0.01423 0.00437 --------------------------------------------------------------------------------------------------- The random part estimates at the student level: Coef. Std. Err. var_Intercept 0.55322 0.01521 cov_Intercept_standlrt -0.01474 0.00640 var_standlrt 0.00066 0.00911 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- > > l2varfn <- mymodel3@RP["RP2_var_Intercept"] + (2 * mymodel3@RP["RP2_cov_Intercept_standlrt"] * mymodel3@data$standlrt) + + (mymodel3@RP["RP2_var_standlrt"] * mymodel3@data$standlrt^2) > > l1varfn <- mymodel3@RP["RP1_var_Intercept"] + (2 * mymodel3@RP["RP1_cov_Intercept_standlrt"] * mymodel3@data$standlrt) + + (mymodel3@RP["RP1_var_standlrt"] * mymodel3@data$standlrt^2) > > varfndata <- as.data.frame(cbind(mymodel3@data$standlrt, l2varfn, l1varfn)[order(mymodel3@data$standlrt), ]) > colnames(varfndata) <- c("standlrt", "l2varfn", "l1varfn") > > if (!require(lattice)) { + warning("package lattice required to run this example") + } else { + xyplot(l2varfn + l1varfn ~ standlrt, data = varfndata, type = "l") + } Loading required package: lattice > > covmatrix <- matrix(, nrow = 3, ncol = 3) > covmatrix[1, 1] <- 1 > covmatrix[2, 1] <- "standlrt" > covmatrix[3, 1] <- "standlrt" > covmatrix[1, 2] <- 1 > covmatrix[2, 2] <- "sexgirl" > covmatrix[3, 2] <- "Intercept" > covmatrix[1, 3] <- 1 > covmatrix[2, 3] <- "standlrt" > covmatrix[3, 3] <- "sexgirl" > > (mymodel4 <- runMLwiN(normexam ~ 1 + standlrt + sex + (1 + standlrt | school) + (1 + standlrt + sex | student), estoptions = list(clre = covmatrix), + 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: -4638.7 Deviance statistic: 9277.4 --------------------------------------------------------------------------------------------------- The model formula: normexam ~ 1 + standlrt + sex + (1 + standlrt | school) + (1 + standlrt + sex | student) Level 2: school Level 1: student --------------------------------------------------------------------------------------------------- The fixed part estimates: Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] Intercept -0.11200 0.04331 -2.59 0.009706 ** -0.19689 -0.02712 standlrt 0.55417 0.01994 27.79 5.388e-170 *** 0.51509 0.59325 sexgirl 0.17580 0.03236 5.43 5.545e-08 *** 0.11238 0.23923 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.08653 0.01722 cov_Intercept_standlrt 0.01957 0.00665 var_standlrt 0.01456 0.00441 --------------------------------------------------------------------------------------------------- The random part estimates at the student level: Coef. Std. Err. var_Intercept 0.58261 0.02082 cov_Intercept_standlrt -0.01281 0.00635 var_sexgirl -0.05412 0.02589 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- > > covmatrix <- matrix(, nrow = 3, ncol = 2) > covmatrix[1, 1] <- 1 > covmatrix[2, 1] <- "standlrt" > covmatrix[3, 1] <- "standlrt" > covmatrix[1, 2] <- 1 > covmatrix[2, 2] <- "sexgirl" > covmatrix[3, 2] <- "Intercept" > > (mymodel5 <- runMLwiN(normexam ~ 1 + standlrt + sex + (1 + standlrt | school) + (1 + standlrt + sex | student), estoptions = list(clre = covmatrix), + 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: -4635.8 Deviance statistic: 9271.6 --------------------------------------------------------------------------------------------------- The model formula: normexam ~ 1 + standlrt + sex + (1 + standlrt | school) + (1 + standlrt + sex | student) Level 2: school Level 1: student --------------------------------------------------------------------------------------------------- The fixed part estimates: Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] Intercept -0.11206 0.04329 -2.59 0.009633 ** -0.19689 -0.02722 standlrt 0.55390 0.02003 27.66 2.237e-168 *** 0.51465 0.59315 sexgirl 0.17532 0.03229 5.43 5.647e-08 *** 0.11203 0.23860 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.08643 0.01720 cov_Intercept_standlrt 0.01957 0.00668 var_standlrt 0.01479 0.00445 --------------------------------------------------------------------------------------------------- The random part estimates at the student level: Coef. Std. Err. var_Intercept 0.58371 0.02090 cov_Intercept_standlrt -0.03368 0.00996 cov_standlrt_sexgirl 0.03209 0.01290 var_sexgirl -0.05827 0.02593 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- > > l2varfn <- mymodel5@RP["RP2_var_Intercept"] + (2 * mymodel5@RP["RP2_cov_Intercept_standlrt"] * mymodel5@data$standlrt) + + (mymodel5@RP["RP2_var_standlrt"] * mymodel5@data$standlrt^2) > > l1varfnboys <- mymodel5@RP["RP1_var_Intercept"] + (2 * mymodel5@RP["RP1_cov_Intercept_standlrt"] * mymodel5@data$standlrt) > > l1varfngirls <- mymodel5@RP["RP1_var_Intercept"] + (2 * mymodel5@RP["RP1_cov_Intercept_standlrt"] * mymodel5@data$standlrt) + + (2 * mymodel5@RP["RP1_cov_standlrt_sexgirl"] * mymodel5@data$standlrt) + mymodel5@RP["RP1_var_sexgirl"] > > varfndata <- as.data.frame(cbind(mymodel5@data$standlrt, l2varfn, l1varfnboys, l1varfngirls)[order(mymodel5@data$standlrt), + ]) > colnames(varfndata) <- c("standlrt", "l2varfn", "l1varfnboys", "l1varfngirls") > > if (!require(lattice)) { + warning("package lattice required to run this example") + } else { + xyplot(l2varfn + l1varfnboys + l1varfngirls ~ standlrt, data = varfndata, type = "l") + } > > # Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . .106 > > ############################################################################ > > proc.time() user system elapsed 3.76 0.51 5.43