[R-sig-ME] differences between lmer and glmmTMB - follow up

[Ben Bolker]

A few things going on here:

1.  It is simply wrong/meaningless to compare REML criteria across 
models with different fixed effects.  lme4 tries to protect you from 
this, to some extent (e.g. the anova.merMod() method refits REML models 
via ML before comparing), but it can't entirely protect you.  So the 
issue in this e-mail is moot.  (In other words, *don't use AIC tables to 
compare models with different fixed effects fitted by REML*.)

2. lme4::lmer() uses REML by default, glmmTMB uses ML by default.  At 
least for the reproducible example below, if we call glmmTMB with 
REML=TRUE we get answers that are the same up to very small 
(floating-point imprecision) differences:

library(lme4)
library(glmmTMB)
m1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
m2 <- glmmTMB(Reaction ~ Days + (Days|Subject), sleepstudy, REML = TRUE)

## results are equivalent

all.equal(logLik(m1), logLik(m2))
coef(summary(m1))
coef(summary(m2))$cond

## compare differences in REML criterion

m3 <- update(m1, . ~ . - (Days|Subject) + (1|Subject))
m4 <- update(m2, . ~ . - (Days|Subject) + (1|Subject))
logLik(m1) - logLik(m3)
logLik(m2) - logLik(m4)

On 2022-11-30 12:26 PM, Don Cohen wrote:
> [Summary: two models that differ only in one added fixed effect
>   the deviance of the one with extra fixed effect coming out higher]
>
> I think this is a related problem:
>
>> mfull <- lmer(data = strict, LN_GCS ~ poly(hour,2) + poly(group_size,2) + poly(age,2) + poly(dom_rank,2) + Pregnancy + takeover + (1|ID) + (1|groupid) + (1|AssayNum))
>> summary(mfull)
> Linear mixed model fit by REML. t-tests use Satterthwaite's method [
> lmerModLmerTest]
> Formula: LN_GCS ~ poly(hour, 2) + poly(group_size, 2) + poly(age, 2) +
>      poly(dom_rank, 2) + Pregnancy + takeover + (1 | ID) + (1 |
>      groupid) + (1 | AssayNum)
>     Data: strict
> REML criterion at convergence: 2088.6
> ...
>> m <- lmer(data = strict, LN_GCS ~ poly(hour,2) + poly(group_size,2) + poly(age,2) + poly(dom_rank,2) + Pregnancy + (1|ID) + (1|groupid) + (1|AssayNum))
>> summary(m)
> Linear mixed model fit by REML. t-tests use Satterthwaite's method [
> lmerModLmerTest]
> Formula: LN_GCS ~ poly(hour, 2) + poly(group_size, 2) + poly(age, 2) +
>      poly(dom_rank, 2) + Pregnancy + (1 | ID) + (1 | groupid) +
>      (1 | AssayNum)
>     Data: strict
> REML criterion at convergence: 2086.4
> ...
>
> The REML criterion is supposeed to be -2 x loglik (deviance), right?
> How can it be lower for m than for mfull when mfull contains all the same
> effects (plus one more) ?
>
>> AICctab(mfull, m, weights = T)
>        dAICc df weight
> m      0.0  15 0.89
> mfull  4.2  16 0.11
>
> If I use glmmTMB I get the expected results - same loglik, AIC
> different by 2
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models