anova() works different for mixed versus lme4 objects

[singmann]

library("afex")
#> Loading required package: lme4
#> Loading required package: Matrix
#> Registered S3 methods overwritten by 'car':
#>   method                          from
#>   influence.merMod                lme4
#>   cooks.distance.influence.merMod lme4
#>   dfbeta.influence.merMod         lme4
#>   dfbetas.influence.merMod        lme4
#> ************
#> Welcome to afex. For support visit: http://afex.singmann.science/
#> - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
#> - Methods for calculating p-values with mixed(): 'KR', 'S', 'LRT', and 'PB'
#> - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
#> - NEWS: library('emmeans') now needs to be called explicitly!
#> - Get and set global package options with: afex_options()
#> - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
#> - For example analyses see: browseVignettes("afex")
#> ************
#> 
#> Attaching package: 'afex'
#> The following object is masked from 'package:lme4':
#> 
#>     lmer

data("Machines", package = "MEMSS") 

# simple model with random-slopes for repeated-measures factor
m1 <- mixed(score ~ Machine + (Machine|Worker), data=Machines)
#> Contrasts set to contr.sum for the following variables: Machine, Worker
#> Fitting one lmer() model. [DONE]
#> Calculating p-values. [DONE]

m2 <- mixed(score ~ Machine + (Machine||Worker), data=Machines, expand_re = TRUE)
#> Contrasts set to contr.sum for the following variables: Machine, Worker
#> Fitting one lmer() model. [DONE]
#> Calculating p-values. [DONE]

anova(m1, m2)
#> Data: data
#> Models:
#> m2: score ~ Machine + (1 + re1.Machine1 + re1.Machine2 || Worker)
#> m1: score ~ Machine + (Machine | Worker)
#>    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
#> m2    7 233.75 247.67 -109.88   219.75                       
#> m1   10 230.51 250.40 -105.25   210.51 9.2424  3    0.02624 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m1$full_model, m2$full_model)
#> refitting model(s) with ML (instead of REML)
#> Data: data
#> Models:
#> m2$full_model: score ~ Machine + (1 + re1.Machine1 + re1.Machine2 || Worker)
#> m1$full_model: score ~ Machine + (Machine | Worker)
#>               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
#> m2$full_model    7 241.49 255.41 -113.75   227.49                       
#> m1$full_model   10 236.42 256.31 -108.21   216.42 11.073  3    0.01134 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2020-12-02 by the reprex package (v0.3.0)}

[Misch007]

Hey @singmann, sorry if this is an obvious question, I am new to R and statistics in general.. I just bumped into this issue, having even more deviating results than in your example (highly significant v.s not significant at all). I was wondering, which result is correct, after all? Best regards, Mischa

[singmann]

The issue above comes from refit = TRUE versus refit = FALSE. For afex mixed objects the default is refit = FALSE, whereas for lme4 objects the default is refit = TRUE (as can be seen by the message "refitting model(s) with ML (instead of REML)").

What is correct, depends on your use case. If the models differ in their fixed effects structure and they were fit with REML (instead of ML) then refit = TRUE has to be the case.

However, in the present example the fixed effects structure is the same, so the anova() results from afex mixed objects are correct and can be replicated with lme4 if refit = FALSE:

anova(m1$full_model, m2$full_model, refit = FALSE)
# Data: data
# Models:
# m2$full_model: score ~ Machine + (1 + re1.Machine1 + re1.Machine2 || Worker)
# m1$full_model: score ~ Machine + (Machine | Worker)
#               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
# m2$full_model    7 233.75 247.67 -109.88   219.75                     
# m1$full_model   10 230.51 250.40 -105.25   210.51 9.2424  3    0.02624
#                
# m2$full_model  
# m1$full_model *
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

So this issue shows actually expected behaviour and should probably be closed.

[singmann]

See also the comments on ?anova.merMod:

refit: logical indicating if objects of class lmerMod should be refitted with ML before comparing models. The default is TRUE to prevent the common mistake of inappropriately comparing REML-fitted models with different fixed effects, whose likelihoods are not directly comparable.

[cjungerius]

Commenting because I just ran into the same 'issue' - I do wonder if it isn't reasonable to match the defaults of anova.merMod, or at least add a warning, either when calling the function or in some documentation, rather than silently using the opposite behaviour. I must admit that when I used the anova call on two mixed objects and it ran without issue I (perhaps foolishly) assumed it was simply forwarding the .$full_model field to anova.merMod, especially because anova.mixed doesn't have any documentation. I only figured out there even was an anova.mixed method by calling methods(anova), and only found out about the refit option by going into the source code in this repo.

Maybe my use case isn't the default though, I understand if you think differently about this!