13.1 Pendientes variables por construcción

Las pendientes variables extienden los interceptos variables permitiendo que el efecto de un predictor difiera entre grupos. El modelo estima una matriz de covarianza de interceptos y pendientes.

13.1.1 El modelo de pendientes variables

En lugar de solo interceptos variables, ahora dejamos que tanto el intercepto como la pendiente varíen por grupo:

$\begin{pmatrix} \gamma_j \\ \beta_j \end{pmatrix} \sim \text{MVNormal}\left(\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_\gamma^2 & \rho\sigma_\gamma\sigma_\beta \\ \rho\sigma_\gamma\sigma_\beta & \sigma_\beta^2 \end{pmatrix}\right)$

La adición clave es la correlación $\rho$ entre interceptos y pendientes — los grupos que empiezan alto pueden también cambiar más rápido o más lento.

# Varying slopes: cafe example
# Each cafe has its own intercept (morning visits) and slope (afternoon effect)

set.seed(42)
n_cafes <- 20

# Population parameters

a_bar <- 3.5 # average morning visits
b_bar <- -0.9 # average afternoon effect (fewer visits)
sigma_a <- 1.0 # between-cafe SD for intercepts
sigma_b <- 0.5 # between-cafe SD for slopes
rho <- -0.6 # correlation: cafes with more morning visits lose more in afternoon

# Draw cafe-specific parameters from MVN

library(MASS)
Sigma <- matrix(c(sigma*a^2, rho * sigma*a * sigma*b,
rho * sigma*a * sigma_b, sigma_b^2), nrow = 2)
cafe_params <- mvrnorm(n_cafes, c(0, 0), Sigma)
a_cafe <- a_bar + cafe_params[, 1]
b_cafe <- b_bar + cafe_params[, 2]

# Simulate observations

n_per <- 10
d <- do.call(rbind, lapply(1:n_cafes, function(j) {
afternoon <- rep(c(0, 1), each = n_per / 2)
y <- a_cafe[j] + b_cafe[j] * afternoon + rnorm(n_per, 0, 0.5)
data.frame(cafe = factor(j), afternoon = afternoon, y = y)
}))

cat("Cafe data simulated:\n")
cat(" Cafes:", n_cafes, "\n")
cat(" Observations per cafe:", n_per, "\n")
cat(" True intercept SD:", sigma_a, "\n")
cat(" True slope SD:", sigma_b, "\n")
cat(" True correlation:", rho, "\n")

# Fit varying slopes model

library(lme4)
m <- lmer(y ~ afternoon + (1 + afternoon | cafe), data = d)

cat("\nVarying slopes model:\n")
cat(" Fixed: intercept =", round(fixef(m)[1], 3),
" afternoon =", round(fixef(m)[2], 3), "\n")
cat(" Random SD: intercept =", round(attr(VarCorr(m)$cafe, "stddev")[1], 3),
  "  slope =", round(attr(VarCorr(m)$cafe, "stddev")[2], 3), "\n")
cat(" Correlation:", round(attr(VarCorr(m)$cafe, "correlation")[1, 2], 3), "\n")

Click Run to execute this code in your browser with WebR

13.1.2 La distribución a priori LKJ para matrices de correlación

La matriz de correlación necesita una distribución a priori. La distribución LKJ es una elección estándar — regulariza hacia correlaciones cero:

$p(\mathbf{R}) \propto |\mathbf{R}|^{\eta - 1}$

Cuando $\eta = 1$ , todas las matrices de correlación son igualmente probables. Cuando $\eta > 1$ , la distribución a priori se concentra hacia $\rho = 0$ .

# Demonstrate LKJ prior for correlation matrices
# Using rethinking or manual simulation

set.seed(42)
n_samples <- 5000

# LKJ(eta = 1): uniform over correlations

# LKJ(eta = 2): favors rho near 0

# LKJ(eta = 4): strongly favors rho near 0

# Simple 2x2 case: sample correlations

# For 2x2, LKJ reduces to: p(rho) proportional to (1 - rho^2)^(eta - 1)

sample_lkj_2x2 <- function(eta, n) {

# Use Beta distribution trick

# rho ~ 2 * Beta(eta, eta) - 1

2 * rbeta(n, eta, eta) - 1
}

rho_1 <- sample_lkj_2x2(1, n_samples)
rho_2 <- sample_lkj_2x2(2, n_samples)
rho_4 <- sample_lkj_2x2(4, n_samples)

cat("LKJ prior on correlation (2x2 case):\n")
cat(" eta = 1: mean rho =", round(mean(rho_1), 3), " SD =", round(sd(rho_1), 3), "\n")
cat(" eta = 2: mean rho =", round(mean(rho_2), 3), " SD =", round(sd(rho_2), 3), "\n")
cat(" eta = 4: mean rho =", round(mean(rho_4), 3), " SD =", round(sd(rho_4), 3), "\n")

cat("\n Higher eta = more regularization toward rho = 0\n")

# Plot

par(mfrow = c(1, 3))
hist(rho_1, breaks = 50, col = "steelblue", main = "LKJ(eta=1)",
xlab = "rho", xlim = c(-1, 1))
hist(rho_2, breaks = 50, col = "steelblue", main = "LKJ(eta=2)",
xlab = "rho", xlim = c(-1, 1))
hist(rho_4, breaks = 50, col = "steelblue", main = "LKJ(eta=4)",
xlab = "rho", xlim = c(-1, 1))
par(mfrow = c(1, 1))

Click Run to execute this code in your browser with WebR

Rethinking

Las pendientes variables son interacciones. Una pendiente variable es solo una interacción entre un predictor y un factor de agrupamiento — pero con agrupamiento parcial. En un modelo de efectos fijos, escribiría y ~ x * group para obtener pendientes específicas por grupo. En un modelo multinivel, escribe (1 + x | group) y las pendientes se extraen de una distribución común. Esto significa que las pendientes se encogen hacia el promedio poblacional, con más encogimiento para grupos con menos datos. La correlación entre interceptos y pendientes es el análogo multinivel de permitir que la interacción covaríe con el efecto principal. Todo está conectado.

Resumen de la Sección 13.1

Las pendientes variables permiten que el efecto de un predictor difiera entre grupos.
El modelo estima una matriz de covarianza con varianzas y una correlación $\rho$ .
Use (1 + x | group) en lme4 para interceptos y pendientes variables.
La distribución a priori LKJ regulariza la matriz de correlación hacia $\rho = 0$ .
Las pendientes variables son interacciones agrupadas — interacciones con encogimiento.

Práctica

Ejemplo: Decisiones de admisión y género