########## R script: bpdBayes ##########

# For Bayesian/MCMC logistic regression 
# analysis of the BPD data.

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

# Read in BPD data from file.

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

# Standardise the continuous predictor birthweight.

sbirthweight <- as.vector(scale(bpd$birthweight))

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

allData <- list(n=223,sbirthweight=sbirthweight,BPD=bpd$BPD)

# Initialise parameters.

parInits <- list(list(beta0=0,beta1=0))

# Obtain MCMC samples for each parameter.

BRugsFit(data=allData,inits=parInits,parametersToSave=c("beta0","beta1"),
         nBurnin=100,nIter=1000,modelFile="bpdBayesModel.txt",numChains=1)

# Display the results.

if (createPDF) pdf("bpdBayes0.pdf")
summariseMCMC("beta0","intercept")
if (createPDF) pdf("bpdBayes1.pdf")
summariseMCMC("beta1","coef. for sbirthweight")

########## End of bpdBayes ##########