Modelos Mixtos

Modelos No Lineales Mixtos (NLMM)

Cuando la función de la media no es lineal en los parámetros: ajuste con nlme, valores iniciales y convergencia.

Datasets: SPRUCE1 (crecimiento de abetos Sitka, Control vs Ozono)

Objetivos de aprendizaje

  • Distinguir un NLMM de un LMM: parámetros que entran de forma no lineal en la función de la media.
  • Interpretar los parámetros del modelo logístico (asíntota, punto de inflexión, tasa de crecimiento).
  • Ajustar un NLMM con nlme::nlme() sobre datos de crecimiento de abetos.
  • Entender por qué los valores iniciales son críticos y cómo elegirlos gráficamente.
  • Comparar la curva poblacional con las curvas individuales predichas.

Del modelo lineal al modelo no lineal: ¿cuándo dar el salto?

En todos los módulos anteriores el modelo tenía la forma general:

yi=Xiβ+Zibi+εi\mathbf{y}_i = \mathbf{X}_i \boldsymbol{\beta} + \mathbf{Z}_i \mathbf{b}_i + \boldsymbol{\varepsilon}_i

La respuesta es una función lineal de los parámetros β\boldsymbol{\beta} y bi\mathbf{b}_i. Esto es conveniente computacionalmente y suficiente para muchas aplicaciones, pero en biología, farmacología y ecología los procesos suelen seguir curvas de crecimiento, saturación o decaimiento que no pueden parametrizarse como combinaciones lineales de las covariables.

Un Modelo No Lineal Mixto (NLMM) generaliza el LMM permitiendo que la función de la media sea no lineal en los parámetros:

yij=f(tij,ϕi)+εijy_{ij} = f(t_{ij}, \boldsymbol{\phi}_i) + \varepsilon_{ij}

donde ff es una función arbitraria (exponencial, logística, Michaelis-Menten, etc.) y los parámetros individuales ϕi\boldsymbol{\phi}_i se descomponen en una parte poblacional (efectos fijos) y una parte individual (efectos aleatorios):

ϕi=Aiβ+Bibi,bi∼N(0,D)\boldsymbol{\phi}_i = \mathbf{A}_i \boldsymbol{\beta} + \mathbf{B}_i \mathbf{b}_i, \qquad \mathbf{b}_i \sim \mathcal{N}(\mathbf{0}, \mathbf{D})
LMM vs NLMM: diferencias clave
AspectoLMMNLMM
Función de la mediaLineal en β: Xβ + ZbNo lineal: f(t, φᵢ)
EstimaciónSolución analítica (GLS/REML)Algoritmo iterativo (GLS linealizado, Laplace)
Valores inicialesNo son críticosCríticos: convergencia sensible a arranque
InterpretaciónPendientes e interceptosParámetros con significado físico/biológico directo
Paquetes Rlmer(), lme()nlme(), nlmer() (experimental)

El modelo logístico de tres parámetros

El modelo que usamos en este módulo —y que es el mismo que Correa y Salazar (2016, Cap. 7) aplican sobre los datos de abetos— es la función logística de tres parámetros:

f(t,ϕ1,ϕ2,ϕ3)=ϕ11+exp⁡ ⁣(−t−ϕ2ϕ3)f(t, \phi_1, \phi_2, \phi_3) = \frac{\phi_1}{1 + \exp\!\left(-\dfrac{t - \phi_2}{\phi_3}\right)}

Los tres parámetros tienen una interpretación geométrica clara:

Interpretación de los parámetros del modelo logístico
ParámetroNombreInterpretación
φ₁ (asíntota)Capacidad máximaValor límite que alcanza la respuesta cuando t → ∞. En los abetos: tamaño máximo del tronco.
φ₂ (punto de inflexión)Tiempo de inflexiónValor de t en que la curva tiene la pendiente máxima, es decir, cuando el crecimiento es más rápido.
φ₃ (escala)Tasa de transiciónControla qué tan 'empinada' es la curva alrededor de φ₂. Valores grandes → curva más suave.

¿Por qué la función logística para crecimiento de árboles?

Los árboles no crecen indefinidamente. Hay una capacidad de carga biológica determinada por luz, nutrientes y espacio. La función logística modela este crecimiento con aceleración inicial, un pico de crecimiento y luego saturación —exactamente el patrón observado en los abetos Sitka.

El dataset SPRUCE1: continuidad con módulos anteriores

Usamos el mismo dataset de abetos (Picea sitchensis) que en los módulos 4 y 5: 24 árboles (12 en condiciones de control, 12 expuestos a ozono), medidos en 4 momentos (días 152, 174, 201 y 227 desde el inicio del experimento). La variable respuesta es el log del tamaño del tronco (log_size).

El modelo logístico mixto que se ajusta por separado para cada grupo (Control / Ozono) es:

log⁡(tronco)ij=β1+b1i1+exp⁡(−tij−β2β3)+εij\log(\text{tronco})_{ij} = \frac{\beta_1 + b_{1i}}{1 + \exp\left(-\dfrac{t_{ij} - \beta_2}{\beta_3}\right)} + \varepsilon_{ij}

donde el intercepto aleatorio actúa sobre la asíntota β1\beta_1: cada árbol tiene su propia capacidad máxima. Los parámetros de curvatura (β2\beta_2 y β3\beta_3) son fijos (poblacionales).

La elección de valores iniciales es el reto principal de los NLMM

A diferencia de los LMM, donde lmer() siempre converge, en los NLMM el algoritmo puede divergir o converger a un mínimo local si los valores iniciales están lejos de la solución. La mejor estrategia es arrancar graficando los datos y estimando visualmente los parámetros: la asíntota ≈ valor máximo observado, el punto de inflexión ≈ tiempo donde el crecimiento se acelera más.

Cómo elegir valores iniciales gráficamente

Antes de ajustar el modelo, ejecuta la celda siguiente para visualizar los datos y estimar los parámetros iniciales de forma visual.

Paso 1: Exploración visual para elegir valores iniciales

Ajuste del NLMM con nlme

Paso 2: Ajuste del NLMM por grupo con nlme()

Interpretación de la salida

Una vez ajustados los modelos, los resultados clave son:

Efectos fijos (Fixed effects): las estimaciones de β^1\hat\beta_1 (asíntota poblacional), β^2\hat\beta_2 (punto de inflexión) y β^3\hat\beta_3 (escala). Correa y Salazar (2016) reportan β^1≈6.63\hat\beta_1 \approx 6.63 para Control y β^1≈8.73\hat\beta_1 \approx 8.73 para Ozono, lo que indica que los árboles expuestos a ozono tienen una asíntota mayor —crecen más, lo cual es consistente con la hipótesis de que el ozono puede estimular el crecimiento a dosis moderadas.

Efectos aleatorios (Random effects: A): la varianza σ^b12\hat\sigma_{b_1}^2 captura cuánto difiere la asíntota individual de la poblacional. Un valor grande indica alta heterogeneidad entre árboles en su tamaño máximo potencial.

Residual (sigma): la varianza σ^2\hat\sigma^2 del error dentro de cada árbol, una vez descontada la curva ajustada.

Curvas ajustadas: populacional e individuales

Paso 3: Visualizar curvas ajustadas (poblacional e individuales)

Impacto de los valores iniciales en la convergencia

La siguiente celda demuestra experimentalmente qué ocurre cuando los valores iniciales están lejos de la solución. Puedes modificar los valores de A_ini, B_ini y C_ini para explorar la sensibilidad.

Paso 4: Sensibilidad a los valores iniciales

Estrategia práctica para valores iniciales

  1. φ₁ (asíntota): toma el valor máximo observado en los datos multiplicado por 1.1.
  2. φ₂ (punto de inflexión): observa gráficamente dónde la pendiente de los perfiles es máxima.
  3. φ₃ (escala): empieza con un valor entre 50 y 200 y ajusta si no converge.
  4. Si el modelo no converge, reduce la tolerancia (pnlsTol) o prueba nlmeControl(pnlsMaxIter = 20, msMaxIter = 100).

Comparación con el LMM: ¿cuándo vale la pena el NLMM?

Comparación NLMM vs LMM para los datos de abetos (grupo Control)

Conceptos clave del NLMM

Repaso rápido: flashcards

Quiz de autoevaluación

Autoevaluación

1. En el NLMM logístico ajustado a los abetos, el grupo Ozono tiene una asíntota (β₁) mayor que el grupo Control. ¿Cuál es la interpretación correcta?

2. ¿Por qué nlme() puede fallar mientras que lmer() nunca lo hace en los módulos anteriores?

3. Se ajusta el NLMM con random = A ~ 1 | ID. ¿Qué implica NO incluir B y C como aleatorios?

Referencias del módulo

  • Correa, J. C. y Salazar, J. C. (2016). Introducción a los modelos mixtos. Universidad Nacional de Colombia, Sede Medellín. [Capítulo 7: Modelo no lineal mixto; pp. 190–206].
  • Pinheiro, J. y Bates, D. (2000). Mixed-Effects Models in S and S-PLUS. Springer. [Capítulos 6–8: modelos no lineales mixtos, estimación y diagnóstico].
  • Davidian, M. y Giltinan, D. M. (1995). Nonlinear Models for Repeated Measurement Data. Chapman & Hall. [Referencia clásica sobre NLMM con aplicaciones en farmacocinética].
  • Pinheiro, J. C. y Bates, D. M. (1995). Approximations to the log-likelihood function in the nonlinear mixed-effects model. Journal of Computational and Graphical Statistics, 4(1), 12–35.