LMM sobre trayectorias de crecimiento - ej. SPRUCE1

Estrategia: construcción progresiva

El principio rector de este módulo es añadir complejidad solo cuando los datos y la teoría lo justifican. Comenzamos con el modelo más simple posible y añadimos componentes uno a uno, verificando en cada paso si la adición mejora el ajuste de forma sustancial. Al final del módulo habremos construido cuatro modelos anidados y podremos decidir cuál es el más apropiado para los datos de SPRUCE1, siguiendo la filosofía de análisis iterativo recomendada por Fitzmaurice, Laird y Ware (2011, cap. 8) y Correa y Salazar (2025, cap. 5).

1. 1

   M0: Referencia: regresión lineal clásica

   Ignora la estructura de árboles (pseudo‑replicación). Útil solo para cuantificar cuánto mejoran los modelos mixtos.

2. 2

   M1: LMM con intercepto aleatorio

   Reconoce que cada árbol tiene su propio nivel basal. Los perfiles son paralelos desplazados verticalmente.

3. 3

   M2: LMM con intercepto y pendiente aleatorios

   Reconoce además que cada árbol crece a su propio ritmo. Los perfiles pueden divergir.

4. 4

   M3: M2 con interacción grupo × tiempo

   Añade la pregunta científica central: ¿el ozono modifica la tasa de crecimiento?

---

Datos y preparación

Cómo leer la salida:

• Estructura del dataset: str() confirma que ID y Group son factores y que time_c es numérica. Verifique que el número de observaciones coincide con n_árboles × n_tiempos (datos balanceados).
• Spaghetti plot: si las líneas son aproximadamente paralelas pero desplazadas verticalmente, la heterogeneidad es principalmente de intercepto (M1 será suficiente). Si las líneas divergen o se cruzan, hay heterogeneidad de pendiente y se justifica M2.
• Boxplot por tiempo: la dispersión constante a lo largo del tiempo sugiere que una varianza residual homogénea ($\sigma^2$) es razonable. Si la dispersión cambia, podría necesitarse una estructura de varianza residual (Módulo 11).
• ICC crudo: un valor alto (> 0.4) indica que una fracción sustancial de la variabilidad se debe a diferencias entre árboles, lo que justifica incluir un intercepto aleatorio (M1). Un ICC cercano a 0 sugeriría que el agrupamiento aporta poco.

Por qué centrar el tiempo:

Sin centrado, el intercepto $\beta_0$ representa el valor esperado de log_growth cuando time1 = 0, que puede estar fuera del rango de los datos y no tiene interpretación práctica. Con time_c = time1 − mean(time1), el intercepto corresponde al punto medio del estudio, mucho más informativo y con mejor comportamiento numérico del optimizador. Correa y Salazar (2025, cap. 4) recomiendan esta práctica para mejorar la interpretabilidad y reducir la correlación entre los estimadores de los efectos fijos.

---

M0: Referencia: regresión lineal clásica

Este modelo es estadísticamente incorrecto para datos con medidas repetidas, pero es necesario como punto de referencia para cuantificar el sesgo que introduce ignorar la estructura de agrupamiento, un paso que Fitzmaurice et al. (2011, p. 187) consideran didácticamente útil.

$$\text{log\_growth}_{ij} = \beta_0 + \beta_1\,\text{time\_c}_{ij} + \beta_2\,\text{Group}_i + \beta_3\,(\text{time\_c} \times \text{Group})_{ij} + \varepsilon_{ij}, \quad \varepsilon_{ij} \overset{iid}{\sim} N(0,\sigma^2)$$

Qué evidencia el diagnóstico:

• La correlación serial de los residuales (correlación entre el residual de un árbol en el tiempo $t$ y en el tiempo $t+1$) debería ser ≈0 si OLS fuera correcto. Si la media es positiva y relevante (> 0.1), confirma que los residuales están correlacionados dentro de cada árbol y que OLS subestima los errores estándar.
• Los valores p del modelo M0 son anticonservadores: el tamaño de muestra efectivo no son todas las observaciones, sino el número de árboles (unidades independientes).

---

M1: LMM con intercepto aleatorio

Qué añade M1 respecto a M0

M1 reconoce que cada árbol tiene su propio nivel basal de crecimiento (intercepto aleatorio $b_{0i}$), independientemente del tiempo o del grupo. Esta heterogeneidad ya era visible en el Módulo 2: en el boxplot temporal, la dispersión dentro de cada tiempo era considerable y relativamente constante.

$$\text{log\_growth}_{ij} = \underbrace{\beta_0 + \beta_1\,\text{time\_c}_{ij} + \beta_2\,\text{Group}_i + \beta_3\,(\text{time\_c} \times \text{Group})_{ij}}_{\text{efectos fijos}} + \underbrace{b_{0i}}_{\text{intercepto aleatorio}} + \varepsilon_{ij}$$ $$b_{0i} \sim N(0, \tau_0^2), \qquad \varepsilon_{ij} \sim N(0,\sigma^2)$$

Lo que resuelve: absorbe la correlación intra‑árbol debida a diferencias basales, produciendo errores estándar de los efectos fijos más honestos.

Lo que supone: todos los árboles crecen al mismo ritmo (misma pendiente $\beta_1 + \beta_3$); solo difieren en su punto de partida.

Correa y Salazar (2025, cap. 5) señalan que el modelo de solo intercepto aleatorio suele ser un excelente punto de partida, y que su simplicidad facilita la comunicación de resultados cuando la heterogeneidad en las pendientes es pequeña.

Código R y SAS

Cómo interpretar el ICC de M1:

• ICC cercano a 0 (< 0.1): las diferencias entre árboles son pequeñas; el modelo mixto aporta poco sobre OLS. Infrecuente en datos longitudinales biológicos.
• ICC entre 0.1 y 0.5: heterogeneidad moderada; el modelo mixto es claramente necesario.
• ICC > 0.5: la mayor parte de la variabilidad está entre árboles, no dentro. La media del grupo es una representación muy imprecisa de cada árbol individual.

---

M2: LMM con intercepto y pendiente aleatorios

Qué añade M2 respecto a M1

M2 reconoce que los árboles no solo difieren en su punto de partida sino también en su tasa de crecimiento ($b_{1i}$). En el Módulo 2, el gráfico de rejilla individual mostraba que las pendientes de las líneas de tendencia de cada árbol eran distintas; este modelo formaliza esa observación, en la línea de lo propuesto por Fitzmaurice et al. (2011, cap. 8).

$$\text{log\_growth}_{ij} = \underbrace{(\beta_0 + b_{0i})}_{\text{intercepto individual}} + \underbrace{(\beta_1 + b_{1i})}_{\text{pendiente individual}}\,\text{time\_c}_{ij} + \beta_2\,\text{Group}_i + \beta_3\,(\text{time\_c} \times \text{Group})_{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)$$

Lo que resuelve: perfiles que divergen con el tiempo. Sin pendiente aleatoria, el modelo M1 subestimaría la variabilidad en los tiempos finales (donde los árboles están más dispersos).

Lo que cuesta: añade dos parámetros ($\tau_1^2$ y $\tau_{01}$), y puede tener problemas de convergencia si hay pocos árboles o pocas mediciones. Correa y Salazar (2025, cap. 5) recomiendan examinar el gráfico de BLUPs de intercepto contra pendiente para juzgar si la correlación es sustantiva y si la Variante B está justificada.

Código R y SAS

Cómo leer la sección de efectos aleatorios de M2

Supongamos que la salida muestra:

consola R Random effects: Groups Name Variance Std.Dev. Corr ID (Intercept) 0.0842 0.290 time_c 0.0003 0.017 0.31 Residual 0.0241 0.155

La correlación $\rho_{01}$ tiene implicación práctica: si es positiva en este estudio, los árboles más grandes en el punto medio también tienen mayor tasa de crecimiento. Si el ozono reduce más el crecimiento en los árboles que ya crecen más rápido, habría una interacción biológica que el término fijo time_c * Group captura a nivel promedio.

---

Comparación de M0, M1 y M2

Cómo leer la comparación:

• Si AIC(M2) < AIC(M1) y el LRT es significativo ($p < 0.05$), la pendiente aleatoria mejora el ajuste y debe incluirse.
• Si BIC(M1) < BIC(M2) pero AIC(M2) < AIC(M1), el BIC favorece parsimonia. En ese caso, consulte la teoría: si hay buenas razones para esperar heterogeneidad de pendientes (el ozono puede afectar de forma diferencial a cada árbol), retenga M2. Esta es la postura defendida por Fitzmaurice et al. (2011, p. 210), quienes insisten en que los criterios estadísticos deben complementarse con el conocimiento del dominio.
• AIC(M1) casi siempre será mucho menor que AIC(M0): esto confirma cuantitativamente que ignorar la estructura de árbol (M0) produce un ajuste peor.

---

M3: Interpretación del modelo completo con interacción

M3 es M2 con el énfasis puesto en la pregunta científica central: ¿el ozono cambia la tasa de crecimiento de los árboles comparado con el control?

Interpretación sistemática de los efectos fijos

Supongamos estos coeficientes (ajuste a los valores reales al ejecutar el código):

---

Visualización del modelo ajustado

Qué leer en este gráfico:

• La separación vertical entre las dos curvas poblacionales cuantifica $\hat\beta_2$ (diferencia de nivel).
• La diferencia de pendiente entre las curvas cuantifica $\hat\beta_3$ (interacción).
• La dispersión de las líneas individuales alrededor de la curva de su grupo refleja $\tau_0^2$ y $\tau_1^2$.
• Si las líneas individuales de un grupo están más dispersas que las del otro, puede haber heterocedasticidad de grupo (un componente que exploraremos en el Módulo 5).

---

Tabla resumen: qué añade y qué cuesta cada componente

---

Verificar la distribución de los efectos aleatorios

Una forma de validar la especificación es verificar que los BLUPs estimados son aproximadamente normales, práctica sugerida por Correa y Salazar (2025, cap. 5) y por la tradición de diagnóstico en modelos mixtos.

Qué buscar:

• QQ-plot del intercepto: los puntos deben seguir la línea diagonal. Desviaciones en las colas indican asimetría en la distribución de niveles basales, lo que podría sugerir que algunos árboles son atípicos.
• QQ-plot de la pendiente: similar. Recuerde que con pocos árboles el QQ‑plot tiene poca potencia; las desviaciones moderadas no son preocupantes.
• Scatter b₀ᵢ vs b₁ᵢ: la pendiente de la recta de regresión visualiza $\hat\rho_{01}$. Si la nube es circular, la covarianza es cercana a cero y podría usarse la Variante C (Módulo 3).

---

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. [Material del curso, disponible en el repositorio].
• Fitzmaurice, G. M., Laird, N. M. y Ware, J. H. (2011). Applied Longitudinal Analysis (2.ª ed.). John Wiley & Sons. [Capítulos 1, 5 y 8: introducción a datos longitudinales, interacciones y construcción progresiva del LMM].