Residual likelihood ratio test for fixed effects in a linear mixed model

Question

I know (but now I have doubts) that "Comparing models that are fitted with REML and differ in their fixed effects never makes sense," just as @BenBolker explains in this answer.

I've been reading the book Ten Projects in Applied Statistics by P. McCullagh and, if I understand correctly, he doesn't agree with the above statement. I quote from Chapter 18, Residual Likelihood, to make sure the math is right.

The optional argument REML=FALSE (in lme4::lmer) is a cop-out, which overrides the default, and reverts to ordinary maximum likelihood instead. This option produces a valid likelihood-ratio statistic, which is not one recommended by Welham and Thompson (1997) or by this author.

The function regress::regress(y ~ X, ~ block + V, kernel = K) has a three-part syntax, permitting greater flexibility, in which the setting for kernel determines the method of estimation. The first part is a standard model-formula for the mean-value subspace $\mathcal{X}$; the second part, which may be empty or missing, is a simple model formula for the covariances. For the third part, the default kernel is $\mathcal{K} = \mathcal{X}$, i.e., REML, not $\mathcal{K} = 0$, i.e., maximum likelihood. (...) For the comparison of mean-values $H_0: \mu \in \mathcal{X}_0$ versus $H_1: \mu \in \mathcal{X}_1 \supset \mathcal{X}_0$ as alternative, residual likelihood may be used in the following manner [comments added]:

> 
> X0 <- model.matrix( ~ mf0); # smaller model matrix
> X1 <- model.matrix( ~ mf1); # larger model matrix (X1 includes X0)
> fit0 <- regress(y ~ mf0, ~ block + V, kernel=X0);
> fit1 <- regress(y ~ mf1, ~ block + V, kernel=X0);
> # likelihood ratio statistic for the X1 fixed effect
> 2 * (fit1$llik - fit0$llik);
> 

There are two part to this question: theory and software.

1. Theory: On one hand, ML gives a valid likelihood-ratio statistic for model comparisons, so there is no need to change favored lmer-based approach. On the other hand, the Welham-Thomposon likelihood-ratio statistic is "more valid" esp. for small samples. Should we prefer the latter over the standard LM approach, and if yes, under what conditions?

2. Software (may not be on-topic on CV but I'm adding it for completeness): We can't get the Welham-Thompson statistic with lme4::lmer. We need to use regress::regress instead.

References:

The Ten Projects book is not freely available online; the rest of the materials are.

• P. McCullagh. Ten Projects in Applied Statistics. Springer Series in Statistics. Springer, 2022.
• S. J. Welham and R. Thompson. Likelihood ratio tests for fixed model terms using residual maximum likelihood. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59(3):701–714, 1997. https://doi.org/10.1111/1467-9868.00092
• D. Clifford and P. McCullagh. The R journal: The regress function. R News, 6(2):R News, 2006. https://journal.r-project.org/articles/RN-2006-008/
• P. McCullagh. REML and residual likelihood. Nelder Lecture, Imperial College, March 8 2012 http://www.stat.uchicago.edu/~pmcc/seminars/nelder/nelder.pdf
• regress() adapted from David Clifford's original function by Peter McCullagh. http://www.stat.uchicago.edu/~pmcc/courses/regress.R

Answer 1

The debate over comparing mixed-effects models fitted by REML (Restricted Maximum Likelihood) against ML (Maximum Likelihood) focuses on the optimal conditions for model comparison, especially when fixed effects are involved.

ML vs REML for Model Comparison:

• ML estimates the parameters in the model by maximising the likelihood of the observed data. It is suitable for comparing nested models that differ in their fixed effects because the likelihood ratio test (LRT) relies on the likelihoods being computed in the same way for both models.

• Restricted Maximum Likelihood (REML): REML modifies the ML estimation to account for the fixed effects degrees of freedom. REML is often used to estimate variance components in models. When comparing models with various fixed effects, however, REML estimates are not directly comparable since the probability of each model is adjusted differently based on its fixed effects.

Welham-Thompson Likelihood-Ratio Statistic:

This statistic refers to a modified version of the LRT that is used to compare models with different fixed effects but are estimated using REML. Because it accounts for the degrees of freedom associated with fixed effects, it is thought to provide a more accurate test, particularly in small samples.

Software lme4 vs. regress:

• The lme4 package defaults to REML estimation.To compare models with different fixed effects, you should refit the models using ML.
• The regress function from the regress package allows greater flexibility in fitting models and can compute the Welham-Thompson statistic, although the version available on CRAN appears not to do so, but apparently there is data a code here.

When to Prefer Welham-Thompson over Standard ML? You might prefer the Welham-Thompson approach over the standard ML approach when:

• You are dealing with small sample sizes where the variance component estimates might be biased.
• You require a more accurate test for the significance of fixed effects under REML estimation.
• You are comparing models that only differ by fixed effects and want to maintain REML estimation for consistency in variance component estimation.