REML Estimation of Variance Components

Why Estimate Variance Components by REML?

A linear mixed model separates systematic mean differences from covariance among phenotypes. In a GWAS setting, the covariance is often represented by a genetic component linked to a genomic relationship matrix and a residual component. Their estimates matter because they determine the fitted covariance used to account for relatedness. Estimating covariance after accounting for fixed effects is therefore not just a computational detail: it changes the likelihood being optimized.

Restricted maximum likelihood, or REML, does this by forming combinations of the observed phenotypes that remove the fixed-effect mean. These combinations retain information about covariance while avoiding treating estimated fixed effects as known. REML is not a method for testing a variant effect; here its role is estimating covariance parameters in the mixed model.

Model and assumptions: Let \(y\) be an \(n\)-vector of phenotypes and \(X\) an \(n\) by \(p\) fixed-effect design matrix of full column rank. The model is \(y=X\beta+u+\varepsilon\), where \(u\sim N(0,\sigma_g^2K)\), \(\varepsilon\sim N(0,\sigma_e^2I)\), and \(u\) and \(\varepsilon\) are independent. The relationship matrix \(K\) is treated as known and positive semidefinite, and \(\sigma_g^2\geq0\), \(\sigma_e^2>0\). Thus \(y\sim N(X\beta,V)\), with \(V=\sigma_g^2K+\sigma_e^2I\). REML relies on this Gaussian covariance model and on the specified \(X\) and \(K\) being appropriate for the data.

Contrasts Remove the Fixed-Effect Mean

Choose an \(n\) by \(m\) matrix \(A\), where \(m=n-p\), whose columns are orthonormal and span the space orthogonal to the columns of \(X\). Thus \(A'X=0\) and \(A'A=I_m\). The contrast vector \(z=A'y\) has mean zero because \(E(z)=A'X\beta=0\), and covariance \(C=A'VA\). If \(C\) is positive definite, \(z\) has a nonsingular multivariate normal distribution. Its likelihood is the restricted likelihood.

Definition: The REML log likelihood, expressed using orthonormal error contrasts, is

\ell_R(\sigma_g^2,\sigma_e^2) =-\frac12\left\{(n-p)\log(2\pi)+\log|A'VA|+(A'y)'(A'VA)^{-1}(A'y)\right\}.

REML estimates the variance components by maximizing this likelihood over their allowable values. The fixed-effect coefficient \(\beta\) does not appear because \(A'X\beta=0\).

This contrast formulation also clarifies the phrase fixed-effect-adjusted data. REML does not simply subtract an ordinary least-squares fit and then pretend the residuals are independent observations. The contrasts are constructed to eliminate the fixed-effect mean, and their covariance \(A'VA\) must still be included. Different valid orthonormal bases for the same contrast space give the same likelihood value.

A commonly used equivalent expression, up to a parameter-independent constant, is

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

The matrix \(P\) removes the fitted fixed-effect space in the covariance-weighted geometry. In particular, \(PX=0\), so \(y'Py\) contains no contribution from adding a vector in the column space of \(X\). The determinant involving \(X'V^{-1}X\) is essential: it accounts for the fixed-effect dimensions removed from the likelihood.

How the REML Score Arises

For a covariance parameter \(\theta_i\), write \(V_i=\partial V/\partial\theta_i\). Differentiating the restricted log likelihood gives the score

\frac{\partial\ell_R}{\partial\theta_i} =\frac12\left\{y'PV_iPy-\operatorname{tr}(PV_i)\right\}.

The first term measures covariance-weighted variation in the observed data in directions affected by \(\theta_i\). The trace term is the corresponding model-based contribution. At an interior maximum, each score is zero; in practice, numerical optimization is often used because the component parameters jointly affect both determinants and the quadratic term.

Specify the mean and covariance.
Choose \(X\), \(K\), and \(V=\sigma_g^2K+\sigma_e^2I\), and state the Gaussian and independence assumptions.

Remove the fixed-effect mean.
Use a contrast basis \(A\) with \(A'X=0\), or the equivalent \(P\)-matrix likelihood.

Optimize over the components.
Maximize \(\ell_R\) subject to the variance constraints, then examine whether the fitted covariance and numerical optimization are credible.

A short justification of the score formula is useful. Since \(P\) is symmetric, the derivative of the quadratic term is \(\partial(y'Py)/\partial\theta_i=-y'PV_iPy\). The derivative of the restricted log-determinant terms is \(\operatorname{tr}(PV_i)\). Combining these derivatives with the factor \(-1/2\) in \(\ell_R\) gives the stated score. The identity for the quadratic derivative follows by differentiating \(P\), using \(\partial V^{-1}/\partial\theta_i=-V^{-1}V_iV^{-1}\), and accounting for the fitted fixed-effect projection in \(P\).

Worked REML Calculations

Worked Example: Estimating Two Components From Contrasts

Consider three hypothetical phenotypes with an intercept as the only fixed effect, so \(X=(1,1,1)'\) and \(n-p=2\). Define the orthonormal contrasts \(a=(1,-1,0)'/\sqrt{2}\) and \(b=(1,1,-2)'/\sqrt{6}\). Let \(K=aa'\), and suppose the observed phenotype vector is \(y=2a+b\), approximately \((1.822,-1.006,-0.816)'\). The contrasts are \(z=A'y=(2,1)'\), for \(A=(a,b)\), and \(A'KA=\operatorname{diag}(1,0)\). The model assumptions are those stated above; here \(K\) is positive semidefinite and \(A'VA\) is positive definite whenever \(\sigma_e^2>0\).

The contrast covariance and restricted log likelihood are

A'VA= \begin{pmatrix} \sigma_g^2+\sigma_e^2&0\\ 0&\sigma_e^2 \end{pmatrix}, \qquad \ell_R=-\frac12\left\{2\log(2\pi)+\log(\sigma_g^2+\sigma_e^2)+\log(\sigma_e^2) +\frac{4}{\sigma_g^2+\sigma_e^2}+\frac{1}{\sigma_e^2}\right\}.

Let \(v_1=\sigma_g^2+\sigma_e^2\) and \(v_2=\sigma_e^2\). For a single observed squared contrast \(q>0\), the variance-dependent log-likelihood term \(-\tfrac12(\log v+q/v)\) is maximized at \(v=q\): its derivative is \(-1/(2v)+q/(2v^2)\), which is zero at \(v=q\) and changes from positive to negative there. Thus the unconstrained contrast-variance estimates are \(\hat v_1=2^2=4\) and \(\hat v_2=1^2=1\). They imply \(\hat\sigma_e^2=1\) and \(\hat\sigma_g^2=4-1=3\), both satisfying the variance constraints.

At these estimates, the determinant is \(4\times1=4\), and the quadratic term is \(2^2/4+1^2/1=2\). Therefore

\ell_R=-\frac12\{2\log(2\pi)+\log(4)+2\} =-\frac12(7.062)\approx-3.531.

This example works because the two contrast directions have different known genetic-covariance contributions. If both directions had the same contribution, these data would not separate the genetic and residual variances in this way.

Worked Example: Why REML Adjusts for the Fixed-Effect Dimension

Suppose four hypothetical independent phenotypes follow \(y_i=\beta+\varepsilon_i\), with \(\varepsilon_i\sim N(0,\sigma^2)\), and the observed values are \(1,2,3,6\). The sample mean is \(3\), so the squared deviations sum to \((1-3)^2+(2-3)^2+(3-3)^2+(6-3)^2=4+1+0+9=14\). Here there is one fixed-effect parameter and three independent error contrasts.

For example, an orthonormal contrast basis has three columns orthogonal to \((1,1,1,1)'\). The sum of squared contrast values is 14, regardless of which such basis is used. The restricted log likelihood, apart from its constant, is

\ell_R=-\frac12\left\{3\log(\sigma^2)+\frac{14}{\sigma^2}\right\}.

Differentiating with respect to \(\sigma^2\) gives \(-3/(2\sigma^2)+14/(2(\sigma^2)^2)=0\), so \(\hat\sigma^2_{\mathrm{REML}}=14/3\approx4.667\) squared phenotype units. The ordinary maximum-likelihood estimate, after estimating the mean, is \(14/4=3.500\) squared phenotype units. The difference comes from REML using the three-dimensional contrast space rather than acting as if the fitted mean had been known in advance. In this simple normal model the REML estimate is unbiased; that exact property should not be assumed for every mixed-model variance component.

Worked Example: Checking a Component Score

Suppose two orthonormal contrasts have genetic-covariance eigenvalues 2 and 0, the residual variance is fixed at 1 squared phenotype unit, and the observed contrast values have squares 8 and 1. The first contrast variance is \(1+2\tau\), where \(\tau=\sigma_g^2\), and the second is 1. The restricted log likelihood as a function of \(\tau\geq0\), omitting constants, is

\ell_R(\tau)=-\frac12\left\{\log(1+2\tau)+\frac{8}{1+2\tau}+1\right\}.

Its derivative is

\frac{d\ell_R}{d\tau} =-\frac{1}{1+2\tau}+\frac{8}{(1+2\tau)^2} =\frac{7-2\tau}{(1+2\tau)^2}.

The score is zero at \(\hat\tau=3.500\) squared phenotype units. It is positive below 3.5 and negative above 3.5, so this is the maximum on the allowed range. This check illustrates how a variance component changes the covariance in a particular contrast direction; the estimate is not obtained by treating the raw phenotype values as independent observations.

Comparisons, Diagnostics, and Limitations

Maximum likelihood and REML both use the Gaussian mixed-model covariance, but they use different likelihoods. Maximum likelihood estimates fixed effects and variance components jointly, then accounts for the fitted fixed effects in its likelihood. REML instead uses the likelihood of fixed-effect-free contrasts. In settings with fixed effects, this distinction often reduces the downward tendency of ML variance estimates, as the preceding simple example shows. REML does not guarantee unbiased estimates for arbitrary covariance structures, and REML likelihood values should not be used as though they were ML likelihood values when comparing models with different fixed-effect designs.

The components are estimable only to the extent that the covariance patterns they induce can be distinguished in the contrast space. If \(K=I\), then \(V=(\sigma_g^2+\sigma_e^2)I\): only the sum is identified, not the two components separately. More generally, if the genetic covariance pattern is nearly indistinguishable from residual variation after removing \(X\), the likelihood may be shallow in some parameter directions. A numerical optimizer can return estimates in this setting, but the estimates may be sensitive to starting values or small changes in the data.

Useful checks include confirming that \(X\) has full column rank, that the relationship matrix has the intended scaling, and that the fitted covariance is positive definite. Inspect optimizer convergence and whether different reasonable starting values lead to similar likelihoods and estimates. A nearly flat profile or very large changes in estimates with little change in likelihood indicate weak information, not necessarily a software failure. Also check that \(K\) and \(V\) align with the assumed phenotype order; a sample-order mismatch can produce a plausible-looking but incorrect fit.

REML depends on the assumed mean and covariance structure. Unmodeled fixed effects can be absorbed into covariance estimates, while a misspecified relationship matrix can make the fitted genetic covariance unrepresentative. Gaussianity and independence of \(u\) and \(\varepsilon\) are modeling assumptions, not facts guaranteed by using REML. Finally, variance parameters are constrained, so numerical optimization must respect those constraints; behavior at constraint boundaries requires separate interpretation.

Common Mistake

Pitfall: “I have subtracted the fitted fixed effects, so I can fit a likelihood to the ordinary residuals as if they were independent observations.” This is generally wrong. Fitted residuals are constrained by the estimation of \(\beta\), and their covariance is not usually \(V\). REML instead uses contrasts \(A'y\) with covariance \(A'VA\), or the equivalent restricted likelihood with the determinant correction. Another error is to interpret a component estimate as well determined merely because an optimizer returned a number; identification and convergence must be checked.

Key takeaway: REML estimates variance components from the \(n-p\) directions that remain after removing the fixed-effect mean. In a genomic mixed model, it balances the observed covariance-weighted variation against the variation expected under \(V=\sigma_g^2K+\sigma_e^2I\).

Check Your Understanding

Work through each item, then compare with the solutions below.

Why does \(A'y\) have mean zero when \(A'X=0\)?
In the first worked example, what would the residual variance estimate be if the observed contrast values were \(z=(2,2)'\) while the covariance directions remained \(A'KA=\operatorname{diag}(1,0)\)?
For the four observations \(1,2,3,6\) with an unknown common mean, calculate the REML estimate of the residual variance and state its units.
Why can genetic and residual variance components not be distinguished when \(K=I\)?
State one reason a converged REML fit could still provide unreliable component estimates.

Solutions. (1) Since \(E(y)=X\beta\), \(E(A'y)=A'X\beta=0\). (2) The contrast variances maximize separately at \(v_1=2^2=4\) and \(v_2=2^2=4\). Hence \(\hat\sigma_e^2=4\) and \(\hat\sigma_g^2=4-4=0\); the genetic and residual estimates lie at the constraint boundary. (3) The squared deviations sum to 14 and the contrast dimension is \(4-1=3\), so \(\hat\sigma^2_{\mathrm{REML}}=14/3\approx4.667\) squared phenotype units. (4) When \(K=I\), \(V=(\sigma_g^2+\sigma_e^2)I\), so the covariance depends only on the sum of the two parameters. (5) For example, if the covariance patterns induced by \(K\) and \(I\) are nearly indistinguishable in the contrast space, the likelihood can be flat and the estimates weakly identified despite numerical convergence.