Statistical Rethinking

4.4 Añadiendo un predictor

Añadir un predictor a la media convierte un modelo Gaussiano simple en una regresión lineal. El valor esperado del resultado se convierte en una función lineal del predictor, mientras que las observaciones se dispersan alrededor de esta expectativa con variación residual.

4.4.1 El modelo lineal

El modelo de regresión lineal para la estatura con el peso como predictor es:

$h_i \sim \text{Normal}(\mu_i, \sigma)$ $\mu_i = \alpha + \beta x_i$ $\alpha \sim \text{Normal}(178, 20)$ $\beta \sim \text{Log-Normal}(0, 1)$ $\sigma \sim \text{Uniform}(0, 50)$

La distribución a priori para $\beta$ es log-normal, lo que asegura que $\beta > 0$ (personas más altas pesan más). La distribución log-normal a priori tiene mediana 1, lo que significa que esperamos aproximadamente 1 cm de estatura por kg de peso.

# Simulate the Howell data with height and weight
set.seed(42)
n <- 352
weight <- rnorm(n, mean = 35, sd = 5)
height <- 100 + 2.5 * weight + rnorm(n, sd = 4)

# Plot the data

plot(weight, height, pch = 16, col = "steelblue",
main = "Howell data: height vs weight",
xlab = "Weight (kg)", ylab = "Height (cm)")

# Fit linear regression

fit <- lm(height ~ weight)
abline(fit, lwd = 2, col = "firebrick")
cat("Intercept (alpha):", round(coef(fit)[1], 2), "\n")
cat("Slope (beta):", round(coef(fit)[2], 2), "\n")
cat("R-squared:", round(summary(fit)$r.squared, 4), "\n")

Click Run to execute this code in your browser with WebR

4.4.2 Simulación predictiva a priori

Antes de ajustar, verificamos qué implican las distribuciones a priori sobre la relación.

# Prior predictive simulation for the linear regression
set.seed(42)

# Priors

n_draws <- 50
alpha_prior <- rnorm(n_draws, 178, 20)
beta_prior <- rlnorm(n_draws, 0, 1) # log-normal: always positive

# Plot prior predictive lines

plot(NULL, xlim = c(0, 100), ylim = c(0, 300),
xlab = "Weight (kg)", ylab = "Height (cm)",
main = "Prior predictive: 50 lines from the prior")
for (i in 1:n_draws) {
abline(a = alpha_prior[i], b = beta_prior[i],
col = rgb(0.27, 0.51, 0.71, 0.2))
}
abline(h = 0, lty = 2, col = "gray50")
abline(h = 272, lty = 2, col = "gray50") # tallest person ever

cat("Prior for alpha: Normal(178, 20)\n")
cat(" Range:", round(range(alpha_prior), 1), "\n")
cat("Prior for beta: Log-Normal(0, 1)\n")
cat(" Range:", round(range(beta_prior), 1), "\n")
cat("Median of beta prior:", round(median(beta_prior), 2), "\n")
cat("\nSome prior lines are unreasonable (e.g., 300 cm at 50 kg)\n")
cat("But the data will constrain the posterior\n")

Click Run to execute this code in your browser with WebR

4.4.3 Ajustando el modelo

Ajustamos el modelo usando aproximación cuadrática. La posterior nos da la distribución conjunta de $\alpha$ , $\beta$ y $\sigma$ .

# Fit the linear regression using lm (quadratic approximation equivalent)
set.seed(42)
n <- 352
weight <- rnorm(n, mean = 35, sd = 5)
height <- 100 + 2.5 * weight + rnorm(n, sd = 4)

fit <- lm(height ~ weight)

# Posterior means and standard deviations

cat("Posterior means (quadratic approximation):\n")
cat("alpha:", round(coef(fit)[1], 3), "\n")
cat("beta:", round(coef(fit)[2], 3), "\n")
cat("sigma:", round(summary(fit)$sigma, 3), "\n")

# 95% credible intervals (approximate)

conf <- confint(fit)
cat("\n95% intervals:\n")
cat("alpha: [", round(conf[1, 1], 3), ",", round(conf[1, 2], 3), "]\n")
cat("beta: [", round(conf[2, 1], 3), ",", round(conf[2, 2], 3), "]\n")

Click Run to execute this code in your browser with WebR

4.4.4 Interpretando la posterior

La posterior para $\beta$ nos dice cuánto cambia la estatura por unidad de peso. La posterior para $\alpha$ nos dice la estatura esperada cuando el peso es cero (lo cual no es significativo — es una extrapolación).

# Interpret the posterior
set.seed(42)
n <- 352
weight <- rnorm(n, mean = 35, sd = 5)
height <- 100 + 2.5 * weight + rnorm(n, sd = 4)

fit <- lm(height ~ weight)

# Posterior samples (using the multivariate normal approximation)

library(MASS)
post_samples <- mvrnorm(1e4, mu = coef(fit), Sigma = vcov(fit))
alpha_samples <- post_samples[, 1]
beta_samples <- post_samples[, 2]
sigma_samples <- summary(fit)$sigma * sqrt((n - 2) / rchisq(1e4, n - 2))

# Posterior summaries

cat("Posterior for beta (cm per kg):\n")
cat(" Mean:", round(mean(beta_samples), 3), "\n")
cat(" 95% interval: [", round(quantile(beta_samples, 0.025), 3), ",",
round(quantile(beta_samples, 0.975), 3), "]\n")

# Plot posterior for beta

hist(beta_samples, breaks = 50, col = "steelblue", border = "white",
freq = FALSE, xlab = "beta (cm per kg)",
main = "Posterior distribution of beta")
abline(v = quantile(beta_samples, c(0.025, 0.975)), lty = 2, col = "firebrick", lwd = 2)

Click Run to execute this code in your browser with WebR

4.4.5 Líneas predictivas a posteriori

Podemos visualizar la incertidumbre en la línea de regresión graficando muchas líneas muestreadas de la posterior.

# Posterior predictive lines
set.seed(42)
n <- 352
weight <- rnorm(n, mean = 35, sd = 5)
height <- 100 + 2.5 * weight + rnorm(n, sd = 4)

fit <- lm(height ~ weight)
library(MASS)
post_samples <- mvrnorm(100, mu = coef(fit), Sigma = vcov(fit))

# Plot data and posterior lines

plot(weight, height, pch = 16, col = rgb(0.27, 0.51, 0.71, 0.5),
main = "Posterior predictive lines",
xlab = "Weight (kg)", ylab = "Height (cm)")
for (i in 1:100) {
abline(a = post_samples[i, 1], b = post_samples[i, 2],
col = rgb(0.8, 0.2, 0.2, 0.1))
}
abline(fit, lwd = 2, col = "darkred")

Click Run to execute this code in your browser with WebR

Rethinking

El intercepto suele carecer de significado. El intercepto $\alpha$ es el valor esperado del resultado cuando el predictor es cero. Pero un peso de cero no es un valor significativo — es una extrapolación muy fuera de los datos. Para hacer que el intercepto sea significativo, centre el predictor: reste el peso medio de cada observación. Entonces $\alpha$ se convierte en la estatura esperada en el peso promedio, que es una cantidad significativa. Centrar también hace que la posterior de $\alpha$ y $\beta$ esté menos correlacionada, lo que mejora el cálculo.

Overthinking

Centrado y estandarización. Centrar el predictor ( $x_i - \bar{x}$ ) hace que el intercepto sea interpretable. La estandarización ( $\frac{x_i - \bar{x}}{s_x}$ ) va más allá: hace que la pendiente sea interpretable como el cambio esperado en $y$ por desviación estándar de $x$ . Esto es útil para comparar los efectos de diferentes predictores medidos en escalas diferentes. En el paquete rethinking, la función standardize() hace esto automáticamente.

Resumen de la Sección 4.4

Añadir un predictor convierte el modelo Gaussiano en una regresión lineal: $\mu_i = \alpha + \beta x_i$ .
La distribución log-normal a priori para $\beta$ asegura una relación positiva (personas más altas pesan más).
La simulación predictiva a priori verifica si las distribuciones a priori implican líneas de regresión razonables.
La posterior para $\beta$ nos dice cuánto cambia el resultado por unidad del predictor.
Centrar el predictor hace que el intercepto sea interpretable y mejora el cálculo.
Las líneas predictivas a posteriori visualizan la incertidumbre en la relación de regresión.

Un modelo gaussiano de altura

Regresión polinomial