Modelos Mixtos

Correlaciones y varianzas heterocedásticas

Estructuras de correlación complejas, heterocedasticidad y criterios de selección para modelos con medidas repetidas

Datasets: Orthodont (nlme), Datos simulados con varianza creciente

Objetivos de aprendizaje

  • Comparar visual y numéricamente estructuras de correlación: simetría compuesta, AR(1), no estructurada y exponencial.
  • Detectar y modelar heterocedasticidad con varIdent y varPower en nlme.
  • Comprender el rol de la matriz de trabajo en modelos marginales con gls y su impacto en la eficiencia.
  • Tomar decisiones informadas usando AIC, BIC, gráficos de calor y gráficos de residuos.

¿Por qué importa la estructura de covarianza?

En estudios longitudinales, las observaciones de un mismo sujeto están correlacionadas. Ignorar esa correlación no siempre sesga los coeficientes estimados, pero sí afecta la eficiencia estadística: los errores estándar pueden ser incorrectos y se pierde potencia para detectar efectos relevantes.

En nlme disponemos de dos herramientas complementarias:

  • Modelos mixtos (lme): los efectos aleatorios inducen una estructura de correlación; las funciones cor* permiten refinarla.
  • Modelos marginales (gls): no incluyen efectos aleatorios, pero modelan directamente la correlación y la varianza.

Ambos enfoques son útiles y comparten las mismas funciones cor* y var*. En este módulo exploramos cuándo y cómo elegir entre las opciones disponibles.

Modelos condicionales vs. marginales en nlme

  • lme(random = ~ 1 | Subject, correlation = corAR1()) → modelo condicional con intercepto aleatorio más correlación AR(1).
  • gls(correlation = corAR1()) → modelo marginal basado en verosimilitud, que modela la correlación sin efectos aleatorios.
    Ambos son válidos; la elección depende de si necesitás los BLUPs de los efectos aleatorios o solo la estimación de los efectos fijos poblacionales.

1. Visualización exploratoria de la correlación y la varianza

Antes de elegir una estructura, conviene inspeccionar los datos. Usamos Orthodont, que contiene medidas repetidas de distancia dental en niños y niñas a los 8, 10, 12 y 14 años.

Spaghetti plot: trayectorias individuales

¿Qué observamos? Las trayectorias tienden a ser crecientes, pero no todas son paralelas. Algunos sujetos parten más arriba y otros más abajo. Esto sugiere que hay variabilidad entre sujetos que podría modelarse con efectos aleatorios o con una estructura de correlación adecuada.

Boxplot: dispersión por edad y sexo

Si el tamaño de las cajas cambia marcadamente entre edades o entre sexos, hay indicios de heterocedasticidad que convendrá modelar con funciones var*.

2. Comparación de estructuras de correlación

Ajustamos varios modelos con lme, variando únicamente la estructura de correlación, y comparamos criterios de información. Usamos method = "ML" para que la comparación por AIC/BIC sea válida.

Ajuste de cuatro estructuras de correlación
Comparación de AIC y BIC

¿Cómo elegir?

  • El menor AIC/BIC sugiere el mejor equilibrio entre ajuste y complejidad.
  • La no estructurada es la más flexible pero con muchos tiempos se vuelve inestable (muchos parámetros que estimar).
  • AR(1) y exponencial asumen que la correlación decae con la distancia temporal: muy razonable en estudios longitudinales.
  • La simetría compuesta asume correlación constante entre cualquier par de mediciones; es simple y útil cuando no hay razón teórica para esperar decaimiento.

3. Visualización de las matrices de correlación estimadas

Extraemos las matrices de cada modelo y las visualizamos en un heatmap. Usamos as.data.frame.table() en lugar de reshape2::melt para mantener la compatibilidad con WebR.

Heatmap de las matrices de correlación

Interpretación de los heatmaps:

  • En simetría compuesta, todas las correlaciones fuera de la diagonal son iguales (un solo valor para todos los pares).
  • En AR(1), la correlación disminuye progresivamente a medida que nos alejamos de la diagonal (más distancia temporal → menos correlación).
  • En no estructurada, cada par de edades tiene un valor único de correlación.
  • En exponencial, el decaimiento es más suave y depende de una función continua.

4. Modelos marginales con gls

Si no necesitamos los BLUPs de los efectos aleatorios, podemos usar gls para modelar directamente la estructura de covarianza. Esto es conceptualmente similar a lo que hacen los modelos marginales (como GEE), pero basado en máxima verosimilitud.

Modelo marginal con gls: comparación de estructuras

¿gls o lme?

  • gls es un modelo marginal: estima los efectos fijos poblacionales modelando directamente la correlación, sin estimar efectos aleatorios.
  • lme es un modelo condicional: incluye efectos aleatorios que inducen correlación.
  • Si solo te interesa la inferencia poblacional y no necesitás predicciones individuales, gls es una alternativa parsimoniosa dentro del ecosistema nlme.
  • GEE (con geepack) es otra alternativa marginal, pero basada en ecuaciones de estimación en lugar de verosimilitud. En WebR usamos gls como aproximación; en tu R local, geepack es el estándar.

5. Diagnóstico y modelado de heterocedasticidad

¿Qué ocurre si la varianza no es constante? Por ejemplo, si la dispersión de la distancia dental cambia con la edad o difiere entre sexos.

Gráfico de residuos del modelo homocedástico
Residuos agrupados por edad

Si los residuos muestran un embudo (dispersión creciente con los valores ajustados o con la edad), conviene modelar la varianza explícitamente.

Modelos con varianza heterocedástica
Parámetros de varianza estimados

No sobreparametrizar sin necesidad

Modelar la varianza añade parámetros adicionales. Si la heterocedasticidad es leve, un modelo más simple puede ser preferible según el principio de parsimonia. Confiá en AIC/BIC, pero también en la inspección visual: ¿el patrón de residuos cambia sustancialmente al modelar la varianza?

6. Validación del modelo final

Una vez elegida la estructura de correlación y de varianza, verificamos que los residuos ya no muestren patrones no deseados.

Residuos del modelo con heterocedasticidad modelada
Q-Q plot de residuos

7. Flujo de decisión recomendado

  1. 1

    Explorar los datos

    Spaghetti plots, boxplots por tiempo y grupo. Detectar visualmente correlación entre mediciones y posibles cambios en la dispersión.

  2. 2

    Elegir estructura de correlación

    Ajustar modelos con CS, AR(1), UN y Exponencial. Comparar AIC/BIC. Visualizar las matrices de correlación estimadas con heatmaps.

  3. 3

    Diagnosticar heterocedasticidad

    Inspeccionar residuos de Pearson vs. ajustados y vs. covariables (edad, grupo). Si hay patrón de embudo, proceder al paso 4.

  4. 4

    Modelar la varianza

    Probar varIdent (grupos), varPower o varExp (covariables continuas). Comparar AIC/BIC con el modelo homocedástico.

  5. 5

    Validar el modelo final

    Verificar que los residuos ya no muestren patrones. Q-Q plot para evaluar normalidad. Reportar la estructura elegida y justificarla.

Flashcard

Pregunta: ¿Qué estructura de varianza usarías si notás que la varianza de la respuesta crece sistemáticamente con el tiempo?

Respuesta: varPower o varExp dentro de lme() o gls(). Ambas permiten modelar la varianza como una función continua de una covariable (como la edad o el tiempo). varPower asume que la desviación estándar es proporcional a una potencia de la covariable; varExp asume que crece o decrece exponencialmente.

Pregunta: En un modelo con corAR1, ¿qué significa un valor de φ (phi) cercano a 0.9?

Respuesta: φ es el parámetro de autocorrelación de AR(1). Un valor de 0.9 indica que las mediciones adyacentes en el tiempo están muy fuertemente correlacionadas (correlación de 0.9), y esa correlación decae a medida que aumenta la distancia temporal: corr(t, t+2) = 0.9² = 0.81, corr(t, t+3) = 0.9³ ≈ 0.73, etc.

Quiz

Código SAS de referencia

Resumen: tabla comparativa de estructuras

EstructuraFunción en R (nlme)Cuándo usarla
Simetría compuestacorCompSymm()Correlación constante entre cualquier par de mediciones. Simple y robusta.
AR(1)corAR1()La correlación decae con la distancia temporal. Ideal para datos longitudinales con tiempo discreto.
No estructuradacorSymm()Cada par de tiempos tiene su propia correlación. Flexible pero costosa en parámetros.
ExponencialcorExp()Decaimiento continuo de la correlación con la distancia temporal. Útil con tiempos irregulares.
Varianza constante(por defecto)Asume homocedasticidad. Punto de partida antes de diagnosticar.
varIdentvarIdent(form = ~ 1 | grupo)Varianza diferente para cada nivel de un factor (ej. sexo, tratamiento).
varPowervarPower(form = ~ covariable)La desviación estándar es proporcional a una potencia de la covariable.

Con estas herramientas, ahora podés elegir y justificar la estructura de covarianza más adecuada para tus datos longitudinales, respaldando cada decisión con criterios numéricos y visualizaciones claras.

Referencias académicas

ReferenciaContribución
Pinheiro, J. & Bates, D. (2000). Mixed-Effects Models in S and S-PLUS. Springer.Texto de referencia para nlme. Cubre exhaustivamente las estructuras de correlación y varianza.
Diggle, P. J., Heagerty, P., Liang, K.-Y. & Zeger, S. L. (2002). Analysis of Longitudinal Data (2nd ed.). Oxford University Press.Texto canónico para datos longitudinales. Discute la elección de estructuras de covarianza.
Fitzmaurice, G. M., Laird, N. M. & Ware, J. H. (2011). Applied Longitudinal Analysis (2nd ed.). Wiley.Aplicaciones en salud con énfasis en diagnóstico y selección de modelos.
Verbeke, G. & Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data. Springer.Referencia completa sobre LMM, incluyendo diagnóstico y estructuras de covarianza.
Zuur, A. F. et al. (2009). Mixed Effects Models and Extensions in Ecology with R. Springer.Excelente guía práctica para selección de estructuras con ejemplos en R.