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

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.

• Does the regress function from the same-named package actually compute the Welham-Thompson statistic ? I don't see any reference to it on cran.r-project.org/web/packages/regress/regress.pdf Commented Nov 13, 2023 at 13:13
• @RobertLong Yes, that's my understanding. regress is what P. McCullagh uses in the Ten Projects book: stat.uchicago.edu/~pmcc/projects (data + R code). In the text he makes an argument for the Welham-Thompson approach & I got interested but don't understand the theory yet. (I know that not everyone has access to the book; I'll look for open source materials.) Commented Nov 13, 2023 at 13:42
• From pp. 20 in Ten Projects: "The Welham-Thompson likelihood-ratio statistic on two degrees of freedom for testing the null hypothesis of additivity can be computed as follows:" fit <- regress(...). Commented Nov 13, 2023 at 13:43

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.
• Thanks, esp. for the final part. I'm not sure about the 2nd point, as I wouldn't in general insist on REML estimation for its own sake but points 1 and 3 do make sense to me. Commented Nov 14, 2023 at 21:30