Convergence and Singularity in Mixed Models

Why Convergence and Singularity Need Separate Checks

A mixed model can produce plausible-looking coefficient estimates and still be unreliable. An optimizer may stop before finding a stable optimum; alternatively, it may find an optimum at the boundary of the allowed variance parameter space. These situations are related but not identical. A convergence warning concerns the numerical search. A singular fit concerns the estimated random-effects covariance structure. Either can occur without the other.

This distinction matters in single-cell differential expression because a model may include donor-level dependence, repeated samples, or cell-level random effects. A warning is not a reason by itself to discard a gene or to remove a random effect. The goal is to determine what the warning says about the fitted model, and whether a defensible simplification still answers the intended fixed-effect question.

Model Setup and Key Terms

For a Gaussian mixed model, write the response vector as \(y=X\beta+Zb+\varepsilon\). Here \(X\) is the fixed-effects design matrix, \(\beta\) contains fixed-effect coefficients, \(Z\) is the random-effects design matrix, and \(b\) is a vector of random effects. A standard model assumes \(b\sim N(0,G(\theta))\), \(\varepsilon\sim N(0,R(\theta))\), and independence between \(b\) and \(\varepsilon\). The variance parameters \(\theta\) must give valid covariance matrices: in particular, \(G(\theta)\) must be positive semidefinite. The model also relies on an appropriate mean structure and covariance structure for the observations.

For a count mixed model, a common specification is \(Y_{ij}\mid b_j\sim\operatorname{NB}(\mu_{ij},\phi)\), with \(\operatorname{Var}(Y_{ij}\mid b_j)=\mu_{ij}+\phi\mu_{ij}^{2}\) and \(\log(\mu_{ij})=\log(s_{ij})+x_{ij}^{T}\beta+z_{ij}^{T}b_j\). The offset \(s_{ij}\) accounts for the chosen library-size or exposure scale. Conditional on the random effects, observations are assumed to follow the specified count distribution; the random effects represent dependence within the relevant grouping unit, such as donor. Diagnostics cannot verify these assumptions automatically, but they can reveal when the specified model is difficult to estimate or poorly supported by the data.

Definition: A fit has converged numerically when the optimizer's stopping conditions indicate that it has reached a stationary point of the objective function, subject to the model's constraints. A fit is singular when an estimated random-effects covariance matrix is rank-deficient, or sufficiently close to rank-deficient that one or more random-effect directions have effectively zero estimated variance. A boundary variance estimate is a variance component estimated at or near zero, the lower boundary of its permitted range.

Convergence is an algorithmic statement, not proof that the model is scientifically appropriate or that a global optimum has been found. Likewise, singularity does not necessarily mean the optimizer failed. A zero variance can be a legitimate maximum-likelihood or restricted maximum-likelihood estimate when the data provide little support for variation in that random-effect direction.

What Singularity Means Geometrically

Consider a donor-specific random intercept and random slope, with covariance matrix \(G\). Its two eigenvalues describe the estimated variance along two perpendicular combinations of intercept and slope. If one eigenvalue is zero, the estimated random effects vary only along one direction; the two-dimensional covariance structure has collapsed to a lower-dimensional one. When the eigenvalue is merely very small, the model is nearly singular and may be numerically fragile.

For a two-dimensional covariance matrix with intercept variance \(\tau_{00}\), slope variance \(\tau_{11}\), and covariance \(\tau_{01}\), its determinant is \(\tau_{00}\tau_{11}-\tau_{01}^{2}\). A zero determinant indicates singularity. If both variances are positive, this occurs when the estimated correlation has magnitude one. A zero variance component also makes the determinant zero. In practice, software may report singularity using a tolerance, so a small positive eigenvalue can be treated as effectively zero.

\det(G)=\tau_{00}\tau_{11}-\tau_{01}^{2} =\tau_{00}\tau_{11}(1-\rho^{2}), \qquad \rho=\frac{\tau_{01}}{\sqrt{\tau_{00}\tau_{11}}}.

A Diagnostic Sequence

Start by recording the exact warning, optimizer, convergence code, and model specification. Then examine the numerical fit and covariance estimates rather than responding to the warning with an automatic model change. A useful sequence is:

Check the numerical stopping evidence.
Inspect the optimizer's convergence code, gradient, and Hessian diagnostics. A small gradient is reassuring only relative to the software's scaling and tolerance; no single universal cutoff applies across implementations.

Check parameter scales and data support.
Look for predictors on very different numerical scales, sparse factor combinations, few groups, or random-effect terms unsupported by the repeated measurements. Centering or rescaling continuous predictors can improve numerical conditioning without changing the underlying contrast when the model is reparameterized correctly.

Inspect the estimated random-effects structure.
Review variance components, correlations, and covariance-matrix eigenvalues where available. A boundary estimate or near-perfect correlation is evidence of weak support for some random-effect dimension, not necessarily a coding error.

Refit and compare defensible specifications.
Use a different optimizer or starting values to see whether the objective and key estimates are stable. If simplifying, state what term is removed and verify that the fixed-effects design and target contrast remain as intended.

A positive-definite Hessian at the reported solution is one useful check of local curvature, but its interpretation can be difficult near a boundary. A flat direction can make the Hessian unstable even when the optimizer has stopped. Conversely, a successful optimizer code does not rule out a singular covariance estimate. Examine both kinds of evidence.

Worked Example: Separating a Gradient Warning From a Singular Fit

Suppose a hypothetical pseudobulk Gaussian model estimates a condition contrast and a donor random intercept. The software reports a scaled gradient norm of 0.018, warns that its stopping tolerance was not met, and estimates the donor-intercept variance as 0.42 log-expression-units squared. The residual variance is 0.30 in the same squared units.

The warning concerns numerical optimization: the reported gradient norm is not zero, and the software says its specified tolerance was exceeded. The positive intercept variance alone does not indicate singularity; with only one random-intercept variance, there is no multi-dimensional covariance matrix whose rank can collapse in the way an intercept-slope matrix can.

A careful response is to refit with a supported alternative optimizer or starting values and compare the objective value, condition coefficient, and standard error. If a second fit has a scaled gradient norm of 0.0004, no warning, and a condition estimate of 0.31 rather than 0.30 log-expression units, the estimates are numerically close, though the original warning should still be documented. The difference is \(0.31-0.30=0.01\) log-expression units. This comparison is evidence of stability, not a universal proof that either model is correct.

Model Simplification and the Target Contrast

A random-effects term should represent a feature of the sampling or dependence structure, not merely be retained or removed to obtain a preferred p-value. A simplification is most interpretable when it changes the covariance model while leaving the scientific contrast explicit. Let the target be \(c^{T}\beta\), where \(c\) specifies the coefficient combination of interest. Keeping the same fixed-effects design \(X\), the same coding, and the same \(c\) preserves the definition of that fixed-effect contrast. It does not guarantee the same estimate, standard error, or validity: changing the random-effects distribution changes the covariance assumptions used for estimation and inference.

For example, removing a random slope while retaining a donor random intercept may preserve a condition coefficient in the linear predictor. But it assumes donors do not vary in their condition slopes, whereas the richer model allowed that variation. Removing a donor intercept as well is a larger change: it removes that modeled source of within-donor dependence. For nonlinear count models, the fixed-effect coefficient is conditional on the random effects; changing the random-effects structure can also affect the relationship between that coefficient and a population-averaged response contrast. State the target and interpretation before simplifying.

Useful simplifications may include removing an unsupported random-effect covariance, using a diagonal covariance structure when the scientific design does not require an estimated correlation, or removing a random slope when the data do not support slope variation and the target model assumptions justify its omission. Do not remove fixed-effect covariates required to adjust the target contrast, or change the contrast coding without translating the target. Preserve grouping terms that are needed to represent the dependence induced by the design unless there is a substantive and statistical justification for changing them.

Worked Example: Near-Singular Random Intercept and Slope

Suppose eight hypothetical donors each contribute observations in two conditions. A fitted model includes a donor random intercept and a donor-specific condition slope. Its estimated covariance matrix, in squared log-expression units, is

\widehat G= \begin{pmatrix} 0.040 & 0.039\\ 0.039 & 0.040 \end{pmatrix}.

The eigenvalues are \(0.040+0.039=0.079\) and \(0.040-0.039=0.001\) squared log-expression units. The smaller eigenvalue is only about \(0.001/0.079=0.0127\), or 1.27%, of the larger one. The estimated correlation is \(0.039/\sqrt{0.040(0.040)}=0.975\). These calculations indicate a nearly singular fit: the covariance structure estimates very little variation along one combination of intercept and slope.

Suppose the target is the fixed condition contrast, estimated as 0.28 log-expression units. A defensible next fit could remove the donor-specific slope while retaining the donor random intercept and the condition fixed effect, provided the analysis is willing to assume a common condition slope across donors. The fixed-effect contrast remains the condition coefficient, but its estimate and uncertainty must be obtained from the refitted model. This diagnosis does not establish that the slope variance is exactly zero, nor does it establish that the simpler model is automatically preferable.

Worked Example: Checking That a Contrast Is Still the Target

Consider a hypothetical count mixed model with a log link, a condition indicator coded 0 for control and 1 for treatment, and a donor random intercept. Suppose the fitted condition coefficient is 0.35 and its standard error is 0.12. The model's conditional treatment-to-control mean ratio is \(\exp(0.35)=1.419\), approximately a 41.9% higher conditional mean for treatment, for otherwise equal modeled terms.

Now suppose a near-zero random-slope variance prompts a refit that removes only the donor-specific condition slope, retains the donor intercept, and keeps the same fixed-effect coding. The refit reports a condition coefficient of 0.32 and standard error 0.10. Its corresponding conditional mean ratio is \(\exp(0.32)=1.377\), approximately 37.7% higher. The target coefficient remains the condition contrast, but the fitted estimate changes by \(0.32-0.35=-0.03\) log units and the standard error by \(0.10-0.12=-0.02\) log units.

This comparison shows what must be reported: the contrast definition, the random-effects change, and the resulting estimate and uncertainty. It does not justify choosing the model with the smaller standard error. The simplification is defensible only if the common-slope assumption is acceptable for the design and question.

Common Failure Modes and Limitations

A boundary estimate is not a license to ignore the design. With few donors, variance components can be imprecise, and a zero estimate may reflect limited information rather than absence of biological heterogeneity. Likewise, a singular fit does not prove that the associated random effect is scientifically irrelevant. It says the specified data and model do not support estimating all covariance dimensions distinctly.

Removing terms solely to silence a warning. First determine whether the warning is numerical, covariance-related, or both; then justify a change in terms of the design and target.
Equating convergence with validity. A converged model may still have a poor mean, dispersion, or dependence structure. Convergence checks complement, rather than replace, residual and influence diagnostics.
Interpreting a zero variance as proof of no heterogeneity. The estimate is constrained to be nonnegative and may land at zero when information is weak. It is not proof that the population variance is exactly zero.
Comparing fits without tracking the estimand. Changing fixed effects, coding, offsets, or random effects can alter what a coefficient means or how it is estimated. Keep the contrast explicit and describe any changed assumptions.
Using a routine chi-square reference for a variance-component test at zero. The null value lies on the boundary of the parameter space, so ordinary interior-parameter likelihood-ratio reference theory need not apply. Use an inference method appropriate to the boundary problem rather than assuming a standard reference distribution.

No diagnostic threshold is universal across software, likelihood approximations, and parameter scales. For cell-level models, computational burden and approximations can add further instability; for pseudobulk models, a small number of independent donors can limit information about random-effect variances. Record the warning, model, software options, and sensitivity results so readers can distinguish a stable simplification from an unexplained change.

Key takeaway: Diagnose numerical convergence and random-effects singularity separately. A simplification can preserve the definition of a fixed-effect contrast while changing its estimate, uncertainty, and model assumptions; justify the change, refit, and report the comparison.

Check Your Understanding

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

A model has a successful optimizer code but estimates a random intercept-slope correlation of 0.999. Is that enough to conclude the fit is nonsingular?
A variance component is estimated as zero. What does this establish, and what does it not establish?
A researcher removes a random slope but keeps the same fixed-effects design and contrast vector. What has been preserved, and what may have changed?
Why should a convergence warning not be treated as equivalent to a singular-fit warning?

Solutions. (1) No. A correlation near one can indicate an almost rank-deficient covariance matrix even when the optimizer reports success; examine the covariance eigenvalues and other diagnostics. (2) It establishes that the fitted estimate is at the permitted lower boundary. It does not prove the true population variance is exactly zero or that the random effect is scientifically unimportant. (3) The definition of the fixed-effect contrast is preserved, but covariance assumptions, estimates, standard errors, and potentially the interpretation of nonlinear-model effects may change. (4) A convergence warning concerns whether numerical stopping criteria were met; singularity concerns the rank of the estimated random-effects covariance. They are different properties and must be checked separately.