4  Semana 4. Pendientes variables e interacciones entre niveles

SP-1653 Modelos Mixtos

4.1 Panorama de la semana

En las semanas anteriores construimos el modelo multinivel de manera progresiva. Primero identificamos la dependencia inducida por el agrupamiento y contrastamos pooling completo, ausencia de pooling y pooling parcial. Después permitimos que el intercepto variara entre grupos, lo que condujo al ICC y al shrinkage. En la semana 3 incorporamos predictores individuales y grupales y separamos explícitamente asociaciones dentro y entre grupos mediante centrado y descomposición.

Hasta ahora, sin embargo, hemos mantenido una restricción fuerte: una vez fijado el predictor individual, la asociación dentro de los grupos se ha supuesto igual para todos ellos. Por ejemplo, en la aplicación de rendimiento matemático escribimos

\[ MathAch_{ij} = \alpha_j +\beta_W SES^W_{ij} +\beta_B \overline{SES}^{B}_j +\varepsilon_{ij}, \]

donde el intercepto \(\alpha_j\) cambia entre escuelas, pero \(\beta_W\) es común a todas.

Esta semana retiramos esa restricción. La pregunta central es:

¿Qué ocurre cuando la asociación entre una respuesta y un predictor cambia entre grupos?

Gelman y Hill presentan los modelos con interceptos y pendientes variables como la extensión natural del modelo de interceptos variables y subrayan que una pendiente variable puede verse como una interacción entre un predictor continuo y los indicadores de grupo (Gelman y Hill 2007, 237-38, 279-83). McElreath desarrolla la misma idea desde una perspectiva generativa: interceptos y pendientes se conciben como realizaciones conjuntas de una población de coeficientes y el pooling ocurre simultáneamente en ambas dimensiones (McElreath 2020, 437-46).

El paso adicional de esta semana será preguntar si parte de la heterogeneidad de pendientes puede explicarse mediante características observadas de los grupos. Esa pregunta conduce a las interacciones entre niveles.

Objetivos de aprendizaje

Al finalizar esta semana, se espera que la persona estudiante pueda:

1. escribir un modelo gaussiano con interceptos y pendientes variables en forma jerárquica y en forma de componentes;
2. distinguir la pendiente poblacional promedio \(\mu_\beta\) de las pendientes específicas \(\beta_j\);
3. interpretar \(\tau_\alpha\), \(\tau_\beta\), \(\rho\) y \(\sigma\) sin confundir variación entre grupos con incertidumbre posterior;
4. representar la matriz de covarianza grupal como \(\Sigma=D R D\) y explicar por qué esa factorización separa escalas y correlaciones;
5. explicar cómo el pooling parcial se extiende de una a dos dimensiones cuando interceptos y pendientes covarían;
6. demostrar por qué la correlación entre interceptos y pendientes depende del punto de referencia del predictor;
7. derivar cómo cambia la varianza entre grupos a lo largo del eje del predictor;
8. explicar por qué, con pendientes variables, la correlación inducida entre observaciones del mismo grupo ya no se resume mediante un único ICC constante;
9. distinguir una interacción convencional, una pendiente variable y una interacción entre niveles;
10. formular una interacción entre niveles como un modelo para la pendiente grupal y obtener la ecuación combinada por sustitución;
11. interpretar una interacción entre niveles como heterogeneidad sistemática de la pendiente, sin confundirla con la heterogeneidad residual \(\tau_\beta\);
12. reconocer cuándo una pendiente variable está débilmente identificada debido a poca variación del predictor dentro de los grupos;
13. justificar por qué normalmente una pendiente variable se acompaña de un intercepto variable y reconocer excepciones sustantivas;
14. especificar previas explícitas para medias poblacionales, desviaciones estándar grupales, correlación y desviación estándar residual;
15. realizar una comprobación predictiva previa sobre familias de líneas grupales;
16. ajustar modelos de pendientes variables en brms usando cmdstanr, diagnosticar el muestreo y realizar comprobaciones predictivas sensibles a la heterogeneidad de pendientes;
17. distinguir predicción para un grupo observado de predicción para un grupo nuevo;
18. proponer, para el proyecto del curso, qué coeficientes deberían variar entre grupos y qué variables grupales podrían explicar parte de esa variación.

4.2 Problema motivador: ¿todas las escuelas tienen la misma pendiente?

Retomemos los datos de estudiantes agrupados en escuelas. En la semana 3 usamos

\[ SES^W_{ij} = SES_{ij}-\overline{SES}_j \]

para representar la posición socioeconómica relativa de un estudiante dentro de su escuela. Un modelo con pendiente común es

\[ MathAch_{ij} \sim \mathcal N(\mu_{ij},\sigma^2), \]

\[ \mu_{ij} = \alpha_j +\beta_W SES^W_{ij} +\beta_B\overline{SES}^{B}_j. \]

La interpretación de \(\beta_W\) es clara: para dos estudiantes de la misma escuela, una unidad de diferencia en SES está asociada con \(\beta_W\) unidades de diferencia esperada en rendimiento.

Pero el modelo hace una afirmación adicional que a veces pasa inadvertida:

\[ \boxed{ \text{la misma pendiente }\beta_W\text{ se aplica a todas las escuelas.} } \]

Esa restricción podría ser razonable. También podría ser demasiado fuerte. El gradiente socioeconómico dentro de una escuela podría cambiar con su organización, composición, recursos, políticas o cualquier otra característica que modifique la relación entre SES individual y rendimiento.

Una primera expansión consiste en escribir

\[ \beta_W \longrightarrow \beta_j, \]

de modo que

\[ \mu_{ij} = \alpha_j+\beta_j SES^W_{ij}+\beta_B\overline{SES}^{B}_j. \]

Ahora cada escuela tiene su propia recta.

ImportanteLa pregunta cambia en dos etapas

Primero preguntamos si existe heterogeneidad de pendientes: ¿cuánto difieren las asociaciones entre grupos?

Después preguntamos si podemos explicar parte de esa heterogeneidad mediante predictores grupales.

Estas son preguntas relacionadas, pero no idénticas.

4.3 De una pendiente común a pendientes variables

Considere inicialmente un único predictor \(x_{ij}\) centrado en una referencia interpretable. El modelo con interceptos y pendientes variables es

\[ y_{ij} \mid \alpha_j,\beta_j,\sigma \sim \mathcal N(\alpha_j+\beta_j x_{ij},\sigma^2). \]

Los pares de coeficientes grupales se modelan conjuntamente:

\[ \begin{pmatrix} \alpha_j\\ \beta_j \end{pmatrix} \sim \mathcal N_2 \left( \begin{pmatrix} \mu_\alpha\\ \mu_\beta \end{pmatrix}, \Sigma \right), \qquad j=1,\ldots,J. \]

Con dos coeficientes variables,

\[ \Sigma = \begin{pmatrix} \tau_\alpha^2 & \rho\tau_\alpha\tau_\beta\\ \rho\tau_\alpha\tau_\beta & \tau_\beta^2 \end{pmatrix}. \]

Esta es la estructura básica de Gelman y Hill en su modelo (13.1) (Gelman y Hill 2007, 279-80). McElreath la presenta como una población conjunta de interceptos y pendientes cuyo modelo regulariza simultáneamente ambos coeficientes y la relación entre ellos (McElreath 2020, 441-43).

Los parámetros tienen papeles diferentes:

Parámetro	Interpretación
\(\alpha_j\)	intercepto del grupo \(j\) en el punto \(x=0\)
\(\beta_j\)	pendiente del predictor en el grupo \(j\)
\(\mu_\alpha\)	intercepto promedio en la población de grupos
\(\mu_\beta\)	pendiente promedio en la población de grupos
\(\tau_\alpha\)	desviación estándar entre interceptos
\(\tau_\beta\)	desviación estándar entre pendientes
\(\rho\)	correlación poblacional entre interceptos y pendientes
\(\sigma\)	desviación estándar residual dentro de los grupos

NotaVariabilidad e incertidumbre siguen siendo conceptos distintos

Una posterior amplia para \(\tau_\beta\) significa que tenemos incertidumbre sobre la magnitud de la heterogeneidad de pendientes. Una posterior concentrada en un valor grande de \(\tau_\beta\) significa que aprendimos con relativa precisión que las pendientes difieren sustancialmente.

4.3.1 Forma de componentes

Es útil separar medias poblacionales y desviaciones grupales:

\[ \alpha_j = \mu_\alpha+u_{0j}, \]

\[ \beta_j = \mu_\beta+u_{1j}, \]

con

\[ \begin{pmatrix} u_{0j}\\ u_{1j} \end{pmatrix} \sim \mathcal N_2 \left[ \begin{pmatrix} 0\\0 \end{pmatrix}, \begin{pmatrix} \tau_\alpha^2 & \rho\tau_\alpha\tau_\beta\\ \rho\tau_\alpha\tau_\beta & \tau_\beta^2 \end{pmatrix} \right]. \]

Sustituyendo,

\[ \begin{aligned} y_{ij} &= \mu_\alpha +\mu_\beta x_{ij} +u_{0j} +u_{1j}x_{ij} +\varepsilon_{ij},\\ \varepsilon_{ij} &\sim \mathcal N(0,\sigma^2). \end{aligned} \]

Esta expresión muestra algo importante: la contribución grupal ya no es una constante \(u_{0j}\). Es

\[ \boxed{ u_{0j}+u_{1j}x_{ij}}, \]

y por tanto cambia con el valor de \(x\).

4.4 La distribución conjunta de coeficientes

Una forma especialmente útil de escribir la matriz de covarianza es

\[ \boxed{ \Sigma=D R D } \]

con

\[ D = \begin{pmatrix} \tau_\alpha&0\\ 0&\tau_\beta \end{pmatrix}, \]

y

\[ R = \begin{pmatrix} 1&\rho\\ \rho&1 \end{pmatrix}. \]

La factorización separa dos preguntas:

1. ¿cuánto varían los interceptos y las pendientes? Esto lo describen \(\tau_\alpha\) y \(\tau_\beta\);
2. ¿cómo se asocian ambos tipos de coeficiente? Esto lo describe \(R\) y, en dos dimensiones, \(\rho\).

Esta descomposición es la que utiliza McElreath para construir el modelo de pendientes variables (McElreath 2020, 441-43) y coincide con la parametrización moderna de coeficientes grupales descrita para brms, donde las matrices de covarianza se expresan mediante desviaciones estándar y una matriz de correlación (Bürkner 2017, 3-4). Bayesian Workflow emplea la misma factorización al desarrollar un modelo de interceptos y pendientes variables para el estudio de privación de sueño (Gelman et al. 2026, 281-84).

4.4.1 ¿Qué significa \(\rho\)?

Si \(\rho>0\), grupos con interceptos relativamente altos tienden a tener pendientes relativamente altas.

Si \(\rho<0\), grupos con interceptos relativamente altos tienden a tener pendientes relativamente bajas.

Si \(\rho\approx0\), conocer el intercepto de un grupo aporta poca información lineal sobre su pendiente dentro de la población modelada.

La palabra tienden es esencial. \(\rho\) describe una distribución de grupos; no impone una relación determinista.

4.4.2 Geometría de la población de líneas

En dos dimensiones, la distribución de \((\alpha_j,\beta_j)\) puede visualizarse como una nube elíptica. Una correlación positiva inclina la nube hacia arriba; una correlación negativa la inclina hacia abajo. La magnitud de \(\tau_\alpha\) y \(\tau_\beta\) determina sus escalas horizontales y verticales.

El siguiente gráfico es puramente geométrico y no depende de un ajuste.

set.seed(1653)

sim_biv <- function(n, mu_alpha, mu_beta, tau_alpha, tau_beta, rho) {
  R <- matrix(c(1, rho, rho, 1), nrow = 2)
  z <- matrix(rnorm(2 * n), ncol = 2)
  z_cor <- z %*% chol(R)

  tibble::tibble(
    alpha = mu_alpha + tau_alpha * z_cor[, 1],
    beta = mu_beta + tau_beta * z_cor[, 2]
  )
}

geom_rho <- purrr::map_dfr(
  c(-0.75, 0, 0.75),
  ~ sim_biv(
    n = 600,
    mu_alpha = 50,
    mu_beta = 3,
    tau_alpha = 8,
    tau_beta = 2,
    rho = .x
  ) |>
    dplyr::mutate(rho = factor(.x))
)

ggplot2::ggplot(geom_rho, ggplot2::aes(x = alpha, y = beta)) +
  ggplot2::geom_point(alpha = 0.25) +
  ggplot2::facet_wrap(~ rho, nrow = 1) +
  ggplot2::labs(
    x = "Intercepto grupal",
    y = "Pendiente grupal",
    title = "Misma variación marginal, distinta correlación"
  ) +
  ggplot2::theme_minimal()

Figura 4.1: Distribuciones bivariadas hipotéticas de interceptos y pendientes para tres valores de correlación. Las escalas marginales son las mismas; solo cambia la asociación entre ambos coeficientes.

4.5 La correlación intercepto–pendiente depende del origen de \(x\)

La interpretación de \(\rho\) exige especial cuidado. Gelman y Hill muestran un ejemplo de regresiones de ingreso sobre estatura donde una correlación extremadamente negativa entre interceptos y pendientes aparece en buena medida porque el intercepto se define en una estatura igual a cero, muy lejos del rango observado. Al reescalar el predictor, la correlación se vuelve mucho más interpretable (Gelman y Hill 2007, 287-89). Hox et al. también muestran que centrar un predictor con pendiente variable cambia la varianza del intercepto, aunque el modelo lineal resultante pueda ser equivalente en ajuste (Hox et al. 2018, 49-50).

La razón puede verse algebraicamente.

Partimos de

\[ y_{ij}=\alpha_j+\beta_j x_{ij}+\varepsilon_{ij}. \]

Defina

\[ x^*_{ij}=x_{ij}-c. \]

Como \(x_{ij}=x^*_{ij}+c\),

\[ \begin{aligned} y_{ij} &=\alpha_j+\beta_j(x^*_{ij}+c)+\varepsilon_{ij}\\ &=(\alpha_j+c\beta_j)+\beta_jx^*_{ij}+\varepsilon_{ij}. \end{aligned} \]

Por tanto,

\[ \alpha^*_j=\alpha_j+c\beta_j, \qquad \beta^*_j=\beta_j. \]

La covarianza transformada es

\[ \begin{aligned} \operatorname{Cov}(\alpha^*_j,\beta_j) &= \operatorname{Cov}(\alpha_j+c\beta_j,\beta_j)\\ &= \operatorname{Cov}(\alpha_j,\beta_j) +c\operatorname{Var}(\beta_j)\\ &= \rho\tau_\alpha\tau_\beta +c\tau_\beta^2. \end{aligned} \]

Además,

\[ \operatorname{Var}(\alpha^*_j) = \tau_\alpha^2 +2c\rho\tau_\alpha\tau_\beta +c^2\tau_\beta^2. \]

Así,

\[ \rho^* = \frac{ \rho\tau_\alpha\tau_\beta+c\tau_\beta^2 }{ \tau_\beta \sqrt{ \tau_\alpha^2+2c\rho\tau_\alpha\tau_\beta+c^2\tau_\beta^2 } }. \]

La misma familia de rectas puede tener una correlación intercepto–pendiente muy diferente solo porque movimos el punto \(x=0\).

ImportanteUna correlación de coeficientes no es invariante al centrado

Antes de interpretar \(\rho\) sustantivamente, pregunte:

• ¿qué significa \(x=0\)?;
• ¿está dentro del rango observado?;
• ¿sería más clara una referencia como la media global, la media del grupo, el inicio de seguimiento o una dosis clínicamente relevante?

Una correlación extrema puede ser una consecuencia geométrica de un origen poco útil, no necesariamente una relación sustantiva extrema entre procesos.

4.6 La heterogeneidad entre grupos cambia a lo largo de \(x\)

Con interceptos variables, la desviación entre dos grupos era paralela en todo el eje \(x\). Con pendientes variables, la separación entre grupos depende de \(x\).

La desviación grupal de la media condicional es

\[ u_{0j}+u_{1j}x. \]

Su varianza poblacional es

\[ \begin{aligned} V_G(x) &= \operatorname{Var}(u_{0j}+u_{1j}x)\\ &= \tau_\alpha^2 +2x\rho\tau_\alpha\tau_\beta +x^2\tau_\beta^2. \end{aligned} \]

Por tanto,

\[ \boxed{ V_G(x) = \tau_\alpha^2 +2x\rho\tau_\alpha\tau_\beta +x^2\tau_\beta^2 } \]

es la varianza entre las medias de grupo en el valor \(x\).

Varias consecuencias son inmediatas:

• \(\tau_\alpha^2\) es la varianza entre grupos en \(x=0\);
• \(\tau_\beta^2\) controla cuánto cambia la heterogeneidad al alejarnos del origen;
• \(\rho\) determina si inicialmente la dispersión entre grupos aumenta o disminuye al movernos en una dirección del eje;
• el punto de centrado del predictor determina dónde se interpreta directamente \(\tau_\alpha\).

Esta es otra razón para no considerar el centrado una operación cosmética.

4.7 ¿Qué ocurre con el ICC?

En la semana 2, con un único intercepto variable, dos observaciones distintas del mismo grupo tenían covarianza \(\tau_\alpha^2\) y el ICC era constante:

\[ ICC = \frac{\tau_\alpha^2}{\tau_\alpha^2+\sigma^2}. \]

Con una pendiente variable, dos observaciones del mismo grupo tomadas en \(x\) y \(x'\) comparten tanto \(u_{0j}\) como \(u_{1j}\). Su covarianza condicional en los valores de los predictores es

\[ \begin{aligned} \operatorname{Cov}(Y(x),Y(x')\mid x,x',\text{mismo grupo}) &= \operatorname{Cov}(u_0+u_1x,\;u_0+u_1x')\\ &= \tau_\alpha^2 +\rho\tau_\alpha\tau_\beta(x+x') +\tau_\beta^2xx'. \end{aligned} \]

La varianza marginal condicional en \(x\) es

\[ \operatorname{Var}(Y(x)\mid x) = V_G(x)+\sigma^2. \]

Entonces la correlación entre dos observaciones del mismo grupo depende de \(x\) y \(x'\):

\[ \operatorname{Corr}(Y(x),Y(x')\mid x,x',\text{mismo grupo}) = \frac{ \tau_\alpha^2 +\rho\tau_\alpha\tau_\beta(x+x') +\tau_\beta^2xx' }{ \sqrt{[V_G(x)+\sigma^2][V_G(x')+\sigma^2]} }. \]

NotaNo siempre existe un único ICC que resuma la dependencia

Una vez que la estructura grupal incluye pendientes variables, la similitud entre observaciones del mismo grupo puede depender de sus valores del predictor. El ICC del modelo nulo sigue siendo útil como punto de partida, pero ya no resume toda la estructura de dependencia del modelo ampliado.

4.8 Pooling parcial en dos dimensiones

En un modelo de interceptos variables, cada \(\alpha_j\) se regulariza hacia \(\mu_\alpha\). Con pendientes variables, el objeto que se regulariza es el vector

\[ \begin{pmatrix} \alpha_j\\ \beta_j \end{pmatrix}. \]

La estimación de un grupo combina:

• información de las observaciones del propio grupo;
• la distribución poblacional de interceptos;
• la distribución poblacional de pendientes;
• la correlación entre ambos tipos de coeficiente.

McElreath muestra este mecanismo como shrinkage en dos dimensiones: una pendiente extrema puede contraerse hacia el centro de la población y, si interceptos y pendientes están correlacionados, esa contracción puede ir acompañada de un desplazamiento del intercepto (McElreath 2020, 444-46). Bayesian Workflow ilustra el mismo fenómeno comparando regresiones independientes por persona con los coeficientes obtenidos por un modelo multinivel conjunto de interceptos y pendientes (Gelman et al. 2026, 283-84).

El pooling no tiene por qué ser igualmente fuerte para intercepto y pendiente. Un grupo puede informar bien su nivel medio pero muy poco su gradiente si el predictor casi no varía dentro del grupo.

4.9 Pendientes variables como interacciones

Gelman y Hill señalan una equivalencia conceptual útil: una pendiente variable es una interacción entre el predictor \(x\) y la identidad del grupo (Gelman y Hill 2007, 237-38).

Si tuviéramos \(J\) grupos y usáramos indicadores \(I(g=j)\), un modelo sin pooling podría escribirse con interacciones

\[ x_i I(g_i=j). \]

Esto produciría una pendiente separada para cada grupo. El modelo multinivel reemplaza esa colección de coeficientes independientes por una distribución poblacional que permite pooling parcial.

Por tanto, los modelos de pendientes variables pueden verse como máquinas de interacciones parcialmente agrupadas: permiten que el efecto de \(x\) cambie con el grupo, pero regularizan esos cambios mediante una estructura común.

4.10 Tres ideas que no deben confundirse

4.10.1 Interacción convencional

Con dos predictores observados \(x\) y \(z\),

\[ y_i = \alpha +\beta_x x_i +\beta_z z_i +\beta_{xz}x_i z_i +\varepsilon_i. \]

La pendiente de \(x\) es

\[ \frac{\partial E(Y\mid x,z)}{\partial x} = \beta_x+\beta_{xz}z. \]

La asociación de \(x\) cambia sistemáticamente con \(z\).

4.10.2 Pendiente variable

En un modelo multinivel,

\[ \beta_j = \mu_\beta+u_{1j}, \]

\[ u_{1j} \sim \mathcal N(0,\tau_\beta^2). \]

Reconocemos que la asociación entre \(x\) y la respuesta puede diferir entre grupos. En este modelo, esas diferencias se representan mediante \(u_{1j}\), pero no se atribuyen todavía a características observadas de los grupos.

4.10.3 Interacción entre niveles

Sea \(x_{ij}\) un predictor individual y \(z_j\) un predictor grupal. Modelamos la pendiente como

\[ \boxed{ \beta_j = \gamma_{10} +\gamma_{11}z_j +u_{1j} } \]

con

\[ u_{1j}\sim\mathcal N(0,\tau_\beta^2). \]

Ahora \(\gamma_{11}\) describe una parte sistemática de la heterogeneidad de pendientes. La parte que no queda explicada por \(z_j\) permanece en \(u_{1j}\).

Hox et al. recomiendan construir el modelo en ecuaciones por nivel y luego sustituir para hacer explícita la interacción entre niveles (Hox et al. 2018, 25-26). Gelman y Hill muestran que un predictor grupal puede entrar en los modelos tanto del intercepto como de la pendiente (Gelman y Hill 2007, 281-83, 379-80).

4.11 Derivación de una interacción entre niveles

Escribamos un modelo de nivel individual:

\[ y_{ij} = \alpha_j +\beta_jx_{ij} +\varepsilon_{ij}. \]

Modelamos intercepto y pendiente mediante un predictor grupal \(z_j\):

\[ \alpha_j = \gamma_{00} +\gamma_{01}z_j +u_{0j}, \]

\[ \beta_j = \gamma_{10} +\gamma_{11}z_j +u_{1j}. \]

Sustituyendo,

\[ \begin{aligned} y_{ij} &= (\gamma_{00}+\gamma_{01}z_j+u_{0j}) + (\gamma_{10}+\gamma_{11}z_j+u_{1j})x_{ij} +\varepsilon_{ij}\\ &= \gamma_{00} +\gamma_{01}z_j +\gamma_{10}x_{ij} +\gamma_{11}x_{ij}z_j +u_{0j} +u_{1j}x_{ij} +\varepsilon_{ij}. \end{aligned} \]

El término

\[ \boxed{ \gamma_{11}x_{ij}z_j } \]

es la interacción entre niveles.

La pendiente esperada en un grupo con característica \(z\) es

\[ E(\beta_j\mid z_j=z) = \gamma_{10}+\gamma_{11}z. \]

La pendiente particular del grupo es

\[ \beta_j = \gamma_{10}+\gamma_{11}z_j+u_{1j}. \]

ImportanteUna interacción entre niveles no elimina automáticamente las pendientes variables

Si incluimos \(x_{ij}z_j\) y mantenemos \(u_{1j}\), estamos diciendo que \(z_j\) explica parte de la heterogeneidad de pendientes y que todavía puede quedar heterogeneidad residual.

Eliminar \(u_{1j}\) impondría la afirmación más fuerte de que, condicional en \(z_j\), todas las diferencias de pendiente entre grupos han desaparecido.

4.12 Cómo interpretar la interacción

En presencia de una interacción, los coeficientes principales son condicionales. Hox et al. enfatizan que los términos que forman una interacción deben interpretarse conjuntamente y que el valor cero de los predictores determina la interpretación de los efectos principales (Hox et al. 2018, 49, 52-53).

En

\[ \mu_{ij} = \gamma_{00} +\gamma_{01}z_j +\gamma_{10}x_{ij} +\gamma_{11}x_{ij}z_j, \]

la interpretación es:

• \(\gamma_{00}\): media esperada cuando \(x=0\) y \(z=0\);
• \(\gamma_{01}\): asociación con \(z\) cuando \(x=0\);
• \(\gamma_{10}\): asociación con \(x\) cuando \(z=0\);
• \(\gamma_{11}\): cambio en la pendiente de \(x\) por una unidad de cambio en \(z\).

Por eso el centrado es especialmente importante. Si \(z_j\) está centrado en su media,

\[ z_j^*=z_j-\bar z, \]

entonces \(\gamma_{10}\) representa la pendiente de \(x\) para un grupo con valor promedio de \(z\).

Si \(x_{ij}\) está centrado por grupo, entonces \(x=0\) representa una unidad situada en la media de su propio grupo, y el término \(\gamma_{01}\) compara grupos para unidades en esa posición relativa.

4.13 ¿Puede variar la pendiente sin variar el intercepto?

Gelman y Hill advierten que, en la mayoría de aplicaciones, si hay razones para permitir que una pendiente cambie entre grupos también suele haber razones para permitir que el intercepto cambie (Gelman y Hill 2007, 283-84). Un modelo

\[ y_{ij}=\alpha+\beta_jx_{ij}+\varepsilon_{ij} \]

obliga a todas las rectas a atravesar exactamente el mismo punto cuando \(x=0\).

Esa restricción puede ser poco plausible.

Sin embargo, no es imposible que tenga sentido. Gelman y Hill dan como ejemplo una colección de experimentos con una condición control verdaderamente común y tratamientos que sí varían entre experimentos. En ese caso, fijar el intercepto y permitir variación en el efecto del tratamiento puede ser una representación razonable (Gelman y Hill 2007, 283-84).

La regla práctica no debe ser “siempre incluya un intercepto variable”, sino:

pregunte qué igualdad sustantiva está imponiendo al fijar el intercepto y si esa igualdad es defendible.

4.14 Información para estimar una pendiente grupal

La estimación de una pendiente depende de disponer de observaciones en distintos valores del predictor. Para el grupo \(j\), una medida básica de la información disponible sobre la pendiente es

\[ S_{xx,j} = \sum_{i=1}^{n_j} (x_{ij}-\bar{x}_j)^2. \]

En una regresión lineal simple con \(\sigma\) conocida, la desviación estándar del estimador de la pendiente es proporcional a

\[ \frac{\sigma}{\sqrt{S_{xx,j}}}. \]

Esta expresión no corresponde a la distribución posterior del modelo jerárquico, pero resume una característica estructural que continúa siendo relevante:

• un mayor número de observaciones puede aumentar la información disponible para estimar la pendiente;
• una mayor dispersión de \(x\) dentro del grupo aumenta la información sobre la pendiente;
• un mayor ruido residual reduce la precisión con que puede estimarse;
• muchas observaciones concentradas en valores muy similares de \(x\) pueden proporcionar bastante información sobre el nivel medio del grupo, pero poca sobre su pendiente.

El pooling parcial utiliza información de la población de grupos para regularizar las estimaciones grupales, pero no puede compensar completamente la ausencia de variación del predictor dentro de un grupo. Si el diseño observado ofrece poco contraste en \(x\), la pendiente específica de ese grupo estará determinada en mayor medida por la estructura jerárquica y presentará mayor incertidumbre.

4.15 Pocos grupos y correlaciones difíciles

\(\tau_\beta\) y, especialmente, \(\rho\) son propiedades de una población de grupos. Si hay pocos grupos, su posterior puede quedar muy influida por la previa y ser amplia. Además, estimar una correlación requiere aprender simultáneamente dos varianzas y una covarianza.

Esto no significa que exista un número mínimo universal de grupos. La información depende del diseño, de los tamaños de grupo, de la variación de los predictores, de la magnitud de las heterogeneidades y de las previas. La simulación generativa es una forma más adecuada de estudiar un diseño concreto que aplicar una regla rígida de tamaño muestral.

AdvertenciaIdentificación débil no es lo mismo que un problema de MCMC

Una posterior amplia de \(\tau_\beta\) o \(\rho\) puede reflejar que los datos contienen poca información. Una divergencia de HMC es, en cambio, una señal computacional sobre la exploración de la posterior.

Ambos problemas pueden aparecer juntos, pero no deben diagnosticarse como si fueran lo mismo. La semana 7 desarrollará esta distinción con mayor detalle.

4.16 Previas para interceptos, pendientes y correlaciones

El modelo introduce nuevos parámetros y, por tanto, nuevas decisiones de escala.

Una especificación genérica es

\[ \mu_\alpha \sim \mathcal N(m_\alpha,s_\alpha^2), \]

\[ \mu_\beta \sim \mathcal N(m_\beta,s_\beta^2), \]

\[ \tau_\alpha \sim \operatorname{HalfNormal}(0,s_{\tau_\alpha}^2), \]

\[ \tau_\beta \sim \operatorname{HalfNormal}(0,s_{\tau_\beta}^2), \]

\[ \sigma \sim \operatorname{HalfNormal}(0,s_\sigma^2), \]

y

\[ R\sim LKJ(\eta). \]

No existe una elección universal de las escalas \(s\). Deben relacionarse con la unidad de la respuesta y con las escalas de los predictores.

4.16.1 Previa LKJ

En dos dimensiones,

\[ R = \begin{pmatrix} 1&\rho\\ \rho&1 \end{pmatrix}. \]

McElreath usa \(LKJ(2)\) como una previa regularizadora que reduce plausibilidad de correlaciones extremas (McElreath 2020, 442-43). Bayesian Workflow señala que \(LKJ(1)\) es uniforme sobre matrices de correlación y, en el caso bidimensional, induce una marginal uniforme sobre \(\rho\); para \(\eta>1\) la densidad se concentra más cerca de la identidad (Gelman et al. 2026, 282-83).

Para una matriz \(2\times2\), la densidad marginal es proporcional a

\[ p(\rho) \propto (1-\rho^2)^{\eta-1}, \qquad -1<\rho<1. \]

Así:

• \(\eta=1\) no favorece un valor particular de \(\rho\) en dos dimensiones;
• \(\eta>1\) regulariza hacia correlaciones moderadas;
• valores muy grandes de \(\eta\) pueden imponer una regularización demasiado fuerte si la correlación es sustantivamente importante.

La semana 6 estudiará las previas de forma sistemática. Aquí nos interesa que las previas produzcan familias de rectas plausibles antes de observar los datos.

4.17 Comprobación predictiva previa para una familia de líneas

En un modelo de pendiente fija, una comprobación previa puede concentrarse en el rango de \(y\). Con pendientes variables debemos revisar también:

• cuánta dispersión hay entre interceptos;
• cuánta dispersión hay entre pendientes;
• si aparecen líneas con gradientes absurdamente grandes;
• si la correlación genera familias de rectas plausibles;
• si las predicciones explotan en los extremos del rango de \(x\).

Bayesian Workflow recomienda evaluar conjuntamente las implicaciones de las previas mediante simulación predictiva, no parámetro por parámetro de manera aislada (Gelman et al. 2026, 89-96, 275-84).

En el laboratorio haremos esta comprobación mediante simulación directa desde las previas antes de ajustar el modelo.

4.18 Laboratorio reproducible en R

El laboratorio sigue una secuencia generativa:

1. simular una población de grupos con interceptos y pendientes correlacionados;
2. observar que los grupos aportan cantidades distintas de información sobre sus pendientes;
3. inspeccionar las implicaciones de las previas antes del ajuste;
4. comparar un modelo con pendiente común con un modelo de pendientes variables;
5. visualizar shrinkage simultáneo en interceptos y pendientes;
6. comprobar empíricamente que la correlación intercepto–pendiente cambia al mover el origen de \(x\);
7. realizar una comprobación predictiva específicamente sensible a la dispersión de pendientes.

4.18.1 Simular una población de grupos

Usaremos \(J=30\) grupos con tamaños diferentes. El proceso generador tendrá

\[ \mu_\alpha=50, \qquad \mu_\beta=3, \]

\[ \tau_\alpha=8, \qquad \tau_\beta=2, \]

\[ \rho=-0.60, \qquad \sigma=6. \]

Además, permitiremos que la dispersión de \(x\) sea distinta entre grupos. Esto hace que algunos grupos informen mucho mejor su pendiente que otros.

sim_bivnorm <- function(n, mu, sd, rho) {
  stopifnot(length(mu) == 2, length(sd) == 2, abs(rho) < 1)

  R <- matrix(c(1, rho, rho, 1), nrow = 2)
  z <- matrix(rnorm(2 * n), ncol = 2)
  z_cor <- z %*% chol(R)

  sweep(
    sweep(z_cor, 2, sd, FUN = "*"),
    2,
    mu,
    FUN = "+"
  )
}

J <- 30
niveles_grupo <- sprintf("G%02d", seq_len(J))

mu_alpha_real <- 50
mu_beta_real <- 3
tau_alpha_real <- 8
tau_beta_real <- 2
rho_real <- -0.60
sigma_real <- 6

coef_reales <- sim_bivnorm(
  n = J,
  mu = c(mu_alpha_real, mu_beta_real),
  sd = c(tau_alpha_real, tau_beta_real),
  rho = rho_real
)

parametros_grupo <- tibble(
  grupo = factor(niveles_grupo, levels = niveles_grupo),
  n = sample(8:38, size = J, replace = TRUE),
  sd_x = runif(J, min = 0.35, max = 1.40),
  alpha_real = coef_reales[, 1],
  beta_real = coef_reales[, 2]
)

parametros_grupo |>
  summarise(
    grupos = n(),
    n_min = min(n),
    n_mediana = median(n),
    n_max = max(n),
    cor_realizada = cor(alpha_real, beta_real)
  )

# A tibble: 1 × 5
  grupos n_min n_mediana n_max cor_realizada
   <int> <int>     <dbl> <int>         <dbl>
1     30    10      20.5    38        -0.763

La correlación observada entre los 30 pares simulados no tiene por qué ser exactamente \(-0.60\): los grupos son una muestra finita de la población generadora.

Generamos ahora las observaciones. Dentro de cada grupo hacemos que la media empírica de \(x\) sea exactamente cero para que el intercepto tenga una interpretación simple: respuesta esperada en la posición media del predictor dentro del grupo.

sim <- parametros_grupo |>
  tidyr::uncount(n, .id = "i") |>
  group_by(grupo) |>
  mutate(
    x = rnorm(n(), mean = 0, sd = first(sd_x)),
    x = x - mean(x),
    y = rnorm(
      n(),
      mean = first(alpha_real) + first(beta_real) * x,
      sd = sigma_real
    )
  ) |>
  ungroup()

stopifnot(
  !anyNA(sim),
  nlevels(sim$grupo) == J,
  max(abs(sim |> group_by(grupo) |> summarise(mx = mean(x)) |> pull(mx))) < 1e-10
)

info_grupo <- sim |>
  group_by(grupo) |>
  summarise(
    n = n(),
    Sxx = sum((x - mean(x))^2),
    rango_x = diff(range(x)),
    alpha_real = first(alpha_real),
    beta_real = first(beta_real),
    .groups = "drop"
  )

info_grupo |>
  arrange(Sxx) |>
  select(grupo, n, Sxx, rango_x) |>
  slice_head(n = 8)

# A tibble: 8 × 4
  grupo     n   Sxx rango_x
  <fct> <int> <dbl>   <dbl>
1 G04      22  2.02    1.00
2 G24      18  3.88    1.91
3 G01      12  4.09    1.89
4 G21      29  4.26    1.54
5 G25      21  4.64    1.59
6 G03      19  5.69    1.88
7 G17      23  6.09    2.14
8 G07      19  6.09    1.92

Dos grupos con el mismo \(n_j\) pueden tener información muy distinta si su dispersión de \(x\) es diferente.

4.18.2 Visualizar las líneas verdaderas y los datos

seleccion <- niveles_grupo[1:12]

ggplot(
  sim |> filter(grupo %in% seleccion),
  aes(x = x, y = y)
) +
  geom_point(alpha = 0.45) +
  geom_abline(
    data = parametros_grupo |> filter(grupo %in% seleccion),
    aes(intercept = alpha_real, slope = beta_real),
    linewidth = 0.8
  ) +
  facet_wrap(~ grupo, ncol = 4) +
  labs(
    x = "Predictor centrado x",
    y = "Respuesta y"
  ) +
  theme_minimal()

Figura 4.2: Datos simulados en una selección de grupos. Por construcción, tanto los interceptos como las pendientes cambian entre grupos.

Las rectas no son paralelas. El modelo de pendiente común sería deliberadamente incorrecto para este proceso generador.

4.18.3 Comprobación predictiva previa por simulación directa

Para la escala simulada propondremos:

\[ \mu_\alpha\sim\mathcal N(50,15^2), \]

\[ \mu_\beta\sim\mathcal N(0,6^2), \]

\[ \tau_\alpha\sim Half\text{-}\mathcal N(0,12^2), \]

\[ \tau_\beta\sim Half\text{-}\mathcal N(0,4^2), \]

\[ \sigma\sim Half\text{-}\mathcal N(0,10^2), \]

\[ R\sim LKJ(2). \]

En el caso bidimensional, si \(R\sim LKJ(\eta)\) puede simularse la correlación mediante

\[ \frac{\rho+1}{2}\sim Beta(\eta,\eta). \]

Usaremos esta identidad únicamente para la comprobación previa \(2\times2\).

simular_familia_prior <- function(draw_id, J = 12) {
  mu_alpha <- rnorm(1, 50, 15)
  mu_beta <- rnorm(1, 0, 6)
  tau_alpha <- abs(rnorm(1, 0, 12))
  tau_beta <- abs(rnorm(1, 0, 4))
  rho <- 2 * rbeta(1, 2, 2) - 1

  B <- sim_bivnorm(
    n = J,
    mu = c(mu_alpha, mu_beta),
    sd = c(tau_alpha, tau_beta),
    rho = rho
  )

  tibble(
    draw = factor(draw_id),
    grupo = factor(seq_len(J)),
    alpha = B[, 1],
    beta = B[, 2],
    rho = rho
  )
}

prior_lineas <- purrr::map_dfr(
  seq_len(9),
  simular_familia_prior
)

x_grid <- seq(-2.5, 2.5, length.out = 60)

prior_lineas |>
  tidyr::crossing(x = x_grid) |>
  mutate(mu = alpha + beta * x) |>
  ggplot(aes(x = x, y = mu, group = grupo)) +
  geom_line(alpha = 0.6) +
  facet_wrap(~ draw, ncol = 3) +
  labs(
    x = "x",
    y = "Media condicional previa"
  ) +
  theme_minimal()

Figura 4.3: Nueve familias de líneas simuladas desde las previas. La comprobación predictiva previa permite evaluar simultáneamente niveles, pendientes, heterogeneidad y correlación antes de usar los datos.

Este gráfico debe leerse como una pregunta científica: ¿son plausibles familias de rectas de esta magnitud en la escala del problema? Si la respuesta fuera no, deberíamos modificar las previas antes de observar la posterior.

4.18.4 Ajustar pendiente común y pendientes variables

Primero ajustamos un modelo que comparte una única pendiente:

\[ \mu_{ij}=\alpha_j+\beta x_{ij}. \]

prior_sim_fija <- c(
  prior(normal(50, 15), class = "Intercept"),
  prior(normal(0, 6), class = "b", coef = "x"),
  prior(normal(0, 12), class = "sd", group = "grupo"),
  prior(normal(0, 10), class = "sigma")
)

fit_sim_fija <- brm(
  y ~ 1 + x + (1 | grupo),
  data = sim,
  family = gaussian(),
  prior = prior_sim_fija,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  seed = 1653,
  backend = "cmdstanr",
  control = list(adapt_delta = 0.95),
  file = "_fits/semana04_sim_pendiente_fija",
  file_refit = "on_change",
  refresh = 0
)

Running MCMC with 4 parallel chains...

Chain 4 finished in 2.5 seconds.
Chain 1 finished in 2.8 seconds.
Chain 2 finished in 2.8 seconds.
Chain 3 finished in 3.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.8 seconds.
Total execution time: 3.3 seconds.

Luego permitimos que la pendiente cambie entre grupos:

prior_sim_var <- c(
  prior(normal(50, 15), class = "Intercept"),
  prior(normal(0, 6), class = "b", coef = "x"),
  prior(
    normal(0, 12),
    class = "sd",
    group = "grupo",
    coef = "Intercept"
  ),
  prior(
    normal(0, 4),
    class = "sd",
    group = "grupo",
    coef = "x"
  ),
  prior(lkj(2), class = "cor", group = "grupo"),
  prior(normal(0, 10), class = "sigma")
)

fit_sim_var <- brm(
  y ~ 1 + x + (1 + x | grupo),
  data = sim,
  family = gaussian(),
  prior = prior_sim_var,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  seed = 1653,
  backend = "cmdstanr",
  control = list(adapt_delta = 0.95),
  file = "_fits/semana04_sim_pendiente_variable",
  file_refit = "on_change",
  refresh = 0
)

Running MCMC with 4 parallel chains...

Chain 4 finished in 6.7 seconds.
Chain 2 finished in 6.9 seconds.
Chain 3 finished in 6.8 seconds.
Chain 1 finished in 7.4 seconds.

All 4 chains finished successfully.
Mean chain execution time: 7.0 seconds.
Total execution time: 7.5 seconds.

La sintaxis

(1 + x | grupo)

indica que brms modelará conjuntamente un intercepto y una pendiente de x por grupo, incluyendo su correlación. Bürkner utiliza la misma estructura de fórmula para representar múltiples coeficientes grupales (Bürkner 2017, 6-7).

4.18.5 Resumen de los parámetros poblacionales

draws_sim <- posterior::as_draws_df(fit_sim_var)

posterior::summarise_draws(
  draws_sim,
  "mean",
  "sd",
  "median",
  ~ quantile(.x, 0.05),
  ~ quantile(.x, 0.95),
  "rhat",
  "ess_bulk",
  "ess_tail"
) |>
  filter(
    variable %in% c(
      "b_Intercept",
      "b_x",
      "sd_grupo__Intercept",
      "sd_grupo__x",
      "cor_grupo__Intercept__x",
      "sigma"
    )
  )

# A tibble: 6 × 9
  variable               mean    sd median   `5%`  `95%`  rhat ess_bulk ess_tail
  <chr>                 <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 b_Intercept          50.6   1.52  50.7   48.0   53.0    1.01     348.     466.
2 b_x                   2.63  0.587  2.62   1.68   3.61   1.00     699.    1087.
3 sd_grupo__Intercept   8.10  1.10   8.00   6.47  10.0    1.00     550.    1040.
4 sd_grupo__x           2.66  0.479  2.62   1.94   3.52   1.00    1271.    2152.
5 cor_grupo__Intercep… -0.534 0.175 -0.549 -0.791 -0.211  1.00    1171.    1613.
6 sigma                 5.69  0.162  5.69   5.43   5.97   1.00    3857.    2798.

Conviene comparar las posteriores con los valores generadores, pero sin esperar recuperación exacta en una sola muestra simulada. La simulación sirve para estudiar qué cantidades están bien informadas y cuáles conservan incertidumbre apreciable.

4.18.6 Shrinkage simultáneo de interceptos y pendientes

Para visualizar el pooling, estimaremos primero una regresión separada en cada grupo. Estas estimaciones no comparten información.

coef_unpooled <- sim |>
  group_by(grupo) |>
  group_modify(
    ~ {
      m <- lm(y ~ x, data = .x)
      tibble(
        alpha_unpooled = unname(coef(m)[1]),
        beta_unpooled = unname(coef(m)[2])
      )
    }
  ) |>
  ungroup()

# coef() es el genérico S3 de stats; brms aporta el método coef.brmsfit.
# Con summary = TRUE, el arreglo tiene dimensiones:
# nivel del grupo x estadístico resumen x coeficiente.
coef_pooled_array <- stats::coef(fit_sim_var, summary = TRUE)$grupo

coef_pooled <- tibble::tibble(
  grupo = factor(
    dimnames(coef_pooled_array)[[1]],
    levels = niveles_grupo
  ),
  alpha_pooled = coef_pooled_array[, "Estimate", "Intercept"],
  beta_pooled = coef_pooled_array[, "Estimate", "x"]
)

shrinkage_sim <- coef_unpooled |>
  left_join(coef_pooled, by = "grupo") |>
  left_join(
    info_grupo |> select(grupo, n, Sxx),
    by = "grupo"
  )

ggplot(shrinkage_sim) +
  geom_segment(
    aes(
      x = alpha_unpooled,
      y = beta_unpooled,
      xend = alpha_pooled,
      yend = beta_pooled
    ),
    alpha = 0.55,
    arrow = grid::arrow(length = grid::unit(0.08, "inches"))
  ) +
  geom_point(
    aes(x = alpha_unpooled, y = beta_unpooled),
    shape = 1,
    size = 2.2
  ) +
  geom_point(
    aes(x = alpha_pooled, y = beta_pooled, size = Sxx),
    alpha = 0.75
  ) +
  labs(
    x = "Intercepto",
    y = "Pendiente",
    size = expression(S[xx]),
    title = "De estimaciones separadas a pooling parcial"
  ) +
  theme_minimal()

Figura 4.4: Shrinkage bidimensional. Cada segmento conecta la estimación separada de un grupo con la estimación parcialmente agrupada del modelo multinivel. La dirección del shrinkage puede no ser horizontal ni vertical cuando interceptos y pendientes están correlacionados.

Los grupos con poca información sobre la pendiente pueden mostrar una contracción importante. El tamaño muestral es parte de esa información, pero \(S_{xx,j}\) deja claro que también importa la cobertura del predictor.

4.18.7 La correlación cambia al mover el origen

No necesitamos volver a ajustar el modelo para estudiar qué ocurriría al redefinir el origen del predictor. Podemos transformar cada draw posterior de \((\tau_\alpha,\tau_\beta,\rho)\) usando las fórmulas de Sección 4.5, sin cambiar las rectas ajustadas.

rho_transformado <- draws_sim |>
  transmute(
    tau_alpha = sd_grupo__Intercept,
    tau_beta = sd_grupo__x,
    rho = cor_grupo__Intercept__x
  ) |>
  tidyr::crossing(c = c(-2, -1, 0, 1, 2)) |>
  mutate(
    cov_star = rho * tau_alpha * tau_beta + c * tau_beta^2,
    tau_alpha_star = sqrt(
      tau_alpha^2 +
        2 * c * rho * tau_alpha * tau_beta +
        c^2 * tau_beta^2
    ),
    rho_star = cov_star / (tau_alpha_star * tau_beta)
  )

ggplot(
  rho_transformado,
  aes(x = rho_star, group = factor(c))
) +
  geom_density(aes(linetype = factor(c)), linewidth = 0.8) +
  labs(
    x = expression(rho^"*"),
    y = "Densidad posterior",
    linetype = "Desplazamiento c"
  ) +
  theme_minimal()

Figura 4.5: Posterior de la correlación intercepto–pendiente bajo distintos desplazamientos del origen de x. Las rectas predichas son las mismas; cambia la parametrización de los coeficientes.

El mensaje del gráfico es más importante que un valor particular: \(\rho\) debe interpretarse junto con la definición del intercepto.

4.18.8 Comprobación predictiva de la heterogeneidad de pendientes

Una superposición global de densidades puede ser casi insensible a que las líneas de los grupos tengan pendientes distintas. Construiremos un estadístico de discrepancia dirigido a la estructura que motivó el modelo.

Como \(x\) tiene media cero dentro de cada grupo, la pendiente OLS descriptiva de un grupo puede calcularse mediante

\[ \widehat b_j = \frac{\sum_i x_{ij}y_{ij}} {\sum_i x_{ij}^2}. \]

Usaremos como estadístico

\[ T(y) = SD_j(\widehat b_j), \]

la dispersión de esas pendientes descriptivas. No la interpretamos como un estimador directo de \(\tau_\beta\), porque contiene ruido muestral. Su función es únicamente comparar el patrón observado con patrones replicados bajo el modelo.

sd_pendientes_observadas <- function(y, x, grupo) {
  indices <- split(seq_along(y), grupo)

  pendientes <- vapply(
    indices,
    function(idx) {
      sum(x[idx] * y[idx]) / sum(x[idx]^2)
    },
    numeric(1)
  )

  sd(pendientes)
}

T_obs <- sd_pendientes_observadas(
  y = sim$y,
  x = sim$x,
  grupo = sim$grupo
)

Generamos réplicas predictivas desde ambos modelos.

set.seed(1653)

yrep_fija <- posterior_predict(fit_sim_fija, ndraws = 80)
yrep_var <- posterior_predict(fit_sim_var, ndraws = 80)

T_fija <- apply(
  yrep_fija,
  1,
  sd_pendientes_observadas,
  x = sim$x,
  grupo = sim$grupo
)

T_var <- apply(
  yrep_var,
  1,
  sd_pendientes_observadas,
  x = sim$x,
  grupo = sim$grupo
)

ppc_T <- tibble(
  modelo = rep(
    c("Pendiente común", "Pendientes variables"),
    each = 80
  ),
  T_rep = c(T_fija, T_var)
)

ggplot(ppc_T, aes(x = T_rep)) +
  geom_density() +
  geom_vline(xintercept = T_obs, linetype = 2) +
  facet_wrap(~ modelo, ncol = 1) +
  labs(
    x = "SD entre pendientes descriptivas por grupo",
    y = "Densidad predictiva posterior"
  ) +
  theme_minimal()

Figura 4.6: Comprobación predictiva dirigida a la heterogeneidad de pendientes. La línea vertical representa la dispersión de pendientes descriptivas observada; las densidades representan réplicas posteriores bajo cada modelo.

Este ejercicio muestra por qué la comprobación predictiva debe diseñarse para el aspecto del modelo que nos interesa. Un modelo puede reproducir razonablemente la distribución marginal de \(y\) y fallar en la dispersión de asociaciones entre grupos.

4.18.9 Diagnóstico MCMC básico del modelo simulado

Antes de interpretar la posterior, revisamos al menos \(\widehat R\), ESS, trazas y divergencias.

posterior::summarise_draws(
  draws_sim,
  "mean",
  "sd",
  "rhat",
  "ess_bulk",
  "ess_tail"
) |>
  filter(
    variable %in% c(
      "b_Intercept",
      "b_x",
      "sd_grupo__Intercept",
      "sd_grupo__x",
      "cor_grupo__Intercept__x",
      "sigma"
    )
  )

# A tibble: 6 × 6
  variable                  mean    sd  rhat ess_bulk ess_tail
  <chr>                    <dbl> <dbl> <dbl>    <dbl>    <dbl>
1 b_Intercept             50.6   1.52   1.01     348.     466.
2 b_x                      2.63  0.587  1.00     699.    1087.
3 sd_grupo__Intercept      8.10  1.10   1.00     550.    1040.
4 sd_grupo__x              2.66  0.479  1.00    1271.    2152.
5 cor_grupo__Intercept__x -0.534 0.175  1.00    1171.    1613.
6 sigma                    5.69  0.162  1.00    3857.    2798.

bayesplot::mcmc_trace(
  as.array(fit_sim_var),
  pars = c(
    "b_x",
    "sd_grupo__Intercept",
    "sd_grupo__x",
    "cor_grupo__Intercept__x",
    "sigma"
  )
)

Trazas de parámetros poblacionales y de covarianza seleccionados del modelo simulado.

np_sim <- brms::nuts_params(fit_sim_var)

tibble(
  divergencias_post_warmup = sum(
    np_sim$Parameter == "divergent__" &
      np_sim$Value == 1
  ),
  alcanzan_treedepth_10 = sum(
    np_sim$Parameter == "treedepth__" &
      np_sim$Value >= 10
  )
)

# A tibble: 1 × 2
  divergencias_post_warmup alcanzan_treedepth_10
                     <int>                 <int>
1                        0                     0

No intentaremos solucionar aquí todos los posibles problemas de geometría. Si aparecen divergencias o ESS muy bajos, deben tratarse como una señal para revisar parametrización, escala, previas e identificación. La parametrización centrada/no centrada y los diagnósticos de HMC se desarrollarán explícitamente en la semana 7. McElreath y Bayesian Workflow discuten que las parametrizaciones matemáticamente equivalentes pueden comportarse de manera muy distinta en HMC, especialmente en modelos jerárquicos complejos (McElreath 2020, 447-54; Gelman et al. 2026, 209-35, 281-83).

4.19 Aplicación: heterogeneidad del gradiente socioeconómico entre escuelas

Retomamos MathAchieve del paquete nlme, usado en la semana 3. Los datos contienen estudiantes agrupados en escuelas y permiten mantener la misma interpretación dentro–entre mientras ampliamos la estructura de pendientes.

Trabajaremos con:

• MathAch: rendimiento matemático;
• SES: nivel socioeconómico individual;
• MEANSES: SES promedio de la escuela;
• School: identificador de escuela.

Hox et al. utilizan estos datos para mostrar que el centrado por grupo separa la asociación dentro de escuelas de la comparación entre escuelas (Hox et al. 2018, 51-52). Esta semana preguntaremos además si la asociación dentro de las escuelas cambia entre ellas.

4.19.1 Preparación de los datos

data("MathAchieve", package = "nlme")

math <- MathAchieve |>
  as_tibble() |>
  mutate(
    School = factor(School)
  )

stopifnot(
  !anyNA(math[c("School", "SES", "MathAch", "MEANSES")]),
  nlevels(math$School) == 160
)

grand_mean_ses <- mean(math$SES)

math <- math |>
  mutate(
    ses_wc = SES - MEANSES,
    meanses_gmc = MEANSES - grand_mean_ses
  )

math |>
  summarise(
    N = n(),
    escuelas = n_distinct(School),
    media_math = mean(MathAch),
    sd_math = sd(MathAch),
    media_ses_wc = mean(ses_wc),
    sd_ses_wc = sd(ses_wc),
    sd_meanses = sd(MEANSES)
  )

# A tibble: 1 × 7
      N escuelas media_math sd_math media_ses_wc sd_ses_wc sd_meanses
  <int>    <int>      <dbl>   <dbl>        <dbl>     <dbl>      <dbl>
1  7185      160       12.7    6.88     -0.00600     0.661      0.414

La variable ses_wc contiene únicamente variación dentro de escuela. Así, permitir que su pendiente varíe responde directamente a la pregunta de si el gradiente socioeconómico intragrupo cambia entre escuelas.

4.19.2 ¿Hay evidencia descriptiva de pendientes diferentes?

No interpretaremos regresiones separadas como estimaciones finales, pero sí pueden servir para explorar la estructura.

pendientes_math_desc <- math |>
  group_by(School) |>
  group_modify(
    ~ {
      m <- lm(MathAch ~ ses_wc, data = .x)
      tibble(
        n = nrow(.x),
        Sxx = sum(.x$ses_wc^2),
        pendiente = unname(coef(m)[2])
      )
    }
  ) |>
  ungroup()

pendientes_math_desc |>
  summarise(
    mediana_pendiente = median(pendiente),
    sd_pendientes = sd(pendiente),
    q10 = quantile(pendiente, 0.10),
    q90 = quantile(pendiente, 0.90)
  )

# A tibble: 1 × 4
  mediana_pendiente sd_pendientes    q10   q90
              <dbl>         <dbl>  <dbl> <dbl>
1              2.31          1.63 0.0847  4.11

ggplot(
  pendientes_math_desc,
  aes(x = Sxx, y = pendiente)
) +
  geom_point(alpha = 0.55) +
  geom_hline(
    yintercept = median(pendientes_math_desc$pendiente),
    linetype = 2
  ) +
  scale_x_log10() +
  labs(
    x = expression(S[xx] ~ "dentro de escuela (escala log)"),
    y = "Pendiente descriptiva separada"
  ) +
  theme_minimal()

Figura 4.7: Pendientes OLS separadas por escuela contra la información Sxx disponible para estimarlas. La gran dispersión en escuelas con poca información ilustra por qué no conviene interpretar estas estimaciones sin pooling.

La variación de estas pendientes combina heterogeneidad real y error de estimación. Precisamente por eso necesitamos un modelo jerárquico.

4.19.3 Modelo 1: pendiente intragrupo variable

Partimos del modelo de la semana 3 y permitimos que la asociación de ses_wc cambie entre escuelas:

\[ MathAch_{ij} \sim \mathcal N(\mu_{ij},\sigma^2), \]

\[ \mu_{ij} = \alpha_j +\beta_j SES^W_{ij} +\beta_B \overline{SES}^{B}_j, \]

\[ \begin{pmatrix} \alpha_j\\ \beta_j \end{pmatrix} \sim \mathcal N_2 \left[ \begin{pmatrix} \mu_\alpha\\ \mu_\beta \end{pmatrix}, \begin{pmatrix} \tau_\alpha^2 & \rho\tau_\alpha\tau_\beta\\ \rho\tau_\alpha\tau_\beta & \tau_\beta^2 \end{pmatrix} \right]. \]

Interpretamos:

• \(\mu_\beta\): asociación socioeconómica intragrupo promedio entre las escuelas;
• \(\tau_\beta\): heterogeneidad residual entre las pendientes intragrupo;
• \(\rho\): asociación entre el rendimiento esperado de un estudiante situado en la media de SES de su escuela y el gradiente socioeconómico dentro de esa escuela.

Esta última interpretación depende de nuestro centrado: como ses_wc = 0 corresponde al SES medio de la escuela, el intercepto está definido justamente en ese punto.

4.19.3.1 Previas para la aplicación

Usaremos previas explícitas en la escala de MathAch:

\[ \mu_\alpha\sim\mathcal N(12,10^2), \]

\[ \mu_\beta,\beta_B\sim\mathcal N(0,5^2), \]

\[ \tau_\alpha\sim Half\text{-}\mathcal N(0,10^2), \]

\[ \tau_\beta\sim Half\text{-}\mathcal N(0,4^2), \]

\[ R\sim LKJ(2), \]

\[ \sigma\sim Half\text{-}\mathcal N(0,10^2). \]

La menor escala de la previa sobre \(\tau_\beta\) refleja que una desviación estándar de varios puntos de rendimiento por una unidad de SES ya implica heterogeneidad considerable entre gradientes. Esta elección es pedagógica y debe revisarse mediante simulación previa en una aplicación sustantiva real.

prior_math_var <- c(
  prior(normal(12, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b", coef = "ses_wc"),
  prior(normal(0, 5), class = "b", coef = "meanses_gmc"),
  prior(
    normal(0, 10),
    class = "sd",
    group = "School",
    coef = "Intercept"
  ),
  prior(
    normal(0, 4),
    class = "sd",
    group = "School",
    coef = "ses_wc"
  ),
  prior(lkj(2), class = "cor", group = "School"),
  prior(normal(0, 10), class = "sigma")
)

fit_math_var <- brm(
  MathAch ~ 1 + ses_wc + meanses_gmc +
    (1 + ses_wc | School),
  data = math,
  family = gaussian(),
  prior = prior_math_var,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  seed = 1653,
  backend = "cmdstanr",
  control = list(adapt_delta = 0.95),
  file = "_fits/semana04_math_pendiente_variable",
  file_refit = "on_change",
  refresh = 0
)

Running MCMC with 4 parallel chains...

Chain 1 finished in 91.9 seconds.
Chain 2 finished in 96.9 seconds.
Chain 4 finished in 97.2 seconds.
Chain 3 finished in 99.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 96.3 seconds.
Total execution time: 99.4 seconds.

4.19.3.2 Resumen de la heterogeneidad de pendientes

draws_math_var <- posterior::as_draws_df(fit_math_var)

posterior::summarise_draws(
  draws_math_var,
  "mean",
  "sd",
  "median",
  ~ quantile(.x, 0.05),
  ~ quantile(.x, 0.95),
  "rhat",
  "ess_bulk",
  "ess_tail"
) |>
  filter(
    variable %in% c(
      "b_Intercept",
      "b_ses_wc",
      "b_meanses_gmc",
      "sd_School__Intercept",
      "sd_School__ses_wc",
      "cor_School__Intercept__ses_wc",
      "sigma"
    )
  )

# A tibble: 7 × 9
  variable              mean     sd median   `5%`  `95%`  rhat ess_bulk ess_tail
  <chr>                <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 b_Intercept         12.7   0.150  12.7   12.4   12.9    1.00    2509.    2854.
2 b_ses_wc             2.19  0.130   2.19   1.98   2.40   1.00    6251.    2993.
3 b_meanses_gmc        5.86  0.367   5.87   5.26   6.46   1.00    2439.    2819.
4 sd_School__Interce…  1.66  0.127   1.66   1.46   1.88   1.00    1913.    2428.
5 sd_School__ses_wc    0.792 0.194   0.802  0.452  1.09   1.00    1241.    1242.
6 cor_School__Interc… -0.182 0.185  -0.182 -0.483  0.123  1.00    3639.    2117.
7 sigma                6.06  0.0523  6.06   5.98   6.15   1.00    7033.    2697.

No resumiremos \(\tau_\beta\) mediante una decisión binaria de “cero/no cero”. Nos interesa cuánto gradiente residual permite la posterior, con qué incertidumbre y si esa heterogeneidad ayuda a reproducir los patrones observados.

4.19.3.3 Líneas parcialmente agrupadas por escuela

# El genérico coef() pertenece a stats; el método para objetos brmsfit
# es despachado automáticamente como coef.brmsfit.
coef_math_array <- stats::coef(fit_math_var, summary = TRUE)$School

coef_math_school <- tibble::tibble(
  School = factor(
    dimnames(coef_math_array)[[1]],
    levels = levels(math$School)
  ),
  intercepto = coef_math_array[, "Estimate", "Intercept"],
  pendiente = coef_math_array[, "Estimate", "ses_wc"]
)

escuelas_muestra <- levels(math$School)[1:12]

math |>
  filter(School %in% escuelas_muestra) |>
  ggplot(aes(x = ses_wc, y = MathAch)) +
  geom_point(alpha = 0.30) +
  geom_abline(
    data = coef_math_school |>
      filter(School %in% escuelas_muestra),
    aes(intercept = intercepto, slope = pendiente),
    linewidth = 0.8
  ) +
  facet_wrap(~ School, ncol = 4) +
  labs(
    x = "SES relativo dentro de escuela",
    y = "Rendimiento matemático"
  ) +
  theme_minimal()

Figura 4.8: Rectas posteriores parcialmente agrupadas para una selección de escuelas. Como ses_wc está centrado por escuela, el intercepto corresponde al rendimiento esperado en el SES medio de cada escuela.

4.19.4 Modelo 2: una interacción entre niveles

El modelo anterior reconoce heterogeneidad, pero no intenta explicarla. Ahora preguntamos si el SES promedio de la escuela está asociado con el gradiente de SES individual.

Definimos

\[ z_j = \overline{SES}_j-\overline{SES}_{\cdot\cdot}. \]

Modelamos

\[ \alpha_j = \gamma_{00} +\gamma_{01}z_j +u_{0j}, \]

\[ \beta_j = \gamma_{10} +\gamma_{11}z_j +u_{1j}. \]

La ecuación combinada es

\[ \mu_{ij} = \gamma_{00} +\gamma_{01}z_j +\gamma_{10}SES^W_{ij} +\gamma_{11}SES^W_{ij}z_j +u_{0j} +u_{1j}SES^W_{ij}. \]

Aquí:

• \(\gamma_{10}\) es la pendiente intragrupo esperada en una escuela cuyo SES medio coincide con la media global;
• \(\gamma_{11}\) es el cambio esperado en esa pendiente por una unidad de aumento del SES medio escolar;
• \(\tau_\beta\) es la heterogeneidad residual de pendientes después de incorporar esta moderación sistemática.

La sintaxis de brms es

MathAch ~ ses_wc * meanses_gmc + (1 + ses_wc | School)

El * agrega ambos efectos principales y la interacción. Mantener (1 + ses_wc | School) permite que todavía exista heterogeneidad no explicada por meanses_gmc.

4.19.4.1 Previa del término de interacción

Una unidad de ses_wc * meanses_gmc combina dos escalas de SES. Usaremos

\[ \gamma_{11}\sim\mathcal N(0,3^2), \]

algo más concentrada que las previas de los efectos principales. La elección no es universal: su adecuación depende del rango sustantivo de ambas variables.

prior_math_cross <- c(
  prior(normal(12, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b", coef = "ses_wc"),
  prior(normal(0, 5), class = "b", coef = "meanses_gmc"),
  prior(
    normal(0, 3),
    class = "b",
    coef = "ses_wc:meanses_gmc"
  ),
  prior(
    normal(0, 10),
    class = "sd",
    group = "School",
    coef = "Intercept"
  ),
  prior(
    normal(0, 4),
    class = "sd",
    group = "School",
    coef = "ses_wc"
  ),
  prior(lkj(2), class = "cor", group = "School"),
  prior(normal(0, 10), class = "sigma")
)

fit_math_cross <- brm(
  MathAch ~ 1 + ses_wc * meanses_gmc +
    (1 + ses_wc | School),
  data = math,
  family = gaussian(),
  prior = prior_math_cross,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  seed = 1653,
  backend = "cmdstanr",
  control = list(adapt_delta = 0.95),
  file = "_fits/semana04_math_interaccion_niveles",
  file_refit = "on_change",
  refresh = 0
)

Running MCMC with 4 parallel chains...

Chain 1 finished in 108.9 seconds.
Chain 3 finished in 109.2 seconds.
Chain 2 finished in 111.7 seconds.
Chain 4 finished in 112.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 110.5 seconds.
Total execution time: 112.2 seconds.

4.19.4.2 Posterior de la interacción y de la heterogeneidad residual

draws_math_cross <- posterior::as_draws_df(fit_math_cross)

cantidades_cross <- draws_math_cross |>
  transmute(
    beta_promedio = b_ses_wc,
    beta_entre = b_meanses_gmc,
    beta_interaccion = .data[["b_ses_wc:meanses_gmc"]],
    tau_pendiente = sd_School__ses_wc,
    rho = cor_School__Intercept__ses_wc,
    sigma = sigma
  )

cantidades_cross |>
  pivot_longer(
    cols = everything(),
    names_to = "cantidad",
    values_to = "valor"
  ) |>
  group_by(cantidad) |>
  summarise(
    mediana = median(valor),
    q05 = quantile(valor, 0.05),
    q95 = quantile(valor, 0.95),
    .groups = "drop"
  )

# A tibble: 6 × 4
  cantidad         mediana    q05   q95
  <chr>              <dbl>  <dbl> <dbl>
1 beta_entre         5.85   5.22  6.43 
2 beta_interaccion   0.285 -0.230 0.800
3 beta_promedio      2.19   1.99  2.40 
4 rho               -0.185 -0.509 0.113
5 sigma              6.06   5.98  6.15 
6 tau_pendiente      0.789  0.426 1.08 

4.19.5 ¿Cuánto cambia la heterogeneidad residual?

Una manera descriptiva de estudiar lo que aporta el predictor grupal es observar cómo cambia la posterior de \(\tau_\beta\) entre los dos modelos. No es una prueba de significancia ni un criterio automático de selección.

tau_compare <- bind_rows(
  tibble(
    modelo = "Sin moderador de pendiente",
    tau_beta = draws_math_var$sd_School__ses_wc
  ),
  tibble(
    modelo = "Con interacción entre niveles",
    tau_beta = draws_math_cross$sd_School__ses_wc
  )
)

ggplot(tau_compare, aes(x = tau_beta)) +
  geom_density() +
  facet_wrap(~ modelo, ncol = 1) +
  labs(
    x = expression(tau[beta]),
    y = "Densidad posterior"
  ) +
  theme_minimal()

Figura 4.9: Posterior de la desviación estándar residual entre pendientes antes y después de incluir la interacción entre SES individual relativo y SES medio escolar.

Si la posterior de \(\tau_\beta\) se desplaza hacia valores menores, meanses_gmc está capturando parte de la heterogeneidad que antes quedaba en las pendientes residuales. Si cambia poco, el predictor grupal no parece explicar mucha de esa heterogeneidad bajo esta especificación. En ambos casos puede seguir existiendo incertidumbre considerable.

4.19.6 Visualizar la interacción en la escala de la respuesta

Los coeficientes de interacción suelen entenderse mejor mediante predicciones. Elegimos tres valores del SES medio escolar: percentiles 20, 50 y 80 de la distribución entre escuelas.

z_vals <- math |>
  distinct(School, meanses_gmc) |>
  summarise(
    bajo = quantile(meanses_gmc, 0.20),
    medio = quantile(meanses_gmc, 0.50),
    alto = quantile(meanses_gmc, 0.80)
  ) |>
  pivot_longer(
    everything(),
    names_to = "nivel",
    values_to = "meanses_gmc"
  )

x_vals <- seq(
  quantile(math$ses_wc, 0.05),
  quantile(math$ses_wc, 0.95),
  length.out = 60
)

nd_cross <- tidyr::crossing(
  ses_wc = x_vals,
  z_vals
)

Calculamos la media poblacional, sin una desviación de escuela particular. Esto permite visualizar la parte sistemática de la interacción.

set.seed(1653)

n_draws_cross <- min(1200L, nrow(draws_math_cross))

draws_cross_small <- draws_math_cross |>
  mutate(.draw = row_number()) |>
  slice_sample(n = n_draws_cross) |>
  transmute(
    .draw,
    b0 = b_Intercept,
    bx = b_ses_wc,
    bz = b_meanses_gmc,
    bxz = .data[["b_ses_wc:meanses_gmc"]]
  )

pred_cross <- tidyr::crossing(
  draws_cross_small,
  nd_cross
) |>
  mutate(
    mu = b0 +
      bx * ses_wc +
      bz * meanses_gmc +
      bxz * ses_wc * meanses_gmc
  ) |>
  group_by(nivel, meanses_gmc, ses_wc) |>
  summarise(
    mediana = median(mu),
    q05 = quantile(mu, 0.05),
    q95 = quantile(mu, 0.95),
    .groups = "drop"
  )

ggplot(
  pred_cross,
  aes(x = ses_wc, y = mediana, linetype = nivel)
) +
  geom_ribbon(
    aes(ymin = q05, ymax = q95, group = nivel),
    alpha = 0.12
  ) +
  geom_line(linewidth = 0.9) +
  labs(
    x = "SES relativo dentro de la escuela",
    y = "Rendimiento esperado",
    linetype = "SES medio escolar"
  ) +
  theme_minimal()

Figura 4.10: Predicciones poblacionales de rendimiento según SES relativo dentro de la escuela para tres niveles de SES medio escolar. Las diferencias entre pendientes representan la interacción entre niveles.

No debemos interpretar el gráfico causalmente sin supuestos adicionales. El modelo describe cómo cambia una asociación condicional entre niveles.

4.19.7 Comprobaciones predictivas de la aplicación

4.19.7.1 Comprobación global

pp_check(
  fit_math_cross,
  type = "dens_overlay",
  ndraws = 60
)

Comprobación predictiva posterior global para el modelo con interacción entre niveles.

Una buena reproducción global no garantiza una buena representación de las pendientes entre escuelas. Por eso repetimos el estadístico dirigido usado en la simulación.

sd_pendientes_math <- function(y, x, grupo) {
  indices <- split(seq_along(y), grupo)

  pendientes <- vapply(
    indices,
    function(idx) {
      denom <- sum(x[idx]^2)
      if (denom <= 1e-12) return(NA_real_)
      sum(x[idx] * y[idx]) / denom
    },
    numeric(1)
  )

  sd(pendientes, na.rm = TRUE)
}

T_math_obs <- sd_pendientes_math(
  y = math$MathAch,
  x = math$ses_wc,
  grupo = math$School
)

set.seed(1653)
yrep_math <- posterior_predict(fit_math_cross, ndraws = 60)

T_math_rep <- apply(
  yrep_math,
  1,
  sd_pendientes_math,
  x = math$ses_wc,
  grupo = math$School
)

ggplot(tibble(T_rep = T_math_rep), aes(x = T_rep)) +
  geom_density() +
  geom_vline(xintercept = T_math_obs, linetype = 2) +
  labs(
    x = "SD entre pendientes descriptivas por escuela",
    y = "Densidad predictiva posterior"
  ) +
  theme_minimal()

Figura 4.11: Comprobación predictiva posterior de la dispersión de pendientes descriptivas entre escuelas.

La discrepancia utilizada es una elaboración pedagógica para esta clase. No sustituye otras comprobaciones: convendría examinar también medias por escuela, colas de la respuesta, residuos respecto de SES y posibles no linealidades.

4.19.8 Diagnóstico MCMC de la aplicación

posterior::summarise_draws(
  draws_math_cross,
  "mean",
  "sd",
  "rhat",
  "ess_bulk",
  "ess_tail"
) |>
  filter(
    variable %in% c(
      "b_Intercept",
      "b_ses_wc",
      "b_meanses_gmc",
      "b_ses_wc:meanses_gmc",
      "sd_School__Intercept",
      "sd_School__ses_wc",
      "cor_School__Intercept__ses_wc",
      "sigma"
    )
  )

# A tibble: 8 × 6
  variable                        mean     sd  rhat ess_bulk ess_tail
  <chr>                          <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 b_Intercept                   12.7   0.150  1.00     2427.    2612.
2 b_ses_wc                       2.19  0.126  1.000    6817.    3544.
3 b_meanses_gmc                  5.84  0.367  1.00     2279.    2830.
4 b_ses_wc:meanses_gmc           0.291 0.314  1.00     6253.    3419.
5 sd_School__Intercept           1.66  0.124  1.00     1718.    2362.
6 sd_School__ses_wc              0.775 0.199  1.00     1180.    1233.
7 cor_School__Intercept__ses_wc -0.187 0.191  1.00     3137.    2127.
8 sigma                          6.06  0.0521 1.00     7295.    2809.

np_math_cross <- brms::nuts_params(fit_math_cross)

tibble(
  divergencias_post_warmup = sum(
    np_math_cross$Parameter == "divergent__" &
      np_math_cross$Value == 1
  ),
  alcanzan_treedepth_10 = sum(
    np_math_cross$Parameter == "treedepth__" &
      np_math_cross$Value >= 10
  )
)

# A tibble: 1 × 2
  divergencias_post_warmup alcanzan_treedepth_10
                     <int>                 <int>
1                        0                     0

bayesplot::mcmc_trace(
  as.array(fit_math_cross),
  pars = c(
    "b_ses_wc",
    "b_meanses_gmc",
    "b_ses_wc:meanses_gmc",
    "sd_School__ses_wc",
    "cor_School__Intercept__ses_wc",
    "sigma"
  )
)

Figura 4.12: Trazas MCMC para parámetros seleccionados del modelo con interacción entre niveles.

Una posterior amplia de cor_School__Intercept__ses_wc no debe “arreglarse” aumentando adapt_delta: si el algoritmo muestrea adecuadamente y la incertidumbre persiste, el problema puede ser simplemente que la correlación está poco identificada por los datos. En cambio, divergencias requieren atención computacional.

4.19.9 Predicción para escuelas existentes y nuevas

La distinción de la semana 2 se vuelve más rica. Para una escuela ya observada, la posterior de su intercepto y su pendiente se ha actualizado con sus propios datos. Para una escuela nueva, debemos generar un par nuevo de coeficientes de la distribución poblacional.

Condicional en los hiperparámetros,

\[ \begin{pmatrix} \alpha_{new}\\ \beta_{new} \end{pmatrix} \sim \mathcal N_2 \left( \begin{pmatrix} \gamma_{00}+\gamma_{01}z_{new}\\ \gamma_{10}+\gamma_{11}z_{new} \end{pmatrix}, \Sigma \right). \]

Por tanto, una predicción para una escuela nueva debe propagar incertidumbre sobre:

• los coeficientes poblacionales;
• \(\tau_\alpha\) y \(\tau_\beta\);
• \(\rho\);
• la realización nueva de intercepto y pendiente;
• y, para una observación futura, también \(\sigma\).

brms puede simular niveles nuevos mediante allow_new_levels = TRUE. Usaremos posterior_epred() para mostrar incertidumbre en la media esperada, sin agregar todavía error residual individual.

escuela_existente <- levels(math$School)[1]
z_existente <- math |>
  filter(School == escuela_existente) |>
  summarise(z = first(meanses_gmc)) |>
  pull(z)

x_pred <- seq(-1.5, 1.5, length.out = 50)

new_existente <- tibble(
  ses_wc = x_pred,
  meanses_gmc = z_existente,
  School = escuela_existente
)

new_nueva <- tibble(
  ses_wc = x_pred,
  meanses_gmc = 0,
  School = "Escuela_nueva"
)

mu_existente <- posterior_epred(
  fit_math_cross,
  newdata = new_existente,
  re_formula = NULL
)

mu_nueva <- posterior_epred(
  fit_math_cross,
  newdata = new_nueva,
  re_formula = NULL,
  allow_new_levels = TRUE,
  sample_new_levels = "gaussian"
)

resumir_epred <- function(mat, x, tipo) {
  tibble(
    ses_wc = x,
    mediana = apply(mat, 2, median),
    q05 = apply(mat, 2, quantile, probs = 0.05),
    q95 = apply(mat, 2, quantile, probs = 0.95),
    tipo = tipo
  )
}

pred_grupos <- bind_rows(
  resumir_epred(
    mu_existente,
    x_pred,
    "Escuela observada"
  ),
  resumir_epred(
    mu_nueva,
    x_pred,
    "Escuela nueva"
  )
)

ggplot(
  pred_grupos,
  aes(x = ses_wc, y = mediana)
) +
  geom_ribbon(
    aes(ymin = q05, ymax = q95),
    alpha = 0.15
  ) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~ tipo, ncol = 1) +
  labs(
    x = "SES relativo dentro de escuela",
    y = "Media esperada de MathAch"
  ) +
  theme_minimal()

Figura 4.13: Predicción de la media esperada para una escuela observada y para una escuela nueva. La escuela nueva requiere simular un nuevo intercepto y una nueva pendiente de la distribución poblacional.

El intervalo para un grupo nuevo puede ser mayor porque no contamos con observaciones de esa escuela para actualizar su par \((\alpha,\beta)\). McElreath destaca esta diferencia entre predicciones para clusters existentes y nuevos en su discusión de predicción multinivel (McElreath 2020, 426-31).

4.20 ¿Cuál modelo debemos preferir?

Durante esta semana hemos construido una secuencia de expansiones:

1. pendiente común;
2. pendiente variable;
3. pendiente variable con un predictor grupal que explica parte de su heterogeneidad.

No las trataremos como una competencia donde un único número elige automáticamente el “mejor” modelo.

Cada expansión debe evaluarse preguntando:

• ¿qué afirmación sustantiva nueva representa?;
• ¿qué cantidad podemos interpretar ahora que antes no existía?;
• ¿la posterior de los nuevos parámetros es suficientemente informativa para la pregunta?;
• ¿el modelo reproduce mejor el aspecto de los datos que motivó la expansión?;
• ¿la complejidad adicional genera problemas de identificación o de cómputo?;
• ¿las conclusiones principales son sensibles a las previas razonables?

La comparación predictiva formal mediante PSIS-LOO se desarrollará en la semana 14. En esta etapa damos prioridad a construcción generativa, comprobaciones predictivas dirigidas e interpretación.

4.21 Perspectiva frecuentista y terminología

En la literatura clásica es común encontrar expresiones como:

• random intercept;
• random slope;
• random coefficient model;
• fixed effect;
• cross-level interaction.

En estas notas preferimos intercepto variable, pendiente variable y coeficiente poblacional, siguiendo la recomendación de hacer explícito que los coeficientes grupales son cantidades modeladas y que podemos describir su incertidumbre (Bürkner 2017, 2-3).

En lme4, una formulación frecuentista conceptualmente paralela sería

lmer(
  MathAch ~ ses_wc * meanses_gmc +
    (1 + ses_wc | School),
  data = math
)

La estructura del predictor lineal y de la matriz de covarianza es muy similar. Cambian el marco inferencial, la forma de incorporar información previa, la representación de incertidumbre y los procedimientos de ajuste.

El objetivo del curso no es organizar la modelación alrededor de pruebas secuenciales de componentes aleatorios. La pregunta de si una pendiente debe variar debe comenzar por el mecanismo científico y el diseño, y luego estudiarse mediante inferencia y comprobación del modelo.

4.22 Errores frecuentes de interpretación

4.22.1 1. Interpretar \(\mu_\beta\) como la pendiente de todos los grupos

\(\mu_\beta\) es el centro de la distribución poblacional de pendientes. Si \(\tau_\beta>0\), las pendientes específicas son

\[ \beta_j=\mu_\beta+u_{1j} \]

y no coinciden exactamente con \(\mu_\beta\).

4.22.2 2. Confundir incertidumbre de \(\beta_j\) con heterogeneidad entre \(\beta_j\)

Un intervalo amplio para una escuela puede deberse a poca información en esa escuela. \(\tau_\beta\) describe la distribución entre pendientes verdaderas bajo el modelo, no la incertidumbre de una escuela particular.

4.22.3 3. Convertir automáticamente todas las pendientes en variables

Cada pendiente variable agrega una nueva desviación estándar grupal y nuevas correlaciones. La estructura debe responder a un mecanismo plausible y a información suficiente en el diseño.

4.22.4 4. Fijar el intercepto sin reconocer la restricción

Una pendiente variable con intercepto común obliga a todas las rectas a coincidir en \(x=0\). Esa igualdad necesita una justificación.

4.22.5 5. Interpretar \(\rho\) sin considerar el centrado

La correlación intercepto–pendiente cambia al mover el origen del predictor. Un valor extremo no debe interpretarse fuera de la parametrización que lo define.

4.22.6 6. Tratar una interacción entre niveles como si explicara toda la heterogeneidad

El término \(\gamma_{11}z_j\) describe heterogeneidad sistemática. Si mantenemos \(u_{1j}\), todavía queda variación residual de pendientes.

4.22.7 7. Eliminar la pendiente variable porque la interacción parece grande

Un predictor grupal puede asociarse fuertemente con la pendiente y aun así dejar heterogeneidad residual importante.

4.22.8 8. Interpretar el coeficiente de interacción de manera aislada

\(\gamma_{11}\) es un cambio de pendiente. Debe interpretarse junto con la escala y el origen de \(x\) y \(z\), idealmente mediante predicciones.

4.22.9 9. Pensar que muchos datos totales garantizan buena información sobre \(\tau_\beta\)

Podemos tener miles de observaciones, pero si hay pocos grupos o escasa variación de \(x\) dentro de ellos, la distribución de pendientes puede seguir estando débilmente informada.

4.22.10 10. Confundir una posterior amplia con una falla computacional

Si las cadenas mezclan bien, no hay divergencias y el ESS es adecuado, una posterior amplia puede ser el resultado correcto de información limitada.

4.22.11 11. Usar una comprobación predictiva global para evaluar heterogeneidad de pendientes

La densidad marginal de \(y\) puede verse bien incluso si el modelo representa mal las relaciones dentro de grupo. La comprobación debe ser sensible a las pendientes.

4.22.12 12. Interpretar asociaciones observacionales como moderación causal

Una interacción estadística entre SES individual y composición escolar no demuestra que cambiar la composición de la escuela causaría un cambio en el gradiente individual.

4.23 Síntesis

El modelo central de la semana es

\[ y_{ij} \sim \mathcal N(\alpha_j+\beta_jx_{ij},\sigma^2), \]

\[ \begin{pmatrix} \alpha_j\\ \beta_j \end{pmatrix} \sim \mathcal N_2 \left[ \begin{pmatrix} \mu_\alpha\\ \mu_\beta \end{pmatrix}, \begin{pmatrix} \tau_\alpha^2 & \rho\tau_\alpha\tau_\beta\\ \rho\tau_\alpha\tau_\beta & \tau_\beta^2 \end{pmatrix} \right]. \]

Su estructura puede entenderse mediante cuatro ideas.

Primera: la heterogeneidad puede afectar relaciones, no solo niveles. \(\tau_\beta\) describe cuánto cambian las pendientes entre grupos.

Segunda: los coeficientes grupales se aprenden conjuntamente. La distribución multivariada permite pooling parcial de interceptos y pendientes, y la correlación puede hacer que información sobre un coeficiente contribuya a regularizar el otro.

Tercera: la parametrización tiene interpretación. El origen de \(x\) define el intercepto y modifica tanto \(\tau_\alpha\) como \(\rho\). Por ello el centrado de la semana 3 es parte integral del modelo de la semana 4.

Cuarta: la heterogeneidad puede modelarse. Si un predictor grupal \(z_j\) explica parte de la variación de pendientes,

\[ \beta_j = \gamma_{10}+\gamma_{11}z_j+u_{1j}, \]

y aparece una interacción entre niveles

\[ \gamma_{11}x_{ij}z_j. \]

Los principios que deben permanecer son:

• una pendiente variable es una distribución de relaciones entre grupos, no una colección de regresiones independientes;
• \(\mu_\beta\) y \(\tau_\beta\) responden preguntas distintas;
• \(\rho\) solo tiene sentido junto con el origen del predictor;
• la variación entre grupos depende de \(x\) cuando las pendientes varían;
• el ICC deja de ser una única constante suficiente para resumir la dependencia;
• una interacción entre niveles explica heterogeneidad sistemática, mientras \(\tau_\beta\) representa heterogeneidad residual;
• la información sobre una pendiente depende de la variación del predictor dentro del grupo, no solo de \(n_j\);
• las previas deben comprobarse en términos de líneas y respuestas plausibles;
• los diagnósticos computacionales y la identificación estadística deben distinguirse;
• las comprobaciones predictivas deben dirigirse a la estructura que motivó la expansión;
• la elección de coeficientes variables debe justificarse por el proceso generador y la pregunta científica.

La próxima semana ampliaremos la estructura de agrupamiento: modelos con más niveles, factores cruzados y ANOVA jerárquico.

4.24 Ejercicios

4.24.1 Ejercicios conceptuales

4.24.1.1 Ejercicio 1. Varianza entre grupos a lo largo del predictor

Considere

\[ y_{ij} = \mu_\alpha+\mu_\beta x_{ij}+u_{0j}+u_{1j}x_{ij}+\varepsilon_{ij}, \]

con

\[ \operatorname{Var}(u_{0j})=\tau_\alpha^2, \qquad \operatorname{Var}(u_{1j})=\tau_\beta^2, \]

\[ \operatorname{Corr}(u_{0j},u_{1j})=\rho. \]

1. Derive \[ V_G(x)=\tau_\alpha^2+2x\rho\tau_\alpha\tau_\beta+x^2\tau_\beta^2. \]
2. Evalúe \(V_G(x)\) para \(x=-2,0,2\) cuando \[ \tau_\alpha=3,\qquad \tau_\beta=1,\qquad \rho=-0.5. \]
3. Explique por qué \(\tau_\alpha\) no debe llamarse simplemente “variabilidad entre grupos” sin mencionar el punto \(x=0\).
4. ¿en qué región del eje esperaría que las líneas estuvieran más próximas?

4.24.1.2 Ejercicio 2. Cambiar el origen cambia \(\rho\)

Suponga

\[ \tau_\alpha=4, \qquad \tau_\beta=1.5, \qquad \rho=-0.7. \]

Defina \(x^*=x-c\).

1. Derive \(\operatorname{Var}(\alpha_j^*)\).
2. Derive \(\operatorname{Cov}(\alpha_j^*,\beta_j)\).
3. Calcule la nueva correlación para \(c=2\).
4. Explique por qué cambiar la correlación no implica cambiar las rectas ajustadas.
5. Dé un ejemplo en que una correlación intercepto–pendiente extrema podría deberse a un origen absurdo del predictor.

4.24.1.3 Ejercicio 3. Interacción, pendiente variable o ambas

Para cada situación, indique si usaría una interacción convencional, una pendiente variable, una interacción entre niveles o una combinación. Justifique.

1. El efecto de una dosis sobre presión arterial puede cambiar con la edad del paciente.
2. La asociación entre horas de estudio y rendimiento puede cambiar entre escuelas, sin un moderador escolar específico propuesto todavía.
3. La asociación entre horas de estudio y rendimiento puede cambiar con el tamaño de la escuela y además quedar heterogeneidad entre escuelas.
4. El efecto de un tratamiento 0/1 cambia entre 25 hospitales.
5. La relación entre temperatura y rendimiento de una parcela cambia entre fincas y se sospecha que la altitud de la finca explica parte de la diferencia.

4.24.1.4 Ejercicio 4. Derivar una interacción entre niveles

Considere

\[ y_{ij}=\alpha_j+\beta_jx_{ij}+\varepsilon_{ij}, \]

\[ \alpha_j=\gamma_{00}+\gamma_{01}z_j+u_{0j}, \]

\[ \beta_j=\gamma_{10}+\gamma_{11}z_j+u_{1j}. \]

1. Sustituya las ecuaciones de segundo nivel en la primera.
2. Identifique el término de interacción.
3. Interprete \(\gamma_{10}\) si \(z\) está centrado en su media.
4. Interprete \(\gamma_{11}\).
5. Interprete \(\tau_\beta\).
6. Explique qué restricción se impondría al eliminar \(u_{1j}\).

4.24.1.5 Ejercicio 5. Información para una pendiente

Dos grupos tienen \(n=30\) observaciones cada uno.

• En A, \(x\) varía aproximadamente entre \(-3\) y \(3\).
• En B, \(x\) varía aproximadamente entre \(-0.2\) y \(0.2\).

Suponga la misma \(\sigma\).

1. ¿qué grupo informará mejor su pendiente?
2. ¿por qué el tamaño de grupo no basta para anticipar la precisión?
3. ¿qué papel desempeña \(S_{xx}\)?
4. explique cómo el pooling parcial afectaría al grupo B.

4.24.1.6 Ejercicio 6. Una pendiente variable sin intercepto variable

Proponga:

1. una situación en la que forzar un intercepto común sea claramente implausible;
2. una situación en la que sí pueda justificarse un intercepto común y una pendiente variable;
3. para cada caso, explique qué representa el punto \(x=0\).

4.24.2 Ejercicios computacionales

4.24.2.1 Ejercicio 7. Cambiar la correlación generadora

Repita la simulación del laboratorio con

\[ \rho=0.7. \]

1. grafique los pares reales \((\alpha_j,\beta_j)\);
2. ajuste el modelo de pendientes variables;
3. grafique el shrinkage bidimensional;
4. compare la dirección de los segmentos de shrinkage con el caso \(\rho=-0.6\);
5. examine la posterior de \(\rho\);
6. repita con \(\rho=0\) y discuta cuánto aprende el modelo con \(J=30\).

4.24.2.2 Ejercicio 8. Poca variación dentro de los grupos

Mantenga \(J\), \(n_j\), \(\tau_\alpha\), \(\tau_\beta\), \(\rho\) y \(\sigma\) del laboratorio, pero cambie la desviación estándar de \(x\) a

\[ SD(x\mid j)=0.15 \]

para todos los grupos.

1. ajuste el mismo modelo;
2. compare la posterior de \(\tau_\beta\) con la del laboratorio original;
3. compare la posterior de \(\rho\);
4. compare ESS y divergencias;
5. distinga qué cambios reflejan falta de información y cuáles, si aparecen, son computacionales.

4.24.2.3 Ejercicio 9. Prior LKJ

Para una matriz \(2\times2\), simule \(10\,000\) valores de \(\rho\) bajo

\[ LKJ(1),\qquad LKJ(2),\qquad LKJ(8). \]

Use

\[ (\rho+1)/2\sim Beta(\eta,\eta). \]

1. grafique las tres densidades;
2. calcule \(P(|\rho|>0.8)\) bajo cada previa;
3. explique el significado de aumentar \(\eta\);
4. discuta por qué LKJ(8) podría ser demasiado restrictiva en una aplicación con una razón sustantiva para esperar correlaciones fuertes.

4.24.2.4 Ejercicio 10. Comprobación predictiva dirigida

Use los datos simulados de la clase.

1. proponga un estadístico predictivo diferente de la desviación estándar de pendientes OLS;
2. explique qué característica del modelo examina;
3. calcúlelo en los datos observados y en 100 réplicas posteriores;
4. compare pendiente común y pendiente variable;
5. explique por qué una comprobación global de densidades podría no detectar la misma discrepancia.

4.24.2.5 Ejercicio 11. Moderador escolar alternativo

Cargue MathAchSchool del paquete nlme y agregue Sector a los datos de estudiantes mediante School.

1. describa el nivel de medición de Sector;
2. proponga una codificación que dé una interpretación clara a la pendiente principal de ses_wc;
3. ajuste un modelo donde Sector modere la pendiente de ses_wc y mantenga una pendiente residual variable por escuela;
4. grafique las pendientes poblacionales por sector;
5. compare la posterior de \(\tau_\beta\) con la del modelo sin Sector;
6. interprete la interacción como asociación, no como efecto causal, salvo que pueda defender supuestos adicionales.

4.24.2.6 Ejercicio 12. Predicción para grupos nuevos

A partir de fit_math_cross:

1. construya una escuela nueva con meanses_gmc en el percentil 20 de la distribución escolar;
2. obtenga predicciones de la media para una secuencia de ses_wc;
3. repita para el percentil 80;
4. compare las distribuciones predictivas de las pendientes nuevas;
5. explique qué fuentes de incertidumbre no aparecerían si utilizara re_formula = NA y predijera solo la media poblacional.

4.24.3 Ejercicio de interpretación de salida

Un modelo produce el siguiente resumen posterior hipotético:

Cantidad	Mediana posterior	Intervalo 90%
\(\mu_\beta\)	2.4	\([1.9,2.9]\)
\(\tau_\beta\)	1.1	\([0.6,1.8]\)
\(\rho\)	-0.55	\([-0.86,0.03]\)
\(\gamma_{11}\)	0.8	\([0.1,1.5]\)
\(\sigma\)	5.9	\([5.6,6.2]\)

El predictor individual \(x\) está centrado por grupo y el predictor grupal \(z\) está centrado globalmente.

Responda:

1. ¿qué representa \(\mu_\beta\)?
2. ¿qué representa \(\tau_\beta\)?
3. ¿cómo interpretaría \(\gamma_{11}\)?
4. ¿puede afirmarse que \(z\) explica toda la heterogeneidad de pendientes?
5. ¿cómo interpretaría prudentemente la posterior de \(\rho\)?
6. si redefinimos \(x^*=x-2\), ¿espera que \(\rho\) sea necesariamente igual?
7. ¿qué comprobación predictiva utilizaría para evaluar si el modelo reproduce la heterogeneidad de pendientes?
8. ¿qué diferencia conceptual existe entre predecir la pendiente de un grupo observado y la de un grupo nuevo?

4.25 Lecturas para profundizar

Para esta semana se recomienda priorizar:

• Gelman y Hill, §13.1: interceptos y pendientes variables, distribución conjunta y predictores grupales para coeficientes (Gelman y Hill 2007, 279-83).
• Gelman y Hill, §13.4: interpretación de la correlación entre interceptos y pendientes y su relación con el reescalamiento del predictor (Gelman y Hill 2007, 287-89).
• McElreath, §14.1: construcción generativa de pendientes variables, covarianza, LKJ y shrinkage bidimensional (McElreath 2020, 437-46).
• McElreath, §13.5: predicciones multinivel para clusters existentes y nuevos (McElreath 2020, 426-31).

Como complemento:

• Hox, Moerbeek y van de Schoot, §§4.2-4.3: centrado de predictores con pendientes variables e interpretación de interacciones, incluidas interacciones entre niveles (Hox et al. 2018, 46-56).
• Bayesian Workflow, §17.2: previas para modelos lineales multinivel con interceptos y pendientes variables y visualización del pooling conjunto (Gelman et al. 2026, 280-85).
• Bürkner: estructura de los parámetros grupales en brms, factorización de matrices de covarianza y sintaxis de múltiples coeficientes por grupo (Bürkner 2017, 2-7).

Bürkner, Paul-Christian. 2017. «brms: An R Package for Bayesian Multilevel Models Using Stan». Journal of Statistical Software 80 (1): 1-28. https://doi.org/10.18637/jss.v080.i01.

Gelman, Andrew, y Jennifer Hill. 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.

Gelman, Andrew, Aki Vehtari, Richard McElreath, et al. 2026. Bayesian Workflow. Chapman & Hall/CRC.

Hox, Joop J., Mirjam Moerbeek, y Rens van de Schoot. 2018. Multilevel Analysis: Techniques and Applications. 3.ª ed. Routledge.

McElreath, Richard. 2020. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. 2.ª ed. Chapman & Hall/CRC.