Especificación de un LMM (Linear Mixed Model)

Del análisis exploratorio a la especificación del modelo

El análisis exploratorio del módulo anterior reveló dos señales claras en sleepstudy:

1. Los interceptos varían entre sujetos: en el día 0, el tiempo de reacción se extiende desde aproximadamente 200 ms hasta más de 310 ms. Cada sujeto tiene un nivel basal propio.
2. Las pendientes también varían: algunos sujetos se deterioran apenas 5 ms por día de privación, mientras que otros empeoran a un ritmo de 25 ms por día.

Estas dos observaciones empíricas se traducen directamente en decisiones de especificación. Como señalan Correa y Salazar (2025) en su texto sobre modelos mixtos, "la estructura de los efectos aleatorios debe reflejar las fuentes de variación identificadas en la etapa exploratoria; no hay una estructura universalmente correcta, sino una que es adecuada para los datos y la pregunta de investigación".

Esta sección formaliza ese razonamiento. Seguimos de cerca la exposición de Fitzmaurice, Laird y Ware (2011, cap. 8), quienes enfatizan que la especificación de un modelo mixto es un proceso iterativo que combina criterios estadísticos, gráficos y teóricos.

---

Las tres variantes fundamentales de la estructura aleatoria

Un LMM con respuesta $y_{ij}$ para el sujeto $i$ en el tiempo $j$ y un predictor continuo $t_{ij}$ admite tres configuraciones comunes de efectos aleatorios. Cada una responde a una pregunta diferente sobre la heterogeneidad de la población. Siguiendo la notación de Correa y Salazar (2025, cap. 5), presentamos las variantes en orden creciente de complejidad.

---

Variante A: solo intercepto aleatorio

Qué dice conceptualmente

Este modelo asume que todos los sujetos comparten la misma pendiente (la misma tasa de cambio a lo largo del tiempo), pero cada uno parte de un nivel basal distinto. La heterogeneidad entre sujetos es exclusivamente de posición: los perfiles individuales son líneas paralelas, desplazadas verticalmente unas de otras.

$$y_{ij} = (\beta_0 + b_{0i}) + \beta_1\, t_{ij} + \varepsilon_{ij}, \qquad b_{0i} \sim N(0, \tau_0^2), \qquad \varepsilon_{ij} \sim N(0, \sigma^2)$$

La varianza total de $y_{ij}$ se descompone en dos fuentes:

• $\tau_0^2$: varianza entre sujetos (cuánto difieren los niveles basales).
• $\sigma^2$: varianza intra-sujeto (ruido residual no explicado).

Una cantidad derivada de gran utilidad es el coeficiente de correlación intraclase (ICC, por sus siglas en inglés), que mide la proporción de la varianza total atribuible a diferencias entre sujetos:

$$\text{ICC} = \frac{\tau_0^2}{\tau_0^2 + \sigma^2}$$

Correa y Salazar (2025, cap. 4) indican que un ICC superior a 0.10 ya justifica la inclusión del efecto aleatorio, pues refleja que al menos el 10 % de la variabilidad observada proviene de diferencias estables entre las unidades. Valores superiores a 0.50 son frecuentes en estudios con diseños cuidadosamente controlados donde la heterogeneidad entre individuos domina sobre el error de medición.

Cuándo está justificada

• La exploración gráfica muestra perfiles aproximadamente paralelos: las líneas del spaghetti plot tienen pendientes similares y mantienen distancias relativas constantes.
• La teoría no sugiere que las unidades deban tener tasas de cambio distintas. Por ejemplo, si el efecto del tiempo es el mismo para todos los pacientes de un grupo homogéneo, una pendiente aleatoria no estaría respaldada conceptualmente, aunque los datos la sugirieran débilmente.

Fitzmaurice et al. (2011, p. 192) advierten que la variante de solo intercepto aleatorio es a menudo un excelente punto de partida, y que conviene resistir la tentación de añadir pendientes aleatorias sin una justificación clara.

Código R y SAS

---

Variante B: intercepto y pendiente aleatorios (correlación libre)

Qué añade conceptualmente

Este modelo reconoce que las unidades difieren tanto en su nivel basal como en su tasa de cambio. Además, permite que ambas características estén relacionadas en la población. La pregunta sustantiva es: ¿los sujetos que parten de un valor más alto tienden también a cambiar más rápido, o por el contrario los que parten más alto cambian menos?

$$y_{ij} = (\beta_0 + b_{0i}) + (\beta_1 + b_{1i})\, t_{ij} + \varepsilon_{ij}$$ $$\begin{pmatrix} b_{0i} \\ b_{1i} \end{pmatrix} \sim N\!\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}\tau_0^2 & \tau_{01} \\ \tau_{01} & \tau_1^2\end{pmatrix} \right)$$

La correlación entre efectos aleatorios, definida como $\rho_{01} = \tau_{01}/(\tau_0\,\tau_1)$, tiene una interpretación sustantiva directa. Un $\rho_{01}$ positivo indica que los sujetos que parten más alto también cambian más rápido (en el caso de sleepstudy, se deterioran más). Un $\rho_{01}$ negativo sugiere convergencia: los que parten más alto cambian menos a lo largo del tiempo, acercándose a los demás.

Correa y Salazar (2025, cap. 5) dedican una sección completa a la interpretación de esta correlación, señalando que en estudios de crecimiento es frecuente encontrar correlaciones negativas (los individuos que parten más altos crecen más lentamente), mientras que en estudios de deterioro clínico pueden aparecer correlaciones positivas (los más afectados al inicio empeoran más rápido).

Cuándo está justificada

• El spaghetti plot muestra trayectorias que divergen o se cruzan visiblemente.
• El análisis exploratorio reveló que la correlación individual entre el predictor temporal y la respuesta varía sustancialmente entre sujetos.
• Hay suficientes mediciones por sujeto para estimar dos parámetros aleatorios más su correlación. La regla práctica citada por Fitzmaurice et al. (2011, p. 200) es disponer de al menos 5 mediciones por unidad.
• La teoría sugiere que el punto de partida y la velocidad de cambio podrían estar relacionados.

Código R y SAS

---

Variante C: intercepto y pendiente aleatorios (correlación forzada a 0)

Por qué forzar la correlación a cero

Estimar la correlación $\rho_{01}$ consume un parámetro adicional y, cuando hay pocos sujetos o pocas mediciones por sujeto, puede generar inestabilidad numérica. El mensaje boundary (singular) fit que a veces emite lmer suele estar relacionado con una matriz $\mathbf{D}$ difícil de estimar. En estas situaciones, forzar $\tau_{01} = 0$ produce un modelo más estable sin renunciar a la heterogeneidad en las pendientes:

$$\mathbf{D} = \begin{pmatrix}\tau_0^2 & 0 \\ 0 & \tau_1^2\end{pmatrix}$$

Esta restricción no niega que los efectos aleatorios existan; simplemente asume que el nivel basal y la tasa de cambio de un sujeto no están asociados en la población bajo estudio.

Código R y SAS

---

Cómo elegir entre variantes: criterios de datos y teoría

La elección de la estructura aleatoria no es un problema puramente estadístico. Fitzmaurice et al. (2011, cap. 8) insisten en que debe combinarse el juicio teórico con la evidencia empírica. Correa y Salazar (2025, cap. 5) proponen un esquema de decisión en tres planos que adaptamos a continuación.

1. 1

   Plano teórico: preguntarse si la heterogeneidad tiene sentido

   ¿Es razonable que cada sujeto tenga su propia tasa de cambio? Si se mide el efecto de un fármaco en pacientes con la misma condición basal, quizás la pendiente aleatoria no está justificada conceptualmente, aunque los datos la sugieran. La teoría debe guiar la especificación; los datos, validarla.

2. 2

   Plano empírico: leer el spaghetti plot y la variabilidad de pendientes

   Perfiles paralelos → solo intercepto aleatorio (Variante A). Perfiles que divergen o se cruzan → considere pendiente aleatoria (Variante B o C). Complemente con el gráfico de pendientes individuales del módulo de exploración.

3. 3

   Plano estadístico: comparar modelos anidados

   Use anova(mod_A, mod_B, refit=FALSE) para una prueba de razón de verosimilitud (LRT) entre variantes anidadas. Complemente con AIC y BIC. Si el modelo más complejo no mejora significativamente, quédese con el más simple por parsimonia.

4. 4

   Plano numérico: verificar convergencia y estabilidad

   Si lmer() emite una advertencia de ajuste singular o de convergencia fallida con la Variante B, intente primero la Variante C (menos parámetros), luego cambie el optimizador a bobyqa, y como último recurso escale los predictores para mejorar el condicionamiento numérico.

Comparación estadística entre variantes

Cómo leer la comparación:

• El LRT (estadístico $\chi^2$ y valor p en la salida de anova()) contrasta la hipótesis nula de que el modelo más simple es suficiente. Un valor $p < 0.05$ favorece la variante más compleja.
  Correa y Salazar (2025, cap. 5) advierten que el LRT para componentes de varianza es conservador cuando el parámetro está en la frontera del espacio paramétrico ($\tau_1^2 = 0$), y recomiendan considerar también los criterios de información.
• El AIC penaliza la verosimilitud por el número de parámetros. El modelo con menor AIC es preferible cuando existe incertidumbre sobre la estructura correcta.
• El BIC penaliza más fuertemente la complejidad y suele favorecer modelos más parsimoniosos que el AIC. En muestras grandes, el BIC tiende a seleccionar el modelo verdadero si está entre los candidatos.
• Si el LRT y el AIC apuntan al modelo más complejo pero el BIC favorece al más simple, la decisión debe apoyarse en la teoría y en los patrones observados en los gráficos exploratorios.

---

Guía completa de lectura del output de summary(lmer())

La salida de summary() de un modelo ajustado con lmer contiene cinco bloques de información. Saber leer cada uno es una competencia fundamental. Fitzmaurice et al. (2011, cap. 8) dedican un apartado completo a la interpretación de salidas de software, y adaptamos aquí sus recomendaciones.

La salida se estructura en cinco bloques. A continuación se detalla cada uno.

Bloque 1 — Fórmula y método de estimación

Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (Days | Subject)
Data: sleepstudy
REML criterion at convergence: 1743.6

REML criterion: la deviance REML, es decir, $-2 \times$ log-verosimilitud restringida. Solo es comparable entre modelos que comparten los mismos efectos fijos. No debe usarse para comparar modelos con distintas variables predictoras. La línea REML criterion at convergence indica que el algoritmo de optimización convergió exitosamente. Si apareciera un mensaje de advertencia aquí, habría que revisar la especificación.

Bloque 2 — Residuos escalados

Scaled residuals:
    Min      1Q  Median      3Q     Max
 -3.9536 -0.4634  0.0231  0.4634  5.1793

Los residuos se presentan escalados (divididos por su desviación estándar estimada). Si el modelo es adecuado:

• La mediana debería estar próxima a 0.
• El primer y tercer cuartil deberían ser aproximadamente simétricos ($\approx \pm 0.67$ para una distribución normal).
• El mínimo y el máximo no deberían alejarse excesivamente de $\pm 3.5$. Valores extremos (como el 5.18 que aparece aquí) pueden indicar la presencia de observaciones atípicas que conviene examinar.

Bloque 3 — Efectos aleatorios ← el corazón del LMM

Random effects:
  Groups Name Variance Std.Dev. Corr
  Subject (Intercept) 612.10 24.741
  Days 35.07 5.922 0.07
  Residual 654.94 25.592
  Number of obs: 180, groups: Subject, 18

La desviación estándar (Std.Dev.) es la raíz cuadrada de la varianza y se expresa en las mismas unidades que la respuesta. Por ejemplo, Std.Dev. (Intercept) = 24.7 ms significa que, aproximadamente, el 95 % de los sujetos tienen su intercepto dentro de $\pm 1.96 \times 24.7 \approx \pm 48$ ms alrededor del intercepto medio poblacional.

Correa y Salazar (2025, cap. 5) recomiendan prestar especial atención a la magnitud relativa de $\tau_1$ respecto a $\beta_1$: si la desviación estándar de las pendientes es grande en comparación con la pendiente media, existe una heterogeneidad sustancial en las trayectorias y el modelo con pendiente aleatoria está plenamente justificado.

Bloque 4 — Efectos fijos

Fixed effects:
  Estimate Std. Error t value
  (Intercept) 251.405 6.825 36.84
  Days 10.467 1.546 6.77

• Intercept (251.4 ms): tiempo de reacción promedio poblacional en el día 0 de privación.
• Days (10.5 ms/día): por cada día adicional de privación de sueño, el tiempo de reacción promedio aumenta 10.5 ms.
• Std. Error: error estándar del estimador, que cuantifica la incertidumbre asociada.
• t value: el cociente entre la estimación y su error estándar. En muestras razonablemente grandes, valores superiores a 2 en valor absoluto indican significación estadística aproximada.

lme4 no reporta valores p para los efectos fijos por defecto. Esto obedece a un debate metodológico sobre los grados de libertad apropiados en modelos con estructura aleatoria compleja o datos desbalanceados. Para obtener valores p, existen tres alternativas:

1. Usar lmerTest::lmer(), que incorpora la aproximación de Satterthwaite.
2. Calcular intervalos de confianza con confint(mod_B, method = "Wald").
3. Usar la aproximación de Kenward-Roger disponible en el paquete pbkrtest.

Bloque 5 — Correlación de los efectos fijos

Correlation of Fixed Effects:
  (Intr)
  Days -0.138

Esta es la correlación entre los estimadores $\hat{\beta}_0$ y $\hat{\beta}_1$, no entre los efectos aleatorios. Un valor cercano a 0 (como −0.14 en este caso) indica que las estimaciones son casi ortogonales entre sí. Una correlación alta en valor absoluto (> 0.7) sugeriría que centrar el predictor temporal (por ejemplo, usar I(Days - mean(Days)) en lugar de Days) mejoraría la interpretación y la estabilidad numérica del modelo.

---

REML vs ML: regla práctica

Fitzmaurice et al. (2011, cap. 4) presentan una discusión accesible sobre la diferencia entre máxima verosimilitud (ML) y máxima verosimilitud restringida (REML). La tabla siguiente resume las recomendaciones prácticas.

---

Valores ajustados y residuales

Una vez ajustado el modelo, conviene extraer los valores ajustados, los residuos y las predicciones de los efectos aleatorios (BLUPs). Fitzmaurice et al. (2011, p. 203) recomiendan examinar gráficamente los residuos para verificar los supuestos del modelo antes de proceder a la inferencia.

Qué observar en este gráfico de ajustes individuales:

• Las rectas no son paralelas: cada sujeto tiene su propio intercepto y su propia pendiente. Esto refleja la estructura de la Variante B y permite apreciar visualmente lo que aporta la inclusión de la pendiente aleatoria.
• Los sujetos cuya recta sube poco (pendiente cercana a 0 o incluso ligeramente negativa) son los más resistentes a la privación de sueño. Aquellos con rectas muy inclinadas son los más vulnerables.
• Compare mentalmente este gráfico con el que se obtendría de la Variante A ((1 | Subject), rectas paralelas) para apreciar la diferencia. En el módulo de exploración gráfica se mostró cómo superponer ambas representaciones.

Correa y Salazar (2025, cap. 5) dedican una sección a la interpretación de los BLUPs, señalando que no deben tratarse como estimaciones insesgadas de los efectos individuales (tienen shrinkage hacia la media poblacional), pero sí como herramientas valiosas para la validación del modelo y la identificación de sujetos con trayectorias atípicas.

---

Conceptos del módulo

Concepto · clic para ver la respuesta

¿Qué diferencia hay entre (Days | Subject) y (1 | Subject) + (0 + Days | Subject)?

Concepto · clic para ver la respuesta

¿Qué es el ICC y qué valor indica que el modelo mixto es necesario?

Concepto · clic para ver la respuesta

¿Por qué lme4 no reporta valores p para los efectos fijos?

Concepto · clic para ver la respuesta

¿Qué son los BLUPs y para qué sirven?

Autoevaluación

1. En la salida de summary(lmer()), ¿qué representa 'Std.Dev. (Days, Subject)'?

El error estándar del coeficiente fijo de DaysLa desviación estándar de los efectos aleatorios de pendiente entre sujetos (√τ₁²)El residual estándar del modeloLa correlación entre los efectos aleatorios

2. ¿Cuándo es correcto usar anova(mod_A, mod_B, refit = FALSE)?

Cuando los modelos tienen distintos efectos fijosCuando los modelos tienen los mismos efectos fijos pero distinta estructura aleatoriaSiempre que se quiera comparar dos modelos mixtosSolo cuando ambos modelos se ajustaron con ML

3. Un spaghetti plot muestra líneas que parten de niveles distintos pero con pendientes muy similares. ¿Qué variante es la más adecuada por parsimonia?

Variante B: intercepto y pendiente aleatorios con correlación libreVariante C: intercepto y pendiente aleatorios sin correlaciónVariante A: solo intercepto aleatorioRegresión clásica sin efectos aleatorios

---

Apéndice A — Bayes empírico y el fenómeno de shrinkage

Los BLUPs que obtuviste con ranef() en las secciones anteriores no son solo "residuales de grupo": son predictores bayesianos empíricos. Esta sección explica por qué, y cómo visualizar el efecto de contracción (shrinkage) hacia la media global.

Los BLUPs como media a posteriori

En el modelo con intercepto aleatorio $b_{0i} \sim \mathcal{N}(0, \tau_0^2)$, la distribución a posteriori de $b_{0i}$ dado $\mathbf{y}_i$ es normal con media:

$$\hat{b}_{0i}^{\text{BLUP}} = \underbrace{\frac{\tau_0^2}{\tau_0^2 + \sigma^2/n_i}}_{\lambda_i} \cdot (\bar{y}_i - \hat{\beta}_0)$$

El factor $\lambda_i \in (0,1)$ es el peso de shrinkage del sujeto $i$:

• Si $\tau_0^2 \gg \sigma^2/n_i$ (mucha variabilidad entre sujetos, o sujeto bien medido): $\lambda_i \to 1$ y el BLUP se acerca a la media cruda del sujeto.
• Si $\tau_0^2 \ll \sigma^2/n_i$ (poca variabilidad entre sujetos, o sujeto poco medido): $\lambda_i \to 0$ y el BLUP se contrae hacia la media global $\hat\beta_0$.

Este es exactamente el estimador de Bayes empírico bajo pérdida cuadrática, y también coincide con las Ecuaciones de Henderson (1950).

---

Apéndice B — El algoritmo EM: convergencia visual paso a paso

El algoritmo EM (Expectation-Maximization, Dempster, Laird y Rubin, 1977) es uno de los métodos clásicos para estimar modelos con datos incompletos o variables latentes. En los LMM, los efectos aleatorios $\mathbf{b}_i$ son las "variables latentes" que el algoritmo EM maneja de manera implícita.

Los dos pasos del algoritmo

Paso E (Expectation): calcula $Q(\theta | \theta^{(t)}) = \mathbb{E}_{\mathbf{b}|\mathbf{y},\theta^{(t)}}[\log L(\theta; \mathbf{y}, \mathbf{b})]$, es decir, la verosimilitud esperada sobre la distribución condicional de los efectos aleatorios dados los datos y los parámetros actuales.

Paso M (Maximization): maximiza $Q$ respecto a $\theta$, produciendo la nueva estimación $\theta^{(t+1)}$.

Para el modelo de intercepto aleatorio balanceado ($n_i = n$ para todo $i$, $N$ sujetos), las actualizaciones cerradas del EM son:

$$\hat\sigma^{2,(t+1)} = \frac{1}{Nn} \sum_{i,j} \left[(y_{ij} - \hat\beta^{(t)} - \hat b_i^{(t)})^2 + \text{Var}(b_i | \mathbf{y})\right]$$ $$\hat\tau^{2,(t+1)} = \frac{1}{N} \sum_i \left[(\hat b_i^{(t)})^2 + \text{Var}(b_i | \mathbf{y})\right]$$

---

Referencias del módulo

• Correa, J. C. y Salazar, J. C. (2025). Modelos Mixtos. Universidad Nacional de Colombia, Sede Medellín, Facultad de Ciencias, Departamento de Estadística. [Disponible en el repositorio del curso].
• Fitzmaurice, G. M., Laird, N. M. y Ware, J. H. (2011). Applied Longitudinal Analysis (2.ª ed.). John Wiley & Sons. [Capítulos 4, 5 y 8: especificación del LMM, estimación REML/ML e interpretación de salidas].
• Pinheiro, J. y Bates, D. (2000). Mixed-Effects Models in S and S-PLUS. Springer. [Texto de referencia del paquete nlme; capítulos 1 y 2 contienen una exposición rigurosa de la notación y las estructuras de covarianza].
• Bates, D., Mächler, M., Bolker, B. y Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. [Documentación técnica de lme4; incluye recomendaciones sobre optimizadores y manejo de advertencias de convergencia].