I'm trying to compare a set of four linear mixed-effect models (run in R, lme4::lmer), and, judging by what I've read, there seems to be some debate as to the best way of doing so. I was hoping someone might provide some much needed input.
I started with two models:
m1 <- lmer(IA_DT_PC ~ CONDITION * IA_LABEL + (1 | PARTICIPANT), data2)
m2 <- lmer(IA_DT_PC ~ GROUP * CONDITION * IA_LABEL + (1 | PARTICIPANT), data2)
My usual approach is to use the anova function, so I began by comparing the first two, which showed no difference.
> anova(m1,m2, refit = FALSE)
Data: data2
Models:
object: IA_DT_PC ~ CONDITION * IA_LABEL + (1 | PARTICIPANT)
..1: IA_DT_PC ~ GROUP * CONDITION * IA_LABEL + (1 | PARTICIPANT)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
object 6 -5670.2 -5634.4 2841.1 -5682.2
..1 10 -5638.6 -5578.8 2829.3 -5658.6 0 4 1
So I figured the between-subject factor GROUP wasn't adding to the overall fit of the model. So I put that aside and fitted a third model with another factor of interest:
m3 <- lmer(IA_DT_PC ~ CONDITION * VERSION * IA_LABEL + (1 | PARTICIPANT), data2)
Now, comparing m1 to m3 did reveal a difference, suggesting that m3 is a better fit:
> anova(m1,m3, refit = FALSE)
Data: data2
Models:
object: IA_DT_PC ~ CONDITION * IA_LABEL + (1 | PARTICIPANT)
..1: IA_DT_PC ~ CONDITION * VERSION * IA_LABEL + (1 | PARTICIPANT)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
object 6 -5670.2 -5634.4 2841.1 -5682.2
..1 10 -5680.1 -5620.3 2850.1 -5700.1 17.882 4 0.001302 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
However, because I'm a fundamentally curious person, I also tried digging that GROUP factor up again, and included it in a fourth model:
m4 <- lmer(IA_DT_PC ~ GROUP * CONDITION * VERSION * IA_LABEL + (1 | PARTICIPANT), data2)
The confusing part (well, to me anyway) is that m4 actually appears to be a better fit than m3 now:
> anova(m3,m4, refit = FALSE)
Data: data2
Models:
object: IA_DT_PC ~ CONDITION * VERSION * IA_LABEL + (1 | PARTICIPANT)
..1: IA_DT_PC ~ GROUP * CONDITION * VERSION * IA_LABEL + (1 | PARTICIPANT)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
object 10 -5680.1 -5620.3 2850.1 -5700.1
..1 18 -5711.6 -5604.0 2873.8 -5747.6 47.502 8 0.000000123 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
So, my question is, why did GROUP not explain any additional variance when added to the m1 structure, but it does when added to the m3 structure?
Somewhat related questions:
1) When is it (in)appropriate to allow anova() to refit the models when comparing them?
2) Is there an altogether better way of comparing models that I might try?
lmer
default. UseREML = FALSE
in thelmer()
call. $\endgroup$refit = FALSE
, the models are refit with ML. I take it that's the correct way in this case, and the ANOVA call will take the non-REML values provided by lmer? $\endgroup$