#APTS Statistical Modelling 2017
#Practical 2 - Bayesian estimation of an LMM

#Using MCMCglmm

hip<-read.table("hip.txt",col.names=c("y","age","sex","subj","time"))
library("MCMCglmm")

#(a)
hip.gibbs1 <- MCMCglmm(y~age+sex+factor(time), random=~subj,data=hip, nitt=10000,
burn=0, thin=1, prior=list(R=list(V=1,nu=0.002), G=list(G1=list(V=1000,nu=0.002))),pr=T)

#(b) trace plots
par(mfrow=c(2,4))
bn<-c("intercept","age","sex","time 2","time 3")
sn<-c("random effect variance","error variance")
for(i in 1:5) plot(1:10000,hip.gibbs1$Sol[,i],type="l",ylab=bn[i],xlab="iteration")
for(i in 1:2) plot(1:10000,hip.gibbs1$VCV[,i],type="l",ylab=sn[i],xlab="iteration")

#with burn in removed
s<-1001:10000
for(i in 1:5) plot(s,hip.gibbs1$Sol[s,i],type="l",ylab=bn[i],xlab="iteration")
for(i in 1:2) plot(s,hip.gibbs1$VCV[s,i],type="l",ylab=sn[i],xlab="iteration")

#(c) posterior summaries
apply(hip.gibbs1$Sol[s,],2,mean)
apply(hip.gibbs1$Sol[s,],2,sd)
apply(hip.gibbs1$VCV[s,],2,mean)
apply(hip.gibbs1$VCV[s,],2,sd)
for(i in 1:5) plot(density(hip.gibbs1$Sol[s,i]),xlab=bn[i],main="")
for(i in 1:2) plot(density(hip.gibbs1$VCV[s,i]),xlab=sn[i],main="")

#(d) prediction
pred<-function(x=c(63,0,3),subj=8)
{
xx<-rep(0,30)
if (subj>0) xx[subj]<-1 #only set subj to 0 for a new subject' o/w, pick out the subject
xx<-c(1,x[1:2],x[3]==2,x[3]==3,xx)
mu<-hip.gibbs1$Sol[s,]%*%xx
if (subj==0) mu<-mu+rnorm(9000,0,sqrt(hip.gibbs1$VCV[s,1]))
y<-rnorm(9000,mu,sqrt(hip.gibbs1$VCV[s,2]))
}
par(mfrow=c(1,3))
p1<-pred(c(63,0,3),8)
mean(p1)
sd(p1)
plot(density(p1),xlab="haematocrit",main="Predictive density for subject 8, time 3")
p2<-pred(c(52,1,1),15)
mean(p2)
sd(p2)
plot(density(p2),xlab="haematocrit",main="Predictive density for subject 15, time 1")
p3<-pred(c(52,1,1),0)
mean(p3)
sd(p3)
plot(density(p3),xlab="haematocrit",main="Predictive density for new female subject, age 52, time 1")

#(e)
hip.gibbs2<-MCMCglmm(y~age+sex+factor(time), random=~subj,data=hip, nitt=10000,
burn=0, thin=1, prior=list(R=list(V=1,nu=0.002), G=list(G1=list(V=1,nu=0.002))),pr=T)
par(mfrow=c(2,4))
for(i in 1:5) plot(1:10000,hip.gibbs2$Sol[,i],type="l",ylab=bn[i],xlab="iteration")
for(i in 1:2) plot(1:10000,hip.gibbs2$VCV[,i],type="l",ylab=sn[i],xlab="iteration")


#(f) partial autocorrelation
par(mfrow=c(2,4))
for(i in 1:6) acf(hip.gibbs1$Sol[s,i],type="partial")
for(i in 1:2) acf(hip.gibbs1$VCV[s,i],type="partial")
#simulation error
apply(hip.gibbs1$Sol[s,1:5],2,mean)
beta.ac<-NULL
for(i in 1:5) beta.ac[i]<-acf(hip.gibbs1$Sol[s,i],type="partial",main=bn[i],plot=F)$acf[1]
sqrt((1+beta.ac)/(1-beta.ac))
sqrt((1+beta.ac)/(1-beta.ac))*apply(hip.gibbs1$Sol[s,1:5],2,sd)/sqrt(9000)

