#APTS Statistical Modelling 2017
#Practical 2 - Bayesian LMMs

#For MCMCglmm, use the solutions in practical2_PartB2_sols.R

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

#(a)
hip.gibbs1<-lmm.gibbs(its=10000)

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

s<-1001:10000
for(i in 1:5) plot(s,hip.gibbs1$beta[s,i],type="l",ylab=bn[i],xlab="iteration")
for(i in 1:2) plot(s,hip.gibbs1$sigma[s,i],type="l",ylab=sn[i],xlab="iteration")
plot(s,hip.gibbs1$b[s,1],type="l",,ylab="b(1)",xlab="iteration")

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

#(d)
pred<-function(x=c(63,0,3),subj=8)
{
xx<-c(1,x[1:2],x[3]==2,x[3]==3)
mu<-hip.gibbs1$beta[s,]%*%xx
if (subj>0) mu<-mu+hip.gibbs1$b[s,subj]
if (subj==0) mu<-mu+rnorm(9000,0,sqrt(hip.gibbs1$sigma[s,2]))
y<-rnorm(9000,mu,sqrt(hip.gibbs1$sigma[s,1]))
}
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<-lmm.gibbs(its=10000,precb.sc=0.001)
for(i in 1:5) plot(1:10000,hip.gibbs2$beta[,i],type="l",ylab=bn[i],xlab="iteration")
for(i in 1:2) plot(1:10000,hip.gibbs2$sigma[,i],type="l",ylab=sn[i],xlab="iteration")
plot(1:10000,hip.gibbs2$b[,1],type="l",,ylab="b(1)",xlab="iteration")

#(f)
par(mfrow=c(2,4))
for(i in 1:5) acf(hip.gibbs1$beta[s,i],type="partial",main=bn[i])
for(i in 1:2) acf(hip.gibbs1$sigma[s,i],type="partial",main=sn[i])
acf(hip.gibbs1$b[s,1],type="partial",main="b(1)")
apply(hip.gibbs1$beta[s,],2,mean)
beta.ac<-NULL
for(i in 1:5) beta.ac[i]<-acf(hip.gibbs2$beta[s,i],type="partial",main=bn[i])$acf[1]
sqrt((1+beta.ac)/(1-beta.ac))
sqrt((1+beta.ac)/(1-beta.ac))*apply(hip.gibbs2$beta[s,],2,sd)/sqrt(9000)

