#These solutions are for the functions in Prac3.txt #For MCMCglmm, use the solutions in Prac3solb.txt hip<-read.table("hip.txt",col.names=c("y","age","sex","subj","time")) library(MASS) hip.gibbs1<-lmm.gibbs(its=10000) 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") 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="") 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") 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) 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")