Estimating Variance Components by Score and Information

Why Score and Information Matter

A Gaussian pedigree model describes related individuals through a covariance matrix whose entries depend on variance components. Maximum likelihood or restricted maximum likelihood can estimate those components, but the likelihood is generally nonlinear in them. The score measures the local direction of increase, while the expected information scales that direction to form a scoring update. Because variance components cannot be negative, the update must also respect the parameter constraints.

The aim here is to derive the covariance-parameter score and expected information, and then show how they guide estimation. A scoring step is an iterative computational device; it is not by itself a guarantee that the resulting estimate is feasible or that an iteration has increased the likelihood.

Gaussian Covariance Model and Score

Let \(y\) be the phenotype vector for the observed members of a pedigree, with \(m\) entries. Use the Gaussian model \(y\sim N(X\beta,V(\theta))\), where \(X\) is a known design matrix, \(\beta\) is a vector of mean parameters, and \(\theta=(\theta_1,\ldots,\theta_q)^\mathsf{T}\) contains covariance parameters. A common variance-component representation is \(V(\theta)=\sum_{k=1}^q\theta_k K_k\), possibly with a known covariance contribution added. The matrices \(K_k\) are known and symmetric; the parameters are constrained to be nonnegative. At every parameter value used for inference, \(V(\theta)\) must be positive definite on the observed data.

Definition: The covariance derivative for component \(k\) is \(V_k=\partial V(\theta)/\partial\theta_k\). The score component \(U_k\) is the derivative of the log likelihood with respect to \(\theta_k\). The expected information entry \(I_{k\ell}\) is the expected negative second derivative of the log likelihood with respect to \(\theta_k\) and \(\theta_\ell\).

For maximum likelihood, write \(r=y-X\hat{\beta}\), where \(\hat{\beta}\) is the generalized least-squares estimate at the current covariance parameters. The covariance score is

U_k^{\mathrm{ML}} =\frac{1}{2}\left\{ r^\mathsf{T}V^{-1}V_kV^{-1}r -\operatorname{tr}(V^{-1}V_k) \right\}.

The first term measures the covariance pattern in the observed residuals, weighted by \(V^{-1}\). The trace term is the model-based adjustment for the amount of variation expected under the current covariance. With known mean parameters, this expression follows directly by differentiating the Gaussian log likelihood. When \(\beta\) is estimated by maximum likelihood, the same score expression applies, evaluated at \(\hat{\beta}\); the derivative through \(\hat{\beta}\) contributes zero because the likelihood score for \(\beta\) is zero there.

For ML, the expected information for covariance parameters is

I_{k\ell}^{\mathrm{ML}} =\frac{1}{2}\operatorname{tr} \left(V^{-1}V_kV^{-1}V_\ell\right).

This is the covariance-parameter block of the expected information, including when the mean parameters are estimated in the Gaussian model. It is an expectation under the current model parameters, not the observed negative Hessian for a particular data set. The latter is used in ordinary Newton updates and need not equal the expected information.

Deriving the Expected Information

A short derivation explains why the trace expression has this form. For clarity, first take the mean as known and put \(z=V^{-1/2}(y-X\beta)\). Under the model, \(z\sim N(0,I)\). The score can be written as a centered quadratic form:

U_k=\frac{1}{2}\left\{ z^\mathsf{T}A_k z-\operatorname{tr}(A_k) \right\}, \qquad A_k=V^{-1/2}V_kV^{-1/2}.

For a standard normal vector \(z\) and symmetric matrices \(A\) and \(B\), the covariance of the centered quadratic forms is \(\operatorname{Cov}(z^\mathsf{T}Az,z^\mathsf{T}Bz)=2\operatorname{tr}(AB)\). This identity follows by expanding both quadratic forms and using the standard normal fourth moments: products involving an index appearing only once have expectation zero, while paired indices contribute the two cross-pairings in addition to the product of the means. Therefore,

\operatorname{E}(U_kU_\ell) =\frac{1}{4}\,2\operatorname{tr}(A_kA_\ell) =\frac{1}{2}\operatorname{tr}(V^{-1}V_kV^{-1}V_\ell).

In a regular Gaussian model, the covariance of the score equals the expected information. This gives the stated formula. The derivation also shows why expected information depends on the covariance model and its derivatives, not on the observed phenotype values.

Restricted Likelihood Version

When fixed mean effects are present, restricted maximum likelihood (REML) bases covariance estimation on contrasts that remove those effects. Define

P=V^{-1}-V^{-1}X(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}.

Assume \(X\) has full column rank. The REML score and expected information are

U_k^{\mathrm{REML}} =\frac{1}{2}\left\{ y^\mathsf{T}PV_kPy-\operatorname{tr}(PV_k) \right\}, \qquad I_{k\ell}^{\mathrm{REML}} =\frac{1}{2}\operatorname{tr}(PV_kPV_\ell).

These are the corresponding formulas for the restricted likelihood, not the ML formulas with a different label. The matrix \(P\) removes the fitted mean directions and is used in both terms. In practice, one selects the likelihood being optimized, evaluates its matching score and information, and does not mix the two sets of expressions.

From Scores to Constrained Updates

At an interior maximum, each free component has score zero. For a Fisher-scoring iteration, collect the score components in \(U\) and the expected information entries in \(I\). The unconstrained update is

\theta^{\mathrm{new}}=\theta^{\mathrm{old}}+I^{-1}U.

The update uses the local score direction, adjusted for the curvature represented by expected information. A Newton update instead uses the observed Hessian, when available. Either method can propose a negative variance component or a covariance matrix that is not positive definite. A practical constrained algorithm therefore uses an active set or another constrained optimization method, and typically checks a step length and the likelihood before accepting an update. Simply replacing negative proposals by zero changes the proposed direction and does not, by itself, establish that the likelihood has improved.

At a boundary maximum, the one-sided condition for a component fixed at zero is that its score in the increasing direction be nonpositive, holding the other components at their constrained optimum. A zero score at a boundary is also possible. For components strictly above zero, the score must be zero at a regular maximum. These local conditions need to be assessed together with feasibility; the sign of the score at zero alone does not establish that zero is the global constrained maximum.

Worked Example: A One-Component ML Estimate

Suppose two hypothetical independent phenotypes have known mean zero and variance \(1+\theta\), where the known residual variance is 1 and the variance component satisfies \(\theta\geq0\). The hypothetical observations are \(0.5\) and \(0.5\). Thus \(n=2\), \(q=\sum_i y_i^2=0.5\), and \(s=1+\theta\geq1\). The model is Gaussian and the variance is positive throughout the feasible range.

Ignoring constants, the log likelihood is \(\ell(\theta)=-\frac{n}{2}\log(s)-\frac{q}{2s}\). Differentiating with respect to \(\theta\) (so \(ds/d\theta=1\)) gives

U(\theta)=\frac{1}{2}\left(\frac{q}{s^2}-\frac{n}{s}\right) =\frac{0.25}{s^2}-\frac{1}{s}.

For every feasible \(s\geq1\), \(0.25/s^2<1/s\), since multiplying by \(s^2>0\) reduces this inequality to \(0.25<s\). Hence \(U(\theta)<0\) throughout the feasible range. The log likelihood is strictly decreasing as \(\theta\) increases, so its constrained maximum is at \(\hat{\theta}=0\). This conclusion follows from the score over the full feasible range, not merely from its value at the boundary.

For reference, the expected information at a general feasible value is \(I(\theta)=n/(2s^2)=1/s^2\). At \(\theta=0\), it is 1, but the score there is \(-0.75\), pointing toward negative \(\theta\), outside the allowed parameter space. The constrained estimate therefore remains at zero.

Worked Example: A Two-Component Pedigree Covariance Update

Consider a hypothetical pair of relatives with known mean zero and covariance \(V=\theta_AK+\theta_EI\), where \(K=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\), \(\theta_A\geq0\), and \(\theta_E>0\). Here \(K\) is a specified relationship-pattern matrix for this example, while \(I\) is the identity matrix. The eigenvalues of \(K\) are 1.5 for the common direction and 0.5 for the contrast direction. The data are \(y=(1,-1)^\mathsf{T}\), so the transformed observations in those directions are \(z_1=0\) and \(z_2=\sqrt{2}\). Evaluate the ML score and expected information at \((\theta_A,\theta_E)=(0.20,1.00)\). The eigenvariances are \(\lambda_1=1.5\theta_A+\theta_E=1.30\) and \(\lambda_2=0.5\theta_A+\theta_E=1.10\), both positive.

In the eigenbasis, for a component with derivative \(d_{jk}=\partial\lambda_j/\partial\theta_k\), the score is \(U_k=\frac{1}{2}\sum_j d_{jk}(z_j^2/\lambda_j^2-1/\lambda_j)\), and the expected information is \(I_{k\ell}=\frac{1}{2}\sum_j d_{jk}d_{j\ell}/\lambda_j^2\). Here \(d_{jA}=(1.5,0.5)\) and \(d_{jE}=(1,1)\). Substitution gives

\begin{aligned} U_A&=\frac{1}{2}\left[ 1.5\left(0-\frac{1}{1.30}\right) +0.5\left(\frac{2}{1.10^2}-\frac{1}{1.10}\right) \right]\approx-0.3910,\\ U_E&=\frac{1}{2}\left[ -\frac{1}{1.30} +\frac{2}{1.10^2}-\frac{1}{1.10} \right]\approx-0.01271,\\ I&\approx \begin{pmatrix} 0.7690&0.6504\\ 0.6504&0.7091 \end{pmatrix}. \end{aligned}

Solving \(I\Delta=U\) gives approximately \(\Delta=(-2.20,2.00)^\mathsf{T}\). The full scoring proposal would put \(\theta_A\) below zero, so it is infeasible. A trial step with length \(0.05\) gives \((0.20,1.00)+0.05\Delta\approx(0.090,1.10)\), which is feasible and retains positive eigenvariances. This is a candidate for likelihood checking, not an automatically accepted update: a constrained algorithm must evaluate its chosen step and reduce or alter it if the objective does not improve. The example illustrates why score and information supply an update direction but do not remove the need to enforce the parameter constraints.

Worked Example: A REML Score for a Family Contrast

Suppose a hypothetical pair has an unknown fixed common mean and covariance \(V=\tau K+I\), with \(\tau\geq0\) and the same \(K\) as above. The residual variance is known to be 1. REML removes the common-mean direction, leaving the contrast \(z=(y_1-y_2)/\sqrt{2}\). Take hypothetical observations \(y=(1,-1)^\mathsf{T}\), giving \(z=\sqrt{2}\). In this contrast direction, the variance is \(\lambda=1+0.5\tau\), which is positive for every feasible \(\tau\).

The restricted log likelihood for this single contrast, up to a constant, is \(-\frac{1}{2}\log(\lambda)-\frac{z^2}{2\lambda}\). Since \(d\lambda/d\tau=0.5\), the REML score and expected information are

U_\tau^{\mathrm{REML}} =\frac{1}{2}(0.5)\left(\frac{z^2}{\lambda^2}-\frac{1}{\lambda}\right), \qquad I_{\tau\tau}^{\mathrm{REML}} =\frac{1}{2}\frac{0.5^2}{\lambda^2}.

At \(\tau=0\), \(\lambda=1\) and \(z^2=2\), so \(U_\tau^{\mathrm{REML}}=0.25\) and \(I_{\tau\tau}^{\mathrm{REML}}=0.125\). The Fisher-scoring proposal is therefore \(\tau^{\mathrm{new}}=0+0.25/0.125=2.00\), a feasible value. In this one-parameter example the restricted likelihood is maximized when \(\lambda=z^2=2\), which indeed corresponds to \(\hat{\tau}=2\). The equality is specific to this simple contrast model; in general, scoring iterations need not reach the maximum in one step.

Practical Checks and Common Pitfalls

Score-based estimation is useful because it turns a complicated covariance likelihood into derivatives that can be evaluated using matrix operations. In a pedigree analysis, the covariance derivatives encode which entries change with each component, while the information indicates how distinguishable those changes are under the current model. If two covariance components have very similar derivative matrices, the information matrix can be poorly conditioned, making updates unstable. This is a numerical warning about the local estimation problem, not evidence by itself that either component is zero.

Using an ML score with REML information. These formulas differ through \(V^{-1}\) versus \(P\). Use the score and information belonging to the likelihood actually being optimized.
Treating a scoring step as the estimate. It is only a proposal. Check nonnegativity, positive definiteness, and the objective value before accepting an iteration.
Inferring a boundary maximum from one score value. A nonpositive score at zero indicates a local direction, not necessarily monotonic decrease over the entire feasible range. To establish a global maximum in a one-dimensional example, inspect the score throughout the range or compare the likelihood at all relevant candidates.
Forgetting the scale of a component. A variance component is expressed in squared phenotype units. Report its estimate in those units and distinguish it from a standard deviation, which is on the phenotype scale.

Key takeaway: The covariance score compares observed residual variation with the model's expected variation, and expected information scales the resulting update. The nonnegative constraints remain part of the estimation problem: feasibility and likelihood improvement must be checked at every proposed step.

Check Your Understanding

Work through these questions using the formulas and constraints described above.

In the ML score, what role is played by the quadratic term, and what role is played by the trace term?
Why does the expected information for two covariance components depend on both covariance derivative matrices?
A Fisher-scoring proposal makes one variance component negative. Why is replacing that value by zero not, by itself, a sufficient update rule?
For a one-component model, what score condition is required at a boundary maximum with the component fixed at zero?
Why must a REML calculation use \(P\) rather than simply substituting fitted residuals into every ML expression?