2  Semana 2. Interceptos variables, ICC y shrinkage

SP-1653 Modelos Mixtos

2.1 Panorama de la semana

En la semana anterior introdujimos tres formas de tratar datos agrupados: pooling completo, ausencia de pooling y pooling parcial. El modelo de interceptos variables apareció entonces como una primera forma de representar probabilísticamente la dependencia entre observaciones del mismo grupo.

Esta semana estudiamos ese modelo con mayor profundidad. La pregunta central es:

¿Cómo se modela e interpreta la heterogeneidad entre grupos?

La respuesta requiere separar tres ideas que a menudo se mezclan:

1. la variación dentro de los grupos, representada por la desviación estándar residual \(\sigma\);
2. la variación entre grupos, representada por la desviación estándar \(\tau_\alpha\) de los interceptos;
3. la incertidumbre sobre ambas cantidades y sobre cada intercepto particular \(\alpha_j\).

Esta separación conduce directamente al coeficiente de correlación intraclase (ICC) y al shrinkage. El ICC resume, bajo el modelo gaussiano de interceptos variables, cuánto de la variación marginal corresponde al nivel grupal y qué correlación induce el grupo entre dos unidades. El shrinkage describe cómo la información de cada grupo y la información de la población de grupos se combinan para aprender \(\alpha_j\).

Gelman y Hill presentan el pooling parcial como un promedio adaptativo entre la información propia del grupo y la distribución poblacional, con mayor contracción cuando los grupos contienen menos información (Gelman y Hill 2007, 251-59). McElreath interpreta el mismo mecanismo como regularización adaptativa y muestra que, en promedio, puede mejorar la predicción especialmente en grupos pequeños (McElreath 2020, 401-13). BDA desarrolla la estructura normal jerárquica que permite ver algebraicamente de dónde surge este promedio ponderado (Gelman et al. 2013, 113-23).

Objetivos de aprendizaje

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

1. escribir el modelo nulo gaussiano de interceptos variables en forma jerárquica y en forma de componentes;
2. distinguir la distribución de las observaciones de la distribución poblacional de los interceptos;
3. derivar la media, la varianza y la covarianza marginal de dos observaciones del mismo grupo;
4. derivar e interpretar el ICC para el modelo nulo gaussiano;
5. obtener una distribución posterior del ICC, en lugar de tratarlo como una única estimación puntual;
6. derivar el promedio posterior condicional de un intercepto grupal y reconocer el factor de shrinkage;
7. explicar cómo \(n_j\), \(\sigma\) y \(\tau_\alpha\) determinan la magnitud de la contracción;
8. relacionar algebraicamente el factor de shrinkage con el ICC;
9. diferenciar predicción para grupos existentes de predicción para grupos nuevos;
10. ajustar un modelo de interceptos variables en brms, realizar diagnósticos MCMC básicos y efectuar comprobaciones predictivas globales y por grupo;
11. reconocer por qué pocos grupos implican mayor incertidumbre en la heterogeneidad entre grupos;
12. resumir en los datos del proyecto el número de grupos, los tamaños grupales y las fuentes preliminares de variación dentro y entre grupos.

2.2 Problema motivador: mismos centros, nueva pregunta

Retomemos la situación de la semana 1. Observamos una respuesta continua en varios centros:

\[ i=1,\ldots,n_j, \qquad j=1,\ldots,J. \]

La semana anterior la pregunta era principalmente estructural: ¿debemos ignorar los centros, estimarlos por separado o compartir información entre ellos?

Ahora suponemos que el modelo de pooling parcial es una representación inicial razonable y hacemos preguntas más específicas:

• ¿cuánto difieren los centros entre sí?
• ¿cuánto varían las personas dentro de un mismo centro?
• ¿qué tan semejantes esperamos que sean dos personas del mismo centro?
• ¿cuánto debe confiar una estimación grupal en su propia media observada?
• ¿cómo cambia esa confianza si un centro aporta 5 observaciones en lugar de 50?
• ¿qué debemos simular si queremos predecir una observación en un centro que nunca estuvo en la muestra?

Estas preguntas apuntan a cantidades diferentes. Conviene no reducir todo el análisis a “estimar un efecto aleatorio”.

ImportanteTres niveles de incertidumbre

En un modelo jerárquico debemos distinguir entre:

• variación de observaciones alrededor de su grupo;
• heterogeneidad real de los parámetros entre grupos;
• incertidumbre posterior acerca de esas cantidades.

Una desviación estándar grupal grande no es lo mismo que una gran incertidumbre acerca de la desviación estándar grupal.

2.3 Modelo nulo de interceptos variables

Comenzamos sin predictores para aislar la estructura de dependencia. El modelo generativo es

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

\[ \alpha_j\mid\mu_\alpha,\tau_\alpha \sim \mathcal N(\mu_\alpha,\tau_\alpha^2). \]

Usaremos la convención de que \(\mathcal N(\mu,\sigma^2)\) se parametriza con varianza en las ecuaciones. En brms y Stan, normal(mu, sigma) usa desviación estándar.

Los parámetros tienen interpretaciones distintas:

Parámetro	Nivel	Interpretación
\(\alpha_j\)	grupo observado	media esperada de la respuesta en el grupo \(j\)
\(\mu_\alpha\)	población de grupos	media de la distribución poblacional de interceptos
\(\tau_\alpha\)	población de grupos	desviación estándar entre los interceptos
\(\sigma\)	observación dentro de grupo	desviación estándar residual alrededor de \(\alpha_j\)
————	————	————

El énfasis en distinguir \(\tau_\alpha\) de \(\sigma\) es central. En el ejemplo lineal multinivel de Bayesian Workflow, la desviación estándar del intercepto variable representa variación entre personas, mientras la desviación estándar residual representa variación dentro de personas una vez incorporados los predictores (Gelman et al. 2026, 280-83).

2.3.1 Forma equivalente con desviaciones grupales

También podemos escribir

\[ \alpha_j=\mu_\alpha+u_j, \qquad u_j\sim\mathcal N(0,\tau_\alpha^2), \]

de modo que

\[ y_{ij} = \mu_\alpha+u_j+\varepsilon_{ij}, \]

con

\[ \varepsilon_{ij}\sim\mathcal N(0,\sigma^2). \]

Esta forma conecta con la notación clásica de modelos mixtos, en la cual \(u_j\) suele denominarse un efecto aleatorio de intercepto. Hox et al. escriben el modelo vacío (nulo) como una suma de un intercepto global, un error de nivel grupal y un error de nivel individual (Hox et al. 2018, 12-13).

En estas notas preferiremos hablar de interceptos variables y de desviaciones grupales, pero señalaremos la terminología tradicional cuando sea útil para leer la literatura.

2.4 De la jerarquía a la distribución marginal

La jerarquía permite derivar propiedades marginales de \(Y_{ij}\). Estas propiedades explican por qué la agrupación induce dependencia incluso cuando los errores \(\varepsilon_{ij}\) son independientes.

2.4.1 Media marginal

Como

\[ Y_{ij} = \mu_\alpha+u_j+\varepsilon_{ij}, \]

y

\[ E(u_j)=0, \qquad E(\varepsilon_{ij})=0, \]

tenemos

\[ E(Y_{ij})=\mu_\alpha. \]

2.4.2 Varianza marginal

Por independencia entre \(u_j\) y \(\varepsilon_{ij}\),

\[ \begin{aligned} \operatorname{Var}(Y_{ij}) &= \operatorname{Var}(u_j+\varepsilon_{ij})\\ &= \operatorname{Var}(u_j) + \operatorname{Var}(\varepsilon_{ij})\\ &= \tau_\alpha^2+\sigma^2. \end{aligned} \]

Por tanto, la variabilidad marginal total se descompone en

\[ \boxed{ \operatorname{Var}(Y_{ij}) = \underbrace{\tau_\alpha^2}_{\text{entre grupos}} + \underbrace{\sigma^2}_{\text{dentro de grupos}} }. \]

2.4.3 Covarianza dentro del mismo grupo

Para dos observaciones distintas \(i\neq i'\) del mismo grupo \(j\),

\[ Y_{ij}=\mu_\alpha+u_j+\varepsilon_{ij}, \]

\[ Y_{i'j}=\mu_\alpha+u_j+\varepsilon_{i'j}. \]

Ambas comparten \(u_j\). Entonces

\[ \begin{aligned} \operatorname{Cov}(Y_{ij},Y_{i'j}) &= \operatorname{Cov} (u_j+\varepsilon_{ij}, u_j+\varepsilon_{i'j})\\ &= \operatorname{Var}(u_j)\\ &= \tau_\alpha^2, \end{aligned} \]

si \(u_j\), \(\varepsilon_{ij}\) y \(\varepsilon_{i'j}\) son mutuamente independientes.

Para observaciones pertenecientes a grupos distintos, \(j\neq k\), suponemos además independencia entre \(u_j\) y \(u_k\), de modo que

\[ \operatorname{Cov}(Y_{ij},Y_{i'k})=0. \]

NotaCondicionalmente independientes, marginalmente correlacionadas

Dado \(\alpha_j\), las observaciones del grupo pueden ser independientes:

\[ p(\mathbf y_j\mid \alpha_j,\sigma) = \prod_i p(y_{ij}\mid\alpha_j,\sigma). \]

Al integrar la incertidumbre sobre \(\alpha_j\), las observaciones del mismo grupo quedan correlacionadas porque comparten el mismo componente \(u_j\).

2.5 Coeficiente de correlación intraclase

El coeficiente de correlación intraclase del modelo nulo gaussiano se obtiene dividiendo la covarianza de dos observaciones del mismo grupo por su desviación estándar marginal:

\[ \rho = \frac{ \operatorname{Cov}(Y_{ij},Y_{i'j}) }{ \sqrt{\operatorname{Var}(Y_{ij}) \operatorname{Var}(Y_{i'j})} }. \]

Como ambas varianzas marginales son \(\tau_\alpha^2+\sigma^2\),

\[ \boxed{ \rho = \frac{\tau_\alpha^2} {\tau_\alpha^2+\sigma^2} }. \tag{2.1}\]

Hox et al. derivan esta expresión a partir del modelo vacío y señalan dos interpretaciones equivalentes en este caso: proporción de la variación total correspondiente al nivel grupal y correlación esperada entre dos unidades tomadas del mismo grupo (Hox et al. 2018, 12-13). Gelman y Hill presentan la misma razón de varianzas al discutir la variación individual y grupal en el modelo de interceptos variables (Gelman y Hill 2007, 258).

2.5.1 Interpretación como proporción de varianza

En el modelo nulo,

\[ \rho = \frac{\text{varianza entre grupos}} {\text{varianza marginal total}}. \]

Por ejemplo, \(\rho=0.30\) significa que, bajo este modelo, 30% de la varianza marginal de la respuesta corresponde a heterogeneidad entre interceptos de grupo y 70% a heterogeneidad dentro de los grupos.

Esta interpretación depende del modelo. El ICC no es una propiedad inmutable del archivo de datos.

2.5.2 Interpretación como correlación

Si elegimos dos unidades distintas del mismo grupo,

\[ \operatorname{Corr}(Y_{ij},Y_{i'j})=\rho. \]

Un ICC cercano a cero implica poca semejanza marginal inducida por el intercepto compartido. Un ICC alto implica que conocer el grupo aporta mucha información acerca del nivel de una observación.

2.5.3 Casos límite

Si \(\tau_\alpha\rightarrow0\),

\[ \rho\rightarrow0, \]

y los grupos dejan de diferir en su intercepto. El modelo se aproxima al pooling completo.

Si \(\sigma\rightarrow0\) mientras \(\tau_\alpha>0\),

\[ \rho\rightarrow1, \]

y las observaciones de un mismo grupo quedan casi completamente determinadas por su intercepto común.

AdvertenciaEl ICC no decide por sí solo si necesitamos un modelo multinivel

Un ICC pequeño no implica automáticamente que el agrupamiento sea irrelevante. Su importancia depende de la pregunta, del tamaño de los grupos, del diseño, de los predictores y del estimando. Tampoco existe un umbral universal de ICC que determine cuándo “usar” o “no usar” un modelo multinivel.

2.5.4 ICC condicional después de incorporar predictores

Si añadimos un predictor,

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

la razón

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

describe la correlación inducida por el intercepto compartido condicionalmente en los predictores incluidos. Las componentes de varianza pueden cambiar al incorporar covariables.

Por eso no conviene interpretar el ICC de un modelo nulo como si fuera una constante física de la población. El modelo determina qué variación se atribuye a cada componente.

2.6 ICC bayesiano: una distribución, no un único número

En un análisis bayesiano, \(\tau_\alpha\) y \(\sigma\) tienen distribuciones posteriores. El ICC es entonces una cantidad derivada:

\[ \rho^{(s)} = \frac{ \left(\tau_\alpha^{(s)}\right)^2 }{ \left(\tau_\alpha^{(s)}\right)^2 + \left(\sigma^{(s)}\right)^2 }, \qquad s=1,\ldots,S. \]

Cada draw posterior produce un valor posible de \(\rho\). La colección

\[ \rho^{(1)},\ldots,\rho^{(S)} \]

aproxima la distribución posterior del ICC.

Esta forma de proceder tiene una ventaja conceptual importante: la incertidumbre acerca de las componentes de varianza se propaga automáticamente a la incertidumbre acerca del ICC.

BDA enfatiza que sustituir componentes jerárquicos desconocidos por estimaciones puntuales ignora incertidumbre que puede ser considerable, especialmente cuando la población de grupos es pequeña (Gelman et al. 2013, 113-23).

2.7 ¿De dónde sale el shrinkage?

La palabra shrinkage describe la contracción de una estimación grupal hacia la distribución poblacional. En el modelo normal, el mecanismo puede verse de manera exacta cuando condicionamos en los hiperparámetros.

Supongamos provisionalmente que \(\mu_\alpha\), \(\tau_\alpha\) y \(\sigma\) son conocidos.

Para el grupo \(j\),

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

La media muestral satisface

\[ \bar y_j\mid\alpha_j,\sigma \sim \mathcal N \left( \alpha_j,\frac{\sigma^2}{n_j} \right). \]

La distribución poblacional del intercepto es

\[ \alpha_j\mid\mu_\alpha,\tau_\alpha \sim \mathcal N(\mu_\alpha,\tau_\alpha^2). \]

Tenemos entonces una verosimilitud normal para \(\bar y_j\) y una distribución normal para \(\alpha_j\). La distribución posterior condicional es normal:

\[ \alpha_j \mid \bar y_j,\mu_\alpha,\tau_\alpha,\sigma \sim \mathcal N(m_j,V_j), \]

donde la precisión posterior es la suma de precisiones,

\[ \frac{1}{V_j} = \frac{n_j}{\sigma^2} + \frac{1}{\tau_\alpha^2}, \]

por lo que

\[ V_j = \left( \frac{n_j}{\sigma^2} + \frac{1}{\tau_\alpha^2} \right)^{-1}. \]

La media posterior es

\[ m_j = V_j \left( \frac{n_j}{\sigma^2}\bar y_j + \frac{1}{\tau_\alpha^2}\mu_\alpha \right). \]

Reorganizando,

\[ \boxed{ m_j = w_j\bar y_j+(1-w_j)\mu_\alpha } \tag{2.2}\]

con

\[ \boxed{ w_j = \frac{n_j\tau_\alpha^2} {n_j\tau_\alpha^2+\sigma^2} } \tag{2.3}\]

y

\[ 1-w_j = \frac{\sigma^2} {n_j\tau_\alpha^2+\sigma^2}. \]

BDA desarrolla esta estructura para el modelo normal jerárquico de medias intercambiables (Gelman et al. 2013, 113-19). Gelman y Hill presentan la estimación multinivel del coeficiente grupal como un promedio ponderado entre la estimación propia del grupo y el promedio poblacional (Gelman y Hill 2007, 258-59).

NotaQué es elaboración pedagógica

La notación \(w_j\) y la derivación paso a paso de Ecuación 2.2 y Ecuación 2.3 se presentan aquí como una elaboración pedagógica del modelo normal conjugado. En el ajuste bayesiano completo no fijamos \(\mu_\alpha\), \(\tau_\alpha\) y \(\sigma\); integramos su incertidumbre mediante draws posteriores.

2.7.1 ¿Qué controla la magnitud de la contracción?

El peso \(w_j\) es el peso de la media observada del grupo. Por tanto, \(1-w_j\) mide cuánta influencia tiene el centro de la distribución poblacional.

2.7.1.1 Tamaño del grupo

Cuando \(n_j\) crece,

\[ w_j\rightarrow1. \]

La media del grupo se estima con más precisión y el grupo depende menos de la distribución poblacional.

Cuando \(n_j\) es pequeño, \(w_j\) disminuye y aumenta el pooling.

2.7.1.2 Variación residual

Si \(\sigma\) es grande, las observaciones dentro de un grupo son ruidosas y

\[ w_j \downarrow. \]

La media observada \(\bar y_j\) contiene menos información precisa acerca de \(\alpha_j\), de modo que la jerarquía tiene mayor influencia.

2.7.1.3 Heterogeneidad entre grupos

Si \(\tau_\alpha\) es grande,

\[ w_j \uparrow. \]

Una población muy heterogénea hace plausible que un grupo se encuentre lejos de \(\mu_\alpha\). La estimación se contrae menos.

Si \(\tau_\alpha\) es pequeña, la población de interceptos está concentrada y la contracción es mayor.

McElreath muestra gráficamente estos patrones: los grupos pequeños se contraen más y la regularización es especialmente útil en términos de error promedio cuando hay poca información por grupo (McElreath 2020, 405-13).

2.8 Relación entre ICC y shrinkage

El ICC y el shrinkage no son la misma cantidad, pero están algebraicamente relacionados.

Partimos de

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

Despejando,

\[ \tau_\alpha^2 = \frac{\rho}{1-\rho}\sigma^2. \]

Sustituimos en Ecuación 2.3:

\[ w_j = \frac{ n_j \left[ \frac{\rho}{1-\rho}\sigma^2 \right] }{ n_j \left[ \frac{\rho}{1-\rho}\sigma^2 \right] + \sigma^2 }. \]

Después de simplificar,

\[ \boxed{ w_j = \frac{n_j\rho} {1+(n_j-1)\rho} }. \tag{2.4}\]

Esta identidad muestra que, para el modelo nulo normal y condicionando en los hiperparámetros, el peso que recibe la información del grupo depende de dos cosas:

• el tamaño del grupo \(n_j\);
• la magnitud de la dependencia intraclase \(\rho\).

Algunas consecuencias:

• si \(\rho\rightarrow0\), entonces \(w_j\rightarrow0\): los datos no apoyan diferencias persistentes entre grupos y el pooling es fuerte;
• si \(\rho\rightarrow1\), entonces \(w_j\rightarrow1\): las diferencias entre grupos dominan la variación residual y hay poca contracción;
• para un \(\rho\) fijo, \(w_j\) aumenta con \(n_j\).

ImportanteShrinkage no equivale a multiplicar por el ICC

El peso de la información grupal no es simplemente \(\rho\). Incluso con el mismo ICC, dos grupos con tamaños distintos pueden experimentar grados muy diferentes de contracción.

2.9 Shrinkage como regularización adaptativa

El shrinkage puede parecer extraño si se compara la mediana posterior de \(\alpha_j\) con la media observada \(\bar y_j\). ¿Por qué no reproducir exactamente la media de cada grupo?

Porque el objetivo no es memorizar la muestra. El modelo intenta aprender parámetros latentes y realizar predicciones.

Los grupos con poca información producen estimaciones no agrupadas de alta varianza. El modelo jerárquico utiliza la población de grupos para estabilizarlas. McElreath describe esta propiedad como regularización adaptativa: la fuerza de la regularización se aprende conjuntamente con el modelo y no tiene que ser igual para todos los grupos (McElreath 2020, 408-13).

Esto produce una distinción importante:

• un modelo sin pooling puede ajustarse más estrechamente a medias grupales observadas;
• un modelo con pooling parcial puede sacrificar parte de ese ajuste dentro de la muestra para reducir error de estimación y mejorar generalización.

El shrinkage no garantiza que cada estimación grupal esté más cerca del verdadero \(\alpha_j\) en cada conjunto de datos. Su beneficio es probabilístico y promedio, no una promesa determinista grupo por grupo (McElreath 2020, 412-13).

2.10 El modelo bayesiano completo

Hasta ahora condicionamos varias derivaciones en hiperparámetros conocidos. En un análisis real, debemos aprenderlos.

Para nuestro ejemplo simulado usaremos

\[ \begin{aligned} y_{ij} &\sim \mathcal N(\alpha_j,\sigma^2),\\ \alpha_j &\sim \mathcal N(\mu_\alpha,\tau_\alpha^2),\\ \mu_\alpha &\sim \mathcal N(50,20^2),\\ \tau_\alpha &\sim \operatorname{Normal}^{+}(0,10^2),\\ \sigma &\sim \operatorname{Exponential}(0.1). \end{aligned} \]

Estas previas corresponden a la escala artificial del ejemplo y no son recomendaciones universales.

La posterior conjunta puede escribirse, salvo constante de proporcionalidad, como

\[ \begin{aligned} p( \boldsymbol\alpha,\mu_\alpha,\tau_\alpha,\sigma \mid \mathbf y ) \propto& \left[ \prod_{j=1}^{J} \prod_{i=1}^{n_j} p(y_{ij}\mid\alpha_j,\sigma) \right]\\ &\times \left[ \prod_{j=1}^{J} p(\alpha_j\mid\mu_\alpha,\tau_\alpha) \right]\\ &\times p(\mu_\alpha) p(\tau_\alpha) p(\sigma). \end{aligned} \]

A diferencia de la derivación conjugada condicional, aquí \(\tau_\alpha\) y \(\sigma\) no se sustituyen por valores puntuales. Sus incertidumbres se propagan a los \(\alpha_j\), al ICC, al grado de shrinkage y a las predicciones.

2.10.1 Sobre las previas para componentes de escala

Las desviaciones estándar son positivas y, con pocos grupos, pueden estar débilmente identificadas por los datos. McElreath advierte que estimar componentes de variación con pocos grupos puede requerir regularización más informativa (McElreath 2020, 406-7). Gelman y Hill también señalan que con pocos grupos la principal dificultad se encuentra en aprender la variación entre grupos (Gelman y Hill 2007, 275-76).

La semana 6 estudiará este problema con detalle. Por ahora conservamos tres principios:

1. expresar la previa en una escala sustantivamente interpretable;
2. evitar previas absurdamente amplias por defecto;
3. examinar qué datos puede generar el modelo antes de condicionarlo en las observaciones.

2.11 Laboratorio reproducible en R

El laboratorio tiene tres objetivos:

1. visualizar la descomposición dentro/entre grupos;
2. estimar la distribución posterior del ICC;
3. conectar el tamaño grupal con la magnitud posterior del shrinkage.

2.11.1 Simular grupos de tamaños muy diferentes

Usaremos doce grupos, con tamaños entre 5 y 60.

J <- 12

n_j <- c(
  5, 6, 8, 10,
  12, 15, 20, 25,
  30, 40, 50, 60
)

mu_alpha_real <- 50
tau_alpha_real <- 7
sigma_real <- 12

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

parametros_grupo <- tibble(
  grupo = factor(niveles_grupo, levels = niveles_grupo),
  n = n_j,
  alpha_real = rnorm(
    J,
    mean = mu_alpha_real,
    sd = tau_alpha_real
  )
)

datos <- parametros_grupo |>
  tidyr::uncount(n, .id = "i") |>
  mutate(
    y = rnorm(
      n(),
      mean = alpha_real,
      sd = sigma_real
    )
  )

stopifnot(
  !anyNA(datos),
  nlevels(datos$grupo) == J,
  all(count(datos, grupo)$n == n_j)
)

datos |>
  count(grupo, name = "n")

# A tibble: 12 × 2
   grupo     n
   <fct> <int>
 1 G01       5
 2 G02       6
 3 G03       8
 4 G04      10
 5 G05      12
 6 G06      15
 7 G07      20
 8 G08      25
 9 G09      30
10 G10      40
11 G11      50
12 G12      60

El proceso generador es

\[ \alpha_j\sim\mathcal N(50,7^2), \]

\[ y_{ij}\sim\mathcal N(\alpha_j,12^2). \]

El ICC verdadero usado en la simulación es

\[ \rho_{\text{real}} = \frac{7^2}{7^2+12^2}. \]

icc_real <- tau_alpha_real^2 /
  (tau_alpha_real^2 + sigma_real^2)

icc_real

[1] 0.253886

2.11.2 Visualizar variación dentro y entre grupos

ggplot(
  datos,
  aes(x = grupo, y = y)
) +
  geom_jitter(
    width = 0.12,
    height = 0,
    alpha = 0.45
  ) +
  geom_point(
    data = parametros_grupo,
    aes(y = alpha_real),
    shape = 4,
    size = 3,
    stroke = 1.1
  ) +
  labs(
    x = "Grupo",
    y = "Respuesta",
    subtitle = "La dispersión vertical dentro de cada grupo refleja σ; la separación entre cruces refleja τα"
  ) +
  theme_minimal(base_size = 12)

Figura 2.1: Datos simulados por grupo. Las cruces indican los interceptos verdaderos del proceso generador y los puntos muestran la variación residual dentro de cada grupo.

La figura permite distinguir dos escalas de variación. No debemos usar el rango de las medias observadas como estimación directa de \(\tau_\alpha\): las medias contienen error muestral, especialmente en grupos pequeños.

2.11.3 Comprobación predictiva previa

Antes de ajustar el modelo, simulamos de las previas propuestas. Hacemos la simulación manualmente para mostrar la jerarquía completa.

set.seed(1654)

S_prior <- 6000

prior_pred <- tibble(
  mu_alpha = rnorm(S_prior, 50, 20),
  tau_alpha = abs(rnorm(S_prior, 0, 10)),
  sigma = rexp(S_prior, rate = 0.1),
  alpha = rnorm(S_prior, mu_alpha, tau_alpha),
  y_rep = rnorm(S_prior, alpha, sigma)
)

prior_pred |>
  summarise(
    q01 = quantile(y_rep, 0.01),
    q50 = quantile(y_rep, 0.50),
    q99 = quantile(y_rep, 0.99)
  )

# A tibble: 1 × 3
    q01   q50   q99
  <dbl> <dbl> <dbl>
1 -13.5  49.8  115.

ggplot(
  prior_pred,
  aes(x = y_rep)
) +
  geom_density() +
  labs(
    x = expression(y^rep),
    y = "Densidad",
    subtitle = "La plausibilidad debe juzgarse en la escala de la respuesta"
  ) +
  theme_minimal(base_size = 12)

Figura 2.2: Distribución predictiva previa para una observación de un grupo genérico bajo las previas pedagógicas del ejemplo.

Esta comprobación no pretende agotar la especificación de previas. Su función es impedir que una previa aparentemente “débil” en la escala de parámetros implique observaciones absurdas en la escala de datos.

2.11.4 Ajustar el modelo nulo de interceptos variables

prior_null <- c(
  prior(normal(50, 20), class = "Intercept"),
  prior(normal(0, 10), class = "sd", group = "grupo"),
  prior(exponential(0.1), class = "sigma")
)

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

La fórmula

y ~ 1 + (1 | grupo)

puede relacionarse con

\[ \alpha_j=\mu_\alpha+u_j, \qquad u_j\sim\mathcal N(0,\tau_\alpha^2). \]

En la parametrización interna de brms, el intercepto de nivel poblacional y la desviación grupal se representan por separado. La distribución normal jerárquica de los coeficientes grupales es parte de la estructura del modelo (Bürkner 2017, 2-4).

2.11.5 Diagnóstico MCMC mínimo

Antes de interpretar el modelo, examinamos algunos parámetros principales.

draws_null <- posterior::as_draws_df(fit_null)

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

# A tibble: 3 × 6
  variable             mean    sd  rhat ess_bulk ess_tail
  <chr>               <dbl> <dbl> <dbl>    <dbl>    <dbl>
1 b_Intercept         50.7  2.57   1.00     808.     879.
2 sd_grupo__Intercept  7.93 1.97   1.00    1118.    1820.
3 sigma               12.1  0.519  1.00    2649.    2623.

np_null <- brms::nuts_params(fit_null)

n_divergencias <- sum(
  np_null$Parameter == "divergent__" &
    np_null$Value == 1
)

tibble(
  divergencias_post_warmup = n_divergencias
)

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

bayesplot::mcmc_trace(
  as.array(fit_null),
  pars = c(
    "b_Intercept",
    "sd_grupo__Intercept",
    "sigma"
  )
)

Figura 2.3: Trazas MCMC para el intercepto poblacional y las dos componentes de escala del modelo nulo.

En esta semana el diagnóstico es deliberadamente básico. La semana 7 desarrollará \(\widehat R\), ESS, divergencias, profundidad del árbol y parametrizaciones centradas/no centradas.

Por ahora exigimos como mínimo:

• \(\widehat R\) muy cercano a 1;
• tamaños efectivos razonables;
• cadenas que exploran regiones semejantes;
• ausencia de divergencias post-warmup.

AdvertenciaConvergencia no implica adecuación

Un algoritmo puede muestrear correctamente de una distribución posterior correspondiente a un modelo sustantivamente inadecuado. Los diagnósticos computacionales y las comprobaciones del modelo responden preguntas diferentes.

2.11.6 Distribución posterior de las componentes de variación

componentes <- draws_null |>
  transmute(
    mu_alpha = b_Intercept,
    tau_alpha = sd_grupo__Intercept,
    sigma = sigma
  )

posterior::summarise_draws(
  posterior::as_draws_df(componentes)
)

# A tibble: 3 × 10
  variable   mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
1 mu_alpha  50.7   50.5  2.57  2.35  46.7   55.0  1.00     804.     877.
2 tau_alpha  7.93   7.63 1.97  1.83   5.27  11.5  1.00    1118.    1813.
3 sigma     12.1   12.1  0.519 0.519 11.3   13.0  1.00    2586.    2622.

componentes |>
  select(tau_alpha, sigma) |>
  pivot_longer(
    everything(),
    names_to = "componente",
    values_to = "valor"
  ) |>
  mutate(
    componente = recode(
      componente,
      tau_alpha = "Entre grupos: τₐ",
      sigma = "Dentro de grupos: σ"
    )
  ) |>
  ggplot(aes(x = valor)) +
  geom_density() +
  facet_wrap(
    ~ componente,
    scales = "free"
  ) +
  labs(
    x = "Desviación estándar",
    y = "Densidad"
  ) +
  theme_minimal(base_size = 12)

Figura 2.4: Distribuciones posteriores de la desviación estándar entre grupos y de la desviación estándar residual.

No debemos comparar únicamente las medianas de \(\tau_\alpha\) y \(\sigma\). La incertidumbre posterior de \(\tau_\alpha\) suele ser mayor porque su información efectiva proviene principalmente del número de grupos \(J\), no del número total de filas.

2.11.7 Posterior del ICC

icc_draws <- componentes |>
  transmute(
    tau_alpha = tau_alpha,
    sigma = sigma,
    icc = tau_alpha^2 /
      (tau_alpha^2 + sigma^2)
  )

icc_resumen <- icc_draws |>
  summarise(
    mediana = median(icc),
    q05 = quantile(icc, 0.05),
    q95 = quantile(icc, 0.95),
    prob_mayor_025 = mean(icc > 0.25)
  )

icc_resumen

# A tibble: 1 × 4
  mediana   q05   q95 prob_mayor_025
    <dbl> <dbl> <dbl>          <dbl>
1   0.285 0.157 0.482          0.639

ggplot(
  icc_draws,
  aes(x = icc)
) +
  geom_density() +
  geom_vline(
    xintercept = icc_real,
    linetype = 2
  ) +
  labs(
    x = "ICC",
    y = "Densidad",
    subtitle = "La incertidumbre de τα y σ se propaga a ρ"
  ) +
  theme_minimal(base_size = 12)

Figura 2.5: Distribución posterior del ICC. La línea vertical marca el ICC verdadero utilizado en la simulación.

En datos reales no conoceremos la línea vertical “verdadera”. La tenemos aquí porque el conjunto de datos fue simulado.

2.11.8 Medias observadas y medias posteriores por grupo

Construimos la distribución posterior de la media esperada de cada grupo.

new_grupos <- tibble(
  grupo = factor(
    niveles_grupo,
    levels = niveles_grupo
  )
)

epred_grupos <- posterior_epred(
  fit_null,
  newdata = new_grupos,
  re_formula = NULL
)

colnames(epred_grupos) <- niveles_grupo

epred_long <- as_tibble(epred_grupos) |>
  mutate(.draw = row_number()) |>
  pivot_longer(
    - .draw,
    names_to = "grupo",
    values_to = "alpha"
  )

resumen_alpha <- epred_long |>
  group_by(grupo) |>
  summarise(
    alpha_post = median(alpha),
    q05 = quantile(alpha, 0.05),
    q95 = quantile(alpha, 0.95),
    .groups = "drop"
  )

medias_grupo <- datos |>
  group_by(grupo) |>
  summarise(
    media_observada = mean(y),
    n = n(),
    .groups = "drop"
  ) |>
  mutate(
    grupo = as.character(grupo)
  )

resumen_alpha <- resumen_alpha |>
  left_join(
    medias_grupo,
    by = "grupo"
  ) |>
  left_join(
    parametros_grupo |>
      transmute(
        grupo = as.character(grupo),
        alpha_real = alpha_real
      ),
    by = "grupo"
  )

resumen_alpha

# A tibble: 12 × 7
   grupo alpha_post   q05   q95 media_observada     n alpha_real
   <chr>      <dbl> <dbl> <dbl>           <dbl> <int>      <dbl>
 1 G01         57.1  49.5  65.1            60.6     5       56.9
 2 G02         48.8  41.7  55.7            48.0     6       47.2
 3 G03         55.5  49.7  61.9            57.3     8       54.7
 4 G04         59.2  53.3  65.0            61.5    10       60.7
 5 G05         42.8  37.4  48.0            41.0    12       39.1
 6 G06         55.7  50.8  60.5            56.6    15       59.3
 7 G07         51.5  47.2  55.8            51.6    20       50.7
 8 G08         49.9  46.2  53.6            49.9    25       47.0
 9 G09         39.2  35.7  42.7            38.1    30       38.8
10 G10         55.6  52.6  58.7            56.0    40       55.8
11 G11         52.6  50.0  55.3            52.7    50       49.6
12 G12         40.0  37.5  42.5            39.6    60       37.0

2.11.9 Figura de shrinkage con flechas

resumen_alpha |>
  arrange(n) |>
  mutate(
    grupo = factor(
      grupo,
      levels = grupo
    )
  ) |>
  ggplot() +
  geom_segment(
    aes(
      x = media_observada,
      xend = alpha_post,
      y = grupo,
      yend = grupo,
      linewidth = n
    ),
    arrow = grid::arrow(
      length = grid::unit(0.12, "cm")
    ),
    alpha = 0.65
  ) +
  geom_point(
    aes(
      x = media_observada,
      y = grupo
    ),
    shape = 1,
    size = 2.5
  ) +
  geom_point(
    aes(
      x = alpha_post,
      y = grupo
    ),
    size = 2.5
  ) +
  scale_linewidth_continuous(
    name = expression(n[j])
  ) +
  labs(
    x = "Media del grupo",
    y = "Grupo",
    subtitle = "Círculo abierto: media observada; punto sólido: mediana posterior"
  ) +
  theme_minimal(base_size = 12)

Figura 2.6: Shrinkage por grupo. Cada flecha parte de la media observada y termina en la mediana posterior de la media grupal. Los grupos se ordenan por tamaño.

No esperamos que la longitud de una flecha dependa únicamente de \(n_j\). También depende de cuán lejos está \(\bar y_j\) del centro poblacional y de los valores plausibles de \(\tau_\alpha\) y \(\sigma\).

2.11.10 Distribución posterior del peso de shrinkage

En vez de sustituir \(\tau_\alpha\) y \(\sigma\) por estimaciones puntuales, podemos calcular un \(w_j\) para cada draw posterior:

\[ w_j^{(s)} = \frac{ n_j \left(\tau_\alpha^{(s)}\right)^2 }{ n_j \left(\tau_\alpha^{(s)}\right)^2 + \left(\sigma^{(s)}\right)^2 }. \]

draws_componentes <- componentes |>
  mutate(
    .draw = row_number()
  )

tabla_n <- medias_grupo |>
  select(grupo, n)

pesos <- tidyr::crossing(
  .draw = draws_componentes$.draw,
  grupo = tabla_n$grupo
) |>
  left_join(
    tabla_n,
    by = "grupo"
  ) |>
  left_join(
    draws_componentes,
    by = ".draw"
  ) |>
  mutate(
    w = n * tau_alpha^2 /
      (n * tau_alpha^2 + sigma^2)
  )

resumen_pesos <- pesos |>
  group_by(grupo, n) |>
  summarise(
    w_mediana = median(w),
    w_q05 = quantile(w, 0.05),
    w_q95 = quantile(w, 0.95),
    .groups = "drop"
  )

resumen_pesos

# A tibble: 12 × 5
   grupo     n w_mediana w_q05 w_q95
   <chr> <int>     <dbl> <dbl> <dbl>
 1 G01       5     0.666 0.482 0.823
 2 G02       6     0.705 0.527 0.848
 3 G03       8     0.761 0.598 0.881
 4 G04      10     0.799 0.650 0.903
 5 G05      12     0.827 0.691 0.918
 6 G06      15     0.857 0.736 0.933
 7 G07      20     0.888 0.788 0.949
 8 G08      25     0.909 0.823 0.959
 9 G09      30     0.923 0.848 0.965
10 G10      40     0.941 0.882 0.974
11 G11      50     0.952 0.903 0.979
12 G12      60     0.960 0.918 0.982

ggplot(
  resumen_pesos,
  aes(x = n, y = w_mediana)
) +
  geom_linerange(
    aes(
      ymin = w_q05,
      ymax = w_q95
    )
  ) +
  geom_point(size = 2.2) +
  labs(
    x = expression(n[j]),
    y = expression(w[j]),
    subtitle = "Valores grandes de wj implican menos contracción hacia μα"
  ) +
  theme_minimal(base_size = 12)

Figura 2.7: Peso posterior de la información propia del grupo, wj, según el tamaño del grupo. Los intervalos reflejan incertidumbre posterior en τα y σ.

La figura hace visible una idea central: el shrinkage también tiene incertidumbre. Con pocos grupos, la posterior de \(\tau_\alpha\) puede ser amplia y, por tanto, la cantidad de pooling compatible con los datos también puede ser incierta.

2.12 Comprobación predictiva posterior

El modelo no debe evaluarse únicamente mediante sus parámetros. Simulamos nuevos conjuntos de datos desde la distribución predictiva posterior y preguntamos si reproducen características relevantes de lo observado.

2.12.1 Comprobación global

pp_check(
  fit_null,
  type = "dens_overlay",
  ndraws = 50
)

Figura 2.8: Comprobación predictiva posterior global del modelo nulo de interceptos variables.

Una buena superposición global no garantiza que el modelo represente adecuadamente la heterogeneidad por grupo.

2.12.2 Comprobación de medias por grupo

pp_check(
  fit_null,
  type = "stat_grouped",
  stat = "mean",
  group = "grupo"
)

Figura 2.9: Comprobación predictiva posterior de las medias grupales.

Esta segunda vista es más sensible a una posible mala representación de la variabilidad entre grupos.

2.12.3 Comprobación de una estadística de heterogeneidad

Podemos comparar la desviación estándar observada de las medias grupales con la misma estadística calculada sobre datos replicados.

yrep <- posterior_predict(
  fit_null,
  ndraws = 500
)

grupo_indice <- as.integer(datos$grupo)

sd_medias <- function(y_vec, grupo_idx) {
  medias <- tapply(
    y_vec,
    grupo_idx,
    mean
  )
  sd(medias)
}

t_obs <- sd_medias(
  datos$y,
  grupo_indice
)

t_rep <- apply(
  yrep,
  1,
  sd_medias,
  grupo_idx = grupo_indice
)

tibble(
  t_rep = t_rep,
  t_obs = t_obs
) |>
  summarise(
    observado = first(t_obs),
    mediana_replicada = median(t_rep),
    q05 = quantile(t_rep, 0.05),
    q95 = quantile(t_rep, 0.95)
  )

# A tibble: 1 × 4
  observado mediana_replicada   q05   q95
      <dbl>             <dbl> <dbl> <dbl>
1      8.03              7.93  5.91  10.3

ggplot(
  tibble(t_rep = t_rep),
  aes(x = t_rep)
) +
  geom_density() +
  geom_vline(
    xintercept = t_obs,
    linetype = 2
  ) +
  labs(
    x = "DE de las medias grupales",
    y = "Densidad"
  ) +
  theme_minimal(base_size = 12)

Figura 2.10: Distribución predictiva posterior de la desviación estándar de las medias grupales. La línea vertical corresponde al valor observado.

Esta comprobación se conecta directamente con el propósito del modelo: representar heterogeneidad entre grupos.

2.13 Predicción: grupo existente versus grupo nuevo

Una de las ventajas de escribir explícitamente la distribución poblacional

\[ \alpha_j \sim \mathcal N(\mu_\alpha,\tau_\alpha^2) \]

es que podemos predecir para grupos que no participaron en el ajuste.

Gelman y Hill distinguen entre nuevas observaciones de un grupo existente y observaciones de grupos nuevos (Gelman y Hill 2007, 272-75). McElreath desarrolla la misma distinción al mostrar que la predicción multinivel requiere decidir si condicionamos en los coeficientes de los grupos observados o generamos un nuevo coeficiente desde la población de grupos (McElreath 2020, 426-31).

Conviene separar cuatro cantidades.

2.13.1 1. Media esperada de un grupo existente

Para un grupo observado \(j\),

\[ E(\tilde Y\mid\alpha_j)=\alpha_j. \]

Usamos la posterior de su \(\alpha_j\).

2.13.2 2. Nueva observación de un grupo existente

Además de la incertidumbre en \(\alpha_j\), aparece variación residual:

\[ \tilde Y_{new,j} \sim \mathcal N(\alpha_j,\sigma^2). \]

2.13.3 3. Media esperada de un grupo nuevo

Debemos generar un nuevo intercepto:

\[ \alpha_{\text{new}} \sim \mathcal N(\mu_\alpha,\tau_\alpha^2). \]

Por tanto, incluso antes de añadir el error residual hay incertidumbre adicional asociada a qué tipo de grupo nuevo obtendremos.

2.13.4 4. Nueva observación de un grupo nuevo

Generamos primero

\[ \alpha_{\text{new}} \sim \mathcal N(\mu_\alpha,\tau_\alpha^2), \]

y luego

\[ \tilde Y_{\text{new}} \sim \mathcal N(\alpha_{\text{new}},\sigma^2). \]

2.13.5 Comparación por simulación posterior

nd_existente <- data.frame(
  grupo = factor(
    "G01",
    levels = niveles_grupo
  )
)

nd_nuevo <- data.frame(
  grupo = "G_nuevo"
)

mu_existente <- posterior_epred(
  fit_null,
  newdata = nd_existente,
  re_formula = NULL
)[, 1]

y_existente <- posterior_predict(
  fit_null,
  newdata = nd_existente,
  re_formula = NULL
)[, 1]

mu_nuevo <- posterior_epred(
  fit_null,
  newdata = nd_nuevo,
  re_formula = NULL,
  allow_new_levels = TRUE,
  sample_new_levels = "gaussian"
)[, 1]

y_nuevo <- posterior_predict(
  fit_null,
  newdata = nd_nuevo,
  re_formula = NULL,
  allow_new_levels = TRUE,
  sample_new_levels = "gaussian"
)[, 1]

predicciones_tipo <- bind_rows(
  tibble(
    valor = mu_existente,
    cantidad = "Media: grupo existente"
  ),
  tibble(
    valor = y_existente,
    cantidad = "Observación: grupo existente"
  ),
  tibble(
    valor = mu_nuevo,
    cantidad = "Media: grupo nuevo"
  ),
  tibble(
    valor = y_nuevo,
    cantidad = "Observación: grupo nuevo"
  )
)

predicciones_tipo |>
  group_by(cantidad) |>
  summarise(
    mediana = median(valor),
    q05 = quantile(valor, 0.05),
    q95 = quantile(valor, 0.95),
    .groups = "drop"
  )

# A tibble: 4 × 4
  cantidad                     mediana   q05   q95
  <chr>                          <dbl> <dbl> <dbl>
1 Media: grupo existente          57.1  49.5  65.1
2 Media: grupo nuevo              50.7  36.4  64.6
3 Observación: grupo existente    57.9  36.3  77.7
4 Observación: grupo nuevo        50.8  25.7  74.7

ggplot(
  predicciones_tipo,
  aes(x = valor)
) +
  geom_density() +
  facet_wrap(
    ~ cantidad,
    scales = "free_y"
  ) +
  labs(
    x = "Respuesta",
    y = "Densidad"
  ) +
  theme_minimal(base_size = 12)

Figura 2.11: Cuatro distribuciones predictivas distintas: media u observación, para un grupo existente o para un grupo nuevo.

ImportanteLa unidad predictiva debe declararse

“Predecir un dato nuevo” es ambiguo en un modelo multinivel. Siempre debemos especificar si el nuevo dato pertenece a un grupo ya observado o a un grupo que también debe ser generado desde la población de grupos.

2.14 Introducir un predictor sin cambiar todavía la lógica

El siguiente paso natural es

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

\[ \alpha_j \sim \mathcal N( \mu_\alpha, \tau_\alpha^2 ). \]

La pendiente \(\beta\) es común a todos los grupos y el intercepto sigue variando.

Ahora:

• \(\sigma\) describe variación residual después de condicionar en \(x\);
• \(\tau_\alpha\) describe heterogeneidad residual entre interceptos después de condicionar en \(x\);
• el ICC calculado con estas componentes es condicional a la estructura del modelo.

No desarrollaremos todavía el significado de centrar \(x\), separar asociaciones dentro/entre grupos ni incorporar predictores grupales. Esos son los temas de la semana 3.

2.15 Aplicación: heterogeneidad entre personas en sleepstudy

Como aplicación corta usaremos sleepstudy, el conjunto de datos empleado también en el estudio de caso de Bayesian Workflow. La respuesta Reaction mide tiempo de reacción y Days representa días de privación de sueño. Hay observaciones repetidas dentro de sujetos.

Esta aplicación se utiliza aquí solo para estudiar interceptos variables. No pretende ser todavía un análisis longitudinal completo. En particular, permitiremos un intercepto por sujeto pero mantendremos una pendiente común para Days; las pendientes variables y la dependencia residual se tratarán más adelante.

2.15.1 Inspección de los datos

data(
  "sleepstudy",
  package = "lme4"
)

sleep <- as_tibble(sleepstudy) |>
  mutate(
    Subject = factor(Subject)
  )

stopifnot(
  !anyNA(sleep),
  nlevels(sleep$Subject) == 18
)

sleep |>
  count(Subject) |>
  summarise(
    sujetos = n(),
    n_min = min(n),
    n_mediana = median(n),
    n_max = max(n)
  )

# A tibble: 1 × 4
  sujetos n_min n_mediana n_max
    <int> <int>     <dbl> <int>
1      18    10        10    10

ggplot(
  sleep,
  aes(
    x = Days,
    y = Reaction,
    group = Subject
  )
) +
  geom_line(alpha = 0.55) +
  geom_point(alpha = 0.65) +
  labs(
    x = "Días de privación de sueño",
    y = "Tiempo de reacción"
  ) +
  theme_minimal(base_size = 12)

Figura 2.12: Trayectorias observadas de tiempo de reacción por sujeto. En esta semana modelaremos diferencias de intercepto, manteniendo una pendiente poblacional común.

2.15.2 Modelo generativo

Usamos

\[ Reaction_{ij} \sim \mathcal N( \alpha_j+\beta\,Days_{ij}, \sigma^2 ), \]

\[ \alpha_j \sim \mathcal N( \mu_\alpha, \tau_\alpha^2 ). \]

Siguiendo la especificación pedagógica de Bayesian Workflow para este modelo, usamos previas explícitas para el intercepto, la pendiente, la desviación estándar entre personas y la desviación estándar residual (Gelman et al. 2026, 280-82):

\[ \begin{aligned} \mu_\alpha &\sim \mathcal N(250,100^2),\\ \beta &\sim \mathcal N(0,20^2),\\ \tau_\alpha &\sim \operatorname{Exponential}(0.02),\\ \sigma &\sim \operatorname{Exponential}(0.02). \end{aligned} \]

prior_sleep <- c(
  prior(
    normal(250, 100),
    class = "Intercept"
  ),
  prior(
    normal(0, 20),
    class = "b",
    coef = "Days"
  ),
  prior(
    exponential(0.02),
    class = "sd",
    group = "Subject"
  ),
  prior(
    exponential(0.02),
    class = "sigma"
  )
)

fit_sleep <- brm(
  Reaction ~ 1 + Days + (1 | Subject),
  data = sleep,
  family = gaussian(),
  prior = prior_sleep,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  seed = 1653,
  backend = "cmdstanr",
  control = list(adapt_delta = 0.95),
  file = "_fits/semana02_sleepstudy",
  file_refit = "on_change",
  refresh = 0
)

2.15.3 Resumen posterior

draws_sleep <- posterior::as_draws_df(
  fit_sleep
)

posterior::summarise_draws(
  draws_sleep,
  "mean",
  "sd",
  "rhat",
  "ess_bulk",
  "ess_tail"
) |>
  filter(
    variable %in% c(
      "b_Intercept",
      "b_Days",
      "sd_Subject__Intercept",
      "sigma"
    )
  )

# A tibble: 4 × 6
  variable               mean     sd  rhat ess_bulk ess_tail
  <chr>                 <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 b_Intercept           252.  10.1    1.00     683.     873.
2 b_Days                 10.4  0.817  1.00    3176.    2721.
3 sd_Subject__Intercept  38.8  7.17   1.01     810.    1614.
4 sigma                  31.2  1.74   1.00    2532.    2634.

La pendiente Days describe el cambio medio esperado por día bajo el supuesto de que todos los sujetos comparten la misma pendiente. La heterogeneidad entre sujetos en su nivel basal queda representada por sd_Subject__Intercept.

2.15.4 ICC condicional en Days

icc_sleep <- draws_sleep |>
  transmute(
    icc = sd_Subject__Intercept^2 /
      (
        sd_Subject__Intercept^2 +
          sigma^2
      )
  )

icc_sleep |>
  summarise(
    mediana = median(icc),
    q05 = quantile(icc, 0.05),
    q95 = quantile(icc, 0.95)
  )

# A tibble: 1 × 3
  mediana   q05   q95
    <dbl> <dbl> <dbl>
1   0.598 0.448 0.739

En este modelo, el ICC describe dependencia residual entre observaciones del mismo sujeto después de incluir una pendiente común para Days.

2.15.5 Comprobación predictiva

pp_check(
  fit_sleep,
  type = "dens_overlay",
  ndraws = 50
)

Figura 2.13: Comprobación predictiva posterior global del modelo de interceptos variables para sleepstudy.

Este PPC global no debe cerrar el análisis. La figura de trayectorias sugiere que la asociación con Days también podría variar entre sujetos. Esa observación motiva naturalmente la semana de pendientes variables.

AdvertenciaNo adelantemos la conclusión longitudinal

El modelo Reaction ~ Days + (1 | Subject) es útil para estudiar heterogeneidad de interceptos, pero no demuestra que una pendiente común ni errores condicionalmente independientes sean suficientes. Más adelante ampliaremos la estructura.

2.16 Perspectiva frecuentista y terminología

En la literatura frecuentista, el modelo

\[ Y_{ij} = \mu_\alpha+u_j+\varepsilon_{ij} \]

se denomina con frecuencia modelo de intercepto aleatorio. Los \(u_j\) pueden predecirse mediante BLUP/EBLUP y esas predicciones también muestran contracción hacia cero, equivalente a contraer los interceptos hacia \(\mu_\alpha\).

La diferencia principal para nuestros propósitos no es que el shrinkage sea exclusivamente bayesiano. También surge en el modelo mixto frecuentista. La ventaja pedagógica del enfoque bayesiano aquí es que podemos tratar de forma unificada:

• incertidumbre de \(\tau_\alpha\) y \(\sigma\);
• incertidumbre de cada \(\alpha_j\);
• posterior del ICC;
• posterior de la cantidad de pooling;
• predicción para nuevos grupos.

El curso seguirá usando la terminología frecuentista cuando facilite la lectura de la literatura, pero organizará la inferencia alrededor del modelo generativo y las distribuciones predictivas.

2.17 ¿Cuántos grupos necesitamos?

No existe un número mágico de grupos a partir del cual un modelo multinivel “se vuelve válido”.

Gelman y Hill señalan que incluso grupos con muy pocas observaciones pueden contribuir información al modelo; la dificultad principal cuando \(J\) es pequeño es aprender la variación entre grupos con precisión (Gelman y Hill 2007, 275-76).

Esto conduce a una regla conceptual importante:

el número total de observaciones y el número de grupos informan componentes diferentes del modelo.

Cientos de observaciones dentro de tres grupos pueden estimar con precisión ciertos aspectos de la variación dentro de grupos, pero siguen proporcionando muy poca información directa para caracterizar una población de interceptos.

Con pocos grupos debemos prestar atención especial a:

• la prior de \(\tau_\alpha\);
• la incertidumbre posterior de \(\tau_\alpha\);
• la sensibilidad del ICC;
• el propósito de inferir o predecir para grupos nuevos;
• la posibilidad de que el modelo sea más ambicioso que la información disponible.

No interpretaremos “pocos grupos” como una prohibición automática. Lo trataremos como un problema de información e identificación.

2.18 Errores frecuentes de interpretación

“El ICC es el porcentaje de observaciones que son iguales dentro de un grupo.”
No. En el modelo nulo gaussiano es una razón de componentes de varianza y también la correlación marginal esperada de dos observaciones distintas del mismo grupo.

“Si el ICC es pequeño, puedo ignorar el agrupamiento.”
No necesariamente. La relevancia depende del estimando, el diseño, los tamaños grupales y la estructura completa del modelo.

“Shrinkage significa sesgar artificialmente todas las medias hacia la media global.”
El shrinkage es consecuencia del modelo jerárquico y combina dos fuentes de información. En escenarios compatibles con la estructura poblacional, funciona como regularización.

“Los grupos pequeños siempre se contraen la misma distancia.”
No. El tamaño afecta el peso \(w_j\), pero la distancia observada a \(\mu_\alpha\), la variación residual, la heterogeneidad entre grupos y la incertidumbre posterior también importan.

“Si dos grupos tienen el mismo \(n_j\), tienen exactamente el mismo shrinkage.”
Condicionalmente en hiperparámetros conocidos tienen el mismo peso \(w_j\), pero no necesariamente la misma distancia de contracción, porque sus medias observadas pueden estar a distancias diferentes de \(\mu_\alpha\).

“La media posterior de un grupo debe coincidir con su media observada si el modelo ajusta bien.”
No. Un modelo jerárquico bien ajustado puede contraer sistemáticamente estimaciones grupales. La reproducción exacta de medias observadas no es el objetivo.

“El ICC es una propiedad fija del conjunto de datos.”
No. Cambia con la familia, los predictores y la estructura del modelo. Debe interpretarse en relación con la especificación utilizada.

“Tener 1000 observaciones compensa tener solamente tres grupos.”
No para todos los parámetros. \(\tau_\alpha\) depende de la variación observada entre grupos, y la información relevante para esa distribución poblacional está limitada por \(J\).

“\(\widehat R=1\) valida la existencia de heterogeneidad entre grupos.”
No. \(\widehat R\) diagnostica el muestreo, no el supuesto sustantivo de una distribución poblacional de interceptos.

2.19 Síntesis

El modelo nulo de interceptos variables

\[ y_{ij}\sim \mathcal N(\alpha_j,\sigma^2), \qquad \alpha_j\sim \mathcal N(\mu_\alpha,\tau_\alpha^2) \]

separa la variación en dos escalas.

Marginalmente,

\[ E(Y_{ij})=\mu_\alpha, \]

\[ \operatorname{Var}(Y_{ij}) = \tau_\alpha^2+\sigma^2, \]

y, para dos observaciones del mismo grupo,

\[ \operatorname{Cov}(Y_{ij},Y_{i'j}) = \tau_\alpha^2. \]

Por tanto,

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

Condicionando en los hiperparámetros, la media posterior del intercepto grupal es

\[ \boxed{ E( \alpha_j \mid \mathbf y_j, \mu_\alpha, \tau_\alpha, \sigma ) = w_j\bar y_j + (1-w_j)\mu_\alpha } \]

con

\[ \boxed{ w_j = \frac{n_j\tau_\alpha^2} {n_j\tau_\alpha^2+\sigma^2} = \frac{n_j\rho} {1+(n_j-1)\rho} }. \]

Estas ecuaciones resumen la lógica de la semana:

• grupos grandes dependen más de su propia información;
• grupos pequeños reciben mayor regularización;
• mayor ruido residual aumenta el pooling;
• mayor heterogeneidad real entre grupos disminuye el pooling;
• el ICC y el shrinkage están relacionados, pero no son equivalentes;
• la incertidumbre sobre las componentes de varianza debe propagarse;
• predecir un grupo nuevo exige generar un nuevo parámetro grupal desde la distribución poblacional.

La próxima semana incorporará predictores individuales y grupales y mostrará por qué una asociación observada entre todos los datos puede mezclar comparaciones dentro y entre grupos.

2.20 Ejercicios

2.20.1 Ejercicios conceptuales

2.20.1.1 Ejercicio 1. Derivación del ICC

Considere

\[ Y_{ij} = \mu+u_j+\varepsilon_{ij}, \]

con

\[ u_j\sim\mathcal N(0,\tau^2), \qquad \varepsilon_{ij}\sim\mathcal N(0,\sigma^2). \]

Suponga independencia entre todos los componentes salvo por el \(u_j\) compartido.

1. Derive \(E(Y_{ij})\).
2. Derive \(\operatorname{Var}(Y_{ij})\).
3. Derive \(\operatorname{Cov}(Y_{ij},Y_{i'j})\) para \(i\neq i'\).
4. Derive \(\operatorname{Cov}(Y_{ij},Y_{i'k})\) para \(j\neq k\).
5. Obtenga \(\operatorname{Corr}(Y_{ij},Y_{i'j})\).
6. Interprete el resultado si \(\tau=4\) y \(\sigma=8\).

2.20.1.2 Ejercicio 2. El peso de shrinkage

Suponga

\[ \tau_\alpha=6, \qquad \sigma=12. \]

Calcule

\[ w_j = \frac{n_j\tau_\alpha^2} {n_j\tau_\alpha^2+\sigma^2} \]

para

\[ n_j\in\{2,5,10,25,100\}. \]

1. Grafique \(w_j\) contra \(n_j\).
2. Explique por qué \(w_j\) no es lineal en \(n_j\).
3. Calcule el ICC.
4. Verifique numéricamente que

\[ w_j = \frac{n_j\rho} {1+(n_j-1)\rho}. \]

2.20.1.3 Ejercicio 3. Mismo peso, distinta distancia

Dos grupos tienen \(n_1=n_2=10\). Suponga

\[ \mu_\alpha=50, \qquad w_1=w_2=0.60. \]

Sus medias observadas son

\[ \bar y_1=52, \qquad \bar y_2=70. \]

1. Calcule la media posterior condicional de cada \(\alpha_j\).
2. Calcule la distancia de contracción \(|\bar y_j-m_j|\).
3. Explique por qué “mismo peso de shrinkage” no significa “misma distancia de shrinkage”.

2.20.1.4 Ejercicio 4. Criticar cuatro afirmaciones

Explique por qué cada afirmación es incompleta o incorrecta.

1. “El ICC de estos datos es 0.20.”
2. “Como el ICC es solo 0.05, no necesitamos un modelo multinivel.”
3. “Shrinkage es un problema porque introduce sesgo en las medias de grupos pequeños.”
4. “Tenemos 800 pacientes, por lo tanto podemos estimar con precisión la variación entre los cuatro hospitales.”

2.20.2 Ejercicios computacionales

2.20.2.1 Ejercicio 5. Recuperación del ICC por simulación

Use el código del laboratorio y considere tres procesos generadores:

\[ (\tau_\alpha,\sigma) \in \{ (2,12), (7,12), (15,12) \}. \]

Para cada escenario:

1. calcule el ICC verdadero;
2. simule los datos con los mismos \(n_j\);
3. ajuste el modelo de interceptos variables;
4. calcule la posterior del ICC;
5. presente mediana e intervalo central del 90%;
6. discuta cómo cambia el shrinkage.

Repita el experimento con una segunda semilla y explique por qué no debe esperarse que el intervalo cubra siempre el valor verdadero en cada simulación particular.

2.20.2.2 Ejercicio 6. Mantener \(N\), cambiar \(J\)

Construya dos diseños con aproximadamente el mismo número total de observaciones:

• Diseño A: pocos grupos grandes;
• Diseño B: muchos grupos pequeños.

Mantenga el mismo proceso generador.

1. Ajuste el mismo modelo a ambos diseños.
2. Compare la incertidumbre posterior de \(\sigma\).
3. Compare la incertidumbre posterior de \(\tau_\alpha\).
4. Compare la incertidumbre posterior del ICC.
5. Explique por qué el número total de filas no resume por sí solo la información jerárquica.

2.20.2.3 Ejercicio 7. Predicción para un grupo nuevo

Usando fit_null:

1. obtenga draws de la media esperada de G01;
2. obtenga draws de una nueva observación de G01;
3. obtenga draws de la media de G_nuevo;
4. obtenga draws de una observación de G_nuevo;
5. compare las cuatro desviaciones estándar posteriores;
6. explique de dónde proviene cada componente adicional de incertidumbre.

2.20.3 Ejercicio de interpretación de salida

Un modelo produce el siguiente resumen hipotético:

Cantidad	Mediana posterior	Intervalo 90%
\(\mu_\alpha\)	72	[68, 76]
\(\tau_\alpha\)	9	[4, 16]
\(\sigma\)	18	[16, 20]
ICC	0.20	[0.05, 0.45]

Responda:

1. ¿qué describe \(\tau_\alpha=9\) que no describe \(\sigma=18\)?
2. ¿por qué el intervalo del ICC puede ser relativamente ancho?
3. ¿sería correcto afirmar que “exactamente 20% de la variación pertenece a los grupos”?
4. ¿qué esperaría sobre el shrinkage de un grupo con \(n_j=4\) comparado con uno de \(n_j=60\)?
5. ¿qué información adicional necesitaría para juzgar si la distribución normal de interceptos es adecuada?

2.21 Lecturas para profundizar

Para esta semana se recomienda priorizar:

• Gelman y Hill, cap. 12, especialmente §§12.2–12.5 y §§12.8–12.9: pooling parcial, componentes de variación, predicción y tamaños de grupo (Gelman y Hill 2007, 251-76).
• McElreath, §§13.1–13.2 y §13.5: shrinkage, regularización adaptativa y predicción multinivel (McElreath 2020, 401-13; 2020, 426-31).
• Gelman et al., BDA3, §§5.3–5.5: análisis bayesiano completo de modelos jerárquicos normales y propagación de incertidumbre (Gelman et al. 2013, 108-23).

Como complemento:

• Hox, Moerbeek y van de Schoot, cap. 2: modelo vacío, componentes de varianza e ICC (Hox et al. 2018, 8-19).
• Gelman et al., Bayesian Workflow, §17.2: previas y especificación en brms de modelos lineales multinivel (Gelman et al. 2026, 280-83).
• Bürkner: estructura general de parámetros poblacionales y grupales en brms (Bürkner 2017, 2-5).

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, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, y Donald B. Rubin. 2013. Bayesian Data Analysis. 3.ª ed. Chapman & Hall/CRC.

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.