Pourquoi lme et aov renvoient-ils des résultats différents pour les mesures répétées ANOVA dans R?

J'essaie de passer de l'utilisation du ezpackage à lmedes mesures répétées ANOVA (car j'espère que je pourrai utiliser des contrastes personnalisés avec lme).

En suivant les conseils de ce billet de blog, j'ai pu configurer le même modèle en utilisant à la fois aov(comme le fait ez, sur demande) et lme. Cependant, alors que dans l'exemple donné dans ce post, les valeurs F sont parfaitement d'accord entre aovet lme(je l'ai vérifié, et elles le font), ce n'est pas le cas pour mes données. Bien que les valeurs F soient similaires, elles ne sont pas identiques.

aovrenvoie une valeur f de 1,3399, lmerenvoie 1,36264. Je suis prêt à accepter le aovrésultat comme étant le "bon" car c'est aussi ce que SPSS renvoie (et c'est ce qui compte pour mon terrain / mon superviseur).

Des questions:

1. Ce serait formidable si quelqu'un pouvait expliquer pourquoi cette différence existe et comment je peux l'utiliser lmepour fournir des résultats crédibles. (Je serais également disposé à utiliser à la lmerplace de lmece type de contenu, s'il donne le résultat "correct". Cependant, je ne l'ai pas utilisé jusqu'à présent.)

2. Après avoir résolu ce problème, je voudrais exécuter une analyse de contraste. En particulier, je serais intéressé par le contraste de la mise en commun des deux premiers niveaux de facteur (c.-à-d. c("MP", "MT")) Et de le comparer avec le troisième niveau de facteur (c.-à "AC"-d.). De plus, tester le troisième par rapport au quatrième niveau de facteur (c.-à-d. Par "AC"rapport à "DA").

Les données:

tau.base <- structure(list(id = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 
9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 
22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 
14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L), .Label = c("A18K", 
"D21C", "F25E", "G25D", "H05M", "H07A", "H08H", "H25C", "H28E", 
"H30D", "J10G", "J22J", "K20U", "M09M", "P20E", "P26G", "P28G", 
"R03C", "U21S", "W08A", "W15V", "W18R"), class = "factor"), factor = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("MP", "MT", "AC", "DA"
), class = "factor"), value = c(0.9648092876, 0.2128662077, 1, 
0.0607615485, 0.9912814024, 3.22e-08, 0.8073856412, 0.1465590332, 
0.9981672618, 1, 1, 1, 0.9794401938, 0.6102546108, 0.428651501, 
1, 0.1710644881, 1, 0.7639763913, 1, 0.5298989196, 1, 1, 0.7162733447, 
0.7871177434, 1, 1, 1, 0.8560509327, 0.3096989662, 1, 8.51e-08, 
0.3278862311, 0.0953598576, 1, 1.38e-08, 1.07e-08, 0.545290432, 
0.1305621416, 2.61e-08, 1, 0.9834051136, 0.8044114935, 0.7938839461, 
0.9910112678, 2.58e-08, 0.5762677121, 0.4750002288, 1e-08, 0.8584252623, 
1, 1, 0.6020385797, 8.51e-08, 0.7964935271, 0.2238374288, 0.263377904, 
1, 1.07e-08, 0.3160751898, 5.8e-08, 0.3460325565, 0.6842217296, 
1.01e-08, 0.9438301877, 0.5578367224, 2.18e-08, 1, 0.9161424562, 
0.2924856039, 1e-08, 0.8672987992, 0.9266688748, 0.8356425464, 
0.9988463913, 0.2960361777, 0.0285680426, 0.0969063841, 0.6947998266, 
0.0138254805, 1, 0.3494775301, 1, 2.61e-08, 1.52e-08, 0.5393467752, 
1, 0.9069223275)), .Names = c("id", "factor", "value"), class = "data.frame", row.names = c(1L, 
6L, 10L, 13L, 16L, 17L, 18L, 22L, 23L, 24L, 27L, 29L, 31L, 33L, 
42L, 43L, 44L, 45L, 54L, 56L, 58L, 61L, 64L, 69L, 73L, 76L, 79L, 
80L, 81L, 85L, 86L, 87L, 90L, 92L, 94L, 96L, 105L, 106L, 107L, 
108L, 117L, 119L, 121L, 124L, 127L, 132L, 136L, 139L, 142L, 143L, 
144L, 148L, 149L, 150L, 153L, 155L, 157L, 159L, 168L, 169L, 170L, 
171L, 180L, 182L, 184L, 187L, 190L, 195L, 199L, 202L, 205L, 206L, 
207L, 211L, 212L, 213L, 216L, 218L, 220L, 222L, 231L, 232L, 233L, 
234L, 243L, 245L, 247L, 250L))

Et le code:

require(nlme)

summary(aov(value ~ factor+Error(id/factor), data = tau.base))

anova(lme(value ~ factor, data = tau.base, random = ~1|id))

Ils sont différents car le modèle lme force la composante de variance de idà être supérieure à zéro. En regardant le tableau anova brut pour tous les termes, nous voyons que l'erreur quadratique moyenne pour id est inférieure à celle pour les résidus.

> anova(lm1 <- lm(value~ factor+id, data=tau.base))

          Df  Sum Sq Mean Sq F value Pr(>F)
factor     3  0.6484 0.21614  1.3399 0.2694
id        21  3.1609 0.15052  0.9331 0.5526
Residuals 63 10.1628 0.16131   

Lorsque nous calculons les composantes de la variance, cela signifie que la variance due à id sera négative. Ma mémoire de la mémoire des carrés moyens attendue est fragile, mais le calcul est quelque chose comme

(0.15052-0.16131)/3 = -0.003597.

Cela peut sembler étrange, mais cela peut arriver. Cela signifie que les moyennes de chaque id sont plus proches les unes des autres que vous ne vous y attendez compte tenu de la quantité de variation résiduelle dans le modèle.

En revanche, l'utilisation de lme force cette variance à être supérieure à zéro.

> summary(lme1 <- lme(value ~ factor, data = tau.base, random = ~1|id))
...
Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev: 3.09076e-05 0.3982667

Cela rapporte les écarts-types, la mise au carré pour obtenir les rendements de variance 9.553e-10pour la variance id et 0.1586164pour la variance résiduelle.

Maintenant, vous devez savoir que l'utilisation aovde mesures répétées n'est appropriée que si vous pensez que la corrélation entre toutes les paires de mesures répétées est identique; c'est ce qu'on appelle la symétrie composée. (Techniquement, sphéricité est nécessaire , mais cela suffit pour l' instant.) L' une des raisons d'utiliser lmeplus aovest qu'il peut gérer différents types de structures de corrélation.

Dans cet ensemble de données particulier, l'estimation de cette corrélation est négative; cela permet d'expliquer comment l'erreur quadratique moyenne pour id était inférieure à l'erreur quadratique résiduelle. Une corrélation négative signifie que si la première mesure d'un individu était inférieure à la moyenne, en moyenne, sa seconde serait supérieure à la moyenne, ce qui rend les moyennes totales des individus moins variables que ce à quoi nous nous attendions s'il y avait une corrélation nulle ou une corrélation positive.

L'utilisation lmeavec un effet aléatoire équivaut à ajuster un modèle de symétrie composé où cette corrélation est forcée d'être non négative; nous pouvons ajuster un modèle où la corrélation peut être négative en utilisant gls:

> anova(gls1 <- gls(value ~ factor, correlation=corCompSymm(form=~1|id),
                    data=tau.base))
Denom. DF: 84 
            numDF   F-value p-value
(Intercept)     1 199.55223  <.0001
factor          3   1.33985   0.267

Cette table ANOVA s'accorde avec la table de l' aovajustement et de l' lmajustement.

Okay, alors quoi? Eh bien, si vous pensez que la variance idet la corrélation entre les observations devraient être non négatives, l' lmeajustement est en fait plus approprié que l'ajustement utilisant aovou lmque son estimation de la variance résiduelle est légèrement meilleure. Cependant, si vous croyez que la corrélation entre les observations pourrait être négative, aovou lmou glsest mieux.

Vous pourriez également être intéressé par une exploration plus approfondie de la structure de corrélation; pour regarder une structure de corrélation générale, vous feriez quelque chose comme

gls2 <- gls(value ~ factor, correlation=corSymm(form=~unclass(factor)|id),
data=tau.base)

Ici, je limite uniquement la sortie à la structure de corrélation. Les valeurs 1 à 4 représentent les quatre niveaux de factor; on voit que le facteur 1 et le facteur 4 ont une corrélation négative assez forte:

> summary(gls2)
...
Correlation Structure: General
 Formula: ~unclass(factor) | id 
 Parameter estimate(s):
 Correlation: 
  1      2      3     
2  0.049              
3 -0.127  0.208       
4 -0.400  0.146 -0.024

Une façon de choisir entre ces modèles consiste à utiliser un test de rapport de vraisemblance; cela montre que le modèle à effets aléatoires et le modèle général de structure de corrélation ne sont pas statistiquement significativement différents; lorsque cela se produit, le modèle le plus simple est généralement préféré.

> anova(lme1, gls2)
     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
lme1     1  6 108.0794 122.6643 -48.03972                        
gls2     2 11 111.9787 138.7177 -44.98936 1 vs 2 6.100725  0.2965

aov()s'adapte au modèle via l' lm()utilisation des moindres carrés, lmes'adapte via le maximum de vraisemblance. Cette différence dans la façon dont les paramètres du modèle linéaire sont estimés explique probablement la (très petite) différence dans vos valeurs f.

En pratique (par exemple pour les tests d'hypothèses), ces estimations sont les mêmes, donc je ne vois pas comment l'une pourrait être considérée comme «plus crédible» que l'autre. Ils proviennent de différents paradigmes d'ajustement de modèle.

Pour les contrastes, vous devez configurer une matrice de contraste pour vos facteurs. Venebles et Ripley montrent comment procéder sur les p 143, p.146 et p.293-294 de la 4e édition.