demo/UserGuide16.R

In R2MLwiN: Running 'MLwiN' from Within 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)
# 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)

hist(height$height)

1 - pnorm((200 - mean(height$height))/sd(height$height))

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))

plot(simdata1$iteration, simdata1$pvar, type = "l")

plot(density(simdata1$pvar))

quantile(simdata1$pvar, c(0.025, 0.975))

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))

plot(simdata2$iteration, simdata2$npvar, type = "l")

plot(density(simdata2$npvar))

quantile(simdata2$npvar, c(0.025, 0.975))

# 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

############################################################################