Maximum Likelihood and REML for Pedigree Models

Why Two Likelihoods?

Phenotypes measured on relatives are correlated, so their joint likelihood must account for the pedigree covariance rather than treating every person as independent. At the same time, the mean phenotype may depend on fixed effects such as an intercept or recorded covariates. Estimating those mean effects uses information in the data. Maximum likelihood (ML) and restricted maximum likelihood (REML) handle that estimation differently, and the difference can matter for estimated variance components, especially in small samples.

The distinction is not that ML ignores fixed effects or that REML assumes they are absent. Both methods allow fixed effects. ML fits the fixed effects and covariance parameters jointly, or profiles the fixed effects out of the likelihood. REML instead constructs a likelihood from linear combinations of the data that do not contain the fixed effects. This accounts for the information used to estimate the mean structure.

Gaussian Pedigree Model and ML Objective

Definition: Let $y$ be the $n$-vector of observed phenotypes, $X$ an $n \times p$ full-column-rank fixed-effect design matrix, and $\beta$ a $p$-vector of fixed-effect coefficients. A Gaussian pedigree model is $$ y \sim N(X\beta,V(\theta)), $$ where $V(\theta)$ is a positive-definite covariance matrix depending on variance-component parameters $\theta$. The pedigree structure contributes to the covariance; for example, a component may be proportional to a known relationship matrix. Variance components are constrained to values that make the covariance model valid.

This likelihood assumes that the observed phenotypes, conditional on the specified mean and covariance, have a joint multivariate normal distribution. The pedigree relationships and any other matrices used to construct $V(\theta)$ are treated as known. The model also assumes that the mean is $X\beta$ and that the covariance has the stated form. Ascertainment or selection that changes the phenotype distribution is not accounted for unless it is explicitly included in the likelihood.

For a fixed value of $\theta$, the log likelihood is

\ell(\beta,\theta) = -\frac{1}{2}\left[ n\log(2\pi)+\log|V(\theta)| +(y-X\beta)^\mathsf{T}V(\theta)^{-1}(y-X\beta) \right].

The determinant term reflects the overall volume of the covariance distribution; the quadratic term measures the residuals relative to that covariance. In particular, a residual in a direction with high modeled variance contributes less to the quadratic than an equally sized residual in a low-variance direction.

For each fixed $\theta$, maximizing over $\beta$ gives the generalized least-squares estimate

\hat{\beta}(\theta) = \left(X^\mathsf{T}V(\theta)^{-1}X\right)^{-1} X^\mathsf{T}V(\theta)^{-1}y.

Substituting this estimate into the likelihood defines the profiled ML objective. The ML estimate $\hat{\theta}_{\mathrm{ML}}$ maximizes that objective over the allowed covariance parameters. The fixed effects are estimated as part of this process, but the profiled log likelihood does not include an extra determinant term for estimating them.

Worked Example: Fitting a Fixed Mean With a Sibling Covariance

Suppose two hypothetical siblings have phenotypes $y=(2,4)^\mathsf{T}$. Consider the model with an intercept $\beta$, known variance component $\tau^2=4$, and covariance $V=4R$, where

R= \begin{pmatrix} 1 & 0.5\\ 0.5 & 1 \end{pmatrix}, \qquad X= \begin{pmatrix} 1\\ 1 \end{pmatrix}.

The matrix $R$ has diagonal entries 1 and sibling covariance 0.5; it is positive definite because its eigenvalues are 1.5 and 0.5. Thus $V$ is also positive definite. This example treats the covariance as known and estimates only the intercept.

Since $R^{-1}=\frac{1}{0.75}\begin{pmatrix}1&-0.5\\-0.5&1\end{pmatrix}$, the generalized least-squares estimate is

\hat{\beta} = \frac{X^\mathsf{T}V^{-1}y}{X^\mathsf{T}V^{-1}X} = \frac{3}{1} =3.

The residual vector is $y-X\hat{\beta}=(-1,1)^\mathsf{T}$. The contrast direction $(-1,1)^\mathsf{T}$ has eigenvalue 0.5 under $R$, so its variance under $V=4R$ is $4(0.5)=2$. Therefore the residual quadratic is $2/2=1$. Also, $|V|=16|R|=16(0.75)=12$. The maximized log likelihood for this fixed covariance is

\ell(\hat{\beta},\theta) = -\frac{1}{2}\left[2\log(2\pi)+\log(12)+1\right] \approx -3.580.

This is the likelihood at the fitted intercept and specified covariance, not a test of whether the siblings have different means. Their dependence is represented in the joint covariance matrix.

REML as a Likelihood of Error Contrasts

An error contrast is a linear combination of the observations that removes the fixed effects. Choose an $n \times (n-p)$ matrix $L$ of full column rank such that $L^\mathsf{T}X=0$. Then $L^\mathsf{T}y$ has mean zero for every value of $\beta$:

E(L^\mathsf{T}y)=L^\mathsf{T}X\beta=0, \qquad \operatorname{Var}(L^\mathsf{T}y)=L^\mathsf{T}V(\theta)L.

Because a linear transformation of a multivariate Gaussian vector is Gaussian, $L^\mathsf{T}y$ has a Gaussian likelihood involving $\theta$ but not $\beta$. Maximizing this likelihood gives the REML estimate. A change of basis among the contrasts changes the log likelihood by a constant that does not depend on $\theta$, so it does not change the REML estimate.

REML objective: Up to an additive constant independent of $\theta$, the restricted log likelihood is $$ \ell_R(\theta) = -\frac{1}{2}\left[ \log|V(\theta)| +\log|X^\mathsf{T}V(\theta)^{-1}X| +y^\mathsf{T}P(\theta)y +(n-p)\log(2\pi) \right], $$ where $$ P(\theta) = V(\theta)^{-1} - V(\theta)^{-1}X \left(X^\mathsf{T}V(\theta)^{-1}X\right)^{-1} X^\mathsf{T}V(\theta)^{-1}. $$

The matrix $P(\theta)$ removes the fitted fixed-effect directions in the covariance-weighted geometry. In particular, $y^\mathsf{T}P(\theta)y$ is the residual quadratic evaluated after fitting the fixed effects. REML maximizes the restricted log likelihood over the allowed covariance parameters; it does not estimate variance components by simply applying an ordinary likelihood to the raw residuals.

Why the REML Determinant Appears

The error-contrast likelihood has covariance $L^\mathsf{T}VL$, so its log likelihood contains $\log|L^\mathsf{T}VL|$. The standard REML expression uses $\log|V|+\log|X^\mathsf{T}V^{-1}X|$ instead. These expressions differ only by a factor independent of $\theta$. The following determinant identity explains the substitution.

Separate fitted and contrast directions.
Choose an orthonormal basis for the columns of $X$ and an orthonormal basis for their orthogonal complement. In this combined basis, the covariance matrix has blocks for the fitted and contrast directions.

Take a block determinant.
The determinant of the full covariance equals the determinant of the contrast covariance times the determinant of the fitted-direction Schur complement. The latter is related to $|X^\mathsf{T}V^{-1}X|^{-1}$.

Discard factors fixed by the design.
For a fixed $X$ and fixed contrast basis, any remaining basis or design factor is constant in $\theta$. Thus the contrast likelihood and the displayed REML objective have the same maximizing covariance parameters.

This derivation also shows why the term $\log|X^\mathsf{T}V^{-1}X|$ belongs to REML. It records the covariance-dependent contribution of estimating the fixed-effect directions when the likelihood is expressed through contrasts. The ML profile likelihood has no such term: ML evaluates the original $n$-dimensional data likelihood at the fitted $\beta$.

Worked Example: REML Removes an Intercept in a Small Gaussian Sample

As a transparent special case, suppose three hypothetical independent observations are $y=(1,2,5)^\mathsf{T}$, with model $y_i=\beta+\varepsilon_i$, where $\varepsilon_i$ are independent $N(0,\sigma^2)$. Here $X$ is an intercept column, $p=1$, and the covariance is $V=\sigma^2 I$, which is positive definite for $\sigma^2>0$.

The fitted mean is $\bar{x}=8/3$. The residual sum of squares is

Q = \sum_{i=1}^3(y_i-\bar{x})^2 = \left(-\frac{5}{3}\right)^2 +\left(-\frac{2}{3}\right)^2 +\left(\frac{7}{3}\right)^2 = \frac{78}{9} = \frac{26}{3}.

The profiled ML log likelihood, apart from constants, is $-\frac{1}{2}[3\log(\sigma^2)+Q/\sigma^2]$. Differentiating with respect to $\sigma^2$, or maximizing this one-parameter expression, gives $\hat{\sigma}^2_{\mathrm{ML}}=Q/3=26/9\approx2.889$.

For REML, two independent contrasts can be formed because $n-p=2$. Their total squared variation is $Q$; the contrast likelihood has two Gaussian dimensions with variance $\sigma^2$. Its scale estimate is therefore $\hat{\sigma}^2_{\mathrm{REML}}=Q/2=13/3\approx4.333$. Both calculations use the same residual quadratic, but ML divides by the number of observations while REML accounts for the one fitted mean by using two error-contrast dimensions. The result illustrates a difference between the estimates; it does not imply that REML is always larger in every model or sample.

Comparing Variance-Component Estimates in a Pedigree Model

A simple pedigree covariance model makes the ML–REML difference concrete. Suppose two hypothetical sibling pairs come from separate families, with no modeled covariance between families. Let the covariance be $V=\tau^2R$, where $R$ has two sibling blocks with diagonal 1 and off-diagonal 0.5. This is a one-component illustration: the unknown $\tau^2>0$ scales a known positive-definite pedigree covariance matrix. The fixed-effect design contains an intercept, so $p=1$.

Worked Example: ML and REML Estimates for Two Sibling Pairs

Let the four hypothetical phenotypes be $y=(1,3,2,4)^\mathsf{T}$, ordered by sibling pair. The covariance matrix $R$ consists of two copies of $\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}$ on its diagonal and zeros between pairs. Each block has eigenvalues 1.5 and 0.5, so $R$ is positive definite. The intercept estimate is the GLS mean, $\hat{\beta}=2.5$, since the two members of each pair have equal marginal weights.

The residuals are $(-1.5,0.5,-0.5,1.5)^\mathsf{T}$. For each pair, the normalized within-pair contrast has residual projection with squared value 2 and covariance eigenvalue 0.5. The two within-pair contributions to $y^\mathsf{T}P y$ are therefore $2/0.5=4$ each. The normalized contrast between pair means has projection $(1+3-2-4)/2=-1$, squared value 1, and eigenvalue 1.5, contributing $1/1.5=2/3$. Thus

Q=y^\mathsf{T}P y = 4+4+\frac{2}{3} = \frac{26}{3}.

For $V=\tau^2R$, the ML log likelihood after fitting the intercept has covariance-scale terms $-\frac{1}{2}[4\log(\tau^2)+Q/\tau^2]$; $\log|R|$ is constant in $\tau^2$. Maximizing gives

\hat{\tau}^2_{\mathrm{ML}} = \frac{Q}{n} = \frac{26/3}{4} = \frac{13}{6} \approx 2.167.

The restricted likelihood uses the three error-contrast dimensions, so its scale-dependent terms are $-\frac{1}{2}[3\log(\tau^2)+Q/\tau^2]$. Hence

\hat{\tau}^2_{\mathrm{REML}} = \frac{Q}{n-p} = \frac{26/3}{3} = \frac{26}{9} \approx 2.889.

Both estimates are in squared phenotype units. Under this specified one-component model, REML estimates a larger covariance scale because it uses three contrast dimensions after fitting the intercept, rather than four dimensions in the ML likelihood. The numerical difference is specific to these data and this model, not a universal ordering rule.

Choosing and Interpreting the Objective

ML and REML answer related but distinct likelihood questions. ML uses the full data likelihood and is appropriate when the likelihoods for competing models need to be compared, including models with different fixed-effect designs. REML is often used for estimating covariance parameters when the fixed-effect design is held constant, because its contrast likelihood accounts for the estimation of those effects. REML likelihood values for models with different fixed-effect designs are generally not directly comparable: the contrasts, and therefore the restricted likelihood, depend on $X$.

Neither method removes the need to check the model. The covariance must be positive definite for the observed data, the mean structure must be specified, and the Gaussian assumption must be reasonable for the intended likelihood analysis. Estimates can also lie at or near allowed boundaries, or covariance components may be difficult to distinguish given the pedigree structures and sample size. ML and REML are likelihood procedures, not guarantees that a component is precisely estimated or that the covariance model is biologically correct.

Common Mistake

Wrong reasoning: “REML is just ML applied to the fitted residuals, so I can compare its likelihood across any two models.”
Correction: REML is the likelihood of error contrasts that eliminate the fixed effects. Its objective includes a covariance-dependent fixed-effect adjustment and depends on the design matrix $X$. Compare REML likelihoods only when the fixed-effect design is the same; use ML when comparing models with different fixed-effect structures.
Wrong reasoning: “REML always gives the larger variance-component estimate.”
Correction: The estimates depend on the covariance structure, data, and fitted mean. REML accounts for the fixed-effect estimation, but there is no universal ordering of the ML and REML estimates across all models and samples.

Key takeaway: ML maximizes the joint Gaussian likelihood after fitting the fixed effects. REML maximizes a likelihood built from contrasts that remove those effects. For a fixed mean design, the REML determinant adjustment and the $n-p$ contrast dimensions can change variance-component estimates; they do not make the two likelihoods interchangeable.

Check Your Understanding

Use the model definitions and objectives above to answer the following questions.

In a Gaussian pedigree model, what roles do $X\beta$ and $V(\theta)$ play?
Why does an error contrast $L^\mathsf{T}y$ have a distribution that does not depend on $\beta$?
In an intercept-only model with $n=6$ observations, how many independent error-contrast dimensions are used by REML?
Why are REML likelihood values generally unsuitable for comparing models with different fixed-effect design matrices?
In the two-sibling-pair example, why are the ML and REML variance-component estimates different even though both use the same residual quadratic?