1.3 Tres herramientas para la ingeniería de golems

Este libro enseña tres herramientas para construir mejores modelos estadísticos: inferencia bayesiana, modelos multinivel y criterios de información. Estas herramientas no son mutuamente excluyentes — se refuerzan mutuamente y forman los cimientos del enfoque del libro hacia el modelado estadístico.

Herramienta 1: Inferencia bayesiana

La inferencia bayesiana es un marco para actualizar creencias bajo incertidumbre. Dada una distribución a priori $P(\theta)$ y una verosimilitud $P(D \mid \theta)$, el teorema de Bayes nos dice cómo calcular la distribución a posteriori:

$\underbrace{p(\theta \mid D)}_{\text{posterior}} \propto \underbrace{p(D \mid \theta)}_{\text{likelihood}} \cdot \underbrace{p(\theta)}_{\text{prior}}$

La distribución a posteriori es la respuesta completa a la pregunta "¿qué sabemos sobre $\theta$ después de ver los datos?" Es una distribución, no una estimación puntual, por lo que captura naturalmente la incertidumbre.

La inferencia bayesiana tiene varias ventajas sobre los métodos frecuentistas:

• Incertidumbre completa: La distribución a posteriori nos da toda la distribución de valores plausibles, no solo una estimación puntual y un error estándar.
• Información a priori: Podemos incorporar formalmente lo que ya sabemos, en lugar de pretender ser ignorantes.
• Coherencia: Las reglas de probabilidad aseguran que nuestras conclusiones sean internamente consistentes.
• Flexibilidad: El mismo marco funciona para modelos simples y modelos jerárquicos complejos.

# Bayesian inference from scratch: grid approximation
# Estimate the proportion of water on a globe
# Data: 6 waters in 9 tosses

# 1. Define the parameter grid

p_grid <- seq(from = 0, to = 1, length.out = 100)

# 2. Define the prior: uniform (uninformative)

prior <- rep(1, 100)

# 3. Compute the likelihood: binomial

likelihood <- dbinom(6, size = 9, prob = p_grid)

# 4. Compute the posterior: multiply and normalize

posterior <- (likelihood _ prior) / sum(likelihood _ prior)

# 5. Summarize the posterior

cat("Posterior mean:", round(sum(p*grid * posterior), 3), "\n")
cat("Posterior std dev:", round(sqrt(sum((p*grid - sum(p_grid * posterior))^2 * posterior)), 3), "\n")
cat("95% credible interval: [",
round(p_grid[which(cumsum(posterior) >= 0.025)[1]], 3), ",",
round(p_grid[which(cumsum(posterior) >= 0.975)[1]], 3), "]\n")

# Plot the posterior

plot(p_grid, posterior, type = "l", lwd = 2, col = "steelblue",
xlab = "Proportion of water (p)",
ylab = "Posterior density",
main = "Bayesian posterior: 6 waters in 9 tosses")
polygon(c(p_grid, rev(p_grid)),
c(posterior, rep(0, length(posterior))),
col = rgb(0.27, 0.51, 0.71, 0.3), border = NA)
abline(v = 6/9, lty = 2, col = "firebrick")
text(0.72, max(posterior) * 0.9, "MLE = 6/9", col = "firebrick")

La distribución a posteriori captura todo lo que sabemos sobre $p$ después de ver los datos. La media, la desviación estándar y el intervalo de credibilidad son todos resúmenes de esta distribución. A diferencia de un valor $p$, la distribución a posteriori nos dice el rango completo de valores plausibles y sus probabilidades relativas.

Herramienta 2: Modelos multinivel

Los modelos multinivel (también llamados modelos jerárquicos o de efectos mixtos) manejan datos con estructura: estudiantes anidados en escuelas, mediciones anidadas en pacientes, ensayos anidados en sujetos. En lugar de asumir que todos los grupos son idénticos (agrupamiento completo) o completamente independientes (sin agrupamiento), los modelos multinivel usan agrupamiento parcial: cada grupo obtiene sus propios parámetros, pero los parámetros de grupo se extraen de una distribución común.

La idea clave es que la distribución poblacional vincula a los grupos. Los grupos con pocos datos son atraídos hacia la media poblacional (contracción), mientras que los grupos con muchos datos se dejan mayormente solos. Esto produce mejores estimaciones para todos los grupos, especialmente aquellos con tamaños de muestra pequeños.

# Demonstrate partial pooling with a simple simulation
# 5 tanks of tadpoles, each with different sample sizes
set.seed(42)
n_tanks <- 5
true_survival <- c(0.8, 0.5, 0.7, 0.3, 0.9)  # true survival per tank
tank_sizes <- c(10, 20, 5, 30, 8)  # different sample sizes

# Simulate data

survived <- rbinom(n_tanks, tank_sizes, true_survival)

# No pooling: estimate each tank independently

no_pool <- survived / tank_sizes

# Complete pooling: ignore tank differences

complete_pool <- rep(sum(survived) / sum(tank_sizes), n_tanks)

# Partial pooling (simplified): shrink toward overall mean

overall_mean <- sum(survived) / sum(tank_sizes)

# Shrinkage factor depends on sample size

shrinkage <- tank*sizes / (tank_sizes + 10) # 10 is a tuning parameter
partial_pool <- shrinkage * no*pool + (1 - shrinkage) * overall_mean

# Compare

comparison <- data.frame(
tank = 1:n_tanks,
true = true_survival,
n = tank_sizes,
no_pooling = round(no_pool, 3),
partial_pooling = round(partial_pool, 3),
complete_pooling = round(complete_pool, 3)
)
print(comparison)

# Plot

plot(1:n_tanks, true_survival, pch = 16, col = "black", cex = 1.5,
ylim = c(0, 1), xlab = "Tank", ylab = "Survival probability",
main = "Partial pooling: shrinkage toward the mean")
points(1:n_tanks, no_pool, pch = 17, col = "firebrick")
points(1:n_tanks, partial_pool, pch = 18, col = "steelblue")
points(1:n_tanks, complete_pool, pch = 15, col = "darkgreen")
legend("bottomleft",
legend = c("True", "No pooling", "Partial pooling", "Complete pooling"),
col = c("black", "firebrick", "steelblue", "darkgreen"),
pch = c(16, 17, 18, 15), cex = 0.8)

Observe cómo el agrupamiento parcial atrae las estimaciones de los tanques pequeños (ej., tanque 3 con $n=5$) hacia la media general, mientras que deja los tanques grandes (ej., tanque 4 con $n=30$) casi sin cambios. Este es el beneficio clave de los modelos multinivel: toman prestada fuerza entre grupos.

Herramienta 3: Criterios de información

Los criterios de información (como WAIC y PSIS) estiman la precisión predictiva fuera de muestra. Nos dicen qué modelo predecirá mejor los datos nuevos, sin necesidad de recolectar datos nuevos. Esta es la cura para el sobreajuste: en lugar de elegir el modelo que mejor se ajusta a la muestra (lo cual siempre favorece a los modelos más complejos), elegimos el modelo que mejor generalizará.

La idea central es la desvianza, que mide qué tan bien el modelo se ajusta a los datos:

$\text{Deviance} = -2 \sum_i \log P(y_i \mid \theta)$

Pero la desvianza de los datos de entrenamiento siempre mejora al añadir parámetros. El truco es estimar la desvianza de datos nuevos, lo cual penaliza la complejidad. WAIC (Criterio de Información Amplio) hace esto calculando la densidad predictiva puntual logarítmica y añadiendo una penalización por el número efectivo de parámetros:

$\text{WAIC} = -2(\text{lppd} - p_\text{WAIC})$

donde $\text{lppd}$ es la densidad predictiva puntual logarítmica y $p_\text{WAIC}$ es el número efectivo de parámetros.

# Demonstrate overfitting: more parameters always fit better
set.seed(42)
n <- 20
x <- runif(n, 0, 10)
y <- 2 * x + rnorm(n, sd = 2)

# Fit polynomials of increasing degree

degrees <- 1:8
train_error <- numeric(length(degrees))

for (i in seq_along(degrees)) {
d <- degrees[i]
fit <- lm(y ~ poly(x, degree = d, raw = TRUE))
train_error[i] <- mean(resid(fit)^2) # training MSE
}

# Plot training error vs model complexity

plot(degrees, train_error, type = "b", pch = 16, lwd = 2,
col = "firebrick", xlab = "Polynomial degree",
ylab = "Training MSE",
main = "Overfitting: training error always decreases")
grid()

# Add a simulated test error curve

test*error <- train_error + 0.3 * degrees^2 _ 0.1
lines(degrees, test_error, type = "b", pch = 17, lwd = 2, col = "steelblue")
legend("topright",
legend = c("Training error (in-sample)", "Test error (out-of-sample)"),
col = c("firebrick", "steelblue"), pch = c(16, 17), lty = 1, lwd = 2)

El error de entrenamiento siempre disminuye al añadir parámetros — el golem ajusta felizmente el ruido. Pero el error de prueba (estimado fuera de muestra) tiene un mínimo: más allá de cierta complejidad, el modelo empieza a sobreajustar y generaliza peor. Los criterios de información estiman este error fuera de muestra sin necesidad de un conjunto de prueba separado.

Cómo funcionan juntas las tres herramientas

Las tres herramientas forman un flujo de trabajo coherente:

1. Especificar: Escribir un modelo científico y traducirlo a un modelo bayesiano (verosimilitud + distribuciones a priori).
2. Ajustar: Usar inferencia bayesiana (aproximación cuadrática, MCMC) para calcular la distribución a posteriori.
3. Criticar: Usar comprobaciones predictivas a posteriori para ver si el modelo puede reproducir los datos.
4. Comparar: Usar criterios de información para comparar modelos y estimar cuál predecirá mejor.

Si el modelo falla la crítica, volvemos al paso 1 y especificamos un modelo mejor. Si múltiples modelos pasan la crítica, los comparamos y posiblemente promediamos sus predicciones.

# The full workflow: specify, fit, critique, compare
# Model 1: Linear relationship
# Model 2: Quadratic relationship

set.seed(42)
x <- runif(50, 0, 10)
y <- 1 + 2 _ x + 0.3 _ x^2 + rnorm(50, sd = 3)

# Step 1 & 2: Specify and fit both models

m1 <- lm(y ~ x) # linear
m2 <- lm(y ~ x + I(x^2)) # quadratic

# Step 3: Critique - compare R-squared

cat("Model 1 (linear): R-squared =", round(summary(m1)$r.squared, 4), "\n")
cat("Model 2 (quadratic): R-squared =", round(summary(m2)$r.squared, 4), "\n")

# Step 4: Compare using AIC (a simple information criterion)

cat("\nAIC (lower is better):\n")
cat("Model 1:", round(AIC(m1), 2), "\n")
cat("Model 2:", round(AIC(m2), 2), "\n")
cat("\nDifference:", round(AIC(m1) - AIC(m2), 2), "\n")
cat("Model 2 is preferred by", round(AIC(m1) - AIC(m2), 2), "AIC units\n")

Resumen de la Sección 1.3

• La inferencia bayesiana actualiza creencias bajo incertidumbre usando el teorema de Bayes: posterior $\propto$ verosimilitud $\times$ a priori.
• Los modelos multinivel manejan datos estructurados mediante agrupamiento parcial: los grupos toman prestada fuerza entre sí.
• Los criterios de información estiman la precisión predictiva fuera de muestra para prevenir el sobreajuste.
• Las tres herramientas forman un flujo de trabajo: especificar, ajustar, criticar, comparar.
• Estas son prácticas de ingeniería, no conjuros mágicos — hacen que el modelado sea transparente e iterativo.