Non-Positive-Definite Random-Effects Covariance

Why a Covariance Estimate Can Stop a Fit

In a population pharmacometric model, random effects describe how individual parameters vary around their population values. Their covariance matrix determines both the amount of variation in each random effect and how the effects vary together. If that matrix is not positive definite, a likelihood calculation that requires its inverse and determinant may fail, or the optimizer may be unable to take a valid step.

A reported failure needs interpretation. A matrix can be indefinite, with at least one negative eigenvalue, or singular, with at least one zero eigenvalue. Both are non-positive-definite, but they do not mean the same thing: singularity can represent a boundary estimate with a zero-variance direction, while indefiniteness is not a valid covariance matrix. The distinction helps identify whether the problem lies in one variance, a covariance or correlation, or a combination of several components.

Definition: Let \(\eta=(\eta_1,\ldots,\eta_q)^\mathsf{T}\) be a vector of random effects with covariance matrix \(\Omega\). Its diagonal entries \(\omega_j^2=\operatorname{Var}(\eta_j)\) are variances, and its off-diagonal entries \(\omega_{jk}=\operatorname{Cov}(\eta_j,\eta_k)\) are covariances. The matrix is positive definite if \(a^\mathsf{T}\Omega a>0\) for every nonzero vector \(a\); it is positive semidefinite if \(a^\mathsf{T}\Omega a\geq 0\) for every \(a\). A covariance matrix must be symmetric and positive semidefinite. A likelihood calculation that uses an ordinary multivariate normal density with an inverse and determinant generally requires it to be positive definite.

When both variances are positive, the correlation for a pair is \(\rho_{jk}=\omega_{jk}/(\omega_j\omega_k)\). It must lie between −1 and 1. Equivalently, the magnitude of the covariance cannot exceed the product of the two standard deviations. This pairwise check is useful, but in a matrix with three or more random effects it does not, by itself, establish that the whole matrix is positive definite.

Conditions and a Two-Component Check

For two random effects, write the covariance matrix using variances \(v_1,v_2\) and covariance \(c\). Positive definiteness requires \(v_1>0\) and a positive determinant. The determinant condition gives a direct connection to the correlation bound.

\Omega= \begin{pmatrix} v_1 & c\\ c & v_2 \end{pmatrix}, \qquad \det(\Omega)=v_1v_2-c^2.

To justify the check, complete the square in the quadratic form. For any \(a_1,a_2\),

\begin{aligned} \begin{pmatrix}a_1&a_2\end{pmatrix} \Omega \begin{pmatrix}a_1\\a_2\end{pmatrix} &=v_1a_1^2+2ca_1a_2+v_2a_2^2\\ &=v_1\left(a_1+\frac{c}{v_1}a_2\right)^2+ \left(v_2-\frac{c^2}{v_1}\right)a_2^2. \end{aligned}

If \(v_1>0\), this expression is positive for every nonzero pair precisely when \(v_2-c^2/v_1>0\), or \(v_1v_2-c^2>0\). Thus positive variances and a positive determinant are sufficient and necessary for positive definiteness in this two-component case. When both variances are positive, dividing \(c^2<v_1v_2\) by \(v_1v_2\) gives \(|\rho_{12}|<1\). Equality gives a singular matrix, not a positive-definite one.

Worked Example: A Pairwise Covariance That Is Too Large

Suppose two hypothetical random effects have estimated variances \(v_1=0.0400\) and \(v_2=0.00250\), in their respective parameter units squared, and estimated covariance \(c=0.0120\), in the product of those units. Check whether these entries form a valid positive-definite covariance matrix.

The standard deviations are \(\sqrt{0.0400}=0.200\) and \(\sqrt{0.00250}=0.0500\), in the random effects' respective units. Their product is \(0.0100\) in the covariance's units. Therefore the implied correlation is

\rho_{12}=\frac{0.0120}{(0.200)(0.0500)}=1.20.

A correlation of 1.20 is outside the allowed range. The determinant confirms the problem:

\det(\Omega)=(0.0400)(0.00250)-(0.0120)^2 =0.000100-0.000144=-0.0000440.

Because the determinant is negative, the matrix is indefinite; it is not a valid covariance matrix. The immediate offending pair is clear: the covariance is too large in magnitude for the estimated variances. This check does not establish why the optimizer produced those entries. The cause could involve an invalid intermediate step, a parameterization or numerical issue, or weak information in the data; the matrix alone does not decide among them.

Use Eigenvalues for the Whole Matrix

For a real symmetric covariance matrix, all eigenvalues are real. The matrix is positive definite exactly when every eigenvalue is positive, and positive semidefinite exactly when every eigenvalue is nonnegative. A negative eigenvalue therefore establishes indefiniteness even when every variance is positive and every pairwise correlation is within its allowed range. An eigenvalue equal to zero indicates a singular direction.

For matrices larger than two by two, eigenvalues provide a whole-matrix check; correlations and selected submatrices help investigate which components are implicated. The most negative eigenvalue signals a violation, but does not alone identify its practical cause. Pay attention to the matrix's scale as well: a very small eigenvalue may be meaningful relative to the other eigenvalues, or it may be a numerical or rounding effect. Inspect estimates at full precision where available, and use tolerances appropriate to the matrix's scale rather than treating every tiny negative display value as equally serious.

Worked Example: Pairwise-Valid Correlations That Fail Together

Consider a hypothetical correlation matrix for three standardized random effects. Each diagonal entry is 1, and each pairwise correlation is −0.600. Every individual correlation is within [−1, 1]. Is the complete matrix positive definite?

R= \begin{pmatrix} 1 & -0.600 & -0.600\\ -0.600 & 1 & -0.600\\ -0.600 & -0.600 & 1 \end{pmatrix}.

The vector \((1,1,1)^\mathsf{T}\) is an eigenvector: multiplying by \(R\) gives \((-0.200,-0.200,-0.200)^\mathsf{T}\). Thus one eigenvalue is −0.200. Any vector whose entries sum to zero, such as \((1,-1,0)^\mathsf{T}\), is multiplied by 1.600, so the other two eigenvalues are 1.600. The eigenvalues are therefore −0.200, 1.600, and 1.600.

The matrix is indefinite despite having individually allowed pairwise correlations. No single pair violates the two-component bound; the three correlations are jointly incompatible with a valid covariance matrix. This points the diagnosis toward the set of components together rather than one obviously excessive pair.

Find Which Components Drive the Failure

A useful diagnostic sequence moves from entries to the whole matrix, then back to implicated components. Start with the variance estimates and units: a negative variance is invalid, while a variance near zero may indicate a boundary or numerical problem. Next calculate standard deviations and correlations for pairs with positive variances. A correlation outside [−1, 1] identifies an impossible pair, but correlations inside the range are only a necessary pairwise check. Then check symmetry and examine the eigenvalues of the full, unrounded matrix. Finally, inspect the pairs and subsets of components associated with the failed matrix, including their estimates at earlier or later optimizer iterations if those are available.

Confirm what failed.
Determine whether the message concerns an intermediate optimizer iterate or the final reported estimate. A transient invalid step and a final invalid estimate call for different interpretations.

Check variances and scale.
Look for negative entries, entries near zero, mismatched units, and values that have been rounded in a report. Covariances have product units; correlations are dimensionless.

Check pairs, then the full matrix.
For positive variances, compute correlations. Then assess symmetry and the eigenvalues: pairwise checks cannot rule out a multivariate failure.

Localize and document.
Review the variances, correlations, and relevant component subsets at full precision. Record the matrix and the fit state so the diagnosis can be reproduced; do not silently replace the estimate with a nearby valid one.

This sequence distinguishes three common patterns. A negative variance is an entry-level failure. A pairwise correlation outside its bounds points to an incompatible covariance and variance for that pair. If all variances are positive and all pairwise correlations appear valid but a full-matrix eigenvalue is negative, the issue is joint compatibility among multiple components. A zero or near-zero eigenvalue instead suggests a singular or nearly singular direction; check whether one or more variances are near zero and whether the result persists at full precision.

Worked Example: A Zero-Variance Boundary

Suppose a hypothetical two-random-effect fit reports

\Omega= \begin{pmatrix} 0.0400 & 0\\ 0 & 0 \end{pmatrix}.

The quadratic form is \(0.0400a_1^2\), which is nonnegative for every \((a_1,a_2)\), but it equals zero for the nonzero vector \((0,1)^\mathsf{T}\). The eigenvalues are 0.0400 and 0. The matrix is positive semidefinite and singular, not positive definite.

The second random effect has estimated variance zero, so its standard deviation is zero and its correlation with the first effect is undefined: the correlation formula would divide by zero. It would be incorrect to report this as an estimated correlation of zero. If the software requires a strictly positive-definite matrix for its likelihood calculation, this boundary matrix can still trigger a failure even though it is positive semidefinite. This pattern is different from an indefinite matrix with a negative eigenvalue.

Practical Interpretation, Alternatives, and Limits

Covariance parameterizations affect how a failure appears. Some implementations optimize unconstrained values and transform them into a covariance matrix; others use a factorization with constraints intended to keep the matrix positive definite. Such parameterizations can help prevent invalid matrices, but they do not guarantee that a fit is numerically stable or that the data support every estimated component. A parameterization may also approach a singular boundary, or reported rounding may make a valid matrix look invalid. Check the implementation's parameterization and the unrounded estimates before concluding that the model itself is defective.

The main diagnostic alternatives answer different questions. Entry and correlation checks are transparent and can quickly identify an offending pair. Eigenvalue checks assess the entire symmetric matrix, including joint failures that pairwise checks miss. A factorization-based check can be convenient when the software provides it, but a failed factorization does not by itself say which estimate caused the problem. These checks diagnose matrix validity; they do not establish that the random-effects structure is supported by the observations, nor do they determine the appropriate scientific model.

A covariance failure can disrupt computations that need a multivariate normal random-effects density, but the failure message alone does not prove that the likelihood is wrong or that a particular variance should be fixed to zero. An estimated variance near zero can arise from limited information, a boundary solution, or numerical instability. Conversely, an indefinite matrix is not a scientifically interpretable covariance estimate. Preserve the estimates and diagnostics, check units and precision, and avoid treating an automatic matrix adjustment as evidence that the underlying fit is resolved.

Common Mistakes

Checking only the variances. Positive diagonal entries do not guarantee a valid covariance matrix. Check the covariances and the full matrix as well.
Checking only pairwise correlations. All pairs can be individually valid while the joint matrix is indefinite. Use a whole-matrix check.
Calling every failure a negative variance. Indefiniteness can arise from covariances, and singularity can arise from a zero-variance direction. Inspect the eigenvalues and entries before naming the problem.
Interpreting a zero-variance correlation. If a standard deviation is zero, its correlation is undefined, not zero.
Relying on rounded output. Rounding can obscure a small negative eigenvalue or make a near-boundary estimate appear exact. Check full-precision values where possible.

Key takeaway: Diagnose covariance failures at three levels: variances, pairwise correlations, and eigenvalues of the complete symmetric matrix. Use the pattern of results to distinguish a bad entry, a pairwise incompatibility, a joint incompatibility, and a singular boundary; then verify precision and fit state before drawing conclusions.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

A two-effect covariance matrix has variances 0.09 and 0.04 and covariance 0.030. Calculate the correlation and determine whether the matrix is positive definite.
What does a negative eigenvalue establish about a symmetric covariance matrix? Does it identify a single offending pair?
For a covariance matrix with eigenvalues 0.12, 0.03, and 0, classify it as positive definite, positive semidefinite, or indefinite.
All three pairwise correlations in a three-effect model lie between −1 and 1, but the full matrix has a negative eigenvalue. Explain why the pairwise checks did not catch the failure.

Solutions. (1) The standard deviations are 0.300 and 0.200, so \(\rho=0.030/[(0.300)(0.200)]=0.500\). The determinant is \((0.09)(0.04)-(0.030)^2=0.00270>0\), and both variances are positive; the matrix is positive definite. (2) It establishes that the matrix is indefinite and therefore cannot be a valid covariance matrix. It does not, by itself, identify one pair; the failure may involve several components jointly. (3) The matrix is positive semidefinite and singular, not positive definite, because its eigenvalues are nonnegative and one is zero. (4) Pairwise bounds are necessary but not sufficient in dimensions three and higher; individually valid correlations can be jointly incompatible.