Why the Estimation Method Matters
A Gaussian mixed model separates the mean trajectory from the covariance among repeated observations. When comparing fixed-effect structures, such as a common time trend versus treatment-specific trends, the covariance model may be the same while the mean model differs. The choice between maximum likelihood (ML) and restricted maximum likelihood (REML) matters because the two methods do not produce likelihoods that are interchangeable for this comparison.
The essential rule is straightforward: fit both models by ML when using a likelihood-ratio test to compare different fixed-effect structures. REML is useful for estimating covariance parameters once the fixed-effect structure is specified, but its likelihood is constructed after accounting for the fixed effects. Changing those effects changes the restricted likelihood being evaluated.
Here, “marginal” means that any random effects have been integrated out. For ML, the marginal log likelihood is evaluated using the observed response vector and the model's mean and covariance. For a given \(\theta\), the generalized least-squares estimate of the fixed effects is \(\hat{\beta}(\theta)=(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}Y\), assuming \(X\) has full column rank and \(V\) is positive definite. ML then chooses \(\theta\) to maximize the likelihood after estimating \(\beta\).
REML instead focuses on contrasts \(K^\mathsf{T}Y\) satisfying \(K^\mathsf{T}X=0\). These contrasts have mean zero, so their distribution does not depend on \(\beta\). The restricted likelihood therefore uses information in a residual space determined by \(X\). In a standard full-rank formulation, its log likelihood includes a term involving \(\log|X^\mathsf{T}V^{-1}X|\), in addition to covariance and residual terms. If \(X\) changes, the contrast space and this term change too. This is the core reason REML likelihoods for different fixed-effect structures are not directly comparable by an ordinary likelihood-ratio test.
ML and REML in a Gaussian Model
For a fixed \(X\), the ML estimate of \(\beta\) at a given covariance parameter value is the generalized least-squares estimate above. Define \(P_V=V^{-1}-V^{-1}X(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}\). Substituting \(\hat{\beta}(\theta)\) into the marginal log likelihood gives the profiled ML log likelihood
REML uses the same residual quadratic form, but also adjusts for estimating the fixed effects. For full-rank \(X\), its log likelihood can be written
where \(p\) is the number of columns of \(X\). The displayed expressions make the distinction precise: REML's likelihood depends on the fixed-effect design, not just on the observed data and covariance model. These formulas assume a nonsingular covariance matrix and a full-rank fixed-effect design; software may report likelihoods using equivalent conventions that differ by constants.
Worked Example: ML and REML Variance Estimates
Suppose four independent hypothetical observations follow \(Y_i=\beta+\epsilon_i\), where \(\epsilon_i\sim N(0,\sigma^2)\). The observations are 1, 2, 4, and 5 units. This is the simplest Gaussian linear model and illustrates the finite-sample difference between ML and REML variance estimation.
The sample mean is \(\bar{Y}=(1+2+4+5)/4=3\). The residual sum of squares is \((1-3)^2+(2-3)^2+(4-3)^2+(5-3)^2=4+1+1+4=10\).
Under ML, the variance estimate divides the residual sum of squares by \(n\): \(\hat{\sigma}^2_{\mathrm{ML}}=10/4=2.50\) squared units. Under REML, the mean has been removed by the contrast information, leaving \(n-p=4-1=3\) independent residual contrasts, so \(\hat{\sigma}^2_{\mathrm{REML}}=10/3=3.333\) squared units.
The estimates differ because ML estimates the mean and variance together, while REML accounts for the degree of freedom used to estimate the mean. This calculation does not imply that REML is always preferable for every inferential task: it illustrates why the variance estimates can differ and why the method used for a model comparison must be considered explicitly.
Likelihood-Ratio Comparisons of Fixed Effects
Suppose a reduced model is nested within a full model: the reduced model is obtained by imposing \(q\) constraints on the full model's fixed-effect coefficients. Both models must use the same response observations and the same covariance-model specification, though the covariance parameters are estimated separately in each fit. Fit each model by ML and let \(\ell_F\) and \(\ell_R\) be their maximized ML log likelihoods.
The statistic measures how much the maximized likelihood improves when the additional fixed effects are allowed. A large value relative to the \(\chi^2_q\) reference distribution is evidence against the reduced-model constraints. The reference distribution is approximate in finite samples; the test does not establish that the full model is scientifically correct, nor does a small p-value measure the size or importance of the effect.
A short justification explains why the degrees of freedom are \(q\). Under the null, the reduced model's constraints are true. In a regular nested model, the log-likelihood difference near that null is approximately a quadratic form in the \(q\) directions released by the full model. Twice that difference has the limiting chi-square distribution with one degree of freedom for each independent constraint. If coefficients are not identifiable, the models are not nested, or relevant nuisance parameters lie at a boundary, this standard justification may fail.
Worked Example: Testing a Two-Group Time Interaction
In a hypothetical longitudinal study, consider two treatment groups and a linear time variable measured in months. The reduced model includes an intercept, treatment indicator, and common time slope. The full model adds treatment-by-time interaction, allowing the groups to have different time slopes. Thus the full model adds one coefficient: the difference between the two slopes.
Suppose both models use the same biomarker observations and the same random-intercept-and-slope covariance specification. Their ML maximized log likelihoods are \(\ell_R=-105.00\) for the reduced model and \(\ell_F=-102.15\) for the full model. The likelihood-ratio statistic is \(T=2[-102.15-(-105.00)]=2(2.85)=5.70\).
There is one added fixed-effect coefficient, so \(q=1\). The p-value is \(P(\chi^2_1\geq5.70)\approx0.01696\). At a 0.05 significance level, this is evidence against the common-slope restriction: the data are inconsistent with equal treatment-group time slopes under the fitted model and the test assumptions. The result is not evidence that treatment caused a difference; causal interpretation requires an appropriate study design and assumptions. The interaction coefficient itself is measured in biomarker units per month.
Counting Parameters and Checking the Comparison
The degrees of freedom count the independent fixed-effect restrictions, not the number of treatment groups or the number of random effects. For example, with three treatment groups and a linear time effect, a common-slope model allows one time slope. A group-specific-slope model has two additional slope differences relative to a reference group, so the interaction contributes \(q=2\). The count must follow from the model matrices or coefficient constraints, not from a guess.
Worked Example: Three Groups and Two Added Slopes
Suppose three hypothetical treatment groups are observed over time. The reduced model has a separate intercept for each group but a common linear time slope. The full model retains those intercepts and adds two time-by-group coefficients, one slope difference for each nonreference group. Therefore \(q=2\). Assume the ML log likelihoods are \(\ell_R=-88.40\) and \(\ell_F=-84.80\), with the same observations and covariance structure in both fits.
The likelihood-ratio statistic is \(T=2[-84.80-(-88.40)]=2(3.60)=7.20\). Under \(H_0\) that both added slope differences are zero, the reference distribution is \(\chi^2_2\). Its upper-tail probability is \(P(\chi^2_2\geq7.20)=e^{-7.20/2}=e^{-3.60}\approx0.0273\). This provides evidence at the 0.05 level that at least one of the two group slopes differs from the reference slope. It does not identify which slope difference is responsible; that requires estimates and suitable follow-up contrasts.
Before interpreting a comparison, verify the following conditions. First, both fits must use the same response records; different missing-data patterns or analysis subsets make the likelihoods refer to different observed data. Second, the covariance structures must represent the same model family, with the same random-effect and residual-covariance terms, so the comparison isolates the intended fixed-effect restrictions. Third, confirm nesting and count the added constraints. Finally, check convergence, parameter identifiability, and whether estimates are on boundaries. A warning, singular Hessian, or variance estimate at a boundary can undermine the usual likelihood-ratio approximation.
Alternatives, Limitations, and Common Pitfalls
| Approach | What it does | Important qualification |
|---|---|---|
| ML likelihood-ratio test | Compares nested fixed-effect models using their maximized ML likelihoods. | Requires nested models and a suitable chi-square approximation. |
| Wald test | Tests a coefficient or group of coefficients using estimates and their estimated covariance matrix. | Can be sensitive to finite-sample approximation and unstable standard errors. |
| REML estimation | Estimates covariance parameters after accounting for fixed effects. | Restricted likelihoods from different fixed-effect designs are not directly comparable in the ordinary likelihood-ratio calculation. |
For the final model, REML is often used to estimate covariance parameters because it adjusts for fixed effects. That common practice does not change the requirement to use ML for a likelihood-ratio comparison of different fixed-effect structures. After selecting a fixed-effect structure, refitting by REML may be appropriate for estimation, provided the model is otherwise suitable. Wald tests may be useful for targeted coefficient hypotheses, but they are not simply another name for the ML likelihood-ratio test.
A common mistake is to fit both fixed-effect models by REML and calculate \(2(\ell_{F,\mathrm{REML}}-\ell_{R,\mathrm{REML}})\), then compare it to a chi-square distribution. The reasoning may sound plausible because REML is designed to estimate covariance parameters well. The correction is that the two restricted likelihoods are based on different fixed-effect contrast spaces; their difference is not the ordinary likelihood-ratio statistic for testing the fixed effects. Fit both models by ML for that comparison.
Another pitfall is assigning the wrong degrees of freedom. With two groups and linear time, the interaction adds one slope difference, not two coefficients. With three groups, it adds two differences. The number of added parameters follows from the particular coding and constraints, but the number of independent restrictions is invariant to a full-rank recoding of the same model.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A two-group model adds one treatment-by-linear-time coefficient to a reduced model. How many degrees of freedom does the usual likelihood-ratio test use?
- Why should two models with different fixed-effect designs not be compared by their ordinary REML likelihood difference?
- For ML log likelihoods \(-54.2\) and \(-57.0\) for full and reduced nested models, respectively, what is the likelihood-ratio statistic?
- Name two checks to make before interpreting a fixed-effect likelihood-ratio comparison.
Solutions. (1) One: the interaction adds one independent slope difference, so \(q=1\). (2) REML uses contrasts determined by the fixed-effect design; changing that design changes the restricted likelihood being evaluated. (3) \(T=2[-54.2-(-57.0)]=2(2.8)=5.60\). Its reference degrees of freedom equal the number of independent added constraints, which must be specified before obtaining a p-value. (4) For example, verify the models use the same observations and covariance structure, confirm nesting and the parameter count, and check convergence and boundary estimates.