Tutorials › Statistical Genetics › Boundary Constraints in Variance Estimation

Linear Mixed Models in Genetics · Tutorial 504 of 989

Boundary Constraints in Variance Estimation

Learn how nonnegative variance constraints shape REML estimation, including how to recognize and interpret a boundary estimate.

Intermediate 10 min read

What You'll Learn

  • State the nonnegative parameter constraints in a genetic mixed model
  • Explain why a variance-component estimate can lie at zero
  • Derive a constrained estimate in a simple REML contrast model
  • Apply boundary conditions to an optimization with linked variance components
  • Interpret a zero estimate without treating it as proof of no genetic variation

Why Variance Estimates Have Boundaries

A variance cannot be negative: it is the expected squared deviation from a mean. In a genetic mixed model, this elementary fact has a practical consequence. Estimation must search only over parameter values that produce valid covariance matrices. The best-fitting model can therefore occur at the edge of the allowed parameter space, such as a genetic variance estimate of zero, rather than at an interior point where the likelihood has an ordinary zero derivative.

A boundary estimate is not automatically an error or a claim that the corresponding biological source of variation is absent. It says that, among the allowed parameter values and given the data and model, the fitted likelihood is best at the boundary. The constraint changes how optimization works and how a zero estimate should be read.

Definition: A variance component is a nonnegative parameter that scales a covariance contribution in a model. For example, in \(y=X\beta+g+e\), suppose \(g\sim N(0,\sigma_g^2K)\), \(e\sim N(0,\sigma_e^2I)\), and the random terms are independent. Then \(V=\operatorname{Var}(y)=\sigma_g^2K+\sigma_e^2I\), with \(\sigma_g^2\geq 0\) and \(\sigma_e^2>0\). Here \(K\) is a known positive semidefinite genetic relationship matrix. The positive residual variance keeps \(V\) positive definite when \(K\) is positive semidefinite.

The constraint on \(\sigma_g^2\) is essential: a negative multiplier of \(K\) need not define a valid covariance matrix. In the usual REML setup, the likelihood is evaluated using error contrasts that remove the fixed effects. If \(C\) has columns spanning the contrast space, so that \(C'X=0\), then \(C'y\) has covariance \(C'VC\). The REML estimate maximizes the likelihood of these contrasts over the permitted variance parameters. Those parameters include their boundaries.

What Constrained Optimization Means

For a single parameter \(\theta\geq 0\), an unconstrained optimization might place its maximum at a negative value. That value is not an admissible variance estimate. The constrained optimizer instead finds the maximum over \([0,\infty)\), which may be at \(\theta=0\) or at an interior value \(\theta>0\). With several components, the allowed region is a set of nonnegative parameter values; a maximum can occur on a face or corner where one or more components equal zero.

The boundary condition can be stated using the negative log likelihood \(Q(\theta)=-\ell(\theta)\), which is minimized. At an interior minimum, the derivative with respect to a component is zero. At a minimum where \(\theta_j=0\), feasible small changes in that component can only increase it. Thus, when \(Q\) is differentiable, its derivative in that coordinate at the boundary must be nonnegative. Equivalently, for the log likelihood, the score in a zero component must be nonpositive. These are necessary local conditions; they do not by themselves guarantee a global optimum.

Key result: A constrained variance estimate is obtained by optimizing the likelihood over the nonnegative parameter space, not by accepting an inadmissible negative estimate. If the best fit over that space occurs when a component is zero, the estimate is on the boundary. Other components can affect where the optimum lies, so simply replacing every negative unconstrained estimate by zero is not generally a valid multicomponent procedure.

A One-Contrast Calculation

A simple likelihood term makes the role of the boundary visible. Suppose a Gaussian contrast has observed value \(z\) and variance \(v\). Apart from a positive constant factor and additive constants, its negative log likelihood is \(Q(v)=\log(v)+z^2/v\), for \(v>0\). Write \(a=z^2\). When \(a>0\),

$$ \frac{dQ}{dv}=\frac{v-a}{v^2}. $$

The derivative is negative for \(0<v<a\) and positive for \(v>a\), so \(Q\) is minimized at \(v=a\). The condition \(a>0\) matters. If \(a=0\), then \(Q(v)=\log(v)\), which has no minimum over \(v>0\): it decreases without bound as \(v\to 0+\). A variance model with a positive known baseline can prevent that limit by imposing a positive lower bound on \(v\).

Worked Example: A Genetic Component Added to a Baseline Variance

Suppose one standardized error contrast has observed value \(z=0.6\), and its variance is \(v=1+\tau\), where the baseline variance is known to be 1 and the genetic component satisfies \(\tau\geq 0\). Assume the contrast is Gaussian and the baseline is fixed. What value of \(\tau\) maximizes the likelihood?

The observed square is \(a=z^2=0.6^2=0.36\). The unrestricted best variance for this one-contrast likelihood term is \(v=a=0.36\), because \(a>0\). But the model permits only \(v=1+\tau\geq 1\), so \(v=0.36\) is inadmissible.

For every feasible \(v\geq 1\), we have \(v>a\), and therefore \(dQ/dv=(v-a)/v^2>0\). The negative log likelihood increases as \(v\) increases from 1, so its minimum over the feasible region is at \(v=1\). Consequently, \(\hat{\tau}=v-1=0\). This is a boundary estimate, not a negative variance.

The conclusion depends on the stipulated baseline: the model never permits \(v\) to approach zero. It is different from optimizing \(\log(v)\) freely over \(v>0\) when the observed square is zero, where the objective has no minimum.

Multiple Contrasts: Deriving a Constrained Estimate

A useful extension is a set of \(m\) independent Gaussian error contrasts, each with variance \(1+\tau\), where the baseline variance is known and \(\tau\geq 0\). This is a simplified model for illustrating the boundary calculation; it is not a general formula for every REML problem. Let the observed contrasts be \(z_1,\ldots,z_m\), and set \(S=\sum_{i=1}^m z_i^2\). Ignoring constants, the negative log likelihood is

$$ Q(\tau)=\frac{1}{2}\left[m\log(1+\tau)+\frac{S}{1+\tau}\right],\qquad \tau\geq 0. $$

Differentiating gives

$$ \frac{dQ}{d\tau} =\frac{1}{2}\frac{m(1+\tau)-S}{(1+\tau)^2}. $$

The denominator is positive. If \(S/m>1\), the derivative changes from negative to positive at \(\tau=S/m-1\), which is feasible and minimizes \(Q\). If \(S/m\leq 1\), the derivative at zero is nonnegative and remains nonnegative for \(\tau\geq 0\), so the minimum is at zero. Thus in this particular model,

$$ \hat{\tau}=\max\left(0,\frac{S}{m}-1\right). $$

This derivation shows explicitly how the constraint truncates the feasible search at zero. It does not imply that REML estimates in more complex genetic models can be calculated by this expression: there, components may interact through the covariance structure and must be optimized jointly.

Worked Example: Three Contrasts Favor the Boundary

Suppose three independent Gaussian error contrasts have variance \(1+\tau\), with known baseline variance 1 and \(\tau\geq 0\). Their hypothetical observed values are \(0.5\), \(-0.7\), and \(0.9\). The assumptions for the calculation are the same as in the contrast model above: independent normal contrasts and a common variance.

The sum of squared contrasts is

$$ S=0.5^2+(-0.7)^2+0.9^2=0.25+0.49+0.81=1.55. $$

Here \(m=3\), so \(S/m=1.55/3=0.5167\), which is below the baseline variance of 1. The constrained estimate is

$$ \hat{\tau}=\max(0,0.5167-1)=0. $$

The derivative at the boundary confirms the direction: \(Q'(0)=\frac{1}{2}(3-1.55)=0.725>0\). Increasing \(\tau\) from zero initially increases the negative log likelihood. The fitted contrast variance is therefore \(1+\hat{\tau}=1\), and the genetic component is estimated at the boundary.

Linked Components: Why Separate Truncation Can Fail

A genetic mixed model commonly has more than one unknown variance component. Even if each component is individually constrained to be nonnegative, the variance contributions may be linked through a shared residual variance. An unconstrained solution that is inadmissible does not generally reveal the constrained solution by itself. The likelihood must be optimized under all constraints together.

Worked Example: Two Linked Variances and a Zero Genetic Estimate

Consider two independent Gaussian error contrasts with covariance matrix \(\operatorname{diag}(\sigma_e^2+\sigma_g^2,\sigma_e^2)\). This is a small contrast-space illustration of a genetic covariance contribution: the first direction has genetic loading 1 and the second has loading 0. The residual variance satisfies \(\sigma_e^2>0\), and \(\sigma_g^2\geq 0\). Suppose the hypothetical squared observed contrasts are \(a=0.4\) in the first direction and \(b=1.6\) in the second.

Write \(v_1=\sigma_e^2+\sigma_g^2\) and \(v_2=\sigma_e^2\). The parameter constraints become \(v_1\geq v_2>0\). Apart from a factor of one half and constants, the negative log likelihood is

$$ Q(v_1,v_2)=\log(v_1)+\frac{0.4}{v_1}+\log(v_2)+\frac{1.6}{v_2}. $$

If the two variances were optimized independently, their minimizing values would be \(v_1=0.4\) and \(v_2=1.6\). Both observed squares are positive, so each scalar minimum follows from the derivative calculation above. But these values violate \(v_1\geq v_2\); they cannot arise from nonnegative \(\sigma_g^2\).

To find the constrained minimum, first fix \(v_2\). The first pair of terms is minimized at \(v_1=0.4\) if that value is feasible, and otherwise at the smallest feasible value \(v_1=v_2\). If \(v_2\leq 0.4\), the best \(v_1\) is 0.4; the second pair of terms decreases as \(v_2\) increases up to 1.6, so the best point in this range is \(v_2=0.4\). If \(v_2\geq 0.4\), the best \(v_1\) is \(v_2\), and the objective reduces to \(2\log(v_2)+2/v_2\). Its derivative is \(2(v_2-1)/v_2^2\), so its minimum on this range is at \(v_2=1\). This value gives a lower objective than the endpoint \(v_2=0.4\), because the reduced objective decreases between 0.4 and 1.

The constrained optimum is therefore \(v_1=v_2=1\). Translating back gives \(\hat{\sigma}_e^2=1\) and \(\hat{\sigma}_g^2=v_1-v_2=0\). This example shows why jointly respecting the covariance constraints matters: separate componentwise estimates do not solve the constrained problem.

How to Interpret a Zero Estimate

A zero estimate means that the fitted likelihood does not improve by assigning a positive value to that component under the specified model and constraints. It is a statement about the estimate, not a proof that the population variance is exactly zero. Limited information, the observed sample, correlations among covariance contributions, or features of the model can all affect where the maximum occurs.

In computation, variance parameters are often optimized directly subject to nonnegativity, or represented in a way that enforces nonnegativity. A boundary-capable method must allow a parameter to equal zero; a method that requires every parameter to be strictly positive can return an arbitrarily small positive value instead. Convergence should be checked against the constraints and the likelihood, especially when a reported component is near zero. A small positive numerical value and an exact boundary estimate can be practically similar, but they are not identical mathematical results.

1
Specify the covariance model.
Write the variance components and the covariance contribution each one scales. State which parameters must be nonnegative and which must be strictly positive for a valid covariance matrix.
2
Optimize within the feasible region.
Use the REML likelihood of the error contrasts and enforce all parameter constraints jointly. Do not treat an inadmissible unconstrained solution as a fitted variance estimate.
3
Check the boundary solution.
If a component is zero, confirm that moving into the feasible region does not improve the likelihood and that the reported covariance matrix remains valid.
4
Report the estimate in context.
Describe a zero as a boundary estimate under the fitted model. Do not interpret it alone as proof that the corresponding source of population variation is absent.

Common Mistakes

  • Allowing a negative variance. A negative fitted component can make the implied covariance invalid. The estimate must respect the covariance model’s permitted parameter space.
  • Assuming every unconstrained estimate can be clipped at zero. With linked components, the constrained optimum may require changing several parameters together. Optimize the likelihood subject to all constraints.
  • Calling a boundary estimate proof of no genetic effect. Zero is the value that maximizes the fitted likelihood among the allowed values; it does not establish that the true population component is exactly zero.
  • Ignoring the zero-observation case in a scalar likelihood argument. For \(a>0\), \(\log(v)+a/v\) is minimized at \(v=a\). For \(a=0\), it is unbounded below as \(v\to0+\) if \(v\) is unrestricted over positive values. A positive baseline changes the feasible range and must be included in the argument.
Key takeaway: Nonnegative variance constraints define the admissible covariance models. A REML optimum can lie at zero, and such a boundary solution must be found by constrained optimization. Its interpretation is limited to the fitted model and data; it is not, by itself, proof that the underlying variance is absent.

Check Your Understanding

Use the covariance constraints and likelihood behavior to answer the following questions.

  1. Why must a genetic variance component multiplying a positive semidefinite relationship matrix be constrained to be nonnegative?
  2. For a single contrast with squared value \(a>0\), explain why \(\log(v)+a/v\) is minimized at \(v=a\). What changes if \(a=0\) and \(v\) can approach zero?
  3. In the three-contrast example, what feature of the data makes the constrained estimate of \(\tau\) equal to zero?
  4. In the two-variance example, why is it invalid to use the independently optimized values \(v_1=0.4\) and \(v_2=1.6\)?
  5. What does a zero variance-component estimate say, and what does it not prove?
Clinical Biostats