Última actualización de este sitio: 28 junio, 2023

————————

Objetivo del reporte

Explicar como se realiza el ajuste de un modelo no-lineal de efectos mixtos.

Contexto

Un modelo de efectos mixtos permite obtener efectos aleatorios que alteran sus parametros (i.e., fijos) para la estructura jerárquica definida Salas et al. (2008). Esto en términos simples, implica que los parametros del modelo tiene una desviación debido a dicha estructura.

Estos efectos aleatorios se predicen al momento de ajustar el modelo de efectos mixtos (generalmente mediante el método de máxima verosimilitud). Es decir, se pueden conocer para los datos usados en el ajuste del modelo, sin embargo, si uno quisiera ocupar el modelo en el contexto de predicción (i.e., para observaciones no ocupadas en el ajuste del modelo), dicho efectos aleatorios no son conocidos. Mayores detalles teóricos y aplicados pueden ser revisados en Gregoire & Schabenberger (1996) y Pinheiro & Bates (2000).

Datos

Se ocupará la dataframe \({\tt orange2}\) del paquete \({\tt datana}\) del libro Análisis de datos con el programa estadístico R: una introducción aplicada de Salas-Eljatib (2021). Esta dataframe contiene mediciones de diámetro de árboles de diferentes especies.

> library(datana)
> library(nlme)
> data(dgorange2)
> df <- dgorange2
> head(df)
> df<-groupedData(dap~tiempo|arbol,data=df)
> plot(dap~tiempo, data=df)

Dado que los datos tiene una estructura temporal, es necesario realizar un gráfico de serie de tiempo, para lo cual Ud. puede emplear una función del autor para dichos fines. Para esto, debe bajar a su computador el siguiente archivo que contiene la función a emplear.

> source("timeSerPlot.R")
> timeSerPlot(df, y="dap", x="tiempo", obs.unit = "arbol")

Otra alternativa es realizando un gráfico de serie de tiempo como es explicado en el libro de Salas-Eljatib (2021), y detallado en el siguiente link, al emplear la librería \({\tt lattice}\), como sigue

> library(lattice)

Attaching package: 'lattice'
The following object is masked _by_ '.GlobalEnv':

    parallel
> xyplot(dap ~ tiempo, group=arbol,
+         data=df, type="b")

El modelo base

Se empleará como modelo de referencia, o base, un modelo no lineal asintotico, expresado como

\[{y_{ij}}= {\alpha} \left(1-e^{-e^{(-\beta t_{ij})}}\right)+ \epsilon_{ij} \tag{1}\] donde \(y_{ij}\) y \(t_{ij}\) son el diámetro a la altura del pecho del \(i\)-ésimo árbol en la \(j\)-ésimo ocasión y la edad del \(i\)-ésimo árbol en la \(j\)-ésimo ocasión; \(\epsilon_{ij}\) es el error aleatorio que sigue una distribucion Gaussiana con valor esperado cero y varianza \(\sigma^2_\epsilon\); mientras que \(\alpha, \beta\) son parámetros a ser estimados.

Definiendo la función no lineal

Representaremos la funcion 1 en R mediante el objeto \({\tt asympOri}\)

> asympOri <- function (b0, b1, x)  {
+   b0*(1-exp(-exp(-b1)*x))
+ }

Ajuste mediante mínimos cuadrados ordinarios no-lineales

16172.    (2.62e+00): par = (70 6)
2487.8    (5.61e-01): par = (60.214 6.8102)
2292.6    (4.60e-01): par = (72.92 7.1722)
2120.2    (3.46e-01): par = (86.604 7.4173)
1975.3    (2.07e-01): par = (109.29 7.7115)
1894.3    (2.71e-03): par = (114.76 7.7262)
1894.3    (5.47e-06): par = (114.75 7.7254)

Formula: dap ~ asympOri(b0, b1, x = tiempo)

Parameters:
   Estimate Std. Error t value Pr(>|t|)
b0  114.752     46.874    2.45     0.02
b1    7.725      0.543   14.22  1.2e-15

Residual standard error: 7.58 on 33 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 5.47e-06

Ajuste modelo de efectos mixtos no-lineales

Emplearemos la función \(\tt{nlme()}\) del paquete \(\tt{nlme}\) Pinheiro et al. (2021)

> m1 <- nlme(dap ~ asympOri(b0, b1, x=tiempo), data=df,
+  fixed = b0+ b1 ~ 1, random = b0 ~ 1,
+  start = c(b0=70,b1=6) )
> 
> summary(m1)
Nonlinear mixed-effects model fit by maximum likelihood
  Model: dap ~ asympOri(b0, b1, x = tiempo) 
  Data: df 
     AIC    BIC  logLik
  206.94 213.16 -99.472

Random effects:
 Formula: b0 ~ 1 | arbol
            b0 Residual
StdDev: 17.109   3.2564

Fixed effects:  b0 + b1 ~ 1 
   Value Std.Error DF t-value p-value
b0 104.6   17.2556 29   6.062       0
b1   7.6    0.2025 29  37.525       0
 Correlation: 
   b0   
b1 0.886

Standardized Within-Group Residuals:
      Min        Q1       Med        Q3       Max 
-2.277806 -0.564852 -0.010944  0.848158  1.485987 

Number of Observations: 35
Number of Groups: 5

Bibliography

Gregoire T.G. & Schabenberger O. (1996). Nonlinear mixed-effects modeling of cumulative bole volume with spatially correlated within-tree data. J. Agric. Biol. Environ. Stat. 1: 107–119.
Pinheiro J.C., Bates D., DebRoy S., Sarkar D. & R Core Team (2021). nlme: Linear and nonlinear mixed effects models. https://CRAN.R-project.org/package=nlme.
Pinheiro J.C. & Bates D.M. (2000). Mixed-effects models in s and splus. Springer-Verlag, New York, USA.
Salas C., Stage A.R. & Robinson A.P. (2008). Modeling effects of overstory density and competing vegetation on tree height growth. Forest Science 54 (1): 107–122. https://doi.org/10.1093/forestscience/54.1.107.
Salas-Eljatib C. (2021). Análisis de datos con el programa estadístico R: una introducción aplicada. Ediciones Universidad Mayor, Santiago, Chile. https://tienda.zigzag.cl/9789566086109-analisis-de-datos-con-el-programa-estadistico-r.html.