Comparing linear mixed-effect models

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?

• If you are testing two mixed models that differ in their fixed effects via some likelihood based score (which you are here), you should be fitting the model with maximum likelihood. REML is the lmer default. Use REML = FALSE in the lmer() call.
– MAB
Jul 5, 2017 at 17:53
• @MAB if using a recent version of lmer that is no longer necessary. lmer will provide non-REML values to ANOVA, c.f. stackoverflow.com/questions/22892063/… Jul 6, 2017 at 4:52
• I'm using v 1.1.13, which I think is the latest version. If I don't use 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? Jul 6, 2017 at 13:31

It probably has to do with your usage of '*'. In your m1 model, you define the following effect: CONDITION * IA_LABEL. It translates to:

Main effects: CONDITION, IA_LABEL.

2-way interactions: CONDITION:IA_LABEL

So, by using '*' you include the two-way interaction and all lower-order effects. In your m2 model, you define the following effect: GROUP * CONDITION * IA_LABEL. It translates to:

Main effects: GROUP, CONDITION, IA_LABEL.

2-way interactions: GROUP:CONDITION, GROUP:IA_LABEL, CONDITION:IA_LABEL

3-way interaction: GROUP:CONDITION:IA_LABEL

So, when comparing m1 to m2, you compare a model that does not include any effect of GROUP to a model that includes the main effect of GROUP and all possible interactions between GROUP and the remaining variables (in addition to the other main effects and interactions).

In your m4 model, you specify the four-way interaction GROUP * CONDITION * VERSION * IA_LABEL, including all lower order effects. It is likely that m4 has a better fit compared to m3 because there is one (or more) interaction that explains variance and involves GROUP and VERSION. This interaction was not specified in the m2 model and may therefore not perform better compared to the m1 model.