############################################################################ # MLwiN MCMC Manual # # 7 Using the WinBUGS Interface in MLwiN . . . . . . . . . . . . . . . .83 # # 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/ ############################################################################ # 7.1 Variance components models in WinBUGS . . . . . . . . . . . . . . . 84 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") ## The highest level comes first, then the second highest and so on ## Uses the results from IGLS to create initial values for bugs ## Fit the model by calling openbugs using the R2WinBUGS package mymodel1 <- runMLwiN(normexam ~ 1 + standlrt + (1 | school) + (1 | student), estoptions = list(EstM = 1, show.file = TRUE), BUGO = c(version = 4, n.chains = 1, debug = FALSE, seed = 1, OpenBugs = TRUE), data = tutorial) summary(mymodel1) summary(mymodel1[, "beta[2]"]) sixway(mymodel1[, "beta[2]", drop = FALSE]) # 7.2 So why have a WinBUGS interface ? . . . . . . . . . . . . . . . . . 92 # 7.3 t distributed school residuals . . . . . . . . . . . . . . . . . . .92 ## Download the model, initial, data files modelfile <- paste0(tempdir(), "/tutorial1_model.txt") download.file("http://www.bristol.ac.uk/cmm/media/r2mlwin/tutorial1_model.txt", modelfile, method = "auto") file.show(modelfile) initfile <- paste0(tempdir(), "/tutorial1_inits.txt") download.file("http://www.bristol.ac.uk/cmm/media/r2mlwin/tutorial1_inits.txt", initfile, method = "auto") file.show(initfile) datafile <- paste0(tempdir(), "/tutorial1_data.txt") download.file("http://www.bristol.ac.uk/cmm/media/r2mlwin/tutorial1_data.txt", datafile, method = "auto") bugEst <- paste0(tempdir(), "/tutorial1_log.txt") chains.bugs1 <- mlwin2bugs(D = "t", levID = c("school", "student"), datafile, initfile, modelfile, bugEst, fact = NULL, addmore = c("sigma2", "df"), n.chains = 1, n.iter = 5500, n.burnin = 500, n.thin = 1, debug = TRUE, bugsWorkingDir = tempdir(), OpenBugs = TRUE) ## Close winbugs manually summary(chains.bugs1) sixway(chains.bugs1[, "df", drop = FALSE]) chains.bugs2 <- mlwin2bugs(D = "t", levID = c("school", "student"), datafile, initfile, modelfile, bugEst, fact = NULL, addmore = c("sigma2", "df"), n.chains = 1, n.iter = 12000, n.burnin = 2000, n.thin = 1, debug = TRUE, bugsWorkingDir = tempdir(), OpenBugs = TRUE) ## Close winbugs manually summary(chains.bugs2) sixway(chains.bugs2[, "df", drop = FALSE]) # Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . . 96 ############################################################################