########## R script: seizure ##########

# For Bayesian Poisson mixed model 
# analysis of the seizure data.

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

# Read in seizure data from file.

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

# Standardise the continuous predictor szcntbas and age.

sszcntbas <- as.vector(scale(seizure$szcntbas))
sage <- as.vector(scale(seizure$age))

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

allData <- list(nObs=236,nSubj=59,szcnt=seizure$szcnt,progabide=seizure$progabide,
                sszcntbas=sszcntbas,sage=sage,idnum=seizure$idnum)

# Initialise parameters.

parInits <- list(list(beta0=0,beta1=0,beta2=0,beta3=0,tauU=1,U=rep(0,59)))

# Obtain MCMC samples for each parameter.

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

# Display the results.

if (createPDF) pdf("seizure0.pdf")
summariseMCMC("beta0","intercept")
if (createPDF) pdf("seizure1.pdf")
summariseMCMC("beta1","coef. for progabide")
if (createPDF) pdf("seizure2.pdf")
summariseMCMC("beta2","coef. for sszcntbas")
if (createPDF) pdf("seizure3.pdf")
summariseMCMC("beta3","coef. for sage")
if (createPDF) pdf("seizure4.pdf")
summariseMCMC("sigU","random intercept std. dev.")

########## End of seizureBayes ##########