Solved – How to interpret the output of lmerTest::ranova in R

anovahypothesis testinglme4-nlmemixed model

I wanted to find the p-value of the random effect in a linear mixed model. I used ranova from lmerTest.
I am not sure how to interpret the result. I thought it would be the same if I compared the two models, one with and one without the random effect. But if I look at AIC and BIC, the two results seem to contradict each other.

Am I missing something?

Image:

  • Top: lmerTest::ranova result.
  • Bottom: anova between two models (simple and mixed).

Top: lmerTest::ranova result. Bottom: anova between two models (simple and mixed).

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:

> library(lmerTest)
Loading required package: lme4
Loading required package: Matrix

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:stats’:

    step

> data("sleepstudy", package="lme4")
> 
> fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy) # Fit model with REML
> ranova(fm1) # REML LR-test
ANOVA-like table for random-effects: Single term deletions

Model:
Reaction ~ Days + (1 | Subject)
              npar  logLik    AIC   LRT Df Pr(>Chisq)    
<none>           4 -893.23 1794.5                        
(1 | Subject)    3 -946.83 1899.7 107.2  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

We can reproduce this test as follows:

> fm2 <- lm(Reaction ~ Days, sleepstudy) 
> 2 * c(logLik(fm1, REML=TRUE) - logLik(fm2, REML=TRUE)) # REML LR statistic
[1] 107.1986

You cannot reproduce this LR-test with anova():

> anova(fm1, fm2, refit=FALSE) # Not sensible
Data: sleepstudy
Models:
fm2: Reaction ~ Days
fm1: Reaction ~ Days + (1 | Subject)
    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
fm2  3 1906.3 1915.9 -950.15   1900.3                             
fm1  4 1794.5 1807.2 -893.23   1786.5 113.83      1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning message:
In anova.merMod(fm1, fm2, refit = FALSE) :
  some models fit with REML = TRUE, some not

If allowed to refit, anova produces the ML LR-test:

> anova(fm1, fm2) # ML LR-test
refitting model(s) with ML (instead of REML)
Data: sleepstudy
Models:
fm2: Reaction ~ Days
fm1: Reaction ~ Days + (1 | Subject)
    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
fm2  3 1906.3 1915.9 -950.15   1900.3                             
fm1  4 1802.1 1814.8 -897.04   1794.1 106.21      1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1'

Which can also be obtained with ranova()

> fm3 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy, REML=FALSE) # Fit model with ML
> ranova(fm3) # ML LR-test
ANOVA-like table for random-effects: Single term deletions

Model:
Reaction ~ Days + (1 | Subject)
              npar  logLik    AIC    LRT Df Pr(>Chisq)    
<none>           4 -897.04 1802.1                         
(1 | Subject)    3 -950.15 1906.3 106.21  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

And also 'by hand':

> 2 * c(logLik(fm3, REML=FALSE) - logLik(fm2, REML=FALSE)) # LR statistic
[1] 106.2144

So if you fit the model with REML=TRUE (default) you get REML likelihood ratio tests - if you use REML=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 using anova() but it can be computed by hand as well as by ranova() 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.

Related Question