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

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

• Faraway's (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." – jvh_ch Sep 25 '14 at 16:56
• Even though AIC is likelihood based, to the best of my knowledge, it was developed for the purposes of prediction. How would one exactly apply a mixed model for prediction? – AdamO Sep 25 '14 at 20:47
• @AdamO, could you be more precise? A fitted mixed model can be used for prediction, either at the population level (predict responses for an unspecified/unknown unit by setting conditional modes/BLUPs to zero) or at the individual level (condition prediction on the estimates of the conditional modes/BLUPs). If you can be more specific, that might make a good new CV question. – Ben Bolker Sep 25 '14 at 22:10
• It was just unclear to me how you were intended to apply this model. Nothing in the problem suggested what kind of prediction, if any, was being done or if it was necessary and if so for what purpose. – AdamO Sep 25 '14 at 23:03

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

• Thanks @janhove, AdamO and Ben Bolker. I also found this link from Aaron to be helpful in answering this question. It says, "The REML likelihood depends on which fixed effects are in the model, and so are not comparable if the fixed effects change. REML is generally considered to give better estimates for the random effects, though, so the usual advice is to fit your best model using REML for your final inference and reporting." – It Figures Sep 26 '14 at 12:04
• @BenBoker following up on your response, would it be acceptable to use the anova (LRT) to compare model these models with different fixed effects if REML = FALSE? – Helio Sep 8 '20 at 12:58

I'll give an example to illustrate why the REML likelihood cannot be used for things like AIC comparisons. Imagine that we a normal mixed effects model. Let $X$ denote the design matrix and assume that this matrix has full rank. We can find a reparametrization of the mean value space, given by the matrix $\tilde{X}$. The two matrices span the same linear subspace of $\mathbb{R}^n$. Thus, the columns of $\tilde{X}$ can be written as linear combinations of the columns of $X$. Therefore, we can find a quadratic matrix, $B$, such that

$\tilde{X} = XB$.

Furthermore, $B$ has full rank (this can be proven by assuming that it didn't; then neither would $X$, a contradiction). This means that $B$ is invertible.

If we start out be using the second parametrization of the mean value space and let $V$ be a covariance matrix then let's consider the REML criterion we should maximize (I'm omitting a constant)

$|V|^{-1/2}|\tilde{X}'V^{-1}\tilde{X}|^{-1/2}\exp((y-\tilde{X}\tilde{\beta})'V^{-1}(y-\tilde{X}\tilde{\beta})/2)$,

over the parameter set, where $\beta = (\tilde{X}V^{-1}\tilde{X})^{-1}y$. Using the fact that $X = \tilde{X}B$, we can realize that this can be rewritten as

$|B||V|^{-1/2}||X'V^{-1}X|^{-1/2}|\exp((y-X\bar{\beta})'V^{-1}(y-X\bar{\beta})/2)$,

where $\bar{\beta} = (XV^{-1}X)^{-1}y$. This is the REML likelihood for the other parametrization times the determinant of $|B|$.

We therefore have an example of two different parametrizations of the same model, giving different likelihood values, assuming that $|B| \neq 1$ (such a matrix can easily be found). The same parameter value will maximize the criterion in both cases but the value of the likelihood will be different. This shows that there is an arbitrary element in the likelihood value and therefore illustrates why one cannot use the value of the likelihood for comparison between models with different fixed effects: you would be able to change the results simply be changing the mean value space parametrization in one of the models.

This an example of why REML should not be used when comparing models with different fixed effects. REML, however, often estimates the random effects parameters better and therefore it is sometimes recommended to use ML for comparisons and REML for estimating a single (perhaps final) model.