########## R script: fossilBayes ##########

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

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

# Read in FOSSIL data from file.

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

# Standardise the continuous predictor age.

sage <- as.vector(scale(fossil$age))
sstrontium.ratio <- as.vector(scale(fossil$strontium.ratio))

# Obtain knots for age variable.

numKnots <- 35
knots <-  quantile(sage,seq(0,1,length=numKnots+2)[-c(1,numKnots+2)])
Z <- outer(sage,knots,"-")
Z <- Z*(Z>0)

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

allData <- list(n=length(sage),numKnots=numKnots,
                sage=sage,sstrontium.ratio=sstrontium.ratio,Z=Z)

# Initialise parameters.

parInits <- list(list(beta0=0,beta1=0,u=rep(0,numKnots),tauU=1,tauEps=1))

# Obtain MCMC samples for each parameter.

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

# Display the fit.

betaMCMC <- rbind(samplesSample("beta0"),samplesSample("beta1"))

uMCMC <- samplesSample("u[1]")
for (k in 2:numKnots)
   uMCMC <- rbind(uMCMC,samplesSample(paste("u[",k,"]")))

ng <- 201
xg <- seq(min(sage),max(sage),length=ng)
Xg <- cbind(rep(1,ng),xg)
Zg <- outer(xg,knots,"-")
Zg <- Zg*(Zg>0)

fhatMCMC <- Xg%*%betaMCMC + Zg%*%uMCMC

fhatg <- apply(fhatMCMC,1,mean)
quant025 <- function(x) return(quantile(x,0.025)) 
quant975 <- function(x) return(quantile(x,0.975)) 
credLower <- apply(fhatMCMC,1,quant025)
credUpper <- apply(fhatMCMC,1,quant975)

par(mfrow=c(1,1))
plot(sage,sstrontium.ratio,xlab="standardised age",
     ylab="standardised strontium ratio")
polygon(c(xg,rev(xg)),c(credLower,rev(credUpper)),col="yellow",border=FALSE)
lines(xg,fhatg,col="forestgreen",lwd=2)
points(sage,sstrontium.ratio,col="blue")

cat("Hit Enter to continue.\n")
ans <- readline()

# Display the posterior distribution of the effective degrees of freedom.

sigUMCMC <- samplesSample("sigU")
sigEpsMCMC <- samplesSample("sigEps")
nMCMC <- length(sigUMCMC)
dfMCMC <- rep(NA,nMCMC)

X <- cbind(rep(1,length(sage)),sage)
CTC <- t(cbind(X,Z))%*%cbind(X,Z)
Dmat <- diag(c(rep(0,2),rep(1,numKnots)))
for (i in 1:nMCMC)
{
   lambda <- (sigEpsMCMC[i]/sigUMCMC[i])^2
   dfMCMC[i] <- sum(diag(solve(CTC+lambda*Dmat,CTC)))
}

library(KernSmooth)
h <- dpik(dfMCMC)
dest <- bkde(dfMCMC,bandwidth=h)
plot(dest,type="l",col="magenta",lwd=3,xlab="posterior degrees of freedom")

########## End of fossilBayes ##########