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
Best Answer
Yes, you are missing the following line from
help(ranova)
: "If the model is fitted with REML the tests are REML-likelihood ratio tests." ;-)Here is a complete example:
We can reproduce this test as follows:
You cannot reproduce this LR-test with
anova()
:If allowed to
refit
,anova
produces the ML LR-test:Which can also be obtained with
ranova()
And also 'by hand':
So if you fit the model with
REML=TRUE
(default) you get REML likelihood ratio tests - if you useREML=FALSE
you get ML likelihood ratio tests. Unfortunately it is not possible to extract the REML likelihood ratio test of the last random-effects term in the model usinganova()
but it can be computed by hand as well as byranova()
as shown above.Also note that the behaviour of
lmerTest::ranova
is deliberate: The most sensible statistical default is to use REML LR tests of random-effects terms if you are fitting the model with REML.