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

Question

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

Answer 1

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

Answer 2

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.