demo/MCMCGuide05.R

In R2MLwiN: Running 'MLwiN' from Within R

############################################################################
#     MLwiN MCMC Manual
#
# 5   Prior Distributions, Starting Values and Random Number Seeds . . . .61
#
#     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/
############################################################################

# 5.1 Prior distributions . . . . . . . . . . . . . . . . . . . . . . . . 61

# 5.2 Uniform on variance scale priors . . . . . . . . . . . . . . . . . .61

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)

# User's input if necessary

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

## IGLS
(mymodel1 <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), data = tutorial))

## Diffuse priors (Gamma priors)
(mymodel2 <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), estoptions = list(EstM = 1), data = tutorial))

## Diffuse priors (Uniform priors)
(mymodel3 <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), estoptions = list(EstM = 1, mcmcMeth = list(priorcode = 0)), 
  data = tutorial))

if (!require(texreg)) {
  warning("texreg package required to use screenreg() function")
} else {
  screenreg(list(mymodel1, mymodel2, mymodel3), custom.model.names=c("IGLS", "Gibbs1", "Gibbs2"), groups = list("Fixed Part" = 1:2, "Level-2" = 3:3, "Level-1" = 4:4),
   stars = numeric(0), include.nobs=FALSE, include.loglik=FALSE, include.deviance=FALSE, include.dbar=FALSE, include.dthetabar=FALSE, include.pd=FALSE, include.dic=FALSE)
}
rm(list = c("mymodel1", "mymodel2", "mymodel3"))

# 5.3 Using informative priors . . . . . . . . . . . . . . . . . . . . . .62

## Informative normal prior for beta_1 Fit the model
(mymodel4 <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), estoptions = list(EstM = 1, mcmcMeth = list(priorParam = list(fixe = list(standlrt = c(1, 
  0.01))))), data = tutorial))

sixway(mymodel4@chains[, "FP_standlrt", drop = FALSE], "beta_1")

## Informative normal prior for beta_1
(mymodel5 <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), estoptions = list(EstM = 1, mcmcMeth = list(priorParam = list(fixe = list(standlrt = c(1, 
  0.1))))), data = tutorial))

sixway(mymodel5@chains[, "FP_standlrt", drop = FALSE], "beta_1")

# 5.4 Specifying an informative prior for a random parameter . . . . . . .65

## Specifies an informative prior for sigma2u
(mymodel6 <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), estoptions = list(EstM = 1, mcmcMeth = list(priorParam = list(rp2 = list(estimate = 0.2, 
  size = 100)))), data = tutorial))

sixway(mymodel6@chains[, "RP2_var_Intercept", drop = FALSE], "sigma^2_u0")

# 5.5 Changing the random number seed and the parameter starting values  .66

## Set starting values for random and fixed parameter estimates
FP.b <- c(-2, 5)
names(FP.b) <- c("FP_Intercept", "FP_standlrt")
RP.b <- c(2, 4)
names(RP.b) <- c("RP2_var_Intercept", "RP1_var_Intercept")
startval <- list(FP.b = FP.b, RP.b = RP.b)

(mymodel7 <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), estoptions = list(EstM = 1, mcmcMeth = list(burnin = 0, 
  iterations = 500), startval = startval), data = tutorial))

rm(list = c("mymodel4", "mymodel5", "mymodel6", "mymodel7"))

## Use different seeds
(mymodel8 <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), estoptions = list(EstM = 1, mcmcMeth = list(nchains=4, seed = 1:4)),
  data = tutorial))

ctable <- round(as.data.frame(lapply(mymodel8@chains[, c(names(mymodel8@FP), names(mymodel8@RP))], colMeans)), 3)
colnames(ctable) <- c("Seed1", "Seed2", "Seed3", "Seed4")
print(ctable)
rm(list = "mymodel8")

# 5.6 Improving the speed of MCMC Estimation . . . . . . . . . . . . . . .69

# Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . . 70





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