Solved – REML or ML to compare two mixed effects models with differing fixed effects, but with the same random effect

fixed-effects-modellme4-nlmemaximum likelihoodmixed modelrandom-effects-model

Background: Note: My data set and R code are included below text

I wish to use AIC to compare two mixed effects models generated using the lme4 package in R. Each model has one fixed effect and one random effect. The fixed effect differs between models, but the random effect remains the same between models. I've found that if I use REML=TRUE, model2 has the lower AIC score, but if I use REML=FALSE, model1 has the lower AIC score.

Support for using ML:

Zuur et al. (2009; p. 122) suggest that "To compare models with nested fixed effects (but with the same random structure), ML estimation must be used and not REML." This indicates to me that I ought to use ML since my random effects are the same in both models, but my fixed effects differ. [Zuur et al. 2009. Mixed Effect Models and Extensions in Ecology with R. Springer.]

Support for using REML:

However, I notice that when I use ML, the residual variance associated with the random effects differs between the two models (model1 = 136.3; model2 = 112.9), but when I use REML, it is the same between models (model1=model2=151.5). This implies to me that I ought instead to use REML so that the random residual variance remains the same between models with the same random variable.

Question:

Doesn't it make more sense to use REML than ML for comparisons of models where the fixed effects change and the random effects remain the same? If not, can you explain why or point me to other literature that explains more?

# Model2 "wins" if REML=TRUE:
REMLmodel1 = lmer(Response ~ Fixed1 + (1|Random1),data,REML = TRUE)
REMLmodel2 = lmer(Response ~ Fixed2 + (1|Random1),data,REML = TRUE)
AIC(REMLmodel1,REMLmodel2)
summary(REMLmodel1)
summary(REMLmodel2)

# Model1 "wins" if REML=FALSE:
MLmodel1 = lmer(Response ~ Fixed1 + (1|Random1),data,REML = FALSE)
MLmodel2 = lmer(Response ~ Fixed2 + (1|Random1),data,REML = FALSE)
AIC(MLmodel1,MLmodel2)
summary(MLmodel1)
summary(MLmodel2)

Dataset:

Response    Fixed1  Fixed2  Random1
5.20    A   A   1
32.50   A   A   1
6.57    A   A   2
24.77   A   B   3
41.69   A   B   3
34.29   A   B   4
1.80    A   B   4
10.00   A   B   5
15.56   A   B   5
4.44    A   C   6
21.65   A   C   6
9.20    A   C   7
4.11    A   C   7
12.52   B   D   8
0.25    B   D   8
27.34   B   D   9
11.54   B   E   10
0.86    B   E   10
0.68    B   E   11
4.00    B   E   11

Best Answer

Zuur et al., and Faraway (from @janhove's comment above) are right; using likelihood-based methods (including AIC) to compare two models with different fixed effects that are fitted by REML will generally lead to nonsense.

Faraway (2006) Extending the linear model with R (p. 156):

The reason is that REML estimates the random effects by considering linear combinations of the data that remove the fixed effects. If these fixed effects are changed, the likelihoods of the two models will not be directly comparable

These two questions discuss the issue further: Allowed comparisons of mixed effects models (random effects primarily) ; REML vs ML stepAIC

Related Question