Uncertainty in Variance-Component Estimates

Why Uncertainty in Variance Components Needs Care

A mixed model can estimate how much variation is associated with genetic similarity and how much remains residual. These estimates describe features of a covariance structure, not directly observed individual effects. Their uncertainty matters: two analyses can produce the same estimated genetic variance but very different precision, and a point estimate alone does not show whether the data tightly constrain the component.

Variance components also differ from ordinary unconstrained regression coefficients. A variance must be nonnegative. Near zero, a symmetric approximation can extend into impossible values, and the estimate may have a sampling distribution that is skewed or concentrated at the boundary. The standard error is still useful, but its interpretation depends on how it was calculated and how far the estimate is from the boundary.

Model and Information for a Variance Component

Consider a Gaussian genetic mixed model with a known genomic relationship matrix $K$:

y=X\beta+u+e,\qquad u\sim N(0,\sigma_g^2K),\qquad e\sim N(0,\sigma_e^2I).

Here $y$ is the vector of outcomes, $X$ is the fixed-effect design matrix, $\beta$ is the fixed-effect coefficient vector, $u$ is the polygenic random effect, and $e$ is residual error. Assume $u$ and $e$ are independent, $K$ is known and positive semidefinite, and the model’s Gaussian covariance structure is appropriate. The marginal covariance matrix is $V=\sigma_g^2K+\sigma_e^2I$. The parameters $\sigma_g^2$ and $\sigma_e^2$ are constrained to be nonnegative.

REML estimates these covariance parameters using likelihood information from contrasts that remove the fixed effects. For uncertainty calculations, collect the parameters in $\theta=(\sigma_g^2,\sigma_e^2)$. The likelihood being differentiated should be the same likelihood used to obtain the estimates—for example, the REML log likelihood if the estimates are REML estimates.

Definition: The observed information matrix at an interior estimate $\hat{\theta}$ is the negative matrix of second derivatives of the log likelihood: $$ J(\hat{\theta})=-\left.\frac{\partial^2\ell(\theta)} {\partial\theta\,\partial\theta^\mathsf{T}}\right|_{\theta=\hat{\theta}}. $$ When the likelihood is locally well approximated by a quadratic and the estimate is away from parameter boundaries, the estimated covariance matrix of $\hat{\theta}$ is $J(\hat{\theta})^{-1}$. The standard error for a component is the square root of its corresponding diagonal entry.

The inverse matrix matters when there is more than one variance component. It accounts for the fact that the estimates can move together: uncertainty in one component may be related to uncertainty in another. Taking the reciprocal of a single diagonal entry of $J$ would generally ignore that dependence.

A Local Approximation and Wald Intervals

Near an interior maximum, a second-order Taylor approximation describes how the log likelihood changes around its maximum. With $\hat{\theta}$ denoting the estimate and $\theta$ a nearby parameter value, the approximation is

\ell(\theta)\approx \ell(\hat{\theta}) -\frac{1}{2}(\theta-\hat{\theta})^\mathsf{T} J(\hat{\theta})(\theta-\hat{\theta}).

This quadratic shape is the reason the inverse information estimates covariance: a sharply curved likelihood corresponds to more precise estimates, while a flatter likelihood corresponds to less precise estimates. Under the usual large-sample, interior-parameter approximation, an estimate has an approximately normal sampling distribution with covariance $J^{-1}$. A two-sided Wald interval for one component is therefore the estimate plus or minus a normal critical value times its standard error.

\hat{\theta}_j\ \pm\ z_{1-\alpha/2}\,SE(\hat{\theta}_j), \qquad SE(\hat{\theta}_j)=\sqrt{\left[J(\hat{\theta})^{-1}\right]_{jj}}.

For a nominal 95% interval, $z_{0.975}$ is approximately 1.96. The interval is an approximation whose repeated-sampling coverage relies on the local normal approximation being adequate. A confidence interval does not assign a probability to the fixed, unknown parameter after the data have been observed. Rather, under the stated model and procedure, intervals constructed this way would cover the parameter at the nominal rate over repeated samples when the approximation is accurate.

Worked Example: Standard Errors From a Two-Component Information Matrix

Suppose a hypothetical REML fit estimates $\sigma_g^2=0.36$ and $\sigma_e^2=1.44$, in squared outcome units. At that fit, suppose the observed information matrix, in the parameter order $(\sigma_g^2,\sigma_e^2)$, is

J= \begin{pmatrix} 400&120\\ 120&100 \end{pmatrix}.

The determinant is $400(100)-120^2=25{,}600$. Thus

J^{-1}=\frac{1}{25{,}600} \begin{pmatrix} 100&-120\\ -120&400 \end{pmatrix} = \begin{pmatrix} 0.00390625&-0.0046875\\ -0.0046875&0.015625 \end{pmatrix}.

The genetic variance-component standard error is $\sqrt{0.00390625}=0.0625$ squared outcome units. The residual variance-component standard error is $\sqrt{0.015625}=0.125$ squared outcome units. The off-diagonal entries show estimated negative covariance between the two component estimates; they are not needed to calculate either marginal standard error, but they are part of their joint uncertainty.

For the genetic component, a nominal 95% Wald interval is $0.36\pm1.96(0.0625)=0.36\pm0.1225$, or approximately $(0.2375,0.4825)$ squared outcome units. For the residual component it is $1.44\pm1.96(0.125)=1.44\pm0.245=(1.195,1.685)$ squared outcome units. Both intervals are positive and the estimates are away from zero, so the main boundary objection to these symmetric intervals is not apparent here. Their coverage still depends on the model and the local normal approximation being adequate.

Profile Likelihood: Following the Likelihood Shape

A Wald interval uses the likelihood’s curvature at one point. A profile-likelihood interval instead examines how well each candidate value of the component can fit after the other parameters have been re-estimated. For a component $\psi$, the profile log likelihood is the maximum log likelihood at that fixed value, allowing nuisance parameters to vary:

\ell_p(\psi)=\max_{\lambda}\ell(\psi,\lambda),

where $\lambda$ denotes the nuisance parameters. A likelihood-ratio confidence set retains candidate values whose profile log likelihood is not too far below its maximum. For a single parameter and a regular interior fit, the familiar approximate 95% rule retains values for which twice the log-likelihood drop is no greater than 3.841, the 95th percentile of a chi-square distribution with one degree of freedom. Unlike a Wald interval, the resulting endpoints need not be equally far from the estimate.

Worked Example: An Asymmetric Profile-Likelihood Interval

Suppose a hypothetical REML analysis estimates a genetic variance component at $0.30$ squared outcome units. Its observed-information standard error is $0.08$ squared outcome units. The 95% Wald interval is $0.30\pm1.96(0.08)=0.30\pm0.1568$, or $(0.1432,0.4568)$ squared outcome units.

Now suppose the profile-likelihood calculation gives the following likelihood-ratio deviances $D(v)=2[\ell_p(0.30)-\ell_p(v)]$ at selected candidate values $v$:

Candidate $v$ (squared outcome units)	Profile deviance $D(v)$
0.12	4.9
0.16	3.6
0.30	0
0.45	2.4
0.52	4.1

Using the interior-parameter cutoff $D(v)\leq3.841$, the lower endpoint lies between 0.12 and 0.16, and the upper endpoint lies between 0.45 and 0.52. For an arithmetic illustration, linearly interpolate between the tabulated values. The lower crossing is $0.12+(4.9-3.841)(0.04)/(4.9-3.6)\approx0.1526$. The upper crossing is $0.45+(3.841-2.4)(0.07)/(4.1-2.4)\approx0.5093$. So this illustrative profile interval is approximately $(0.1526,0.5093)$ squared outcome units. In an actual analysis, a numerical profile curve would be optimized and the threshold crossing found from that curve; interpolation of a few points is only an approximation.

The profile interval is wider on the upper side than on the lower side, unlike the symmetric Wald interval. This is not a contradiction: the methods summarize different features of the likelihood. The profile follows its shape across candidate values, while the Wald calculation uses only local curvature at the estimate. The chi-square cutoff has its usual approximate coverage justification for a regular interior parameter; it should not automatically be treated as calibrated at a boundary.

What Changes Near a Boundary?

A variance component cannot be negative. If an estimate is close to zero, a symmetric normal approximation may assign substantial probability to impossible negative values. At the boundary itself, there is no two-sided neighborhood of admissible parameter values. The usual argument that a smooth likelihood maximum behaves like an unconstrained normal estimate therefore does not apply in its standard form.

Boundary proximity does not guarantee that the absolute standard error will increase in every model. Instead, it makes the usual uncertainty summary less dependable: the estimator’s sampling distribution can be skewed or pile up at zero, and the local curvature may be weak or unstable. A component that is poorly distinguished from other covariance contributions can also have a flat profile and broad uncertainty. These issues motivate checking the likelihood shape and using methods with appropriate boundary calibration when nominal coverage is important.

Worked Example: A Component Estimated at Zero

Suppose a hypothetical constrained REML fit estimates $\sigma_g^2=0$. The estimate is at the nonnegative boundary, so a conventional Wald interval cannot be interpreted as if the estimate were an ordinary interior normal estimate. For example, if software reports a local standard error of $0.12$ squared outcome units, the mechanically calculated symmetric interval is $0\pm1.96(0.12)=(-0.2352,0.2352)$. Negative genetic variance is not an admissible model value; simply clipping this interval at zero does not, by itself, establish correct 95% coverage.

For illustration, suppose the profile deviance for candidate genetic variances $v\geq0$ is $D(v)=(v/0.12)^2$. Applying the conventional interior cutoff $D(v)\leq3.841$ gives $v\leq0.12\sqrt{3.841}\approx0.2352$. The resulting set is $[0,0.2352]$ squared outcome units. This calculation demonstrates how a profile set respects the parameter space and is one-sided at the boundary. However, because the estimate is on the boundary, the conventional cutoff does not automatically guarantee 95% coverage. A boundary-appropriate calibration is needed before describing the set as a calibrated 95% confidence interval.

Practical Interpretation and Common Mistakes

A standard error for a variance component is measured in the same units as that variance component—for example, squared outcome units—not in the original outcome units. It summarizes estimated sampling variability under the fitted model and the approximation used. It is not the standard deviation of the observed outcomes, nor is it a direct measure of biological importance.

The observed-information and profile-likelihood approaches answer related but different questions. Observed information is efficient to compute and communicates local precision through a standard error. Profile likelihood can reveal asymmetry, a nearby boundary, or a likelihood that is flatter on one side. Neither method repairs a misspecified covariance model, and neither turns an estimated variance into a causal effect.

Mistake: treating every standard error as proof of a normal sampling distribution. The correction is to check whether the estimate is interior and whether a local quadratic approximation is plausible.
Mistake: reporting a negative endpoint as a possible variance. A negative endpoint signals that a symmetric interval is a poor representation of uncertainty under the constraint; it is not an admissible variance value.
Mistake: truncating a Wald interval at zero and assuming its coverage is unchanged. The correction is to recognize that truncation changes the interval procedure and does not automatically restore nominal coverage.
Mistake: treating a conventional profile-likelihood cutoff as universally valid. Its standard chi-square justification is for regular interior settings; boundary cases need appropriate calibration.
Mistake: comparing uncertainty without checking units or parameterization. State whether estimates and standard errors are for variances or standard deviations, and give the units in context.

Key takeaway: Use the inverse observed information for local standard errors when the fit is regular and interior; use the profile to examine likelihood shape. Near a variance boundary, normal and conventional likelihood-ratio approximations can be unreliable, so a boundary-aware uncertainty procedure may be required.

Check Your Understanding

Answer each question using the model assumptions and the distinction between local and profile-based uncertainty.

Why does the covariance matrix of multiple variance-component estimates generally come from the full inverse information matrix rather than reciprocals of the information matrix’s diagonal entries?
A variance-component estimate is $0.20$ squared outcome units with standard error $0.05$ squared outcome units. Calculate its approximate 95% Wald interval and state the units.
In what way can a profile-likelihood interval differ from a Wald interval even when both are computed from the same fitted model?
Why does an interval with a negative lower endpoint require caution when the parameter is a variance?
Why is it not enough to clip a Wald interval at zero to claim a calibrated 95% confidence interval at a boundary estimate?