############### # R code for 2-stage Bayesian GLM (Lindley and Smith, 1972) # See Example 4.1 in Carlin and Louis (2008) ############### X <-cbind(rep(1,n),lgage) mu <- c(0,0) R <- matrix(c(0.001,0,0,0.001), nrow=2) tau <- c(1,100) postfn <- function(Y, X, mu, R, tau){ p <- dim(X)[2] D <- solve( tau*t(X)%*%X + R) d <- tau*t(X)%*%Y + R%*%diag(rep(1,p))%*%mu postmean <- D%*%d postsd <- sqrt(diag(D)) return(list(mean= D%*%d, sd=sqrt(diag(D)) )) } post1 <- postfn(Y, X, mu, R, tau[1]) # posterior for tau=1 beta <- post1$mean[2] postsd <- post1$sd[2] CI1 <- c(beta+qnorm(0.025)*postsd , beta+qnorm(0.975)*postsd) post2 <- postfn(Y, X, mu, R, tau[2]) # posterior for tau=100 beta <- post1$mean[2] postsd <- post1$sd[2] CI2 <- c(beta+qnorm(0.025)*postsd , beta+qnorm(0.975)*postsd) theta <- seq(-0.2,0.6,length.out=100) dens1<-dnorm(theta, mean=post1$mean[2], sd=post1$sd[2] ) dens2<-dnorm(theta, mean=post2$mean[2], sd=post2$sd[2] ) plot(theta, dens1, xlab=expression(beta), xlim=c(-0.2, 0.6), ylim=c(0,20), ylab="posterior density", type="n") lines(theta, dens1, lty=1) # posterior density for tau=1 lines(theta, dens2, lty=2) # posterior density for tau=100 postdraw <- rnorm(2000, mean=post2$mean[2], sd=post2$sd[2] ) r1<-hist(postdraw,freq=F,breaks=20, plot=F) lines(r1,lty=3, freq=F, col="gray90") #posterior draws, tau=100 legend(-0.2, 20, legend=c(expression(paste(tau,"=1", sep="")), expression(paste(tau,"=100", sep=""))), lty=c(1,2), ncol=1)