########## R script: pigwtsBayes ##########

# For Bayesian/MCMC analysis of the 
# pigweights data.

library(BRugs)
source("summariseMCMC.r")
createPDF <- FALSE

# Read in pigweights data from file.

pigweights <- read.table("pigweights.txt",header=TRUE)

# Standardise the continuous predictor weeks.

sweeks <- as.vector(scale(pigweights$weeks))

# Store data in a list in format required by BRugsFit.

allData <- list(nObs=432,nSubj=48,sweeks=sweeks,weight=pigweights$weight,
                idnum=pigweights$idnum)

# Initialise parameters.

parInits <- list(list(beta0=0,beta1=0,tauU=1,tauEps=1,U=rep(0,48)))

# Obtain MCMC samples for each parameter.

BRugsFit(data=allData,inits=parInits,parametersToSave=c("beta0","beta1","sigU","sigEps"),
         nBurnin=5000,nIter=5000,nThin=5,modelFile="pigwtsBayesModel.txt",numChains=1)

# Display the results.

if (createPDF) pdf("pigwtBayes0.pdf")
summariseMCMC("beta0","intercept")
if (createPDF) pdf("pigwtBayes1.pdf")
summariseMCMC("beta1","coef. for sweeks")
if (createPDF) pdf("pigwtBayes2.pdf")
summariseMCMC("sigU","random intercept std. dev.")
if (createPDF) pdf("pigwtBayes3.pdf")
summariseMCMC("sigEps","error std. dev.")

########## End of pigwtsBayes ##########