########## R script: pigwtsBayesRS ##########

# For Bayesian/MCMC analysis of the 
# pigweights data, with a random slope
# included.

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,tauEps=1,
                      UV=matrix(0,48,2),Omega=diag(2)))

# Obtain MCMC samples for each parameter.

BRugsFit(data=allData,inits=parInits,
         parametersToSave=c("beta0","beta1","sigU","sigV","sigEps"),
         nBurnin=5000,nIter=5000,nThin=5,modelFile="pigwtsBayesModelRS.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("sigV","random slope std. dev.")
if (createPDF) pdf("pigwtBayes4.pdf")
summariseMCMC("sigEps","error std. dev.")

# Check against lme() answers.

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

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