Modelos mixtos en R: ¿lme, lmer o ambos? ¿cuál es relevante para mis datos y por qué?

Question

Dispongo de un dato obtenido a través de un inventario forestal realizado anualmente (1994-2015) en un país de África Occidental. Se seleccionaron 10 parcelas de igual tamaño (1 ha cada una) de bosque natural no gestionado y luego se identificaron y contaron las especies de árboles y arbustos. Se calcularon índices de biodiversidad como la abundancia, Shannon y Simpson. He elegido sólo 9 años en los que se recogieron datos en las 10 parcelas y he descartado los años incompletos y he considerado el "Año" como factor.

Los datos están estructurados como:

str(BIData)
'data.frame':   90 obs. of  9 variables:
 $ Year          : Factor w/ 9 levels "1994","1995",..: 1 1 1 1 1 1 1 1 1 1 ...
     $ Plot          : Factor w/ 10 levels "Bas Kolel","Bougou",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Richness      : int  8 21 13 14 8 10 6 10 8 20 ...
     $ Abundance     : int  286 1471 1121 466 242 97 250 790 208 2015 ...
 $ Shannon       : num  1.33 1.79 1.55 1.68 1.44 1.71 1.35 1.27 1.27 1.86 ...
     $ Simpson       : num  0.656 0.71 0.682 0.694 0.665 0.714 0.66 0.647 0.649 0.718 ...
 $ InverseSimpson: num  2.91 3.45 3.14 3.28 2.99 3.52 2.95 2.83 2.86 3.54 ...
     $ Topography    : Factor w/ 3 levels "Plateau","Slope",..: 3 1 1 3 3 2 2 2 3 1 ...
 $ Land_use      : Factor w/ 2 levels "Cultivated","Pasture": 1 2 2 2 1 1 2 2 1 2 ...

Además, las parcelas están situadas en diferentes topografías (ladera, valle, meseta) y usos del suelo (cultivo, pastos).

Tengo los siguientes dos modelos en lmer y lme:

model=lmer(Abundance~Year+Topography+Land_use+(1|Plot), method="ML", data=BIData)
model=lme(Abundance~Year+Topography+Land_use, random=~1|Plot, method="ML", data=BIData)

Obtuve resultados totalmente diferentes: ¿Mis preguntas?

No soy un experto pero he encontrado que lme proporciona una especie de resultados "bonitos" con valores p. Puedo ver muchos factores significativos como los años, la topografía y el uso del suelo, mientras que en lmer sólo los valores t sin los valores p. No sé cuál es el correcto para mis datos. En ambos casos, muestra buenos y aceptables gráficos de residuos.

Por favor, ayúdenme a entender cuál es el correcto para mis datos.

---

Gracias @fcoppens. No, no he probado ese parámetro. Aquí está la salida de ambos lme y lmer .

lmer

model=lmer(Abundance~Year+Topography+Land_use+(1|Plot), method="ML", data=BIData)
summary(model)
Linear mixed model fit by REML ['lmerMod']
Formula: Abundance ~ Year + Topography + Land_use + (1 | Plot)
   Data: BIData

REML criterion at convergence: 1106.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.5754 -0.5024 -0.0186  0.4015  3.4341 

Random effects:
 Groups   Name        Variance Std.Dev.
 Plot     (Intercept) 51753    227.5   
 Residual             48592    220.4   
Number of obs: 90, groups:  Plot, 10

Fixed effects:
                 Estimate Std. Error t value
(Intercept)       1073.15     252.41   4.252
Year1995             0.40      98.58   0.004
Year1996           -32.70      98.58  -0.332
Year1998          -198.10      98.58  -2.010
Year1999          -341.90      98.58  -3.468
Year2002          -295.80      98.58  -3.001
Year2004          -324.90      98.58  -3.296
Year2010          -291.60      98.58  -2.958
Year2015          -371.00      98.58  -3.763
TopographySlope   -756.87     206.36  -3.668
TopographyValley  -645.82     236.71  -2.728
Land_usePasture    178.07     200.85   0.887

lme

model=lme(Abundance~Year+Topography+Land_use, random=~1|Plot, method="ML", data=BIData)
summary(model)
Linear mixed-effects model fit by maximum likelihood
 Data: BIData 
       AIC      BIC    logLik
  1264.675 1299.673 -618.3377

Random effects:
 Formula: ~1 | Plot
        (Intercept) Residual
StdDev:    171.5578 209.1232

Fixed effects: Abundance ~ Year + Topography + Land_use 
                     Value Std.Error DF   t-value p-value
(Intercept)      1073.1495  213.5506 72  5.025271  0.0000
Year1995            0.4000  100.4595 72  0.003982  0.9968
Year1996          -32.7000  100.4595 72 -0.325504  0.7457
Year1998         -198.1000  100.4595 72 -1.971938  0.0525
Year1999         -341.9000  100.4595 72 -3.403360  0.0011
Year2002         -295.8000  100.4595 72 -2.944469  0.0044
Year2004         -324.9000  100.4595 72 -3.234138  0.0018
Year2010         -291.6000  100.4595 72 -2.902661  0.0049
Year2015         -371.0000  100.4595 72 -3.693029  0.0004
TopographySlope  -756.8671  171.7008  6 -4.408058  0.0045
TopographyValley -645.8214  196.9543  6 -3.279041  0.0168
Land_usePasture   178.0654  167.1213  6  1.065486  0.3276
Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.6851599 -0.5159528 -0.0222693  0.4401886  3.6493837 

Number of Observations: 90
Number of Groups: 10

Answer 1

No parece que se trate de "resultados totalmente diferentes"; los coeficientes de efectos fijos son los mismos entre los dos resultados, con estimaciones ligeramente diferentes para los errores estándar y los efectos aleatorios. Esto no es sorprendente, ya que el lmer utilizó la máxima verosimilitud restringida (REML) mientras que el lme ajuste utilizó la máxima verosimilitud (ML). Véase esta página con validación cruzada para una introducción a las diferencias entre estos enfoques.

La falta de p -valores en la salida de lmer es una elección consciente de los autores del paquete, como se comenta en la documentación del paquete y en esta página con validación cruzada . El argumento es que p -Los valores que se suelen calcular para los coeficientes en estos análisis son engañosos. Cuando los lme4 está cargado el comando help("pvalues") proporciona una guía para las formas de proceder, como se indica en la página 35 de la actual viñeta .

Consulta de seguimiento:

Gracias por las respuestas. Puedo ver que en los efectos fijos los coeficientes son los mismos en ambas salidas. Los valores de t son casi los mismos. Aunque no me interesan los efectos parcelarios (aleatorios) los resultados no son los mismos (residuos). Mis preguntas: (1) En el paquete lme4 (lmer) el autor omitió intencionadamente los valores p con algunas advertencias, así que ¿qué pasa con el paquete nlme? (2) ¿Son fiables los valores p resultantes de lme? Estos valores p coinciden con las tendencias de mis datos. (3) ¿Hay alguna razón para que mis datos sean más relevantes para usar lmer que lme? ¿Por favor?

Respuesta:

Hay dos cuestiones en las que hay que pensar. La primera es la elección entre los enfoques REML y ML para el ajuste del modelo. (El argumento "método" que dio a lmer no se reconoce; para conseguir un ajuste con ML en lmer tiene que especificar "REML=FALSE", como ha señalado @f_coppens. Puede obtener un ajuste con REML en lme con el argumento 'method="REML"' o simplemente omitir el argumento en lme ya que REML es la opción por defecto para lme ). El REML lmer frente al ML lme es casi seguro que las diferencias en los efectos aleatorios estimados, las diferencias en los errores estimados de los coeficientes y las diferencias resultantes en t -valores. Para elegir entre REML y ML, véase esta página con validación cruzada . Para un estudio relativamente pequeño como éste (al menos en relación con el número de parámetros que se intenta estimar), probablemente se quiera utilizar el método REML.

La segunda cuestión es cómo calcular p -valores. Un problema subyacente importante es cómo estimar los grados de libertad que quedan en el modelo después de tener en cuenta el número de parámetros estimados; es necesario conocer los grados de libertad para traducir t -valores a p -valores. Supongo que no hay nada malo en informar de la p -valores proporcionados por lme siempre que se especifique que se ha utilizado lme para conseguirlos. Esos p -son probablemente demasiado bajos, pero un lector podrá al menos saber cómo pensar en ellos y sus limitaciones. Sin embargo, es mejor que utilice los enfoques recomendados por help("pvalues") bajo la lme4 o proporcionada por paquetes accesorios como lmerTest . En cualquier caso, es importante reflexionar sobre las cuestiones subyacentes en lugar de limitarse a aceptar el resultado estándar del software estadístico.

Consulta de seguimiento:

Muchas gracias. Todas tus sugerencias han funcionado. Ejecuté los modelos tanto en lme como en lmer en la forma predeterminada sin especificar un método (REML es el predeterminado) y encontré los resultados casi idénticos excepto por una ligera diferencia en los valores p asociados con la diferencia en los DF. Si le interesa, publicaré los resultados?

model1=lme(Abundancia~Año+Topografía+Uso del suelo, random=~1|Plot, data=BIData) model2=lmer(Abundancia~Año+Topografía+uso del suelo+(1|Plan), data=BIData)