¿Cuál es la diferencia entre usar aov () y lme () al analizar un conjunto de datos longitudinal?

¿Alguien puede decirme la diferencia entre usar aov()y lme()analizar datos longitudinales y cómo interpretar los resultados de estos dos métodos?

A continuación, se analiza el mismo conjunto de datos usando aov()y lme()y tengo 2 resultados diferentes. Con aov(), obtuve un resultado significativo en la interacción tiempo por tratamiento, pero ajustando un modelo lineal mixto, la interacción tiempo por tratamiento es insignificante.

> UOP.kg.aov <- aov(UOP.kg~time*treat+Error(id), raw3.42)
> summary(UOP.kg.aov)

Error: id
          Df  Sum Sq Mean Sq F value Pr(>F)
treat      1   0.142  0.1421  0.0377 0.8471
Residuals 39 147.129  3.7725               

Error: Within
            Df  Sum Sq Mean Sq  F value  Pr(>F)    
time         1 194.087 194.087 534.3542 < 2e-16 ***
time:treat   1   2.077   2.077   5.7197 0.01792 *  
Residuals  162  58.841   0.363                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> UOP.kg.lme <- lme(UOP.kg~time*treat, random=list(id=pdDiag(~time)), 
                    na.action=na.omit, raw3.42)
> summary(UOP.kg.lme)
Linear mixed-effects model fit by REML
 Data: raw3.42 
       AIC      BIC    logLik
  225.7806 248.9037 -105.8903

Random effects:
 Formula: ~time | id
 Structure: Diagonal
        (Intercept)      time  Residual
StdDev:   0.6817425 0.5121545 0.1780466

Fixed effects: UOP.kg ~ time + treat + time:treat 
                 Value Std.Error  DF   t-value p-value
(Intercept)  0.5901420 0.1480515 162  3.986059  0.0001
time         0.8623864 0.1104533 162  7.807701  0.0000
treat       -0.2144487 0.2174843  39 -0.986042  0.3302
time:treat   0.1979578 0.1622534 162  1.220053  0.2242
 Correlation: 
           (Intr) time   treat 
time       -0.023              
treat      -0.681  0.016       
time:treat  0.016 -0.681 -0.023

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.198315285 -0.384858426  0.002705899  0.404637305  2.049705655 

Number of Observations: 205
Number of Groups: 41 

Según su descripción, parece que tiene un modelo de medidas repetidas con un solo factor de tratamiento. Como no tengo acceso al conjunto de datos ( raw3.42), utilizaré los datos de Orthodont del nlmepaquete para ilustrar lo que está sucediendo aquí. La estructura de datos es la misma (mediciones repetidas para dos grupos diferentes, en este caso, hombres y mujeres).

Si ejecuta el siguiente código:

library(nlme)
data(Orthodont)

res <- lme(distance ~ age*Sex, random = ~ 1 | Subject, data = Orthodont)
anova(res)

obtendrá los siguientes resultados:

            numDF denDF  F-value p-value
(Intercept)     1    79 4123.156  <.0001
age             1    79  122.450  <.0001
Sex             1    25    9.292  0.0054
age:Sex         1    79    6.303  0.0141

Si tu corres:

res <- aov(distance ~ age*Sex + Error(Subject), data = Orthodont)
summary(res)

conseguirás:

Error: Subject
          Df Sum Sq Mean Sq F value   Pr(>F)   
Sex        1 140.46 140.465  9.2921 0.005375 **
Residuals 25 377.91  15.117                    

Error: Within
          Df  Sum Sq Mean Sq  F value  Pr(>F)    
age        1 235.356 235.356 122.4502 < 2e-16 ***
age:Sex    1  12.114  12.114   6.3027 0.01410 *  
Residuals 79 151.842   1.922                     

Tenga en cuenta que las pruebas F son exactamente idénticas.

Para lme(), que utilizó list(id=pdDiag(~time)), que no sólo añade una intersección aleatoria al modelo, sino también una pendiente aleatoria. Además, al usar pdDiag, está configurando la correlación entre la intersección aleatoria y la pendiente a cero. Este es un modelo diferente al que especificó aov()y, por lo tanto, obtiene resultados diferentes.

Solo agregaría que es posible que desee instalar el carpaquete y usar el Anova()que proporciona este paquete en lugar de anova()porque aov()y para lm()objetos, la vainilla anova()usa una suma secuencial de cuadrados, lo que da el resultado incorrecto para tamaños de muestra desiguales, mientras que para lme()el tipo -I o la suma de cuadrados de tipo III dependiendo del typeargumento, pero la suma de cuadrados de tipo III viola la marginalidad, es decir, trata las interacciones de manera diferente a los efectos principales.

La lista de R-help no tiene nada bueno que decir sobre las sumas de cuadrados de tipo I y tipo III, ¡y sin embargo, estas son las únicas opciones! Imagínate.

Editar: En realidad, parece que el tipo II no es válido si hay un término de interacción significativo, y parece que lo mejor que alguien puede hacer es usar el tipo III cuando hay interacciones. Me dieron una respuesta a una de mis propias preguntas que a su vez me apuntó a esta publicación .

Me parece que tiene varias medidas para cada identificación en cada momento. Debe agregar estos para el aov porque aumenta injustamente el poder en ese análisis. No estoy diciendo que hacer el agregado hará que los resultados sean los mismos, pero debería hacerlos más similares.

dat.agg <- aggregate(UOP.kg ~ time + treat + id, raw3.42, mean)

Luego ejecute su modelo aov como antes de reemplazar los datos con dat.agg.

Además, creo que anova (lme) es más lo que quieres hacer para comparar los resultados. La dirección y la magnitud de un efecto no es lo mismo que la relación entre la varianza del modelo y el error.

(Por cierto, si haces el análisis de lme en los datos agregados, que no deberías, y compruebas anova (lme) obtendrás casi los mismos resultados que aov)