R version 3.6.3 (2020-02-29) -- "Holding the Windsock"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ############################################################################
> #     MLwiN MCMC Manual
> #
> # 1   Introduction to MCMC Estimation and Bayesian Modelling . . . . . . . 1
> #
> #     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/
> ############################################################################
> 
> # 1.1 Bayesian modelling using Markov Chain Monte Carlo methods . . . . . .1
> 
> # 1.2 MCMC methods and Bayesian modelling . . . . . . . . . . . . . . . . .2
> 
> # 1.3 Default prior distributions . . . . . . . . . . . . . . . . . . . . .4
> 
> # 1.4 MCMC estimation . . . . . . . . . . . . . . . . . . . . . . . . . . .5
> 
> # 1.5 Gibbs sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
> 
> # 1.6 Metropolis Hastings sampling . . . . . . . . . . . . . . . . . . . . 8
> 
> # 1.7 Running macros to perform Gibbs sampling and Metropolis
> #     Hastings sampling on the simple linear regression model . . . . . . 10
> 
> library(R2MLwiN)

R2MLwiN: A package to run models implemented in MLwiN from R
Copyright 2013-2017 Zhengzheng Zhang, Christopher M. J. Charlton, Richard M. A. Parker, William J. Browne and George Leckie

Support provided by the Economic and Social Research Council (ESRC)
(Grants RES-149-25-1084, RES-576-25-0032 and ES/K007246/1)

To cite R2MLwiN in publications use:

  Zhengzheng Zhang, Richard M. A. Parker, Christopher M. J. Charlton,
  George Leckie, William J. Browne (2016). R2MLwiN: A Package to Run
  MLwiN from within R. Journal of Statistical Software, 72(10), 1-43.
  doi:10.18637/jss.v072.i10

A BibTeX entry for LaTeX users is

  @Article{,
    title = {{R2MLwiN}: A Package to Run {MLwiN} from within {R}},
    author = {Zhengzheng Zhang and Richard M. A. Parker and Christopher M. J. Charlton and George Leckie and William J. Browne},
    journal = {Journal of Statistical Software},
    year = {2016},
    volume = {72},
    number = {10},
    pages = {1--43},
    doi = {10.18637/jss.v072.i10},
  }

The MLwiN_path option is currently set to C:/Program Files/MLwiN v3.05/
To change this use: options(MLwiN_path="<path to MLwiN>")
> # 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")
> 
> set.seed(1)
> 
> ## Set variables
> y <- tutorial$normexam
> x <- tutorial$standlrt
> N <- length(y)
> xsq <- x^2
> xy <- x * y
> sumy <- sum(y)
> sumx <- sum(x)
> sumxsq <- sum(xsq)
> sumxy <- sum(xy)
> 
> ## Starting values for parameter estimates
> beta0 <- 0
> beta1 <- 0
> sigma2e <- 1
> epsilon <- 0.001
> burnin <- 500
> chain <- 5000
> totaliterations <- burnin + chain
> thinning <- 1
> estimates <- matrix(, nrow = floor(chain/thinning), 3)
> rownames(estimates) <- which((1:chain)%%thinning == 0)
> colnames(estimates) <- c("beta0", "beta1", "sigma2e")
> 
> j <- 1
> for (i in 1:totaliterations) {
+   beta0 <- rnorm(1, (sumy - beta1 * sumx)/N, sqrt(sigma2e/N))
+   beta1 <- rnorm(1, (sumxy - beta0 * sumx)/sumxsq, sqrt(sigma2e/sumxsq))
+   
+   e2i <- (y - (beta0 + beta1 * x))^2
+   sume2i <- sum(e2i)
+   
+   sigma2e <- 1/rgamma(1, epsilon + N/2, epsilon + sume2i/2)
+   
+   if ((i%%thinning == 0) && (i > burnin)) {
+     estimates[j, ] <- round(c(beta0, beta1, sigma2e), 3)
+     j <- j + 1
+   }
+ }
> 
> sumstat <- round(rbind(colMeans(estimates), apply(estimates, 2, sd)), 4)
> rownames(sumstat) <- c("mean", "sd")
> print(sumstat)
       beta0  beta1 sigma2e
mean -0.0012 0.5954  0.6492
sd    0.0128 0.0127  0.0144
> 
> # 1.8 Dynamic traces for MCMC . . . . . . . . . . . . . . . . . . . . . . 12
> 
> par(mfrow = c(3, 1), mar = c(4, 4.5, 2, 2))
> plot(1:nrow(estimates), estimates[, "beta0"], xlab = "iteration", ylab = expression(paste("Est. of ", beta[0])), type = "l")
> plot(1:nrow(estimates), estimates[, "beta1"], xlab = "iteration", ylab = expression(paste("Est. of ", beta[1])), type = "l")
> plot(1:nrow(estimates), estimates[, "sigma2e"], xlab = "iteration", ylab = expression(paste("Est. of ", sigma[e]^2)), 
+   type = "l")
> 
> # 1.9 Macro to run a hybrid Metropolis and Gibbs sampling method
> #     for a linear regression example . . . . . . . . . . . . . . . . . . 15
> set.seed(1)
> 
> ## Set variables
> y <- tutorial$normexam
> x <- tutorial$standlrt
> N <- length(y)
> xsq <- x^2
> xy <- x * y
> sumy <- sum(y)
> sumx <- sum(x)
> sumxsq <- sum(xsq)
> sumxy <- sum(xy)
> 
> ## Starting values for paramter estimates
> beta0 <- 0
> beta1 <- 0
> sigma2e <- 1
> epsilon <- 0.001
> burnin <- 500
> chain <- 5000
> beta0sd <- 0.01
> beta1sd <- 0.01
> beta0accept <- 0
> beta1accept <- 0
> totaliterations <- burnin + chain
> thinning <- 1
> estimates <- matrix(, nrow = floor(chain/thinning), 3)
> rownames(estimates) <- which((1:chain)%%thinning == 0)
> colnames(estimates) <- c("beta0", "beta1", "sigma2e")
> 
> j <- 1
> for (i in 1:totaliterations) {
+   # Update beta0 Propose a new beta0
+   beta0prop <- rnorm(1, beta0, beta0sd)
+   beta0logpostdiff <- -1 * (2 * (beta0 - beta0prop) * (sumy - beta1 * sumx) + N * (beta0prop^2 - beta0^2))/(2 * 
+     sigma2e)
+   if (beta0logpostdiff > 0) {
+     # Definitely accept as higher posterior
+     beta0 <- beta0prop
+     beta0accept <- beta0accept + 1
+   } else {
+     # Only sometimes accept
+     if (runif(1) < exp(beta0logpostdiff)) {
+       beta0 <- beta0prop
+       beta0accept <- beta0accept + 1
+     }
+   }
+   
+   # Update beta1
+   beta1prop <- rnorm(1, beta1, beta1sd)
+   beta1logpostdiff <- -1 * (2 * (beta1 - beta1prop) * (sumxy - beta0 * sumx) + sumxsq * (beta1prop^2 - beta1^2))/(2 * 
+     sigma2e)
+   if (beta1logpostdiff > 0) {
+     # Definitely accept as higher posterior
+     beta1 <- beta1prop
+     beta1accept <- beta1accept + 1
+   } else {
+     # Only sometimes accept
+     if (runif(1) < exp(beta1logpostdiff)) {
+       beta1 <- beta1prop
+       beta1accept <- beta1accept + 1
+     }
+   }
+   
+   e2i <- (y - (beta0 + beta1 * x))^2
+   sume2i <- sum(e2i)
+   
+   sigma2e <- 1/rgamma(1, epsilon + N/2, epsilon + sume2i/2)
+   
+   if ((i%%thinning == 0) && (i > burnin)) {
+     estimates[j, ] <- round(c(beta0, beta1, sigma2e), 3)
+     j <- j + 1
+   }
+ }
> 
> sumstat <- round(rbind(colMeans(estimates), apply(estimates, 2, sd)), 4)
> rownames(sumstat) <- c("mean", "sd")
> print(sumstat)
       beta0  beta1 sigma2e
mean -0.0009 0.5954  0.6489
sd    0.0123 0.0119  0.0143
> cat(paste("The acceptance rate of beta0: ", round(beta0accept/totaliterations, 3), "\n"))
The acceptance rate of beta0:  0.763 
> cat(paste("The acceptance rate of beta1: ", round(beta1accept/totaliterations, 3), "\n"))
The acceptance rate of beta1:  0.756 
> 
> # 1.10 MCMC estimation of multilevel models in MLwiN . . . . . . . . . . .18
> 
> # Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . . 19
> 
> ## reinstate par settings
> par(mypar)
> 
> 
> 
> ############################################################################
> 
> proc.time()
   user  system elapsed 
   3.57    0.25    4.04