Profile Likelihood for Variance Components

Why Profile a Variance Component?

A variance component describes one source of variation in a phenotype, such as variation associated with a pedigree covariance structure. Its estimate is rarely determined in isolation: the likelihood also depends on nuisance parameters, including the mean and other variance components. An interval for the component should reflect the best fit available at each candidate value, not an arbitrary choice of those other parameters.

Profile likelihood does this by fixing the parameter of interest in turn and maximizing the likelihood over the nuisance parameters. The resulting likelihood curve can be asymmetric, particularly when the variance component is constrained to be nonnegative and its estimate is near zero. A symmetric estimate-plus-or-minus-standard-error summary can then be misleading.

Definition: Let $\tau$ be the variance component of interest and let $\eta$ denote all nuisance parameters. If $\ell(\tau,\eta)$ is the log likelihood, the profile log likelihood is $$ \ell_p(\tau)=\max_{\eta\in\mathcal{H}(\tau)}\ell(\tau,\eta), $$ where $\mathcal{H}(\tau)$ is the allowed set of nuisance-parameter values when $\tau$ is fixed. The profile likelihood retains the best-fitting nuisance parameters separately for every candidate $\tau$.

A Gaussian Pedigree Model

For a simple quantitative-trait model, let $y$ be the vector of phenotypes for $n$ relatives. Suppose the mean is $X\beta$, and the covariance is $$ V(\tau,\sigma_e^2)=\tau K+\sigma_e^2 I, $$ where $K$ is a known pedigree covariance matrix, $\tau\geq 0$ is the variance component of interest, and $\sigma_e^2\geq 0$ is a residual variance component. The matrix $I$ is the identity matrix. Conditional on the model parameters, assume $y\sim N(X\beta,V)$, with $V$ positive definite at the parameter values used in the likelihood.

The Gaussian log likelihood is $$ \ell(\beta,\tau,\sigma_e^2)=-\frac{1}{2}\left[n\log(2\pi)+\log|V|+(y-X\beta)^{T}V^{-1}(y-X\beta)\right]. $$

For a fixed candidate $\tau$, profiling means maximizing this expression over $\beta$ and $\sigma_e^2$, subject to the model constraints. The resulting value $\ell_p(\tau)$ is then compared with its maximum over all allowed $\tau$. This tutorial uses maximum likelihood; if an analysis is based on a restricted likelihood, the restricted log likelihood must be profiled instead. The chosen likelihood and its assumptions should be kept consistent throughout the calculation.

Key distinction: A profile value at $\tau_0$ is the likelihood after nuisance parameters have been reoptimized with $\tau=\tau_0$. It is not generally the likelihood obtained by fixing nuisance parameters at their estimates from the unrestricted fit.

Worked Example: Setting Up a Small Pedigree Likelihood

Consider four hypothetical relatives with phenotypes $y=(0,1,0,1)^T$, measured in units $u$. Use an intercept-only mean, so $X$ is a column of four ones and $\beta$ has units $u$. For illustration, take the known dimensionless covariance matrix $$ K=\begin{pmatrix} 1&0.5&0.5&0.5\\ 0.5&1&0.5&0.5\\ 0.5&0.5&1&0.5\\ 0.5&0.5&0.5&1 \end{pmatrix}. $$ The model is $V=\tau K+\sigma_e^2 I$; both variance parameters have units $u^2$.

The sample mean is $\bar{x}=0.50\,u$, where here the four entries are the observed phenotypes. The residual vector about that mean is $(-0.5,0.5,-0.5,0.5)^T u$. It sums to zero, has squared length $1\,u^2$, and lies in the three-dimensional subspace orthogonal to the constant vector. For this $K$, the eigenvalue in the constant-vector direction is $2.5$, and each of the three orthogonal directions has eigenvalue $0.5$.

Consequently, at any candidate $\tau$ and residual variance $s=\sigma_e^2$, the covariance eigenvalues are $s+2.5\tau$ once and $s+0.5\tau$ three times. The intercept estimate is $\hat{\beta}(\tau,s)=0.50\,u$, and the residual quadratic form is $1/(s+0.5\tau)$. Thus the log likelihood after fitting the intercept is $$ \ell(\hat{\beta},\tau,s)=-\frac{1}{2}\left[4\log(2\pi)+3\log(s+0.5\tau)+\log(s+2.5\tau)+\frac{1}{s+0.5\tau}\right]. $$ This expression will let us check each profile calculation directly.

Optimizing the Nuisance Variance

At a fixed $\tau$, the intercept and residual variance must be fitted for that candidate. In the example, define $a=s+0.5\tau$ and $b=s+2.5\tau$. Differentiating the fitted log likelihood with respect to $s$ gives $$ \frac{\partial\ell}{\partial s}=-\frac{1}{2}\left[\frac{3}{a}+\frac{1}{b}-\frac{1}{a^2}\right]. $$

An interior optimum for $s$ makes this derivative zero. If the unconstrained solution would have $s<0$, the constrained maximum must instead be checked at the boundary $s=0$. This is an example of why the allowed parameter space is part of the profile calculation, not an afterthought.

Worked Example: Profiling at a Fixed Candidate Value

For the hypothetical data above, fix $\tau=0.10\,u^2$. Then $a=s+0.05$ and $b=s+0.25=a+0.20$, with all these quantities in $u^2$. The interior score equation for $s$ is $$ \frac{3}{a}+\frac{1}{a+0.20}-\frac{1}{a^2}=0. $$ Multiplying by $a^2(a+0.20)$ gives $4a^2-0.4a-0.2=0$. Its positive root is $a=(1+\sqrt{21})/20=0.2791\,u^2$, so the nuisance-variance estimate is $\hat{s}=a-0.05=0.2291\,u^2$, which satisfies the constraint. The fitted intercept remains $0.50\,u$.

Substitute those optimized values into the log likelihood: $$ \ell_p(0.10)=-\frac{1}{2}\left[4\log(2\pi)+3\log(0.2791)+\log(0.4791)+\frac{1}{0.2791}\right]\approx -3.185. $$ The profile value is therefore approximately $-3.185$, not the value obtained by inserting an arbitrary residual variance.

For comparison, if $s$ is held at $0.40\,u^2$ while the intercept is fitted, the eigenvalues are $0.45\,u^2$ three times and $0.65\,u^2$ once. The log likelihood is $$ -\frac{1}{2}\left[4\log(2\pi)+3\log(0.45)+\log(0.65)+\frac{1}{0.45}\right]\approx -3.374. $$ This is a valid likelihood at that parameter choice, but it is lower than the profile value because $s=0.40\,u^2$ is not the maximizing nuisance variance at $\tau=0.10\,u^2$.

From a Profile Curve to an Interval

Let $\hat{\tau}$ maximize $\ell_p(\tau)$ over the allowed values. A commonly used likelihood-ratio interval retains candidate values whose profile log likelihood is not too far below its maximum. With a one-degree-of-freedom chi-square reference, the nominal 95% set is $$ \left\{\tau\geq 0: 2\left[\ell_p(\hat{\tau})-\ell_p(\tau)\right]\leq 3.841\right\}. $$ The threshold $3.841$ is the 95th percentile of a chi-square distribution with one degree of freedom. This familiar calibration is an approximation that relies on regular likelihood theory. When a variance component is on its boundary, regularity can fail, so the resulting interval should not be presented as having guaranteed 95% coverage without suitable calibration.

Unlike a Wald interval, which is often centered on an estimate and scaled by a standard error, a profile interval follows the likelihood itself. The likelihood need not fall at the same rate on both sides of its maximum. Also, because $\tau\geq0$, a profile interval cannot extend into negative values. If the maximum is at zero, an interval may have a zero lower endpoint and a positive upper endpoint rather than being symmetric around the estimate.

Worked Example: An Asymmetric Interval Near Zero

Continue with the four hypothetical phenotypes. Maximizing the profile likelihood over $\tau\geq0$ gives $\hat{\tau}=0$, $\hat{s}=0.25\,u^2$, and $\ell_p(0)\approx-2.903$. The boundary constraint is respected. At this fit, the covariance is $0.25I$, and the fitted residual quadratic form is $1/0.25=4$.

Now consider the candidate $\tau=1.80\,u^2$. The nuisance-variance maximum is at $s=0$. To check the constraint, at $s=0$ we have $a=0.90$ and $b=4.50$, so the derivative with respect to $s$ is $$ -\frac{1}{2}\left[\frac{3}{0.90}+\frac{1}{4.50}-\frac{1}{0.90^2}\right]\approx -1.160<0. $$ The likelihood decreases as $s$ increases from zero, confirming the constrained maximum at the boundary. Substitution gives $$ \ell_p(1.80)=-\frac{1}{2}\left[4\log(2\pi)+3\log(0.90)+\log(4.50)+\frac{1}{0.90}\right]\approx-4.825. $$

The likelihood-ratio statistic comparing this candidate with the maximum is $$ 2\left[\ell_p(0)-\ell_p(1.80)\right]=2[-2.903-(-4.825)]\approx3.844. $$ This is approximately the $3.841$ cutoff, so the conventional chi-square-one profile interval has an upper endpoint near $1.80\,u^2$. For this example it is approximately $[0,1.80]\,u^2$. It is strongly asymmetric: the estimate is at zero and the interval extends upward. Because the maximum is on the boundary, regard this as an approximate likelihood interval, not an automatically exact 95% confidence statement.

Why Profiling Matters in Pedigree Analysis

The likelihood for a pedigree phenotype vector depends jointly on its mean structure and covariance. A candidate variance component changes the covariance matrix, and the best-fitting residual variance may change with it. Reoptimizing nuisance parameters makes the comparison fair: each candidate is evaluated using its own best fit under the stated model.

Profile intervals are also useful because they display information that a single standard error can hide. A curve that falls gradually in one direction and sharply in another gives different amounts of support across the two sides of an estimate. Near zero, the nonnegative constraint can truncate one side entirely. Reporting the estimate, the interval, the units, and the boundary status helps convey that shape without implying that uncertainty is symmetric.

Specify the likelihood and constraints.
State the phenotype distribution, mean, covariance, and allowed parameter ranges before computing a profile.

Choose a candidate value.
Fix the variance component of interest at $\tau_0$, keeping its units and constraint explicit.

Refit nuisance parameters.
Maximize the likelihood over all nuisance parameters permitted at $\tau_0$, checking boundary values as needed.

Repeat and compare.
Evaluate the profile over candidate values, find its maximum, and use a stated likelihood-ratio reference to form an interval.

Common Mistake

A frequent error is to fix the nuisance variance at its unrestricted estimate, or at a convenient round value, and call the resulting likelihood a profile likelihood. That is not profiling: at each candidate $\tau$, nuisance parameters must be reoptimized subject to the model constraints. Another error is to report a symmetric estimate-plus-or-minus-standard-error interval when the estimate is at zero. The correction is to inspect and use the profile curve, respect the nonnegative parameter space, and state the calibration limitation when the estimate is on a boundary.

Key takeaway: A profile likelihood interval for a pedigree variance component is built by fixing that component and maximizing over nuisance parameters at every candidate value. This produces an interval shaped by the likelihood and constraints, often asymmetric near zero.

Check Your Understanding

Answer each question using the model and definitions in this tutorial.

In the Gaussian pedigree model, which parameters are nuisance parameters when $\tau$ is the target?
Why does fixing $s$ at $0.40\,u^2$ when $\tau=0.10\,u^2$ not produce the profile likelihood?
What constraint must be checked when optimizing a variance component and residual variance?
For the worked data, why is the profile interval asymmetric around the estimate?
Why should the usual chi-square-one calibration be described cautiously when the estimate is at zero?