12.4 Multilevel posterior predictions

Posterior predictions from multilevel models must account for group-level uncertainty. There are three types of predictions:

Existing groups: Use the posterior of the group's specific effect.
New groups: Sample a new group effect from the population distribution.
Population average: Integrate over the population distribution.

12.4.1 Prediction types

$\tilde{y}_j \sim \int p(\tilde{y}_j \mid \theta_j) \, p(\theta_j \mid D) \, d\theta_j \quad \text{(existing group)}$

$\tilde{y}_{\text{new}} \sim \int p(\tilde{y} \mid \theta) \, p(\theta \mid \alpha, \sigma_\gamma) \, p(\alpha, \sigma_\gamma \mid D) \, d\theta \, d\alpha \, d\sigma_\gamma \quad \text{(new group)}$

# Multilevel posterior predictions: three types
set.seed(42)
n_tanks <- 20
density <- sample(c(10, 25, 35), n_tanks, replace = TRUE)

# Simulate data

a_pop <- 0.5
sigma_tank <- 1.5
gamma_tank <- rnorm(n_tanks, 0, sigma_tank)
p_tank <- 1 / (1 + exp(-(a_pop + gamma_tank)))
survival <- rbinom(n_tanks, density, p_tank)
d <- data.frame(tank = factor(1:n_tanks), density = density, survival = survival)

# Fit multilevel model

library(lme4)
m <- glmer(cbind(survival, density - survival) ~ (1 | tank),
family = binomial, data = d)

# Extract posterior-like samples via simulation

nsim <- 1000
sim <- simulate(m, nsim = nsim)

# Type 1: Predictions for existing tanks

p_existing <- plogis(fixef(m)[1] + ranef(m)$tank[, 1])

# Type 2: Predictions for a NEW tank

# Sample new tank effect from Normal(0, sigma_tank)

sigma_est <- attr(VarCorr(m)$tank, "stddev")
new_gamma <- rnorm(nsim, 0, sigma_est)
p_new <- plogis(fixef(m)[1] + new_gamma)

# Type 3: Population average

# Integrate over the population distribution

pop_gamma <- rnorm(nsim, 0, sigma_est)
p_pop <- plogis(fixef(m)[1] + pop_gamma)

cat("Three types of predictions:\n")
cat(" Type 1 - Existing tanks (first 5):\n")
for (i in 1:5) {
cat(sprintf(" Tank %d: p = %.3f (observed: %.3f)\n",
i, p_existing[i], survival[i] / density[i]))
}

cat("\n Type 2 - New tank:\n")
cat(" Mean p:", round(mean(p_new), 3), "\n")
cat(" 95% interval: [", round(quantile(p_new, 0.025), 3),
",", round(quantile(p_new, 0.975), 3), "]\n")

cat("\n Type 3 - Population average:\n")
cat(" Mean p:", round(mean(p_pop), 3), "\n")
cat(" 95% interval: [", round(quantile(p_pop, 0.025), 3),
",", round(quantile(p_pop, 0.975), 3), "]\n")

cat("\n The new-tank prediction has MORE uncertainty\n")
cat(" because it includes the between-tank variation.\n")

Click Run to execute this code in your browser with WebR

12.4.2 Posterior predictive checks for multilevel models

# Posterior predictive check for multilevel model
set.seed(42)
n_tanks <- 20
density <- sample(c(10, 25, 35), n_tanks, replace = TRUE)
a_pop <- 0.5
sigma_tank <- 1.5
gamma_tank <- rnorm(n_tanks, 0, sigma_tank)
p_tank <- 1 / (1 + exp(-(a_pop + gamma_tank)))
survival <- rbinom(n_tanks, density, p_tank)
d <- data.frame(tank = factor(1:n_tanks), density = density, survival = survival)

library(lme4)
m <- glmer(cbind(survival, density - survival) ~ (1 | tank),
family = binomial, data = d)

# Simulate posterior predictive data

nsim <- 500
sim_survival <- matrix(NA, nrow = nsim, ncol = n_tanks)
for (s in 1:nsim) {

# Sample new tank effects

gamma_sim <- rnorm(n_tanks, 0, attr(VarCorr(m)$tank, "stddev"))
p_sim <- plogis(fixef(m)[1] + gamma_sim)
sim_survival[s, ] <- rbinom(n_tanks, density, p_sim)
}

# Compare observed vs simulated

obs_prop <- survival / density
sim_prop <- sim_survival / matrix(density, nrow = nsim, ncol = n_tanks, byrow = TRUE)

cat("Posterior predictive check:\n")
cat(" Observed mean proportion:", round(mean(obs_prop), 3), "\n")
cat(" Simulated mean proportion:", round(mean(sim_prop), 3), "\n")
cat(" Observed SD across tanks:", round(sd(obs_prop), 3), "\n")
cat(" Simulated SD across tanks:", round(mean(apply(sim_prop, 1, sd)), 3), "\n")

# Plot

boxplot(sim_prop, col = adjustcolor("steelblue", 0.3),
main = "Posterior predictive check",
xlab = "Tank", ylab = "Survival proportion")
points(1:n_tanks, obs_prop, pch = 16, col = "firebrick", cex = 1.2)
legend("topright", legend = c("Observed", "Simulated"),
col = c("firebrick", "steelblue"), pch = c(16, 15))

Click Run to execute this code in your browser with WebR

Rethinking

New groups have more uncertainty. When predicting for an existing group, you use that group's estimated effect, which is informed by its data. When predicting for a new group, you must sample from the population distribution — you do not know the new group's effect, so the prediction includes the full between-group variation. This is why the new-group prediction interval is always wider than the existing-group interval. This distinction is crucial for scientific communication: "we predict for this specific tank" is very different from "we predict for a new, unseen tank."

Summary of Section 12.4

Existing group: Use the group's posterior — narrower uncertainty.
New group: Sample from the population distribution — wider uncertainty.
Population average: Integrate over all group effects.
Posterior predictive checks verify the model reproduces the observed variation.
The distinction between prediction types matters for scientific communication.

More than one type of cluster

Summary