Why Restricted Likelihood Is Used
A bivariate genomic mixed model separates covariance between two traits into contributions associated with genetic relatedness and residual variation. The model also includes fixed effects, such as trait means or covariates. Restricted maximum likelihood, usually abbreviated REML, estimates covariance components while accounting for the loss of information involved in estimating those fixed effects. It does so by forming linear combinations of the observations that remove the fixed-effect means.
The focus here is the construction and use of that restricted likelihood. Its value depends on the assumed genomic relationship matrix, the fixed-effect design, and the covariance components. A REML estimate is obtained by maximizing the restricted likelihood over admissible covariance parameters; it is not, in general, obtained by separately estimating each covariance from an ordinary sample covariance matrix.
The Bivariate Genomic Mixed Model
Suppose \(n\) individuals have measurements on two traits. Stack the trait-specific vectors in trait-major order as \(y=(y_1',y_2')'\), so \(y\) has length \(N=2n\). Let \(X\) be the \(N\times p\) fixed-effect design matrix and \(\beta\) its coefficient vector. Let \(g\) collect the genetic effects for both traits and \(e\) the residuals. The model is
where \(K\) is a known \(n\times n\) genomic relationship matrix, \(G\) is the \(2\times2\) genetic covariance matrix across traits, and \(R\) is the \(2\times2\) residual covariance matrix across traits. Assume \(g\) and \(e\) are independent, and that \(K\), \(G\), and \(R\) are symmetric, with \(K\), \(G\), and \(R\) positive semidefinite. The \((a,b)\) entry \(G_{ab}\) is the genetic covariance component on the scale defined by \(K\); for individual \(i\), the genetic covariance between traits \(a\) and \(b\) is \(G_{ab}K_{ii}\), and between individuals \(i\) and \(j\) it is \(G_{ab}K_{ij}\). The residual covariance matrix \(R\) describes within-individual residual covariance and, under this model, residuals from different individuals have zero covariance.
After integrating over the random effects, the observations have marginal distribution \(y\sim N(X\beta,V)\), where
For the ordinary likelihood and REML formulas below, assume \(X\) has full column rank \(p\) and \(V\) is positive definite. Positive-semidefinite \(G\) and \(R\) alone do not guarantee that \(V\) is invertible. A positive-definite \(R\), together with positive-semidefinite \(K\) and \(G\), is one sufficient condition for \(V\) to be positive definite. Singular covariance cases require separate treatment rather than direct use of formulas involving \(V^{-1}\) and \(|V|\).
Because \(V\) is positive definite and \(C\) has full column rank, \(C'VC\) is positive definite. Consequently, \(z\) has an ordinary multivariate normal density. Its log likelihood, viewed as a function of the covariance parameters, is
The transformation removes the fixed effects without requiring their unknown values. Different full-rank choices of contrasts can give different numerical normalizing constants, but those constants do not depend on the covariance parameters; they therefore yield the same REML estimates.
Worked Example: Building a Full Set of Error Contrasts
Suppose three individuals have measurements on each of two traits, with a separate fixed mean for each trait and no other fixed effects. Then \(N=6\), \(p=2\), and \(X\) has one indicator column for each trait. An orthonormal pair of contrasts across the three individuals is
For trait-major stacking, set \(C=\operatorname{blockdiag}(Q,Q)\). This is a \(6\times4\) matrix. Its four columns are linearly independent because the columns of \(Q\) are orthonormal. Each column sums to zero within its own trait, so \(C'X=0\). Since \(N-p=4\), the four columns span the full null space of \(X'\): they are a complete set of error contrasts, not just a subset. The transformed data contain two contrasts for each trait, with both trait means removed.
The Standard REML Criterion
When \(V\) is positive definite and \(X\) has full column rank, an equivalent expression for the restricted log likelihood, up to constants independent of the covariance parameters, is
where
The matrix \(P\) removes the part of the data explained by the fixed-effect design in the \(V^{-1}\) metric. In particular, \(PX=0\), so the quadratic form \(y'Py\) does not depend on \(\beta\). The determinant terms adjust for covariance volume and for the fixed-effect information removed from the likelihood.
A short justification connects this expression to error contrasts. Since \(C\) spans the null space of \(X'\), the matrix \(C(C'VC)^{-1}C'\) is the covariance-weighted residual precision matrix; its range lies in that null space, and it equals \(P\) (it is not generally idempotent). Therefore the contrast quadratic form is \(y'Py\). The determinant identity relating the two likelihoods is
where \(c\) depends on the fixed matrices \(X\) and \(C\), but not on \(V\). Substituting these two facts into the contrast log likelihood gives the displayed REML criterion, with only a covariance-parameter-independent constant difference. Thus maximizing either form gives the same covariance-component estimates.
Worked Example: Evaluating a Bivariate Restricted Likelihood
Consider three individuals with a separate fixed mean for each trait, as above. Let the hypothetical relationship matrix be \(K=0.5I_3+0.5J_3\), where \(J_3\) is the \(3\times3\) matrix of ones. Every contrast across individuals has relationship-matrix eigenvalue \(0.5\). Suppose the covariance components are
Both matrices are positive definite: their determinants are \(0.84\) and \(0.24\), respectively. The covariance of the two traits for either orthonormal error-contrast direction is therefore
Suppose the two observed contrast vectors are \(z_1=(0.5,0.2)'\) and \(z_2=(-0.1,0.4)'\). For a vector \((a,b)'\), the quadratic form under \(S\) is \((a^2+b^2-0.6ab)/0.91\). Thus
There are four scalar contrast observations, so their log likelihood is
This is the restricted log likelihood at the stated covariance components and for the specified orthonormal contrasts. It is a likelihood value, not by itself an estimate: estimation requires comparing it with values at other admissible component matrices.
How REML Estimates Covariance Components
Write the distinct entries of \(G\) and \(R\) as a parameter vector \(\theta\). For each candidate \(\theta\), construct \(V(\theta)\), check that it is positive definite, and evaluate the restricted likelihood. REML chooses an admissible value maximizing that likelihood. In computation, this is usually done by numerical optimization, with parameterization or constraints ensuring that \(G\) and \(R\) remain positive semidefinite. The likelihood may have a maximum at a boundary, such as a zero variance component; a component estimate at zero is not automatically a computational error.
The covariance components are coupled. For example, the cross-trait covariance in \(V\) contains both the genetic contribution \(G_{12}K\) and the residual contribution \(R_{12}I_n\). Observed cross-trait covariance is therefore not assigned entirely to genetic covariance. The restricted likelihood assesses candidate components jointly, accounting for the known pattern of relatedness in \(K\) and for the fixed effects removed by the contrasts.
Worked Example: Maximizing a One-Parameter REML Profile
A small hypothetical calculation illustrates how a covariance component can be estimated while the other components are held fixed. Use \(K=0.5I_3+0.5J_3\), trait-specific fixed means, and two orthonormal contrast directions. Set \(R=0\) and consider the genetic covariance matrices
This is a one-parameter profile, not an estimate of every component in a general model. The matrix \(G(\tau)\) is positive semidefinite for \(|\tau|\leq1\). Since \(K\) is positive definite, \(V=G(\tau)\otimes K\) is positive definite for \(|\tau|<1\). Each error-contrast direction has covariance \(0.5G(\tau)=\begin{pmatrix}1&\tau\\\tau&1\end{pmatrix}\). Suppose the observed contrast vectors are \(z_1=(1,1)'\) and \(z_2=(0,1)'\).
Let \(d=1-\tau^2\). The determinant for each contrast vector is \(d\), and the sum of their quadratic forms is
There are two bivariate contrast vectors, hence four scalar observations. Apart from the constant \(-\tfrac12[4\log(2\pi)]\), maximizing the restricted log likelihood is equivalent to minimizing
Differentiation gives
The polynomial \(h(\tau)=2\tau^3-\tau^2+\tau-1\) is strictly increasing because \(h'(\tau)=6\tau^2-2\tau+1>0\) for every \(\tau\); the quadratic has negative discriminant and positive leading coefficient. Its unique root in the admissible interval is approximately \(\hat{\tau}=0.73898\). Since \(f'\) has the sign of \(h\), \(f\) decreases before this root and increases after it. Also, the quadratic term in \(f\) diverges to positive infinity at either boundary for these data. The root therefore gives the profile maximum of the restricted log likelihood. Substitution gives \(f(\hat{\tau})\approx1.7735\), and
At \(\tau=0\), \(f(0)=3\), giving \(\ell_R(0)\approx-5.1758\), which is lower. In this deliberately restricted profile, the estimated genetic cross-trait covariance is \(G_{12}=2\hat{\tau}\approx1.478\) in the model's trait-product units. The calculation illustrates likelihood maximization; it does not imply that the other covariance components can be ignored in a full analysis.
Common Pitfalls and Practical Checks
- Using too few or rank-deficient contrasts. The conditions \(C'X=0\) and the stated dimensions are not enough on their own. \(C\) must have full column rank \(N-p\) and its columns must span the null space of \(X'\); otherwise \(C'VC\) can be singular or the likelihood can omit information.
- Using the inverse-based criterion when \(V\) is singular. The displayed criterion requires \(V\) positive definite and \(X\) full rank. Positive-semidefinite component matrices do not by themselves ensure this condition.
- Interpreting a likelihood value as a component estimate. A single value such as \(-3.814\) only evaluates a candidate. Estimation requires maximizing over the allowed parameter space.
- Assigning all cross-trait covariance to the genetic component. Both \(G_{12}\) and \(R_{12}\) contribute to the marginal covariance. The model estimates them as distinct components using their different covariance structures.
- Ignoring the fixed-effect design. Changing \(X\) changes the space of error contrasts and thus the restricted likelihood. The fixed effects must represent the intended mean structure and \(X\) must have full column rank for the standard formula.
REML matters because covariance-component estimates feed into subsequent descriptions of genetic architecture across traits. Its central operation is precise: remove fixed-effect directions through a complete set of error contrasts, then fit the covariance structure to what remains. Correct interpretation depends on respecting both parts of that construction.
Check Your Understanding
Try these questions using the definitions and likelihood construction above.
- For a model with \(N=40\) observations and a full-rank fixed-effect design with \(p=5\) columns, how many columns must a complete error-contrast matrix \(C\) have?
- Why does \(C'X=0\) make the mean of \(C'y\) independent of \(\beta\), and what additional rank condition is required of \(C\)?
- In the bivariate model, which two covariance components contribute to the marginal cross-trait covariance, and how does the relationship matrix enter the genetic contribution?
- What conditions are needed to use the standard inverse-and-determinant REML criterion involving \(V^{-1}\)?
- In the one-parameter profile example, why is the root of the score equation a maximum of the restricted likelihood rather than a minimum?