demo/MCMCGuide02.R

In R2MLwiN: Running 'MLwiN' from Within R

############################################################################
#     MLwiN MCMC Manual
#
# 2   Single Level Normal Response Modelling . . . . . . . . . . . . . . .21
#
#     Browne, W.J. (2009) MCMC Estimation in MLwiN, v2.13. 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)

## save current par settings
mypar <- par(no.readonly = TRUE)

## Read tutorial data
data(tutorial, package = "R2MLwiN")

## Choose IGLS algoritm for estimation
(mymodel1 <- runMLwiN(normexam ~ 1 + standlrt + (1 | student), data = tutorial))

# 2.1 Running the Gibbs Sampler . . . . . . . . . . . . . . . . . . . . . 26

## Choose MCMC algoritm for estimation
(mymodel2 <- runMLwiN(normexam ~ 1 + standlrt + (1 | student), estoptions = list(EstM = 1), data = tutorial))


if (!require(coda)) {
  warning("package coda required to run this example")
} else {
  par(mfrow = c(2, 2))
  estimates <- mymodel2@chains
  plot(1:niter(estimates), estimates[, "deviance"], xlab = "iteration", ylab = expression(paste("Est. of deviance")), 
    type = "l")
  plot(1:niter(estimates), estimates[, "FP_Intercept"], xlab = "iteration", ylab = expression(paste("Est. of ", beta[0])), 
    type = "l")
  plot(1:niter(estimates), estimates[, "FP_standlrt"], xlab = "iteration", ylab = expression(paste("Est. of ", beta[1])), 
    type = "l")
  plot(1:niter(estimates), estimates[, "RP1_var_Intercept"], xlab = "iteration", ylab = expression(paste("Est. of ", 
    sigma[e0]^2)), type = "l")
  ## reinstate par settings
  par(mypar)
}

# 2.2 Deviance statistic and the DIC diagnostic . . . . . . . . . . . . . 28

# 2.3 Adding more predictors . . . . . . . . . . . . . . . . . . . . . . .29

## Choose IGLS algoritm for estimation
(mymodel3 <- runMLwiN(normexam ~ 1 + standlrt + sex + schgend + (1 | student), data = tutorial))

## Choose MCMC algoritm for estimation
(mymodel4 <- runMLwiN(normexam ~ 1 + standlrt + sex + schgend + (1 | student), estoptions = list(EstM = 1), data = tutorial))

# 2.4 Fitting school effects as fixed parameters . . . . . . . . . . . . .32

tutorial$school <- as.factor(tutorial$school)

## Choose MCMC algoritm for estimation (IGLS will be used to obtain starting values for MCMC)
(mymodel5 <- runMLwiN(normexam ~ 1 + standlrt + sex + school + (1 | student), estoptions = list(EstM = 1), data = tutorial))

# Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . . 33




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