Detecting Singular Random-Effects Covariance Matrices

Why a Covariance Matrix Can Be Singular

A GLMM may assign each group a vector of random effects, such as a random intercept and one or more random slopes. The covariance matrix describes the population variation and joint dependence of those group-level effects. When the fitted matrix is singular, at least one linear combination of the modeled random effects has estimated variance zero. The estimate therefore lies on the boundary of the set of positive-semidefinite covariance matrices.

This is a matrix-level diagnostic. Looking only at the variances on the diagonal can miss a weakly varying combination involving several random effects. Conversely, a small eigenvalue does not by itself establish that a particular named slope has no population variation. The associated eigenvector identifies the combination whose variation is small.

Definition: Let \(G\) be the \(q \times q\) covariance matrix of a group’s \(q\)-dimensional random-effects vector \(b\), where \(b \sim N(0,G)\). The matrix is singular if its rank is less than \(q\), equivalently if at least one of its eigenvalues is zero. A fitted covariance estimate is numerically near-singular when its smallest eigenvalue is positive but small relative to the largest eigenvalue under a stated numerical tolerance.

A covariance matrix is symmetric and positive semidefinite: for every vector \(a\), \(a^\mathsf{T}Ga \geq 0\). Its eigenvalues are therefore nonnegative in exact arithmetic. In a GLMM with multiple grouping factors, apply the diagnostic to each estimated random-effects covariance block separately; a singular block identifies a rank deficiency in that block.

Eigenvalues Identify Weak Directions

The spectral decomposition of a symmetric covariance matrix is \(G=U\Lambda U^\mathsf{T}\), where the columns \(u_j\) of \(U\) are orthonormal eigenvectors and \(\Lambda\) is diagonal with eigenvalues \(\lambda_j\). Each eigenvalue measures variance in the direction of its eigenvector:

\operatorname{Var}(u_j^\mathsf{T}b) = u_j^\mathsf{T}G u_j = \lambda_j.

Thus, if \(\lambda_j=0\), the linear combination \(u_j^\mathsf{T}b\) has zero variance under the fitted distribution. If \(\lambda_j\) is very small relative to the other eigenvalues, that direction has little estimated variation. The estimated covariance has exact rank equal to the number of strictly positive eigenvalues; for a numerical estimate, rank is assessed relative to a tolerance.

Key diagnostic: Examine both the eigenvalues and their eigenvectors. The eigenvalues indicate how much variation is estimated in each principal direction; the corresponding eigenvectors indicate which combination of random-effect coordinates defines that direction.

In computation, an exact zero is not required to flag a practical problem. One common diagnostic compares each eigenvalue with the largest:

r_j=\frac{\lambda_j}{\lambda_{\max}}, \qquad \text{flag direction }j\text{ if }r_j \leq \tau,

where \(\tau\) is a chosen relative tolerance. This is a diagnostic convention, not a universal inferential cutoff. The tolerance should be stated, and conclusions should be checked for sensitivity to reasonable choices. Eigenvalues have units of variance, and their magnitudes depend on the scale of the random-effect coordinates. Rescaling a predictor can therefore change the numerical eigenvalue ratios even though an invertible change of coordinates preserves exact matrix rank.

Worked Diagnostics

Worked Example: An Exactly Singular Random-Intercept and Random-Slope Matrix

Suppose a hypothetical fitted GLMM has a random intercept and a random slope for a centered, scaled predictor. Its estimated covariance matrix, in a common coefficient scale, is

\widehat G= \begin{pmatrix} 0.64 & 0.48\\ 0.48 & 0.36 \end{pmatrix}.

The trace is \(0.64+0.36=1.00\), and the determinant is \[ 0.64(0.36)-(0.48)^2=0.2304-0.2304=0. \] The characteristic equation is therefore \(\lambda^2-\lambda=0\), giving eigenvalues \(1\) and \(0\). The estimated matrix has rank 1, not 2, so it is singular.

An eigenvector for eigenvalue 1 is \((0.8,0.6)^\mathsf{T}\). An eigenvector for eigenvalue 0 is \((-0.6,0.8)^\mathsf{T}\). Accordingly, the fitted variance of \(0.8b_0+0.6b_1\) is 1, while the fitted variance of \(-0.6b_0+0.8b_1\) is zero. The supported and unsupported directions are combinations of the intercept and slope coordinates; the calculation does not say that either coordinate alone has zero variance.

Worked Example: A Full-Rank but Near-Singular Three-Dimensional Estimate

In a hypothetical fit, suppose three random-effect coordinates are on comparable scales and the estimated covariance matrix is

\widehat G= \begin{pmatrix} 2 & 1 & 0\\ 1 & 2 & 0\\ 0 & 0 & 0.0008 \end{pmatrix}.

The first two coordinates have eigenvectors proportional to \((1,1,0)^\mathsf{T}\) and \((1,-1,0)^\mathsf{T}\), with eigenvalues 3 and 1. The third coordinate is an eigenvector with eigenvalue 0.0008. Thus the eigenvalues are \(3,1,0.0008\). All are positive, so the matrix has exact rank 3; nevertheless, the smallest eigenvalue is small relative to the largest.

For illustration, use the relative threshold \(\tau=0.001\). The ratio for the smallest eigenvalue is \(0.0008/3=0.000267\), which is below the threshold. This diagnostic flags the third direction as near-singular, while the other ratios are \(1/3\) and 1. The estimate does not have an exact zero eigenvalue, and another reasonable tolerance could classify it differently. The interpretation is that the fitted model estimates very little variation in the third direction relative to the others, not that the population variance is known to be exactly zero.

Worked Example: Why Scale and Tolerance Must Be Reported

Consider a hypothetical diagonal covariance estimate with eigenvalues \(100\), \(1\), and \(0.02\). The matrix is full rank because all three values are positive. Under a relative threshold of \(0.001\), the largest eigenvalue is 100, so the smallest ratio is \(0.02/100=0.0002\). This flags the third direction as near-singular.

Under an absolute cutoff of \(0.001\) on the eigenvalues, none would be flagged, since \(0.02>0.001\). The diagnostics differ because an absolute cutoff does not account for the overall scale of the matrix. A relative cutoff is more scale-aware, but it is not invariant to arbitrary changes in the scales of different random-effect coordinates. A responsible report would give the eigenvalues, the ratio or cutoff used, the predictor scaling, and the resulting random-effect direction. The threshold alone is not a scientific conclusion.

A Practical Diagnostic Sequence

A singular fit and a convergence failure are related only in the sense that both can occur in difficult models; they are not the same diagnosis. The earlier material on convergence messages, Hessian checks, and covariance parameterization helps assess the optimizer’s behavior. A converged optimizer can return a singular covariance estimate, and a convergence warning does not establish that the covariance estimate is singular.

Identify the covariance block.
Record which grouping factor and which random-effect coordinates the matrix describes. Confirm the predictor coding and scaling used to define those coordinates.

Check the matrix and its eigenvalues.
Verify that the estimated matrix is symmetric and positive semidefinite up to numerical precision. Obtain all \(q\) eigenvalues; do not infer matrix rank from diagonal variances alone.

Assess relative size.
Compare the smallest eigenvalues with \(\lambda_{\max}\), state the tolerance, and note whether the conclusion changes under nearby tolerances. Tiny negative values can arise from rounding or numerical computation; substantial negative values call for checking the fit and reported matrix.

Interpret the eigenvectors.
Translate each flagged direction back into the modeled random-effect coordinates. Describe it as a weakly supported combination unless the eigenvector aligns clearly with one coordinate.

Check the broader fit.
Read optimizer and Hessian diagnostics separately, and examine whether the fitted matrix and substantive estimates are stable under appropriate changes in scaling or model specification. A diagnostic flags a limitation of the fitted structure; it does not by itself determine the remedy.

The eigenvalue check is more informative than several common alternatives. Diagonal variances show variation in named coordinates but can miss a weak linear combination. The determinant is zero for an exactly singular matrix, but a determinant can be numerically tiny merely because of the units and scale of the matrix, so it is a poor standalone near-singularity measure. Covariance entries alone likewise do not reveal the principal directions. Eigenvalues and eigenvectors jointly expose both the amount and the direction of estimated variation.

What the Diagnostic Does—and Does Not—Establish

A singular estimate says that the fitted covariance matrix has fewer varying dimensions than the model allows, either exactly or to the chosen numerical tolerance. It does not prove that the true population covariance is singular, that a random effect is scientifically absent, or that the optimizer failed. It may reflect limited information in the observed data, a model structure that is too ambitious for the available information, a boundary estimate, or numerical sensitivity. The diagnostic alone cannot distinguish these explanations.

Interpretation also depends on parameterization. Under an invertible coordinate transformation \(b^*=Ab\), the covariance becomes \(G^*=AGA^\mathsf{T}\), which has the same exact rank as \(G\). However, its eigenvalues and eigenvectors generally change. This is why exact rank deficiency is a coordinate-invariant property, while a numerical near-singularity judgment based on eigenvalue ratios requires attention to scaling and coding.

A common mistake is to conclude, “The smallest eigenvalue is nearly zero, so the random slope has no variation.” That conclusion is justified only if the corresponding eigenvector points almost entirely in the slope coordinate, and even then it describes the fitted estimate rather than proving a population fact. The correction is to report the eigenvalue ratio and its tolerance, inspect the eigenvector, and state precisely which fitted combination has little estimated variation.

Another mistake is to treat an optimizer’s convergence message as a singularity test. Convergence concerns the numerical search for an optimum; singularity concerns the rank of the covariance matrix at the reported estimate. Check and report these issues separately. A covariance diagnostic is most useful when it is tied to the fitted model’s coordinates and considered alongside, rather than substituted for, other fit diagnostics.

Key takeaway: Diagnose singularity from the rank structure of the fitted covariance matrix. Use eigenvalues to quantify weak directions and eigenvectors to identify the random-effect combinations involved; report the scaling and tolerance, and do not overinterpret a near-zero estimate as proof about the population.

Check Your Understanding

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

A fitted covariance matrix has four eigenvalues \(2.0, 0.7, 0.1,\) and \(0\). What is its rank, and what does the zero eigenvalue imply?
A matrix has eigenvalues \(5, 1,\) and \(0.002\). Under a relative threshold of \(0.001\), is it flagged as near-singular?
Why can diagonal variance estimates alone fail to identify a weakly varying random-effect direction?
Does a near-zero eigenvalue establish that the corresponding population random-effect variance is exactly zero? Explain.

Solutions. (1) Its rank is 3 because three eigenvalues are positive. The zero eigenvalue means one linear combination of the four random-effect coordinates has zero variance under the fitted covariance matrix. (2) The smallest ratio is \(0.002/5=0.0004\), which is below \(0.001\), so it is flagged under that stated rule. (3) A weak direction may be a combination of coordinates even when each diagonal entry is positive; eigenvectors identify such combinations. (4) No. The near-zero value describes the fitted estimate relative to the stated numerical rule. It is not proof that the true population variance is exactly zero.