REML as Likelihood of Error Contrasts

Why Restrict the Likelihood?

A linear mixed model represents a response as a fixed mean structure plus random variation. In genetics, that random variation can include covariance associated with genetic similarity. Maximum likelihood estimates the model parameters by using the likelihood of the full response vector. But that vector also contains information used to estimate the fixed effects. Restricted maximum likelihood (REML) takes a different route: it forms contrasts that remove the fixed effects, then constructs a likelihood from those contrasts.

This perspective gives REML its name and explains its adjustment for the mean structure. The likelihood is restricted to combinations of observations that cannot contain information about the fixed-effect coefficients. The resulting contrast distribution still depends on the covariance parameters, which are the target of variance-component estimation.

Model and Error Contrasts

Assume the response vector \(y\) follows the Gaussian linear mixed model \(y=X\beta+\varepsilon\), where \(y\) has \(n\) entries, \(X\) is a known \(n\)-by-\(p\) design matrix of full column rank, and \(\beta\) is an unknown vector of fixed-effect coefficients. The random vector \(\varepsilon\) has mean zero and covariance matrix \(V(\theta)\), which is positive definite and depends on covariance parameters \(\theta\). Thus \(y\sim N(X\beta,V(\theta))\). The design matrix is treated as fixed, and the model assumes the stated Gaussian distribution and covariance structure.

For example, a genetic mixed model might have \(V(\theta)=\sigma_g^2G+\sigma_e^2I\), where \(G\) is a specified genetic relationship matrix, \(I\) is the identity matrix, and the nonnegative parameters \(\sigma_g^2\) and \(\sigma_e^2\) describe genetic and residual variation. The key construction below applies to any positive-definite \(V(\theta)\); it does not require this particular covariance form.

Definition: An error contrast is a linear combination of the observations \(a'y\) whose coefficient vector \(a\) satisfies \(a'X=0\). More generally, a full set of error contrasts is \(K'y\), where \(K\) is an \(n\)-by-\((n-p)\) matrix of rank \(n-p\) satisfying \(K'X=0\). These conditions make the contrasts eliminate every fixed-effect contribution.

Since \(K'X=0\), the contrasts have mean \(K'X\beta=0\), whatever the value of \(\beta\). Their covariance is \(K'V(\theta)K\), so their distribution is \(K'y\sim N(0,K'V(\theta)K)\). Because \(V(\theta)\) is positive definite and \(K\) has full column rank, this contrast covariance is positive definite. The contrast likelihood therefore depends on \(\theta\), but not on \(\beta\).

Likelihood of the Contrasts

Write \(m=n-p\) for the number of linearly independent error contrasts. The Gaussian likelihood for the \(m\)-dimensional vector \(K'y\), with \(K\) held fixed, is

L_R(\theta;y)\ \propto\ |K'V(\theta)K|^{-1/2} \exp\left\{-\frac{1}{2}(K'y)'[K'V(\theta)K]^{-1}(K'y)\right\}.

Taking its logarithm gives the restricted log-likelihood, up to constants that do not depend on \(\theta\):

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

This is a likelihood for the contrasts, not for the original observations treated as if the fixed effects were known. If a different full-rank contrast basis is used, it describes the same contrast space. The numerical log-likelihood can shift by a constant independent of \(\theta\), but the covariance-parameter estimates are unchanged.

Worked Example: Removing an Intercept from Three Observations

Suppose three hypothetical responses are \(y=(2,4,6)'\), the only fixed effect is an intercept, and the covariance model is \(V=\sigma^2I\), with \(\sigma^2>0\). Here \(X=(1,1,1)'\), so a valid orthonormal contrast matrix is \(K=\big((1,-1,0)'/\sqrt{2},\ (1,1,-2)'/\sqrt{6}\big)\). Each column sums to zero, and the columns are linearly independent.

The contrast values are \(K'y=(-2/\sqrt{2},-6/\sqrt{6})'=(-\sqrt{2},-\sqrt{6})'\). Their squared length is \(2+6=8\). Since \(K'VK=\sigma^2I_2\), the restricted log-likelihood is

\ell_R(\sigma^2) =-\frac{1}{2}\left\{2\log(2\pi)+2\log(\sigma^2)+\frac{8}{\sigma^2}\right\}.

For comparison, the full-data maximum-likelihood estimate of the intercept is \(\bar{x}=4\), and the residual sum of squares is \((2-4)^2+(4-4)^2+(6-4)^2=8\). The full-data likelihood estimates \(\sigma^2\) as \(8/3\), because it divides that sum by \(n=3\). The contrast likelihood estimates it as \(8/2=4\), because it has \(n-p=2\) independent contrasts. This simple case shows how accounting for the estimated intercept changes the variance estimate; it does not imply that REML is unbiased in every model.

From Contrasts to the Standard REML Expression

The contrast definition is useful conceptually, but computations often use an expression involving the full covariance matrix. Let \(P=V^{-1}-V^{-1}X(X'V^{-1}X)^{-1}X'V^{-1}\), with the dependence on \(\theta\) understood. The full-rank condition on \(X\) and positive definiteness of \(V\) ensure that \(X'V^{-1}X\) is invertible. For any full-rank \(K\) satisfying \(K'X=0\), the following identities hold:

y'Py=(K'y)'(K'VK)^{-1}(K'y), \qquad |K'VK|=\frac{|K'K|}{|X'X|}|V|\,|X'V^{-1}X|.

The determinant factor \(|K'K|/|X'X|\) depends on the chosen fixed design and contrast basis, but not on \(\theta\). Substituting these identities into the contrast log-likelihood and dropping only terms independent of \(\theta\) gives the standard form:

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

Why the Fixed-Effect Adjustment Appears

A short argument using the special case \(V=\sigma^2I\) shows the role of the adjustment. Let \(H=X(X'X)^{-1}X'\) be the projection onto the column space of \(X\). The residual sum of squares after fitting the fixed effects is \(y'(I-H)y\). Under the Gaussian model, it is the squared length of \(n-p\) independent error contrasts, each with variance \(\sigma^2\). Its expectation is therefore

E\{y'(I-H)y\}=(n-p)\sigma^2.

The full-data likelihood, after fitting \(\beta\), divides this residual sum of squares by \(n\) when estimating \(\sigma^2\); the contrast likelihood divides it by \(n-p\). The contrast count reflects that \(p\) dimensions of the data have been used to estimate the fixed-effect mean structure. In a general mixed model, the adjustment is not merely a change from \(n\) to \(n-p\): the terms \(\log|V|\), \(\log|X'V^{-1}X|\), and \(y'Py\) together express the likelihood and its covariance geometry.

In particular, the term \(\log|X'V^{-1}X|\) is part of the restricted-likelihood adjustment. It accounts for how the covariance model weights the fixed-effect design. Omitting it while retaining the other terms does not give the REML likelihood.

Worked Example: A Genetic Covariance and a Fixed Genotype Trend

Consider four hypothetical individuals with responses \(y=(2,5,7,8)'\) and a fixed intercept plus a genotype trend \(g=(0,0,1,1)'\). Thus \(X\) has columns \((1,1,1,1)'\) and \(g\), and \(p=2\). Let \(V=\sigma_g^2G+\sigma_e^2I\), where \(\sigma_g^2,\sigma_e^2>0\) and

G= \begin{pmatrix} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{pmatrix}.

This particular relationship matrix makes \(V=(\sigma_g^2+\sigma_e^2)I\); it is chosen to keep the contrast calculation transparent. The two vectors \(k_1=(1,-1,0,0)'/\sqrt{2}\) and \(k_2=(0,0,1,-1)'/\sqrt{2}\) each have zero sum within their genotype group. Consequently, both are orthogonal to the intercept and to \(g\), and they form a full-rank error-contrast basis.

The observed contrasts are \(k_1'y=(2-5)/\sqrt{2}=-3/\sqrt{2}\) and \(k_2'y=(7-8)/\sqrt{2}=-1/\sqrt{2}\). Their squared length is \(9/2+1/2=5\). With \(s=\sigma_g^2+\sigma_e^2\), their covariance is \(K'VK=sI_2\), so their log-likelihood is

\ell_R =-\frac{1}{2}\left\{2\log(2\pi)+2\log(s)+\frac{5}{s}\right\}.

The likelihood uses within-genotype contrasts, so neither the unknown intercept nor the unknown genotype trend appears in their mean. This example does not distinguish the two variance components: because \(G=I\), the contrast likelihood depends on them only through their sum. The example illustrates the contrast construction, not a general claim that every genetic covariance model has this limitation.

Worked Example: Changing the Contrast Basis

Return to the three-observation intercept model with \(y=(2,4,6)'\). Let \(K\) be the orthonormal basis used earlier, and form a new contrast matrix \(\widetilde K=KA\), where \(A\) is the diagonal matrix with entries \(2\) and \(1\). This is still a full-rank basis for the same contrast space and satisfies \(\widetilde K'X=0\).

The new contrast vector is \(\widetilde K'y=A'K'y=(-2\sqrt{2},-\sqrt{6})'\), and its covariance is \(\widetilde K'V\widetilde K=\sigma^2A'A\). Its determinant is \(4\sigma^4\), while the quadratic term is \[ (\widetilde K'y)'(\sigma^2A'A)^{-1}(\widetilde K'y)=\frac{8}{\sigma^2}. \] Thus the new log-likelihood has the same \(\sigma^2\)-dependent terms as before, with the determinant term increased by \(\log(4)\). That shift is constant in \(\sigma^2\), so maximizing either contrast likelihood gives the same covariance-parameter estimate.

Interpretation and a Common Pitfall

REML matters when covariance parameters are estimated in the presence of fixed effects, as they are in genetic mixed models with covariates or other mean terms. It bases the likelihood on information that is free of those fixed-effect coefficients and includes an adjustment for the dimension and weighting of the fixed-effect design. This is why the restricted likelihood is not simply the ordinary likelihood with fitted fixed effects substituted into it.

Common mistake: “Once the fixed effects have been estimated, they can be treated as known, so ordinary maximum likelihood already accounts for them.” The correction is that estimating the fixed effects uses information from the data. Treating their fitted values as known does not account for that use. REML instead uses \(n-p\) independent contrasts and the corresponding restricted-likelihood adjustment. This does not guarantee an unbiased estimate in every covariance model, nor does it remove the need for the model assumptions to be reasonable.

Key takeaway: REML is the likelihood of linear combinations \(K'y\) satisfying \(K'X=0\). Those contrasts eliminate the unknown fixed effects while retaining covariance information; the standard REML expression includes a fixed-effect adjustment through \(\log|X'V^{-1}X|\) and the \(n-p\) contrast dimension.

Check Your Understanding

Answer each question using the model and contrast definitions above.

In a model with \(n=12\) observations and a full-rank fixed-effect design with \(p=3\) columns, how many independent error contrasts are required?
Why does \(K'y\) have a mean that does not depend on \(\beta\) when \(K'X=0\)?
Under \(V=\sigma^2I\), why does the restricted likelihood use \(n-p\), rather than \(n\), in its variance information?
What changes in the contrast log-likelihood if the contrast basis is replaced by another full-rank basis for the same space?
Why is the fixed-effect adjustment in REML relevant to genetic variance-component estimation?