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 > # > # 16 An Introduction to Simulation Methods of Estimation . . . . . . . .241 > # > # 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) > > > # 16.1 An illustration of parameter estimation with Normally distributed . . > # . . .data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .242 > > data(height, package = "R2MLwiN") > summary(height) height Min. :154.0 1st Qu.:168.0 Median :175.0 Mean :175.3 3rd Qu.:181.0 Max. :204.0 > > hist(height$height) > > 1 - pnorm((200 - mean(height$height))/sd(height$height)) [1] 0.006861973 > > heightsim1 <- function() { + heightsim <- 175.35 + 10.002 * qnorm(runif(100)) + c(pmean = mean(heightsim), pvar = var(heightsim)) + } > > set.seed(1) > > # Note: To obtain estimates as close as possible to the manual, increase the > # number of reps to 10000. > > simdata1 <- as.data.frame(t(replicate(1000, heightsim1()))) > simdata1$iteration <- 1:nrow(simdata1) > > plot(simdata1$iteration, simdata1$pmean, type = "l") > > plot(density(simdata1$pmean)) > > quantile(simdata1$pmean, c(0.025, 0.975)) 2.5% 97.5% 173.3495 177.1498 > > plot(simdata1$iteration, simdata1$pvar, type = "l") > > plot(density(simdata1$pvar)) > > quantile(simdata1$pvar, c(0.025, 0.975)) 2.5% 97.5% 75.53613 129.95900 > > heightsim2 <- function(variable) { + samp <- sample(variable, replace = TRUE) + c(npmean = mean(samp), npvar = var(samp)) + } > > simdata2 <- as.data.frame(t(replicate(1000, heightsim2(height$height)))) > simdata2$iteration <- 1:nrow(simdata2) > > plot(simdata2$iteration, simdata2$npmean, type = "l") > > plot(density(simdata2$npmean)) > > quantile(simdata2$npmean, c(0.025, 0.975)) 2.5% 97.5% 173.53 177.38 > > plot(simdata2$iteration, simdata2$npvar, type = "l") > > plot(density(simdata2$npvar)) > > quantile(simdata2$npvar, c(0.025, 0.975)) 2.5% 97.5% 73.26895 128.29689 > > # 16.2 Generating random numbers in MLwiN . . . . . . . . . . . . . . . .249 > > female <- as.integer(runif(100) <= 0.6) > > height2 <- (1 - female) * (175 + 10 * qnorm(runif(100))) + female * (160 + 8 * qnorm(runif(100))) > > hist(height2[female == 0]) > hist(height2[female == 1]) > > # Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . .253 > > ############################################################################ > > proc.time() user system elapsed 3.10 0.29 3.57