Potential error with model.comparison() for mixed models?

[tonycarlsen]

Hi Dustin,
I've noticed the following issue when using the model.comparison() function with mixed models whereby if the model is not specified with REML = F it provides different output compared to when you do specify REML=F - even though model.comparison() suggests that it is refitting the models with ML instead of REML:

model <- lmer(Power_Output ~ Training_Hours + Diet + (Training_Hours | Athlete_ID), data = athlete_data)
model2 <- lmer(Power_Output ~ Training_Hours * Diet + (Training_Hours | Athlete_ID), data = athlete_data)
model.comparison(model, model2)

refitting model(s) with ML (instead of REML)
$statistics
aic bic bayes.factor p
model 2435.166 2461.093 0.009 <2e-16
model2 2421.948 2451.579 116.403

$predicted_differences
0% 25% 50% 75% 100%
0.001 0.131 0.328 0.600 2.013

$r_squared_change
(Intercept) Training_Hours Residual
0.0002032730 0.2423298854 0.0006849655

anova(model, model2)

refitting model(s) with ML (instead of REML)
Data: athlete_data
Models:
model: Power_Output ~ Training_Hours + Diet + (Training_Hours | Athlete_ID)
model2: Power_Output ~ Training_Hours * Diet + (Training_Hours | Athlete_ID)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model 7 2446.8 2472.7 -1216.4 2432.8
model2 8 2435.2 2464.8 -1209.6 2419.2 13.618 1 0.000224

Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

#---AIC/BIC appear simular to anova after re-specifying with REML=F

model <- lmer(Power_Output ~ Training_Hours + Diet + (Training_Hours | Athlete_ID), data = athlete_data, REML = F)
model2 <- lmer(Power_Output ~ Training_Hours * Diet + (Training_Hours | Athlete_ID), data = athlete_data, REML = F)
model.comparison(model, model2)

$statistics
aic bic bayes.factor p
model 2446.771 2472.697 0.019 <2e-16
model2 2435.153 2464.783 52.313

$predicted_differences
0% 25% 50% 75% 100%
0.002 0.117 0.334 0.625 2.136

$r_squared_change
(Intercept) Training_Hours Residual
-0.0000533067 0.2605783618 0.0006725911