# 8 schools example of Section 5.5
# Hierarchical normal model with known variances
# Key in data
rm(list=ls()); set.seed(20210209);
y <- c(28, 8, -3, 7, -1, 1, 18, 12) J <- length(y) sigma <- c(15, 10, 16, 11, 9, 11, 10, 18) Sim.size <- 1000 # Separate (unpooled) estimation Ints <- y + qnorm(.975) * sigma %*% t(c(-1,1)) plot((1:J)-.05, y, xlim=c(.5,J+.5), ylim=range(Ints), pch=19, xlab="School", ylab="Effect", xaxp=c(1,J,J-1), col="blue") for(j in 1:J){ lines(c(j,j)-.05, Ints[j,], col="blue") } # Fully pooled estimation post.var <- 1 / sum(1/(sigma^2)) post.mean <- sum(y/(sigma^2)) * post.var Int <- post.mean + c(-1,1) * qnorm(.975) * sqrt(post.var) points((1:J)+.05, rep(post.mean,J), pch=19, col="red") for(j in 1:J){ lines(c(j,j)+.05, Int, col="red") } # Posterior distribution for hierarchical normal model log.post.tau.fun <- function(tau, y, sigma) # (5.21 on p. 117) { V.mu <- 1 / sum( 1/(sigma^2 + tau^2) ) mu.hat <- V.mu * sum( y / (sigma^2 + tau^2) ) log.post <- ( log(V.mu) - sum(log(sigma^2 + tau^2) ) ) / 2 log.post <- log.post - 0.5 * sum( (y-mu.hat)^2 / (sigma^2 + tau^2) ) log.post } log.tau <- seq(-3.4, 3.4, .01); T <- length(log.tau); T; log.post.tau <- rep(NA, T) for(t in 1:T) { log.post.tau[t] <- log.post.tau.fun(tau=exp(log.tau[t]), y, sigma) } maxie <- max(log.post.tau) log.post.tau <- log.post.tau - maxie; rm(maxie); post.tau <- exp(log.post.tau) # Figure 5.5 on page 121 plot(exp(log.tau), post.tau, type="l", xlab="tau", ylab="p(tau|y)") # Calculate E(theta|tau,y), i.e., the "estimated treatment effects" ETE <- matrix(NA, T, J) for(t in 1:T) { tau <- exp(log.tau[t]) V.mu <- 1 / sum( 1 / (sigma^2 + tau^2) ) mu.hat <- V.mu * sum( y / (sigma^2 + tau^2) ) ETE[t,] <- (y/sigma^2 + mu.hat/tau^2) / (1/sigma^2 + 1/tau^2) } # Figure 5.6 on page 122 matplot(exp(log.tau), ETE, type="l", xlab="tau", ylab="E(theta|tau,y)", ylim=range(ETE) + c(-5,5)) legend("topleft", inset=.05, legend=paste("School", 1:J, sep=" "), lty=1:5, col=1:6, cex=.9) # Calculate posterior SD of treatment effects, sd(theta,tau,y) VTE <- matrix(NA, T, J) for(t in 1:T) { tau <- exp(log.tau[t]) V <- 1 / (1/sigma^2 + 1/tau^2) V.mu <- 1 / sum( 1 / (sigma^2 + tau^2) ) VTE[t,] <- V + V^2 * V.mu / tau^4 } # Figure 5.7 on page 122 matplot(exp(log.tau), sqrt(VTE), type="l", xlab="tau", ylab="sd(theta|tau,y)", ylim=range(sqrt(VTE)) + c(-4,4)) legend("topleft", inset=.05, legend=paste("School", 1:J, sep=" "), lty=1:5, col=1:6, cex=.9) # Sampling from the posterior p(theta|y) S <- Sim.size theta.sim <- matrix(NA, J, S) ytilde.sim <- matrix(NA, J, S) delta <- (log.tau[2] - log.tau[1]) / 2 post.log.tau <- exp(log.tau) * post.tau for(s in 1:S) { t <- sample(T, 1, prob=post.log.tau) tau <- exp(log.tau[t] + runif(1, -delta, delta)) V.mu <- 1 / sum( 1 / (sigma^2 + tau^2) ) mu.hat <- V.mu * sum( y / (sigma^2 + tau^2) ) mu <- rnorm(1, mean=mu.hat, sd=sqrt(V.mu)) V <- 1 / (1/sigma^2 + 1/tau^2) theta.hat <- V * (y/sigma^2 + mu/tau^2) theta.sim[,s] <- rnorm(J, mean=theta.hat, sd=sqrt(V)) ytilde.sim[,s] <- rnorm(J, mean=theta.sim[,s], sd=sigma) } quant <- function(x){ quantile(x, probs=c(.025,.25,.50,.75,.975)) } Results <- t(apply(theta.sim, 1, quant)) rownames(Results) <- paste("School", 1:J, sep="_") round(Results) # Table 5.3 on page 123 # Two quantities of interest: theta_1 and max{theta_j} theta1 <- theta.sim[1,]; theta.max <- apply(theta.sim, 2, max); hist.range <- range(c(theta1,theta.max)) + c(-5,5) op <- par(mfrow=c(1,2)) hist(theta1, freq=F, right=F, breaks=30, xlim=hist.range) lines(density(theta1), col="red") hist(theta.max, freq=F, right=F, breaks=30, xlim=hist.range) lines(density(theta.max), col="red") par(op) # Figure 5.8 on page 123