#APTS Statistical Modelling 2013 #Practical 3 - Bayesian LMMs #These solutions are for MCMCglmm #For the functions in Prac3.txt, use the solutions in Prac3sola.txt hip<-read.table("hip.txt",col.names=c("y","age","sex","subj","time")) library("MCMCglmm") 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) 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:6) 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") s<-1001:10000 for(i in 1:6) 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") 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:6) 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="") pred<-function(x=c(63,0,3),subj=8) { xx<-rep(0,30) if (subj>0) xx[subj]<-1 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") 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") 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) 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:6) 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")