# Bernd Beber # Example code for MCMC implementations for multilevel data ############################################################ # Gibbs sampler for normal model # generalized from Gelman et al., Bayesian Data Analysis, 2nd edition, pp. 601-602 # generate some data for this example set.seed(38894) # normally distributed data, with group-specific means that are also distributed normally # number of groups j J <- 100 # group means, with mu=10 and tau=3 theta.j <- rnorm(J, 10, 3) # number of observations in each group n.j <- round(1 + abs(rnorm(J, 5, 50))) # total number of observations n <- sum(n.j) # group indicator group <- rep(1:100, n.j) # generate data, with sigma=5 y <- rnorm(n, rep(theta.j, n.j), 5) # suppose we observe data y and group indicators # number of observations n <- length(y) # number of observations in each group n.j <- as.vector(table(group)) # number of groups J <- length(n.j) # empirical group means y.j <- as.vector(by(y, group, mean)) # update functions theta.update <- function() { theta.hat <- (mu / tau^2 + y.j / (sigma^2 / n.j)) / (1 / tau^2 + 1 / (sigma^2 / n.j)) V.theta <- 1 / (1 / tau^2 + 1 / (sigma^2 / n.j)) rnorm (J, theta.hat, sqrt(V.theta)) } sigma.update <- function() { sqrt(sum((y - rep(theta, n.j))^2) / rchisq(1, n-1)) } mu.update <- function() { rnorm (1, mean(theta), tau / sqrt(J)) } tau.update <- function() { sqrt(sum((theta - mu)^2) / rchisq(1, J-1)) } # set number of iterations and chains n.chains <- 3 n.iter <- 150 # set up array to save output sims <- array(NA, c(n.iter, n.chains, J+3)) dimnames (sims) <- list(NULL, NULL, c(paste ("theta[", 1:J, "]", sep=""), "sigma", "mu", "tau")) # run the sampler for (m in 1:n.chains) { # initialize sigma <- runif (1, 0, sd(y)) mu <- rnorm(1, mean(y.j), sd(y.j)) tau <- runif(1, 0, sd(y.j)) for (t in 1:n.iter) { theta <- theta.update() sigma <- sigma.update() mu <- mu.update() tau <- tau.update() sims[t, m, ] <- c(theta, sigma, mu, tau) } } # plot chains to examine convergence names <- c(paste ("theta[", 1:J, "]", sep=""), "sigma", "mu", "tau") chain.plot <- function(parameter, xlab="Iterations", ylab=paste("Updated", names[parameter]), main=paste("Time series plot", names[parameter]), ...) { plot(c(1:t), sims[, 1, parameter], type="l", xlab=xlab, ylab=ylab, main=main, ...) lines(c(1:t), sims[, 2, parameter], lty=2) lines(c(1:t), sims[, 3, parameter], lty=3) } # for example, tau chain.plot(103) # show for all estimated parameters (i.e. theta of length J, sigma, mu, tau) n.par <- length(names) rows <- sqrt(n.par) par(mfrow=c(ceiling(rows), round(rows))) par(mar=rep(0, 4)) for (i in 1:n.par) chain.plot(i, xlab="", ylab="", main="", xaxt="n", yaxt="n") # discard first half of iterations out <- sims[round(n.iter/2):n.iter, , ] # get estimates mean <- apply(out, 3, mean) sd <- apply(out, 3, function(x) sd(c(x))) # show "shrinkage", i.e. posterior inferences are less dispersed than sample averages theta.hat <- mean[1:J] par(mfrow=c(1, 1)) par(mar=c(5, 4, 4, 2) + 0.1) plot(c(theta.hat, y.j), c(rep(1, J), rep(0, J)), main="comparing posterior inferences and sample averages", yaxt="n", ylab="", xlab="value of estimator") for (i in 1:J) { lines(c(theta.hat[i], y.j[i]), c(1,0)) } axis(2, c(0,1), c("sample average", "posterior inference")) # show that groups with large samples tend to have less shrinkage than groups with small samples plot(n.j, abs(theta.hat - y.j), ylab="distance between group's theta and its sample average", xlab="sample size of group", main="shrinkage for groups of different sample size")