########## R script: bpdBayesFanc ##########

# For Bayesian/MCMC logistic regression 
# analysis of the BPD data. In this `fancy'
# version higher-quality starting values are
# used, there is a longer burn-in perdiod,
# thinning is employed, the continuous predictor 
# is scaled for MCMC, and the results are then 
# converted to the original scale.

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

# Read in BPD data from file.

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

# Standardise the continuous predictor birweight.

meanBirthweight <- mean(bpd$birthweight)
sdBirthweight <- sqrt(var(bpd$birthweight))
sbirthweight <- (bpd$birthweight-meanBirthweight)/sdBirthweight

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

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

# Initialise parameters.

glmFit <- glm(BPD~sbirthweight,family=binomial,data=allData)
parInits <- list(list(sbeta0=glmFit$coef[1],sbeta1=glmFit$coef[2]))

# Obtain MCMC samples for each parameter.

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

# Display the results.

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

# Check against results from maximum likelihood fit.

MLfit <- glm(BPD~birthweight,family=binomial,data=bpd)
print(summary(MLfit))

########## End of bpdBayesFanc ##########