########## R script: pigwtsBayesHC ##########

# For Bayesian/MCMC nalysis of the 
# pigweights data. This version uses
# hierarchical centring.

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

# Read in pigweights data from file.

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

# Standardise the continuous predictor birweight.

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,gamma=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="pigwtsBayesModelHC.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.")

# Check against lme() answers.

library(nlme)
gD <- groupedData(weight~sweeks|idnum,data=pigweights)
fit <- lme(weight~sweeks,random=~1,data=gD)
print(intervals(fit))

########## End of pigwtsBayesHC ##########