12.4 Predicciones posteriores multinivel

Las predicciones a posteriori de modelos multinivel deben contabilizar la incertidumbre de nivel de grupo. Hay tres tipos de predicciones:

Grupos existentes: Use la posterior del efecto específico del grupo.
Grupos nuevos: Muestree un nuevo efecto de grupo de la distribución poblacional.
Promedio poblacional: Integre sobre la distribución poblacional.

12.4.1 Tipos de predicción

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

$\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{(grupo nuevo)}$

# 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 Comprobaciones predictivas a posteriori para modelos multinivel

# 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

Los grupos nuevos tienen más incertidumbre. Al predecir para un grupo existente, se usa el efecto estimado de ese grupo, que está informado por sus datos. Al predecir para un grupo nuevo, se debe muestrear de la distribución poblacional — no se conoce el efecto del grupo nuevo, así que la predicción incluye la variación completa entre grupos. Por esto el intervalo de predicción para grupos nuevos siempre es más amplio que para grupos existentes. Esta distinción es crucial para la comunicación científica: "predecimos para este tanque específico" es muy diferente de "predecimos para un tanque nuevo, no observado".

Resumen de la Sección 12.4

Grupo existente: Use la posterior del grupo — incertidumbre más estrecha.
Grupo nuevo: Muestree de la distribución poblacional — incertidumbre más amplia.
Promedio poblacional: Integre sobre todos los efectos de grupo.
Las comprobaciones predictivas a posteriori verifican que el modelo reproduce la variación observada.
La distinción entre tipos de predicción importa para la comunicación científica.

Más de un tipo de cluster

Resumen