After finishing the model specification, I went back to the R environment and stored data from the data frame, tab14.3, into several independent vectors.
# Load the data:
y = tab14.3$Score
x1 = tab14.3$Wk
S = tab14.3$Subj
x1names = unique(tab14.3$Wk)
Snames = unique(tab14.3$Subj)
Ntotal = length(y)
Nx1Lvl = length(unique(x1))
NSLvl = length(unique(S))
The above codes extract dependent variable (y) and independent variables (x1 and S) and store them into separate vectors. The label for each level of the two independent variables were assigned to x1names and Snames, respectively. For the BUGS code, i and j1, I also created the variables, Ntotal, Nx1Lv1 and NSLv1. Later, they were sent as a data list to BUGS.
ySDorig = sd(y)
yMorig = mean(y)
z = ( y - yMorig ) / ySDorig
dataList = list(
y = z ,
x1 = x1 ,
S = S ,
Ntotal = Ntotal ,
Nx1Lvl = Nx1Lvl ,
NSLvl = NSLvl
)
This step converted the predicted variable (number of headache attack) into z score to speed up the computation. The data were then stored into a data list, dataList, so that BUGS can process it. In total, six data vectors were sent to BUGS.
#------------------------------------------------------------------------------
# INTIALIZE THE CHAINS.
theData = data.frame( y=dataList$y , x1=factor(x1,labels=x1names) ,
S=factor(S,labels=Snames) )
a0 = mean( theData$y )
a1 = aggregate( theData$y , list( theData$x1 ) , mean )[,2] - a0
aS = aggregate( theData$y , list( theData$S ) , mean )[,2] - a0
initsList = list(
a0 = a0 ,
a1 = a1 ,
aS = aS ,
sigma = sd(theData$y)/2 , # there may be a better way to specify the sigma
a1SDunabs = sd(a1) , # to accelerate the computation
aSSDunabs = sd(aS) )
The standardised data were again stored into a separate data frame, theData. Then, the aggregate function calculated the deflection from the grand mean in the two independent variables (a1 and subject). The aggregate works very much like tapply. The will-be estimated parameters are then stored into a list, initsList. Therefore, the sample "standardised" grand mean was the initial a0. The initial a1s were the differences of the "standardised" treatment means from each week and the "standardised" grand mean as well as the subject factors. The standard deviation of each observation unit was initialised as half of the standard deviation from the sample data.
#------------------------------------------------------------------------------
# RUN THE CHAINS
parameters = c( "a0", "a1", "aS", "sigma", "a1SD", "aSSD" )
adaptSteps = 1000 # dont quite understand what adaptSteps function?
burnInSteps = 2000
nChains = 5
numSavedSteps=50000
thinSteps=1
nIter = ceiling( ( numSavedSteps * thinSteps ) / nChains )
# Create, initialise, and adapt the model:
jagsModel = jags.model( "modelonewayRM.txt" , data=dataList , inits=initsList ,
n.chains=nChains , n.adapt=adaptSteps )
# Burn-in: cat( "Burning in the MCMC chain... " )
update( jagsModel , n.iter=burnInSteps )
# The saved MCMC chain: cat( "Sampling final MCMC chain... " )
codaSamples = coda.samples( jagsModel , variable.names=parameters ,
n.iter=nIter , thin=thinSteps )
The first 7 lines just specify some necessary options that need to pass to the jags.model function, including (1) the parameters that you want to estimate, (2) the first few steps that you want to throw away because they may not have had appropriately converged, (3) how many chains you want to run, (4) how many step you wish the R to save for later examination, (5) thinning step to decrease possible autocorrelation, (6) lastly how many iteration you want.
Then the jags.model function loaded the model file, data list, initial (parameter) list, the chain number and the adapt step to start the simulation.
Note that in the burn-in step, the n.iter equals to the burnInSteps specified before. The burn-in and MCMC iteration are separated into two steps. The MCMC iteration is run in the coda.samples function. This finished the BUGS process and next step will start to examine the result. To be continued.