############################################################################################################################
# In the following examples, we use the BRugs, an open-source version of the BUGS embedded in R, to illustrate the Bayesian  
# mediation analysis (Yuan and MacKinnon, 2009). You need to install the package BRugs from the menu "Packages" in R before using it. 
# Alternatively, the stand-alone software WinBUGS can also be used to conduct the Bayesian mediation analysis. The WinBUGS has 
# a better user interface and detailed tutorials. It can be downloaded from http://www.mrc-bsu.cam.ac.uk/bugs/.
############################################################################################################################
          

# load BRRugs
library(BRugs) 
set.seed(6)

## Now setting the working directory to where data and BUGS code are located
DIR = c("Z:\\paper\\Mediation\\Bayesian\\data_example\\Release\\")
setwd(DIR)



######################### Bayesian single-level mediation analysis #########################


# read in the firefighters data provided by Elliot, Goldberg, Kuehl, Moe, Berger and Pickring (2007)
# the independent variable x represents the randomized exposure to the intervention 
# the mediator m is the change from baseline (m1) to followup (m2) in knowledge of the benefits of eating fruits and vegetables
# the dependent variable y is the change of the reported eating of fruits and vegetables from baseline (y1) to followup (y2)
firefighter = read.table("firefighters.txt",  col.name=c("x", "m1", "m2", "y1", "y2"))

# generate varaibles for mediation analysis
m = firefighter$m2-firefighter$m1
y = firefighter$y2-firefighter$y1
x = firefighter$x

# center covariates to improve the convergence of the Markov Chain Monte Carlo (MCMC)
x.c = x-mean(x)
m.c = m-mean(m)
  
# output data in the BUGS format
bugsData(list(y=y, m=m.c, x=x.c, N=length(y)), "simplemed_data.txt")
 
# initialize the parameters for the MCMC and ouput the initial values in the BUGS format
bugsData(list(beta3=0, beta=0, tau.prime=0, prec.y=1, beta2=0, alpha=0, prec.m=1), "simplemed_init.txt")  

# Using the BRugs package in R to do MCMC. Alternatively, the WinBUGS can be used without typing commands
modelCheck("simplemed.txt") # check model file
modelData("simplemed_data.txt") # read data file
modelCompile(numChains=1) # compile model with 1 chain
modelInits("simplemed_init.txt") # read init data file

# burn in 1000 iterations
modelUpdate(1000) 

# to record posterior samples of parameters of interest
samplesSet(c("alpha", "beta", "theta"))    # theta=alpha*beta

# 10000 more iterations to obtain posteriors 
modelUpdate(10000) # 10000 more iterations ....
  
# get summary statistics for the parameters of interest
samplesStats(c("alpha", "beta", "theta"), beg=1001, end=11000)





######################### Bayesian multilevel mediation analysis #########################


# read in the multilevel data provided by Kenny, Korchmaros and Bolger (2003)
# subjid is the subject id, meas is the measure time
# x is the independent variable, m is the mediating variable, and y is dependent variable
kenny <- read.table("kennydata.txt", col.names=c("subjid","meas", "x", "m","y", "const"))

# mediating variable
m = matrix(kenny$m, ncol=10, byrow=T)
  
# independent variable
x = matrix(kenny$x, ncol=10, byrow=T)

# dependent variable
y = matrix(kenny$y, ncol=10, byrow=T)

# the number of observations for the first- and second-level model
N2 = dim(y)[1]
N1 = dim(y)[2]

# output the data in the BUGS format
bugsData(list(y=y, m=m, x=x, N1=N1, N2=N2), "mlmed_data.txt")

# initialize the parameters for Markov Chain Monte Carlo (MCMC)
AlphaBeta = matrix(rnorm(2*N2), ncol=2)
init = list( Beta2=rnorm(N2), Beta3=rnorm(N2), Tau.p=rnorm(N2), AlphaBeta=AlphaBeta, beta2=1, 
beta3=1, tau.p=1, alphabeta=c(1, 1), prec.y=1, prec.m=1, sigma.beta2=1, sigma.beta3=1, sigma.taup=1, 
prec.ab=matrix(rep(10, 4), ncol=2)) 
 
# output the initial values in the BUGS format
bugsData(init, "mlmed_init.txt") 

# Using the BRugs package in R to do MCMC. Alternatively, the WinBUGS can be used without typing commands
modelCheck("mlmed.txt")        # check model file
modelData("mlmed_data.txt")    # read data file
modelCompile(numChains=1)      # compile model with 1 chain
modelInits("mlmed_init.txt")   # read init data file

# burn in with 1000 iterations
modelUpdate(1000) 
# list of parameters of interest
monitvar = c("ab", "c", "r", "alphabeta", "tau.p", "var.ab", "sigma2.taup")
# record posterior samples of the parameters
samplesSet(monitvar)
# 10000 more iterations to obtain posteriors 
modelUpdate(10000) 

# get summary statistics for the parameters of interest
samplesStats(monitvar, beg=1001, end=11000)