Restricted Likelihood and Error-Contrast Derivation

Why restrict the likelihood?

In a Gaussian mixed model, the fixed effects describe the mean while variance components describe the covariance. Estimating both from the same data creates a finite-sample issue: fitting the mean uses information that would otherwise contribute to estimating variation. Restricted maximum likelihood (REML) addresses this by constructing likelihood from linear combinations of the observations that contain no fixed-effect contribution. These combinations are called error contrasts.

The central idea is precise: remove the directions in the data space that can be explained by the fixed-effects design, then estimate covariance parameters from the remaining directions. The method does not guarantee that every variance-component estimate is unbiased, but it corrects the particular loss of residual degrees of freedom that makes the ordinary maximum-likelihood estimate of a Gaussian residual variance too small.

Model and error contrasts

Let \(y\) be an \(n\)-vector of longitudinal measurements, and consider the marginal Gaussian model

y \sim N(X\beta,V(\theta)).

Here \(X\) is a known \(n\times p\) fixed-effects design matrix of full column rank, \(\beta\) is an unknown \(p\)-vector, and \(V(\theta)\) is a positive-definite covariance matrix indexed by variance parameters \(\theta\). A mixed model with Gaussian random effects has this form after the random effects are integrated out. Assume \(n>p\), so there are observations left after removing the fixed-effect directions.

Definition: An error-contrast matrix is an \(n\times(n-p)\) full-column-rank matrix \(K\) such that \(K'X=0\). The vector \(K'y\) is an error contrast: its mean is zero because \(K'X\beta=0\). Its distribution is \(K'y\sim N(0,K'V(\theta)K)\).

The columns of \(K\) span the part of the data space orthogonal to the columns of \(X\). Since \(K\) has \(n-p\) columns, its transformed data retain \(n-p\) independent directions. Different valid choices of \(K\) give different coordinates for the same contrast space; they lead to likelihoods that differ at most by a factor independent of \(\theta\), so they give the same REML estimates.

Worked Example: Removing an Intercept

Suppose three hypothetical biomarker measurements from one subject follow \(y_i=\beta+\varepsilon_i\), with independent errors \(\varepsilon_i\sim N(0,\sigma^2)\). The observed values are \(y=(2,4,6)'\). Here \(X=(1,1,1)'\), \(n=3\), and \(p=1\).

Choose the orthonormal contrast columns

K= \begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{6}\\ -1/\sqrt{2} & 1/\sqrt{6}\\ 0 & -2/\sqrt{6} \end{pmatrix}.

Each column sums to zero, so \(K'X=0\), and \(K'K=I_2\). The two contrasts are

K'y= \begin{pmatrix} (2-4)/\sqrt{2}\\ (2+4-2(6))/\sqrt{6} \end{pmatrix} = \begin{pmatrix} -\sqrt{2}\\ -6/\sqrt{6} \end{pmatrix}.

Their squared length is \((- \sqrt{2})^2+(-6/\sqrt{6})^2=2+6=8\). Because \(V=\sigma^2 I_3\), the contrasts have covariance \(K'VK=\sigma^2 I_2\). Their likelihood therefore estimates \(\sigma^2\) as \(8/2=4\). The contrast construction has removed the unknown intercept without first treating its fitted value as known.

Deriving the restricted likelihood

Because \(K'y\) is Gaussian with mean zero and covariance \(K'V(\theta)K\), its log likelihood, apart from no omitted parameter-dependent terms, is

\ell_R(\theta) = -\frac{1}{2}\left[ (n-p)\log(2\pi) +\log|K'V(\theta)K| +y'K\{K'V(\theta)K\}^{-1}K'y \right].

This is the error-contrast form of the restricted log likelihood. The quadratic term measures the contrast data against their covariance; the log determinant accounts for the scale and shape of that covariance. Maximizing this likelihood over \(\theta\) gives the REML estimates.

A frequently used equivalent expression makes the removal of fixed effects explicit without displaying \(K\). Define

P=V^{-1}-V^{-1}X(X'V^{-1}X)^{-1}X'V^{-1}.

For full-rank \(X\) and positive-definite \(V\), the REML log likelihood can be written, up to a constant independent of \(\theta\), as

\ell_R(\theta) = -\frac{1}{2}\left[ (n-p)\log(2\pi) +\log|V| +\log|X'V^{-1}X| +y'Py \right].

The matrix \(P\) removes the fitted fixed-effect contribution in the covariance-weighted geometry. In particular, \(PX=0\), so \(y'Py\) depends on the component of \(y\) not explained by \(X\beta\). The determinant involving \(X'V^{-1}X\) is essential: REML is not obtained merely by replacing \(\beta\) with an estimate in the ordinary likelihood.

Derivation: Choose \(Q_1\) with orthonormal columns spanning the column space of \(X\), and \(K=Q_2\) with orthonormal columns spanning its orthogonal complement. Then \([Q_1\ Q_2]\) is orthogonal, and \(X=Q_1R\) for an invertible \(p\times p\) matrix \(R\). The block determinant identity gives \(|K'VK|=|V|\,|Q_1'V^{-1}Q_1|\). Also, \(|X'V^{-1}X|=|R|^2|Q_1'V^{-1}Q_1|\) and \(|X'X|=|R|^2\). Thus \(\log|K'VK|=\log|V|+\log|X'V^{-1}X|-\log|X'X|\). The last term depends only on the fixed design, not on \(\theta\). The contrast quadratic is \(y'Py\), yielding the stated criterion up to that design-only constant.

Worked calculations with fixed-effect directions

Worked Example: Contrasts for an Intercept and Time Trend

Consider four hypothetical measurements at times \(0,1,2,3\), with mean \(\beta_0+\beta_1t_i\). Thus the rows of \(X\) are \((1,0),(1,1),(1,2),(1,3)\). For the observations \(y=(5,7,8,12)'\), use the contrast columns

k_1=\frac{1}{2}(1,-1,-1,1)',\qquad k_2=\frac{1}{\sqrt{20}}(1,-3,3,-1)'.

Both columns have entries summing to zero. Their time-weighted sums are \(0-1-2+3=0\) and \(0-3+6-3=0\), so each is orthogonal to both columns of \(X\). Their lengths are one, and their inner product is \((1+3-3-1)/(2\sqrt{20})=0\). Therefore \(K=(k_1,k_2)\) is an orthonormal basis for the two-dimensional error-contrast space.

The observed contrasts are \(k_1'y=(5-7-8+12)/2=1\) and \(k_2'y=(5-21+24-12)/\sqrt{20}=-4/\sqrt{20}\), approximately \(-0.8944\). If the observation errors are independent with variance \(\sigma^2\), then \(K'VK=\sigma^2I_2\). The contrasts contain no intercept or linear time trend, but they do retain information about departures from that mean trajectory. In this balanced, independent-error example, the contrast sum of squares is \(1^2+(-4/\sqrt{20})^2=1.8\).

Why the adjustment can reduce variance-estimation bias

The simplest case shows the source of the adjustment exactly. Suppose \(y\sim N(X\beta,\sigma^2 I)\), with \(X\) of rank \(p\). The residual sum of squares after least squares fitting is \(RSS=y'(I-H)y\), where \(H=X(X'X)^{-1}X'\). The matrix \(I-H\) is an orthogonal projection of rank \(n-p\). Under the model, \(RSS/\sigma^2\) has a chi-square distribution with \(n-p\) degrees of freedom; consequently \(E(RSS)=(n-p)\sigma^2\).

The ordinary Gaussian maximum-likelihood estimate divides by \(n\), whereas the restricted likelihood divides by the number of residual directions, \(n-p\). The difference and its expectation are shown directly below.

Worked Example: ML and REML for a Residual Variance

Suppose a hypothetical Gaussian model has \(n=8\) observations and \(p=3\) fixed-effect coefficients. Assume its residual variance is \(\sigma^2=4\), so \(E(RSS)=(8-3)4=20\). The ML estimator is \(\hat{\sigma}^2_{\mathrm{ML}}=RSS/n\), and hence \(E(\hat{\sigma}^2_{\mathrm{ML}})=20/8=2.5\). Its bias is \(2.5-4=-1.5\) squared measurement units.

The REML estimator in this model is \(\hat{\sigma}^2_{\mathrm{REML}}=RSS/(n-p)\). Its expectation is \(20/5=4\), so it is unbiased under these assumptions. The adjustment is not an arbitrary inflation: it accounts for the three fitted fixed-effect directions and uses the five remaining directions to estimate residual variation.

In more general mixed models, REML uses the full covariance \(K'V(\theta)K\), not simply a residual sum of squares divided by a degrees-of-freedom count. The contrast likelihood adjusts for fixed effects, but it does not make every component estimate exactly unbiased in every design or covariance model. Finite-sample bias, weak identification, and boundary estimates can still occur.

Conditions, comparisons, and failure modes

The derivation requires a correctly specified Gaussian marginal model, a known fixed-effects design of full column rank, \(n>p\), and a positive-definite covariance matrix for the likelihood as written. The covariance model must also contain enough information in the contrast space to estimate the parameters of interest. These conditions concern the validity and identifiability of REML; they do not ensure that the chosen mean and covariance structures are scientifically appropriate.

Approach	What is used	Main implication
ML	The full likelihood for \(y\), estimating fixed effects and covariance parameters together	Allows likelihood comparisons of models with different fixed-effect structures; the usual residual variance estimate can be downward biased in finite samples.
REML	The likelihood of \(n-p\) error contrasts	Adjusts covariance estimation for fitted fixed-effect directions; restricted likelihoods generally should not be compared across different fixed-effect designs.

For a random-intercept model with one subject and an intercept already included as a fixed effect, every error contrast sums to zero. The random intercept contributes a covariance proportional to \(11'\), and \(K'11'K=0\). Thus those contrasts contain no information about that random-intercept variance. This is not a numerical defect in REML: with only one subject, the subject-level random-intercept variation cannot be separated from the overall intercept using these data.

A practical failure mode is a rank-deficient \(X\), often caused by redundant columns or an empty factor level. Then \(p\) is not the stated number of columns, the inverse in the criterion is unavailable, and the contrast dimension must be based on the actual rank. Another warning sign is a variance parameter that has little effect on \(K'V(\theta)K\); it may be weakly identified even if software returns an estimate. Check the design rank, the number of independent contrasts, and whether the covariance parameter changes their covariance before interpreting a fitted REML value.

Common mistake

Wrong reasoning: “REML is ML with the fitted fixed effects substituted, so its only change is to the residual sum of squares.” This misses the likelihood normalization and the determinant adjustment for uncertainty in estimating \(\beta\). The correct construction is the likelihood of \(K'y\), whose distribution has covariance \(K'V(\theta)K\); equivalently, the criterion includes both \(\log|V|\) and \(\log|X'V^{-1}X|\). Also, REML is not a universal cure for bias: its exact unbiasedness in the simple residual-variance example does not extend automatically to each variance component in a complex model.

Key takeaway: REML estimates covariance parameters from linear combinations that eliminate the fixed-effect mean. The lost \(p\) directions appear explicitly in the restricted likelihood, explaining the residual degrees-of-freedom adjustment and why REML and ML can yield different variance estimates.

Check Your Understanding

Answer each question, then compare with the solutions.

If \(X\) has rank \(p\) and \(n=12\), how many columns does a full error-contrast matrix have?
Why does \(K'y\) have mean zero under the model \(y\sim N(X\beta,V)\)?
In the independent-error model with \(n=10\), \(p=2\), and \(\sigma^2=3\), what is \(E(RSS)\), and what is the expected ML residual-variance estimate?
Why should REML likelihoods generally not be used to compare two models with different fixed-effect designs?

Solutions. (1) It has \(n-p=12-p\) columns. (2) Its mean is \(K'X\beta=0\), by the defining condition \(K'X=0\). (3) \(E(RSS)=(n-p)\sigma^2=8(3)=24\); therefore \(E(RSS/n)=24/10=2.4\), below the true variance of 3. (4) The contrast spaces, and thus the data being assigned restricted likelihoods, change with the fixed-effect design; the resulting likelihood values are not directly comparable as if they were based on the same restricted data.