Quelle est la différence entre l'utilisation de aov () et lme () dans l'analyse d'un ensemble de données longitudinales?

Quelqu'un peut-il me dire la différence entre l'utilisation aov()et l' lme()analyse de données longitudinales et comment interpréter les résultats de ces deux méthodes?

Ci-dessous, j'analyse le même ensemble de données en utilisant aov()et lme()et j'ai obtenu 2 résultats différents. Avec aov(), j'ai obtenu un résultat significatif dans l'interaction temps par traitement, mais pour un modèle mixte linéaire, l'interaction temps par traitement est insignifiante.

> 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 

D'après votre description, il semble que vous ayez un modèle à mesures répétées avec un seul facteur de traitement. Comme je n'ai pas accès à l'ensemble de données ( raw3.42), j'utiliserai les données Orthodont du nlmepackage pour illustrer ce qui se passe ici. La structure des données est la même (mesures répétées pour deux groupes différents - dans ce cas, les hommes et les femmes).

Si vous exécutez le code suivant:

library(nlme)
data(Orthodont)

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

vous obtiendrez les résultats suivants:

            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 vous exécutez:

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

tu auras:

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                     

Notez que les tests F sont exactement identiques.

Pour lme(), vous avez utilisé list(id=pdDiag(~time)), ce qui ajoute non seulement une interception aléatoire au modèle, mais également une pente aléatoire. De plus, en utilisant pdDiag, vous définissez la corrélation entre l'interception aléatoire et la pente sur zéro. Il s'agit d'un modèle différent de celui que vous avez spécifié via aov()et vous obtenez donc des résultats différents.

J'ajouterais simplement que vous souhaiterez peut-être installer le carpackage et l'utiliser Anova()que ce package fournit à la place de et anova()pour les objets, la vanille utilise une somme séquentielle de carrés, ce qui donne le mauvais résultat pour des tailles d'échantillons inégales tandis qu'elle utilise soit le type -I ou la somme des carrés de type III selon l' argument, mais la somme des carrés de type III viole la marginalité - c'est-à-dire qu'elle ne traite pas les interactions différemment des effets principaux.aov()lm()anova()lme()type

La liste R-help n'a rien de bon à dire sur les sommes des carrés de type I et de type III, et pourtant ce sont les seules options! Allez comprendre.

Edit: En fait, il semble que le type II ne soit pas valide s'il existe un terme d'interaction significatif, et il semble que le mieux que tout le monde puisse faire est d'utiliser le type III lorsqu'il y a des interactions. J'ai été averti par une réponse à l'une de mes propres questions qui, à son tour, m'a dirigé vers ce poste .

Il me semble que vous avez à chaque fois plusieurs mesures pour chaque identifiant. Vous devez les regrouper pour l'aov, car cela augmente injustement la puissance dans cette analyse. Je ne dis pas que faire l'agrégat rendra les résultats les mêmes, mais cela devrait les rendre plus similaires.

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

Exécutez ensuite votre modèle aov comme avant de remplacer les données par dat.agg.

De plus, je crois que l'anova (lme) est plus ce que vous voulez faire pour comparer les résultats. La direction et l'ampleur d'un effet ne sont pas les mêmes que le rapport de la variance du modèle à l'erreur.

(BTW, si vous faites l'analyse lme sur les données agrégées, ce que vous ne devriez pas faire, et vérifiez anova (lme), vous obtiendrez presque les mêmes résultats que aov)