Solved – Are degrees of freedom in lmerTest::anova correct? They are very different from RM-ANOVA

anovadegrees of freedomlme4-nlmemixed modelrepeated measures

I am analyzing the results of a reaction time experiment in R.

I ran a repeated measures ANOVA (1 within-subject factor with 2 levels and 1 between-subject factor with 2 levels). I ran a similar linear mixed model and I wanted to summarize the lmer results in the form of ANOVA table using lmerTest::anova.

Don't get me wrong: I did not expect the identical results, however I am not sure about the degrees of freedom in lmerTest::anova results. ~~It seems to me it rather reflects an ANOVA with no aggregation on subject-level.~~

I am aware of the fact that calculating degrees of freedom in mixed-effect models is tricky, but lmerTest::anova is mentioned as one possible solution in the updated ?pvalues topic (lme4 package).

Is this calculation correct? Do the results of lmerTest::anova correctly reflect the specified model?

Update: I made the individual differences larger. The degrees of freedom in lmerTest::anova are more different from simple anova, but I am still not sure, why they are so large for the within-subject factor/interaction.

# 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)

Results of the code above [updated]:

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(model)

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

Best Answer

I think that lmerTest is getting it right and ezanova is getting it wrong in this case.

the results from lmerTest agree with my intuition/understanding
two different computations in lmerTest (Satterthwaite and Kenward-Roger) agree
they also agree with nlme::lme
when I run it, ezanova gives a warning, which I don't entirely understand, but which should not be disregarded ...

Re-running example:

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)]

Figure out experimental design

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

So it looks like individuals (subnum) are placed in either control or treatment groups, and each is tested for both directions -- i.e. direction can be tested within individuals (denominator df is large), but group and group:direction can only be tested among individuals

(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

Here I get 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. The denominator DF look a little funky (all equal to 18): I think they should be larger for direction and group:direction, which can be tested independently (but would be smaller if you added (direction|subnum) to the model)?

# 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

the Df column here refers to the numerator df, Denom (second-to-last) gives the estimated denominator df; they agree with the classical intuition. More important, we also get different answers for the F values ...

We can also double-check with Kenward-Roger (very slow because it involves refitting the model several times)

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

The results are identical.

For this example lme (from the nlme package) actually does a perfectly good job guessing the appropriate denominator df (the F and p-values are very slightly different):

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

If I fit an interaction between direction and subnum the df for direction and group:direction are much smaller (I would have thought they would be 18, but maybe I'm getting something wrong):

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

Related Solutions

Solved – How to assign degrees of freedom for two-way ANOVA with two within-subjects factors

I'm not sure I understand the question exactly, but if you are asking about the df for the two-way, factorial, within-subjects ANOVA, here they are:

A = a - 1, where a = number of levels of A
B = b - 1, where b = number of levels of B
A x B = (a - 1)(b - 1)
S = n - 1, where s = number of levels of S (i.e., number of subjects)
A x S = (a - 1)(n - 1)
B x S = (b - 1)(n - 1)
A x B x S = (a - 1)(b - 1)(n - 1)

E.g.:

A = cond (a = 3); B = rnd (b = 6); S (s = 44)
- df_A = 2
- df_B = 5
- df_{A x B} = 10
- df_S = 43
- df_{A x S} = 86
- df_{B x S} = 215
- df_{A x B x S} = 430

Solved – Lme4 and lmertest: Different degrees of freedom from same dataset

Before jumping into a (partial, non-rigorous) explanation, let me point out that all of your denominator dfs are so large that there should be no practical difference between them. For any df greater than about 50, there's little difference between the $t$ distribution and the corresponding Normal distribution, or between the $F$ distribution and the corresponding (scaled) $\chi^2$ distribution -- unless you are doing something like estimating the 99.999% confidence intervals. For example, qnorm(0.75) (the upper tail cutoff of a symmetric 95% interval for a standard Normal) is 1.959964, while qt(0.975,df=1000) is 1.962339; changing the df to 10000 will hardly affect the result. (The comparison for $F$ vs. $\chi^2$ is almost identical, but slightly more annoying to compute because of the difference in scaling.)

You can play with lme to get the classical estimate of the degrees of freedom (be aware that lme's df-calculation heuristic can fail badly for random-slope models ...)

dd <- expand.grid(subj=1:67,stim=1:50,dur=1:3)
dd$y <- rnorm(nrow(dd))
library(nlme)
anova(lme(y~1,random=~1|subj/dur/stim,data=dd))
##            numDF denDF  F-value p-value
## (Intercept)     1  9849 2.520658  0.1124

The df value here is 67*3*49 ... one important thing you haven't told us is what fixed effects you are actually testing, and at what level they vary. I would guess that they don't actually vary the lowest level, so that in the classical case you would essentially be testing over a denominator based on the subject $\times$ duration SS.

I think the fact that the lowest-level variance (cd:dur:subj) is zero in the high-df model is actually a good clue. No variance among stimuli within duration essentially means that the responses are more less independent, so you get closer to the full possible df for the model ...

It can't be a coincidence that the other df you are getting are near 67*27=1809, but I'm not sure where the 27 would come from ...

I would welcome more rigorous/better-informed answers, especially about what the df would be in the classical case ...

Best Answer

Related Solutions

Solved – How to assign degrees of freedom for two-way ANOVA with two within-subjects factors

Solved – Lme4 and lmertest: Different degrees of freedom from same dataset

Related Question