Nothing
############################################################################
# MLwiN MCMC Manual
#
# 6 Random Slopes Regression Models . . . . . . . . . . . . . . . . . . 71
#
# 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)
## Read tutorial data
data(tutorial, package = "R2MLwiN")
## Choose MCMC algoritm for estimation (IGLS will be used to obtain starting values for MCMC)
(mymodel <- runMLwiN(normexam ~ 1 + standlrt + school + school:standlrt + (1 | student), estoptions = list(EstM = 1),
data = tutorial))
## Define the model Choose IGLS algoritm for estimation Fit the model
(mymodel0a <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student), data = tutorial))
## Choose MCMC algoritm for estimation (IGLS will be used to obtain starting values for MCMC)
(mymodel0 <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student), estoptions = list(EstM = 1),
data = tutorial))
# 6.1 Prediction intervals for a random slopes regression model . . . . . 75
## Save level 2 residual chains
(mymodel <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student), estoptions = list(EstM = 1,
mcmcMeth = list(iterations = 5001), resi.store.levs = 2), data = tutorial))
predLines(mymodel, xname = "standlrt", lev = 2, selected = NULL, probs = c(0.025, 0.975), legend.space = "right",
legend.ncol = 2)
dev.new()
predLines(mymodel, xname = "standlrt", lev = 2, selected = c(30, 44, 53, 59), probs = c(0.025, 0.975))
# 6.2 Alternative priors for variance matrices . . . . . . . . . . . . . .78
# 6.3 WinBUGS priors (Prior 2) . . . . . . . . . . . . . . . . . . . . . .78
## Change the starting values for Level 2 variance matrix to .1 on diagonal 0 otherwise.
RP.b <- c(0.1, 0, 0.1, 0.554)
names(RP.b) <- c("RP2_var_Intercept", "RP2_cov_Intercept_standlrt", "RP2_var_standlrt", "RP1_var_Intercept")
## Fit the model
(mymodel1 <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student), estoptions = list(EstM = 1,
startval = list(RP.b = RP.b)), data = tutorial))
# 6.4 Uniform prior . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
## Diffuse priors (Uniform priors)
(mymodel2 <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student), estoptions = list(EstM = 1,
mcmcMeth = list(priorcode = 0)), data = tutorial))
# 6.5 Informative prior . . . . . . . . . . . . . . . . . . . . . . . . . 80
## Informative normal prior for Sigma_u
(mymodel3 <- runMLwiN(normexam ~ 1 + standlrt + (1 + standlrt | school) + (1 | student), estoptions = list(EstM = 1,
mcmcMeth = list(priorParam = list(rp2 = list(estimate = matrix(c(0.09, 0.018, 0.018, 0.015), 2, 2), size = 65)))),
data = tutorial))
# 6.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
if (!require(texreg)) {
warning("texreg package required to use screenreg() function")
} else {
screenreg(list(mymodel0a, mymodel0, mymodel1, mymodel2, mymodel3), custom.model.names=c("IGLS", "default", "prior 2", "uniform", "prior 4"), groups = list("Fixed Part" = 1:2, "Level-2" = 3:5, "Level-1" = 6:6),
stars = numeric(0), include.nobs=FALSE, include.loglik=FALSE, include.deviance=FALSE, include.dbar=FALSE, include.dthetabar=FALSE, include.pd=FALSE, include.dic=FALSE)
}
# Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . . 81
############################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.