demo/MCMCGuide13.R

In R2MLwiN: Running 'MLwiN' from Within R

############################################################################
#     MLwiN MCMC Manual
#
# 13  Ordered Categorical Responses . . . . . . . . . . . . . . . . . . .181
#
#     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/
############################################################################

# 13.1 A level chemistry dataset . . . . . . . . . . . . . . . . . . . . 181

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 alevchem data
data(alevchem, package = "R2MLwiN")

alevchem$gcseav <- alevchem$gcse_tot/alevchem$gcse_no - 6

hist(alevchem$gcseav, breaks = 20)

# 13.2 Normal response models . . . . . . . . . . . . . . . . . . . . . .184

(mymodel <- runMLwiN(a_point ~ 1 + (1 | pupil), estoptions = list(EstM = 1), data = alevchem))

(mymodel <- runMLwiN(a_point ~ 1 + gcseav + I(gcseav^2) + I(gcseav^3) + gender + (1 | pupil), estoptions = list(EstM = 1,
  resi.store = TRUE), data = alevchem))

predCurves(mymodel, xname = "gcseav", group = "genderfemale")

# 13.3 Ordered multinomial modelling . . . . . . . . . . . . . . . . . . 186

##IGLS
(mymodel <- runMLwiN(logit(a_point, cons, 6) ~ 1, D = "Ordered Multinomial", data = alevchem))

##MCMC
(mymodel <- runMLwiN(logit(a_point, cons, 6) ~ 1, D = "Ordered Multinomial", estoptions = list(EstM = 1), data = alevchem))

# 13.4 Adding predictor variables . . . . . . . . . . . . . . . . . . . .191
##MCMC
(mymodel <- runMLwiN(logit(a_point, cons, 6) ~ 1 + gcseav[1:5] + I(gcseav^2)[1:5] + I(gcseav^3)[1:5] + gender[1:5],
  D = "Ordered Multinomial", estoptions = list(EstM = 1), data = alevchem))

# 13.5 Multilevel ordered response modelling . . . . . . . . . . . . . . 192

# Note: Establishment codes on their own do not uniquely identify schools.
# Schools are instead uniquely identified by LEA code, establishment ID
# combination. Thus, here we generated a unique school ID.

alevchem$school <- as.numeric(factor(paste0(alevchem$lea, alevchem$estab)))

##MCMC
(mymodel <- runMLwiN(logit(a_point, cons, 6) ~ 1 + gcseav[1:5] + I(gcseav^2)[1:5] + gender[1:5] + (1[1:5] | school),
  D = "Ordered Multinomial", estoptions = list(EstM = 1), data = alevchem))

##MCMC
(mymodel <- runMLwiN(logit(a_point, cons, 6) ~ 1 + gcseav[1:5] + I(gcseav^2)[1:5] + gender[1:5] + (1[1:5] + gcseav[1:5] | school),
 D = "Ordered Multinomial", estoptions = list(EstM = 1), data = alevchem))
sixway(mymodel@chains[, "RP2_var_Intercept_12345", drop = FALSE], acf.maxlag = 300, "sigma2v6")

##Increases iterations to 50,000
(mymodel <- runMLwiN(logit(a_point, cons, 6) ~ 1 + gcseav[1:5] + I(gcseav^2)[1:5] + gender[1:5] + (1[1:5] + gcseav[1:5] | school),
 D = "Ordered Multinomial", estoptions = list(EstM = 1, mcmcMeth = list(iterations = 50000)), data = alevchem))
sixway(mymodel@chains[, "RP2_var_Intercept_12345", drop = FALSE], acf.maxlag = 300, "sigma2v6")

# Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . .128





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