Estimating Variance Components by REML

Why Estimate Variance Components by REML?

In a genetic mixed model, the covariance among observations is partly attributed to a genomic relationship matrix and partly to residual variation. Estimating these sources of variation requires accounting for fixed effects, such as an intercept or covariates. Restricted maximum likelihood (REML) does this by basing its likelihood on contrasts that eliminate the fixed-effect mean. This tutorial develops the estimation steps for one genetic covariance kernel and shows how to assess an iterative fit.

The One-Kernel Model and Its REML Likelihood

Let \(y\) be an \(n\)-vector of quantitative outcomes, \(X\) an \(n \times p\) fixed-effect design matrix of full column rank, and \(K\) a known positive semidefinite genomic relationship matrix. Use the model

where \(g\) and \(e\) are independent, \(\beta\) is an unknown fixed-effect vector, and \(\sigma_g^2\) and \(\sigma_e^2\) are nonnegative genetic and residual variance components. Thus \(y\sim N(X\beta,V)\), with \(V=\sigma_g^2K+\sigma_e^2I_n\). The model assumes the covariance structure \(K\) is specified and that the residuals are independent with a common variance. For positive \(\sigma_e^2\), \(V\) is positive definite.

Choose an \(n \times m\) matrix \(Q\) whose columns are orthonormal and span the space orthogonal to the columns of \(X\), where \(m=n-p\). Then \(Q'X=0\), so the contrast vector \(z=Q'y\) does not depend on \(\beta\). Its covariance is

Let \(Q'KQ=U\operatorname{diag}(\lambda_1,\ldots,\lambda_m)U'\), where \(U\) is orthogonal and each \(\lambda_j\geq 0\). Define \(w=U'z\). The transformed contrasts are independent under the model, with \(w_j\sim N(0,v_j)\), where \(v_j=\lambda_j\sigma_g^2+\sigma_e^2\). The REML log likelihood, omitting a constant that does not depend on the variance components, is therefore

Scores and the Average-Information Update

Write \(\theta_g=\sigma_g^2\), \(\theta_e=\sigma_e^2\), and \(d_{jg}=\lambda_j,\ d_{je}=1\). Differentiating the restricted log likelihood gives the score component for either variance parameter \(\theta_i\):

The two terms inside the sum compare the observed squared contrast with its model variance. A positive score indicates that a small increase in the corresponding component would increase the log likelihood locally, holding the other component fixed.

The average-information (AI) matrix averages observed and expected information. In the diagonal contrast representation its entries are

An AI step solves \(AI\,\Delta=U\), where \(U=(U_g,U_e)'\), and proposes \(\theta^{new}=\theta+\Delta\). This is a local update, not a guarantee that the proposed point is admissible or increases the likelihood. In practice, an implementation must keep variance estimates nonnegative and check the likelihood after proposed steps.

Expectation-Maximization Updates

An EM update treats the genetic and residual parts of each transformed contrast as latent components. Write \(w_j=a_j+e_j\), where \(a_j\sim N(0,\sigma_g^2\lambda_j)\) and \(e_j\sim N(0,\sigma_e^2)\) are independent. For \(\lambda_j>0\), the conditional means and variances given \(w_j\) are

Also, \(E(e_j\mid w_j)=\sigma_e^2w_j/v_j\) and \(\operatorname{Var}(e_j\mid w_j)=\sigma_g^2\lambda_j\sigma_e^2/v_j\). The expected complete-data variance estimates yield the following simultaneous EM updates, with \(r\) the number of positive eigenvalues:

A zero eigenvalue contributes no latent genetic variance; its contrast is entirely residual. EM updates are often stable and, under the stated model with exact updates, do not decrease the restricted likelihood. They can converge slowly near an optimum. AI commonly takes larger, faster local steps, but its proposed step needs admissibility and likelihood checks.

Convergence Checks in a One-Kernel Fit

A practical fit iterates updates until parameter changes and the score are small, while tracking the restricted likelihood. The score should be assessed relative to the scale of the problem; a small absolute score alone can be misleading when parameters or information differ greatly in scale. Also check that the current variances are admissible and that the AI matrix can be solved reliably. Multiple starting values can help reveal whether an apparent stopping point depends on initialization.

Common Pitfall and Practical Interpretation

A frequent error is to treat REML as if it used all \(n\) outcomes without accounting for the fixed effects. The restricted likelihood is built from \(m=n-p\) contrasts, and the projected kernel \(Q'KQ\)—not the unprojected \(K\) alone—determines the contrast eigenvalues. Ignoring this adjustment changes the likelihood being optimized and can lead to incorrect variance-component updates.

AI and EM are alternative numerical routes to estimate the same model parameters; neither changes what \(\sigma_g^2\) and \(\sigma_e^2\) mean under the assumed covariance model. A converged estimate is conditional on the specified fixed effects, kernel, and residual assumptions. It describes a fitted covariance decomposition, not a causal effect of genetic variation.