J'analyse les résultats d'une expérience de temps de réaction chez R.

J'ai effectué une ANOVA à mesures répétées (1 facteur intra-sujet avec 2 niveaux et 1 facteur inter-sujet avec 2 niveaux). J'ai exécuté un modèle mixte linéaire similaire et je voulais résumer les résultats de lmer sous la forme d'une table ANOVA en utilisant lmerTest::anova.

Ne vous méprenez pas: je ne m'attendais pas à des résultats identiques, mais je ne suis pas sûr du degré de liberté des lmerTest::anovarésultats. Il me semble qu'elle reflète plutôt une ANOVA sans agrégation au niveau du sujet.

Je suis conscient du fait que le calcul des degrés de liberté dans les modèles à effets mixtes est délicat, mais lmerTest::anovaest mentionné comme une solution possible dans le ?pvaluessujet mis à jour ( lme4package).

Ce calcul est-il correct? Les résultats de lmerTest::anovareflètent-ils correctement le modèle spécifié?

Mise à jour: j'ai agrandi les différences individuelles. Les degrés de liberté lmerTest::anovasont plus différents de la simple anova, mais je ne sais toujours pas pourquoi ils sont si grands pour le facteur / interaction intra-sujet.

# mini example with ANT dataset from ez package
library(ez); library(lme4); library(lmerTest)

# repeated measures ANOVA with ez package
data(ANT)
ANT.2 <- subset(ANT, !error)
# update: make individual differences larger
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]

anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum), 
  within = .(direction), between = .(group))
anova.ez

# similarly with lmer and lmerTest::anova
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
lmerTest::anova(model)

# simple ANOVA on all available data
m <- lm(rt ~ group * direction, data = ANT.2)
anova(m)

Résultats du code ci-dessus [ mis à jour ]:

anova.ez

$ ANOVA

           Effect DFn DFd         F          p p<.05          ges
2           group   1  18 2.6854464 0.11862957       0.1294475137
3       direction   1  18 0.9160571 0.35119193       0.0001690471
4 group:direction   1  18 4.9169156 0.03970473     * 0.0009066868

lmerTest :: anova (modèle)

Analysis of Variance Table of type 3  with  Satterthwaite 
approximation for degrees of freedom
                Df Sum Sq Mean Sq F value Denom Pr(>F)
group            1  13293   13293  2.6830    18 0.1188
direction        1   1946    1946  0.3935  5169 0.5305
group:direction  1  11563   11563  2.3321  5169 0.1268

anova (m)

Analysis of Variance Table

Response: rt
                  Df   Sum Sq Mean Sq  F value Pr(>F)    
group              1  1791568 1791568 242.3094 <2e-16 ***
direction          1      728     728   0.0985 0.7537    
group:direction    1    12024   12024   1.6262 0.2023    
Residuals       5187 38351225    7394                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Je pense que lmerTestc'est bien et ezanovamal dans ce cas.

• les résultats de d' lmerTestaccord avec mon intuition / compréhension
• deux calculs différents dans lmerTest(Satterthwaite et Kenward-Roger) concordent
• ils sont également d'accord avec nlme::lme
• quand je l'exécute, ezanovadonne un avertissement, que je ne comprends pas tout à fait, mais qui ne doit pas être ignoré ...

Exemple de réexécution:

library(ez); library(lmerTest); library(nlme)
data(ANT)
ANT.2 <- subset(ANT, !error)
set.seed(101)  ## for reproducibility
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]

Comprendre la conception expérimentale

with(ANT.2,table(subnum,group,direction))

Il semble donc que les individus ( subnum) soient placés dans des groupes de contrôle ou de traitement, et chacun est testé pour les deux directions - c'est-à-dire que la direction peut être testée au sein des individus (le dénominateur df est grand), mais le groupe et le groupe: la direction ne peut être testée que parmi personnes

(anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum), 
    within = .(direction), between = .(group)))
## $ANOVA
##            Effect DFn DFd         F          p p<.05          ges
## 2           group   1  18 2.4290721 0.13651174       0.1183150147
## 3       direction   1  18 0.9160571 0.35119193       0.0002852171
## 4 group:direction   1  18 4.9169156 0.03970473     * 0.0015289914

Ici, je reçois Warning: collapsing data to cell means. *IF* the requested effects are a subset of the full design, you must use the "within_full" argument, else results may be inaccurate. le dénominateur DF un peu génial (tous égaux à 18): je pense qu'ils devraient être plus grands pour la direction et le groupe: direction, qui peuvent être testés indépendamment (mais seraient plus petits si vous les ajoutiez (direction|subnum)au modèle)?

# similarly with lmer and lmerTest::anova
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
lmerTest::anova(model)
##                 Df  Sum Sq Mean Sq F value Denom Pr(>F)
## group            1 12065.7 12065.7  2.4310    18 0.1364
## direction        1  1952.2  1952.2  0.3948  5169 0.5298
## group:direction  1 11552.2 11552.2  2.3299  5169 0.1270

la Dfcolonne se réfère ici au numérateur df, Denom(avant-dernier) donne le dénominateur estimé df; ils sont d'accord avec l'intuition classique. Plus important encore, nous obtenons également des réponses différentes pour les valeurs F ...

On peut aussi revérifier avec Kenward-Roger ( très lent car cela implique de remonter plusieurs fois le modèle)

lmerTest::anova(model,ddf="Kenward-Roger")

Les résultats sont identiques.

Pour cet exemple lme(du nlmepackage) fait en fait un très bon travail en devinant le dénominateur approprié df (les valeurs F et p sont très légèrement différentes):

model3 <- lme(rt ~ group * direction, random=~1|subnum, data = ANT.2)
anova(model3)[-1,]
##                 numDF denDF   F-value p-value
## group               1    18 2.4334314  0.1362
## direction           1  5169 0.3937316  0.5304
## group:direction     1  5169 2.3298847  0.1270

Si je correspond à une interaction entre directionet subnumle df pour directionet group:directionsont beaucoup plus petits (j'aurais pensé qu'ils auraient 18 ans, mais peut-être que je me trompe):

model2 <- lmer(rt ~ group * direction + (direction | subnum), data = ANT.2)
lmerTest::anova(model2)
##                 Df  Sum Sq Mean Sq F value   Denom Pr(>F)
## group            1 20334.7 20334.7  2.4302  17.995 0.1364
## direction        1  1804.3  1804.3  0.3649 124.784 0.5469
## group:direction  1 10616.6 10616.6  2.1418 124.784 0.1459

Je suis généralement d'accord avec l'analyse de Ben, mais permettez-moi d'ajouter quelques remarques et un peu d'intuition.

Premièrement, les résultats globaux:

1. Les résultats de lmerTest utilisant la méthode Satterthwaite sont corrects
2. La méthode Kenward-Roger est également correcte et en accord avec Satterthwaite

Ben décrit la conception dans laquelle subnumest imbriqué dans grouptout direction et group:directionsont croisées avec subnum. Cela signifie que le terme d'erreur naturelle (c'est-à-dire la "strate d'erreur englobante") pour groupest subnumalors que la strate d'erreur englobante pour les autres termes (y compris subnum) est les résidus.

Cette structure peut être représentée dans un diagramme dit de structure factorielle:

names <- c(expression("[I]"[5169]^{5191}),
           expression("[subnum]"[18]^{20}), expression(grp:dir[1]^{4}),
           expression(dir[1]^{2}), expression(grp[1]^{2}), expression(0[1]^{1}))
x <- c(2, 4, 4, 6, 6, 8)
y <- c(5, 7, 5, 3, 7, 5)
plot(NA, NA, xlim=c(2, 8), ylim=c(2, 8), type="n", axes=F, xlab="", ylab="")
text(x, y, names) # Add text according to ’names’ vector
# Define coordinates for start (x0, y0) and end (x1, y1) of arrows:
x0 <- c(1.8, 1.8, 4.2, 4.2, 4.2, 6, 6) + .5
y0 <- c(5, 5, 7, 5, 5, 3, 7)
x1 <- c(2.7, 2.7, 5, 5, 5, 7.2, 7.2) + .5
y1 <- c(5, 7, 7, 3, 7, 5, 5)
arrows(x0, y0, x1, y1, length=0.1)

Ici, les termes aléatoires sont mis entre parenthèses, 0représentent la moyenne globale (ou interception), [I]représentent le terme d'erreur, les numéros de super-script sont le nombre de niveaux et les numéros de sous-script sont le nombre de degrés de liberté en supposant une conception équilibrée. Le diagramme indique que le terme d'erreur naturelle (englobant la strate d'erreur) pour groupest subnumet que le numérateur df pour subnum, qui est égal au dénominateur df pour group, est 18: 20 moins 1 df pour groupet 1 df pour la moyenne globale. Une introduction plus complète aux diagrammes de structure factorielle est disponible dans le chapitre 2 ici: https://02429.compute.dtu.dk/eBook .

Si les données étaient exactement équilibrées, nous serions en mesure de construire les tests F à partir d'une décomposition SSQ comme fourni par anova.lm. Étant donné que l'ensemble de données est très étroitement équilibré, nous pouvons obtenir des tests F approximatifs comme suit:

ANT.2 <- subset(ANT, !error)
set.seed(101)
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]
fm <- lm(rt ~ group * direction + subnum, data=ANT.2)
(an <- anova(fm))
Analysis of Variance Table

Response: rt
                  Df   Sum Sq Mean Sq  F value Pr(>F)    
group              1   994365  994365 200.5461 <2e-16 ***
direction          1     1568    1568   0.3163 0.5739    
subnum            18  7576606  420923  84.8927 <2e-16 ***
group:direction    1    11561   11561   2.3316 0.1268    
Residuals       5169 25629383    4958                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ici, toutes les valeurs F et p sont calculées en supposant que tous les termes ont les résidus comme strate d'erreur englobante, et cela est vrai pour tous sauf «groupe». Le test F «équilibré-correct» pour le groupe est plutôt:

F_group <- an["group", "Mean Sq"] / an["subnum", "Mean Sq"]
c(Fvalue=F_group, pvalue=pf(F_group, 1, 18, lower.tail = FALSE))
   Fvalue    pvalue 
2.3623466 0.1416875 

où nous utilisons le subnumMS au lieu du ResidualsMS dans le dénominateur de valeur F.

Notez que ces valeurs correspondent assez bien aux résultats de Satterthwaite:

model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
anova(model, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group           12065.3 12065.3     1    18  2.4334 0.1362
direction        1951.8  1951.8     1  5169  0.3936 0.5304
group:direction 11552.2 11552.2     1  5169  2.3299 0.1270

Les différences restantes sont dues au fait que les données ne sont pas exactement équilibrées.

L'OP se compare anova.lmà anova.lmerModLmerTest, ce qui est correct, mais pour comparer comme avec, nous devons utiliser les mêmes contrastes. Dans ce cas, il y a une différence entre anova.lmet anova.lmerModLmerTestpuisqu'ils produisent respectivement des tests de type I et III, et pour cet ensemble de données, il y a une (petite) différence entre les contrastes de type I et III:

show_tests(anova(model, type=1))$group
               (Intercept) groupTreatment directionright groupTreatment:directionright
groupTreatment           0              1    0.005202759                     0.5013477

show_tests(anova(model, type=3))$group # type=3 is default
               (Intercept) groupTreatment directionright groupTreatment:directionright
groupTreatment           0              1              0                           0.5

Si l'ensemble de données avait été complètement équilibré, les contrastes de type I auraient été les mêmes que les contrastes de type III (qui ne sont pas affectés par le nombre d'échantillons observé).

Une dernière remarque est que la `` lenteur '' de la méthode de Kenward-Roger n'est pas due au réajustement du modèle, mais parce qu'elle implique des calculs avec la matrice marginale de variance-covariance des observations / résidus (5191x5191 dans ce cas) qui n'est pas le cas de la méthode de Satterthwaite.

Concernant model2

Quant au modèle 2, la situation devient plus complexe et je pense qu'il est plus facile de commencer la discussion avec un autre modèle où j'ai inclus l'interaction «classique» entre subnumet direction:

model3 <- lmer(rt ~ group * direction + (1 | subnum) +
                 (1 | subnum:direction), data = ANT.2)
VarCorr(model3)
 Groups           Name        Std.Dev.  
 subnum:direction (Intercept) 1.7008e-06
 subnum           (Intercept) 4.0100e+01
 Residual                     7.0415e+01

Étant donné que la variance associée à l'interaction est essentiellement nulle (en présence de l' subnumeffet principal aléatoire), le terme d'interaction n'a aucun effet sur le calcul des degrés de liberté du dénominateur, des valeurs F et des valeurs p :

anova(model3, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group           12065.3 12065.3     1    18  2.4334 0.1362
direction        1951.8  1951.8     1  5169  0.3936 0.5304
group:direction 11552.2 11552.2     1  5169  2.3299 0.1270

Cependant, subnum:directionest la strate d'erreur englobante pour subnumdonc si nous supprimons subnumtous les SSQ associés retombe danssubnum:direction

model4 <- lmer(rt ~ group * direction +
                 (1 | subnum:direction), data = ANT.2)

Maintenant, le terme d'erreur naturelle pour group, directionet group:directionest subnum:directionet avec nlevels(with(ANT.2, subnum:direction))= 40 et quatre paramètres, les degrés de liberté du dénominateur pour ces termes devraient être d'environ 36:

anova(model4, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
group           24004.5 24004.5     1 35.994  4.8325 0.03444 *
direction          50.6    50.6     1 35.994  0.0102 0.92020  
group:direction   273.4   273.4     1 35.994  0.0551 0.81583  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ces tests F peuvent également être approximés par les tests F «corrects équilibrés» :

an4 <- anova(lm(rt ~ group*direction + subnum:direction, data=ANT.2))
an4[1:3, "F value"] <- an4[1:3, "Mean Sq"] / an4[4, "Mean Sq"]
an4[1:3, "Pr(>F)"] <- pf(an4[1:3, "F value"], 1, 36, lower.tail = FALSE)
an4
Analysis of Variance Table

Response: rt
                   Df   Sum Sq Mean Sq F value Pr(>F)    
group               1   994365  994365  4.6976 0.0369 *  
direction           1     1568    1568  0.0074 0.9319    
group:direction     1    10795   10795  0.0510 0.8226    
direction:subnum   36  7620271  211674 42.6137 <2e-16 ***
Residuals        5151 25586484    4967                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Passons maintenant au modèle 2:

model2 <- lmer(rt ~ group * direction + (direction | subnum), data = ANT.2)

Ce modèle décrit une structure de covariance à effet aléatoire assez compliquée avec une matrice de variance-covariance 2x2. La paramétrisation par défaut n'est pas facile à gérer et on préfère une re-paramétrisation du modèle:

model2 <- lmer(rt ~ group * direction + (0 + direction | subnum), data = ANT.2)

Si l' on compare model2à model4, ils ont aussi de nombreux effets aléatoires; 2 pour chacun subnum, soit 2 * 20 = 40 au total. Bien qu'il model4stipule un seul paramètre de variance pour les 40 effets aléatoires,model2 stipule que chaque subnumpaire d'effets aléatoires a une distribution normale bi-variée avec une matrice de variance-covariance 2x2 dont les paramètres sont donnés par

VarCorr(model2)
 Groups   Name           Std.Dev. Corr 
 subnum   directionleft  38.880        
          directionright 41.324   1.000
 Residual                70.405        

Cela indique un ajustement excessif, mais gardons cela pour un autre jour. Le point important est que model4est un cas particulier de model2 et qui modelest aussi un cas particulier de model2. Parler librement (et intuitivement) (direction | subnum)contient ou capture la variation associée à l'effet principal subnum ainsi que l'interactiondirection:subnum . En termes d'effets aléatoires, nous pouvons considérer ces deux effets ou structures comme capturant la variation entre les lignes et les lignes par colonnes respectivement:

head(ranef(model2)$subnum)
  directionleft directionright
1    -25.453576     -27.053697
2     16.446105      17.479977
3    -47.828568     -50.835277
4     -1.980433      -2.104932
5      5.647213       6.002221
6     41.493591      44.102056

Dans ce cas, ces estimations de l'effet aléatoire ainsi que les estimations des paramètres de variance indiquent toutes deux que nous n'avons vraiment qu'un effet principal aléatoire de subnum(variation entre les lignes) présent ici. Tout cela conduit à ce que les degrés de liberté du dénominateur Satterthwaite dans

anova(model2, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF   DenDF F value Pr(>F)
group           12059.8 12059.8     1  17.998  2.4329 0.1362
direction        1803.6  1803.6     1 125.135  0.3638 0.5475
group:direction 10616.6 10616.6     1 125.136  2.1418 0.1458

est un compromis entre ces structures d'effet principal et d'interaction: le groupe DenDF reste à 18 (imbriqué subnumpar conception) mais le directionet group:directionDenDF sont des compromis entre 36 ( model4) et 5169 ( model).

Je ne pense pas que quoi que ce soit ici indique que l'approximation de Satterthwaite (ou son implémentation dans lmerTest ) est défectueuse.

Le tableau équivalent avec la méthode de Kenward-Roger donne

anova(model2, type=1, ddf="Ken")
Type I Analysis of Variance Table with Kenward-Roger's method
                 Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
group           12059.8 12059.8     1 18.000  2.4329 0.1362
direction        1803.2  1803.2     1 17.987  0.3638 0.5539
group:direction 10614.7 10614.7     1 17.987  2.1414 0.1606

Il n'est pas surprenant que KR et Satterthwaite puissent différer, mais à toutes fins pratiques, la différence de valeurs p est minime. Mon analyse ci-dessus indique que le DenDFfor directionet group:directionne devrait pas être inférieur à ~ 36 et probablement plus grand que celui étant donné que nous n'avons essentiellement que l'effet principal aléatoire de directionpresent, donc si je pense que c'est une indication que la méthode KR devient DenDFtrop faible dans ce cas. Mais gardez à l'esprit que les données ne prennent pas vraiment en charge la (group | direction)structure, la comparaison est donc un peu artificielle - il serait plus intéressant si le modèle était réellement pris en charge.