Restricted Likelihood for Structured Covariance Components

Why Restrict the Likelihood?

In a Gaussian linear model, covariance parameters are estimated while the mean includes unknown fixed effects. Ordinary maximum likelihood estimates those effects from the same observations used to estimate covariance, and this can lead to downward bias in variance estimation. Restricted maximum likelihood (REML) addresses the issue by forming likelihood from linear combinations of the data that eliminate the fixed effects. These combinations retain information about covariance while excluding the fitted mean directions.

The construction is especially useful for structured covariance models, where $V(\theta)$ may be built from variance components. REML does not make every component identifiable, nor does it guarantee that estimates are unbiased or interior. It provides a likelihood for covariance parameters after accounting for the dimension of the fixed-effect space.

Definition: Suppose $y\in\mathbb{R}^n$ follows $N(X\beta,V(\theta))$, where $X$ is a known $n\times p$ matrix of full column rank, $V(\theta)$ is positive definite, and $\beta$ is an unknown fixed-effect vector. An error contrast matrix is a matrix $K$ with $n-p$ rows, full row rank, and $KX=0$. The vector $Ky$ is an error contrast: its distribution has mean zero and covariance $KV(\theta)K'$.

The model assumptions matter. The Gaussian assumption gives an exact likelihood for $Ky$. Full column rank of $X$ makes the fixed-effect space have dimension $p$, while full row rank of $K$ makes the contrasts span its $(n-p)$-dimensional orthogonal complement. We assume the parameter values under consideration make $V(\theta)$ positive definite. For covariance components, this typically restricts the parameter space.

Deriving the Restricted Likelihood

Since $KX=0$, $Ky=K\varepsilon$ when $y=X\beta+\varepsilon$. Thus $Ky$ has a nonsingular multivariate normal distribution with mean zero and covariance $KVK'$, where $V=V(\theta)$. Its density gives a likelihood for $\theta$ that does not involve $\beta$.

L_R(\theta;Ky) =(2\pi)^{-(n-p)/2}|KV(\theta)K'|^{-1/2} \exp\!\left\{-\frac12 y'K'[KV(\theta)K']^{-1}Ky\right\}.

Different full-rank contrast matrices can describe the same fixed-effect-free information. Replacing $K$ by $AK$, with $A$ nonsingular, changes the density by a factor that depends on $A$, but not on $\theta$. Therefore the likelihood for covariance parameters is defined up to a parameter-independent factor. A standard form, omitting such constants, is

L_R(\theta;y)\ \propto\ |V(\theta)|^{-1/2} |X'V(\theta)^{-1}X|^{-1/2} \exp\!\left\{-\frac12 y'P(\theta)y\right\},

where

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

The matrix $P$ is symmetric and satisfies $PX=0$. The following derivation establishes both the contrast representation and the determinant terms in this expression.

Theorem: For any error contrast matrix $K$, the quadratic form and determinant in the Gaussian likelihood of $Ky$ satisfy $$ K'(KVK')^{-1}K=P,\qquad |KVK'|=|KK'|\frac{|V|\,|X'V^{-1}X|}{|X'X|}. $$ Consequently, the error-contrast likelihood and the standard restricted likelihood above have the same dependence on $\theta$.

Proof. To establish the matrix identity, let $V^{1/2}$ be the symmetric positive-definite square root of $V$, and define $A=KV^{1/2}$ and $B=V^{-1/2}X$. Because $AB=KX=0$, the row space of $A$ is orthogonal to the column space of $B$. Their dimensions are $n-p$ and $p$, respectively, so they are orthogonal complements in $\mathbb{R}^n$. The matrix $A'(AA')^{-1}A$ is the orthogonal projection onto the row space of $A$. It is therefore $I-B(B'B)^{-1}B'$. Multiplying on both sides by $V^{-1/2}$ gives

K'(KVK')^{-1}K =V^{-1/2}\{I-B(B'B)^{-1}B'\}V^{-1/2} =P.

For the determinant identity, first take $Q=[Q_1\ Q_2]$ orthogonal, where the columns of $Q_1$ span the columns of $X$ and those of $Q_2$ span the orthogonal complement. Write $X=Q_1R$ for an invertible $p\times p$ matrix $R$. The block determinant formula, applied to $Q'VQ$ and its inverse, gives

|Q_2'VQ_2| =|V|\,|Q_1'V^{-1}Q_1| =\frac{|V|\,|X'V^{-1}X|}{|X'X|}.

Every full-rank error contrast matrix can be written $K=A_0Q_2'$ for some nonsingular $(n-p)\times(n-p)$ matrix $A_0$. Hence $|KVK'|=|A_0|^2|Q_2'VQ_2|$ and $|KK'|=|A_0|^2$. Substitution yields the stated determinant identity. Inserting both identities into the density of $Ky$ leaves only factors independent of $\theta$ beyond the standard restricted likelihood. $\square$

Worked Example: An Intercept and an Independent Error Variance

Suppose $y_i=\beta+\varepsilon_i$, $i=1,\ldots,4$, with independent $\varepsilon_i\sim N(0,\sigma^2)$. The observations are $2,4,5,7$. An error contrast removes the unknown intercept; for example, the three differences from the first observation give $K y=(2,3,5)'$, with rows of $K$ equal to $(-1,1,0,0)$, $(-1,0,1,0)$, and $(-1,0,0,1)$.

Here $V=\sigma^2 I$. The quadratic form in the restricted likelihood is $y'Py=\sum_i(y_i-\bar y)^2$, because fitting the intercept leaves deviations from the sample mean. The mean is $\bar y=(2+4+5+7)/4=4.5$, so

$$ y'Py=(-2.5)^2+(-0.5)^2+(0.5)^2+(2.5)^2=6.25+0.25+0.25+6.25=13. $$

There are $n-p=3$ error contrasts. Up to a factor independent of $\sigma^2$, the restricted likelihood is

L_R(\sigma^2)\propto(\sigma^2)^{-3/2}\exp\!\left(-\frac{13}{2\sigma^2}\right).

Its log, omitting constants, is $-\tfrac32\log(\sigma^2)-13/(2\sigma^2)$. Differentiating with respect to $\sigma^2$ and setting the derivative to zero gives $-3/(2\sigma^2)+13/(2(\sigma^2)^2)=0$, hence $\hat{\sigma}^2_R=13/3$. The divisor is the number of residual degrees of freedom, not the sample size.

Estimating Covariance Components

In a structured covariance model, write $V(\theta)$ as a differentiable function of covariance-component parameters. For example, a linear variance-component model has $V(\theta)=\sum_j\theta_jA_j$, with known symmetric component matrices $A_j$ and parameter restrictions chosen to keep $V$ positive definite. The restricted log-likelihood, up to constants independent of $\theta$, is

\ell_R(\theta)=-\frac12\left\{\log|V|+\log|X'V^{-1}X|+y'Py\right\}.

The score has a useful form. Let $V_j=\partial V/\partial\theta_j$. The result below applies at an interior parameter value where $V$ is positive definite and the derivatives exist.

Theorem: The restricted-likelihood score for $\theta_j$ is $$ \frac{\partial\ell_R}{\partial\theta_j} =\frac12\left\{y'PV_jPy-\operatorname{tr}(PV_j)\right\}. $$ Thus an interior REML estimate solves these score equations for all covariance parameters.

Proof. For a differentiable positive-definite matrix, $d\log|V|=\operatorname{tr}(V^{-1}dV)$, and $dV^{-1}=-V^{-1}(dV)V^{-1}$. Set $C=X'V^{-1}X$. Then

d\log|C| =\operatorname{tr}(C^{-1}dC) =-\operatorname{tr}\{C^{-1}X'V^{-1}(dV)V^{-1}X\}.

Combining this with $d\log|V|$, and using invariance of trace under cyclic rearrangement, gives

d\{\log|V|+\log|C|\}=\operatorname{tr}(P\,dV).

For fixed $K$, the contrast representation $P=K'(KVK')^{-1}K$ and the inverse-matrix derivative give $dP=-P(dV)P$. It follows that $d(y'Py)=-y'P(dV)Py$. Substituting $dV=V_j\,d\theta_j$ into the differential of $\ell_R$ yields the stated score. $\square$

The score also clarifies why REML accounts for the fixed effects. Under the model, $E[y]=X\beta$ and $\operatorname{Cov}(y)=V$. Since $PX=0$ and $PVP=P$, the expectation of the quadratic term is

E[y'PV_jPy] =\operatorname{tr}(PV_jPV)+\beta'X'PV_jPX =\operatorname{tr}(PV_j).

The fixed-effect contribution vanishes, so the score has expectation zero at the true covariance parameter, subject to the stated model and regularity conditions. This does not imply that each REML estimate is unbiased, or that a score equation has a unique solution.

Worked Example: Two Random-Intercept Groups

Consider four observations in two groups of size two, with covariance $V=\tau^2ZZ'+\sigma^2I$, where $Z$ indicates group membership, $\tau^2\geq0$, and $\sigma^2>0$. Include an unknown overall fixed intercept, so $X=\mathbf{1}_4$. The observed values are $y=(0,2,4,8)'$.

An orthonormal set of three error contrasts consists of $u_1=(-1,1,0,0)'/\sqrt2$, $u_2=(0,0,-1,1)'/\sqrt2$, and $u_3=(1,1,-1,-1)'/2$. Each is orthogonal to $\mathbf{1}_4$. The first two are within-group contrasts, so the random group intercept cancels; their variances are $\sigma^2$. The third contrasts the group means, so its variance is $\sigma^2+2\tau^2$. The observed contrast values are

u_1'y=\sqrt2,\qquad u_2'y=2\sqrt2,\qquad u_3'y=-5.

Their squared values are $2,8,25$. The contrast likelihood is the product of three independent normal densities with the variances just calculated. For fixed $\sigma^2$, its first two factors are maximized at $\sigma^2=(2+8)/2=5$. Its third factor is maximized when $\sigma^2+2\tau^2=25$, giving $\hat{\tau}^2=(25-5)/2=10$. Both estimates satisfy the parameter restrictions. This example shows how within-group contrasts inform the residual variance while the between-group contrast supplies information about the group component.

Worked Example: A Single Structured Covariance Scale

Let $y\sim N(\beta\mathbf{1}_3,\theta A)$, where $A=\operatorname{diag}(1,2,1)$ is known and $\theta>0$. Suppose $y=(1,2,5)'$. There are $n-p=2$ error-contrast dimensions. The scale form $V=\theta A$ implies $P_V=\theta^{-1}P_A$, where

P_A=A^{-1}-A^{-1}\mathbf{1}_3 (\mathbf{1}_3'A^{-1}\mathbf{1}_3)^{-1}\mathbf{1}_3'A^{-1}.

Here $A^{-1}=\operatorname{diag}(1,1/2,1)$, $\mathbf{1}_3'A^{-1}\mathbf{1}_3=5/2$, $y'A^{-1}y=28$, and $\mathbf{1}_3'A^{-1}y=7$. Therefore

y'P_Ay=28-\frac{7^2}{5/2}=28-\frac{98}{5}=\frac{42}{5}=8.4.

Ignoring constants, the restricted log-likelihood as a function of $\theta$ is $-\tfrac{2}{2}\log\theta-\tfrac{8.4}{2\theta}$. Differentiation gives $-1/\theta+4.2/\theta^2=0$, so the restricted-likelihood estimate is $\hat{\theta}_R=4.2$. The known matrix $A$ determines the relative covariance pattern; the two error contrasts determine the scale estimate.

Identifiability and a Common Pitfall

A component is estimable by REML only if it affects the distribution of the error contrasts. For example, with a single group and $V=\tau^2\mathbf{1}\mathbf{1}'+\sigma^2I$ while $X=\mathbf{1}$, every error contrast $K$ satisfies $K\mathbf{1}=0$. Consequently $KVK'=\sigma^2KK'$, which contains no $\tau^2$. The random-intercept component is not identifiable from these contrasts when the same direction is already an unrestricted fixed intercept. The restriction does not create information that the model and design do not contain.

Common mistake: Counting all $n$ observations as independent information about the residual variance, then using a divisor of $n$ after fitting fixed effects. The fixed-effect directions have been estimated from the data, so the restricted likelihood is based on $n-p$ contrast dimensions. In the intercept example, this is why the estimate is $13/3$, not $13/4$. The correction is to form the restricted likelihood—or its equivalent projection-matrix expression—before optimizing covariance parameters.

REML is useful when covariance estimation should reflect the fixed-effect design, especially in variance-component models. Its likelihood is tied to the fixed-effect space: changing $X$ changes which directions are removed. Thus restricted likelihoods from different fixed-effect designs are not automatically comparable as though they were likelihoods for the same set of contrasts. REML estimates also remain subject to parameter constraints and possible boundary solutions; the score equation describes interior stationary points, not a guarantee that the optimum lies in the interior.

Check Your Understanding

Use the definitions and derivations above to answer these questions.

What conditions must an error contrast matrix $K$ satisfy, and why does $Ky$ have mean zero?
Explain why the restricted likelihood has $n-p$ rather than $n$ normal contrast dimensions.
In the two-group example, which contrasts contain information about $\tau^2$, and why?
For a covariance component with derivative $V_j$, write the restricted-likelihood score and identify its quadratic and trace terms.
Why can a random-intercept covariance component be unidentifiable when the model also includes an unrestricted fixed intercept in the same direction?

Restricted Likelihood for Structured Covariance Components

What You'll Learn