Why Random-Effects Complexity Needs More Than a Convergence Check
A random-effects structure can be more elaborate than the observed grouping pattern can support. For instance, a model may assign each group a random intercept and several random slopes, along with an unrestricted covariance matrix among them. The fitted model may return coefficient estimates and a convergence message, yet still have covariance directions that the data barely identify. In that situation, estimates can be unstable or effectively lie on a boundary.
The central task here is to read three kinds of evidence together: the estimated random-effects covariance matrix, the optimizer’s gradient, and the curvature of the objective near its optimum. These diagnostics answer related but different questions. A singular covariance estimate points to a redundant or unsupported random-effect direction; a large gradient suggests the optimizer may not have reached a stationary point; and weak curvature indicates that the objective changes little in some parameter direction. None, by itself, proves that a random-effects structure is wrong.
For clarity, let \(Y_{ij}\) be the response for observation \(i\) in group \(j\), \(x_{ij}\) the fixed-effects predictor vector, and \(z_{ij}\) the random-effects predictor vector. A GLMM specifies a conditional response distribution from an exponential-family model, with conditional mean \(\mu_{ij}=E(Y_{ij}\mid b_j)\) and link relation \(g(\mu_{ij})=x_{ij}^{T}\beta+z_{ij}^{T}b_j\). The group random effects satisfy \(b_j\sim N(0,G)\), independently across groups; independence from covariates may also be assumed by the model. Responses are assumed conditionally independent given the random effects and predictors. The diagnostics below presume that this response model, link, predictor structure, and grouping structure are scientifically defensible. They assess support for the specified random-effects covariance, not the correctness of every part of the GLMM.
Three Diagnostics, Three Questions
If \(G\) is positive definite, all of its eigenvalues are positive. A zero eigenvalue means that some linear combination of random effects has estimated variance zero; a very small eigenvalue relative to the others indicates a nearly unsupported direction. In practice, software often calls a fit singular when an eigenvalue is zero or falls below a numerical tolerance. The numerical threshold is implementation-dependent, so the eigenvalues, scale, and tolerance should be examined rather than treating a software label as a universal scientific cutoff.
The eigenvectors matter as much as the eigenvalues: they show which combinations of random effects correspond to weakly estimated directions. With two random effects, a variance near zero or a correlation near \(+1\) or \(-1\) can signal near-singularity. With three or more, a weak direction can be a combination of effects even when no individual variance appears especially small and no pairwise correlation is exactly extreme. Eigenvalues also depend on units. Comparing them is most interpretable when random-effect predictors have been scaled sensibly, or when the covariance is expressed on a comparable standardized scale.
The optimizer provides a different check. Let \(\theta\) denote the parameter vector in the coordinates used by the optimizer, including its chosen representation of covariance parameters, and let \(Q(\theta)\) be the minimized objective, such as a negative log likelihood or its approximation. The gradient \(g=\nabla Q(\hat{\theta})\) records local slopes. A small gradient is evidence that the objective is near stationary in those coordinates; a large component can indicate incomplete optimization. But a gradient norm has no universal interpretation independent of parameterization, units, and the optimizer’s convergence tolerance.
Curvature is described locally by the Hessian \(H=\nabla^2 Q(\hat{\theta})\). At a well-behaved interior minimum, the Hessian should be positive definite. An eigenvalue near zero indicates a locally flat direction: moving along it changes the objective very little, so the corresponding parameter combination is weakly determined. A negative eigenvalue can indicate that the point is not a local minimum, or that numerical Hessian estimation is unreliable. Near a boundary, such as a variance of zero, interior conditions do not apply directly; boundary constraints and transformations can alter how gradients should be interpreted.
Why Gradient and Curvature Must Be Read Together
The local quadratic approximation to the objective at \(\hat{\theta}\), for a displacement \(\delta\), is
This approximation explains why a small gradient alone is not enough. If the gradient is small but the Hessian has a near-zero eigenvalue, the objective can still be almost flat over a substantial distance in that eigenvector direction. Conversely, an apparently small gradient may be concerning if it is large relative to the local curvature or the convergence tolerance used by the optimizer.
To see the role of curvature explicitly, take a unit eigenvector \(v\) of \(H\) with eigenvalue \(\lambda\), and move by \(\delta=t v\). Since \(Hv=\lambda v\), the quadratic approximation becomes
When \(\lambda\) is very small, the quadratic contribution is small even for a noticeable displacement \(t\). If \(g^{T}v\) is also small, the objective provides little information about movement in that direction. This is a justification for treating weak Hessian curvature as an identifiability warning, not a standalone test of model complexity. The approximation is local, and the conclusion depends on the coordinates and numerical accuracy of \(H\).
Worked Example: A Nearly Singular Intercept–Slope Covariance
Suppose a hypothetical binary-outcome study has 24 clinics and 12–30 patients per clinic. The fitted model includes a clinic random intercept and a random slope for a centered, standardized treatment-intensity predictor. Assume the GLMM’s response distribution and conditional-independence assumptions are appropriate. The estimated covariance matrix, in the stated predictor units, is
The estimated standard deviations are \(\sqrt{0.6400}=0.8000\) for the intercept and \(\sqrt{0.0900}=0.3000\) for the slope. Their estimated correlation is
The determinant is \(0.6400(0.0900)-(0.2394)^2=0.057600-0.05731236=0.00028764\), a small positive value. The two eigenvalues are approximately \(0.729606\) and \(0.000394\), so their ratio is \(0.000394/0.729606\approx0.000540\). The covariance is technically positive definite, but one direction is extremely weak relative to the other. Together with the estimated correlation of \(0.9975\), this is strong evidence of near-singularity. It indicates that the fitted data provide little support for separately estimating both random-effect directions as freely varying. It does not establish that the population correlation is exactly one or that the scientific slope effect is absent.
The eigenvalue ratio is meaningful here because the predictor was standardized and the random effects are in compatible model-scale units. If the slope predictor were measured in a very different scale, raw covariance eigenvalues could be misleading; one should first consider the parameter scaling and corresponding standardized covariance.
Worked Example: A Small Gradient Does Not Rule Out a Flat Direction
Consider a hypothetical fit to count outcomes from 18 hospitals, with a random intercept and random slope for a standardized workload predictor. Suppose the optimizer works in two scaled covariance-parameter coordinates, and the estimated negative log-likelihood has gradient \(g=(0.0008,-0.0004)^T\) and Hessian \(H=\begin{pmatrix}0.02&0\\0&8\end{pmatrix}\). Assume the optimizer reports a small gradient under its declared tolerance and the Hessian estimate is numerically stable. These conditions make the local calculation informative, but do not by themselves verify the response model or the design’s identifiability.
For a one-unit displacement along the first coordinate, \(\delta=(1,0)^T\), the gradient contribution is \(g^T\delta=0.0008\) and the curvature contribution is \(\frac12\delta^TH\delta=\frac12(0.02)=0.0100\). Thus the full quadratic change is approximately \(0.0008+0.0100=0.0108\). Along the second coordinate, \(\delta=(0,1)^T\), the contributions are \(-0.0004\) and \(\frac12(8)=4.000\), giving a full approximate change of \(-0.0004+4.000=3.9996\).
The first direction is much flatter locally than the second, even though the gradient components are both small. The displayed objective changes are in negative-log-likelihood units and are local approximations, not confidence intervals or exact changes. The result warrants checking whether the weak direction corresponds to a covariance combination and whether the Hessian calculation is reliable. It does not alone justify dropping a term: coordinates, parameter scaling, boundary constraints, and the optimizer’s tolerance must be considered.
Worked Example: Reading Diagnostics Across a Three-Slope Structure
Suppose a hypothetical continuous-outcome study has 16 schools, each with 20–35 pupils, and fits school-specific random intercepts and slopes for two standardized predictors: instructional time and class size. Assume the conditional response model is suitable, predictors vary sufficiently within schools for the intended slopes to be considered, and the covariance estimate is computed using a stable fitting procedure. The three eigenvalues of the estimated covariance matrix, on a common standardized scale, are \(0.810\), \(0.180\), and \(0.000600\).
The smallest-to-largest eigenvalue ratio is \(0.000600/0.810=0.000741\). Thus one of the three random-effect combinations has variance tiny relative to the largest direction. The eigenvector associated with \(0.000600\) should be inspected to determine which combination of intercept and slopes is weakly supported. A pairwise correlation table alone may not identify that combination. Suppose also that the optimizer reports a gradient norm of \(0.0007\), within its stated scaled tolerance, while the Hessian eigenvalues are \(0.015\), \(1.6\), and \(9.2\), all positive. The small positive Hessian eigenvalue is additional evidence of a flat parameter direction, whereas the small gradient suggests that the optimizer may have located a stationary point.
Taken together, these outputs support the conclusion that the fitted three-effect covariance contains a weakly supported direction, not that the optimization necessarily failed. The conclusion remains conditional on the Hessian being computed accurately and on reasonable scaling. With only 16 schools, the available group-level information is limited; this context reinforces caution but does not create a universal group-count threshold or prove that the structure is scientifically inappropriate.
How to Weigh the Evidence
| Diagnostic | Main question | Important limitation |
|---|---|---|
| Estimated covariance eigenvalues and eigenvectors | Is a random-effect variance direction estimated near zero, and which combination is involved? | Thresholds depend on tolerance and scale; estimates do not prove a population variance is exactly zero. |
| Gradient at the fitted parameters | Is the objective locally close to stationary in the optimizer’s coordinates? | Size depends on parameterization, scaling, and convergence criteria; a small gradient cannot rule out flat curvature. |
| Hessian eigenvalues | Is the objective sharply or weakly curved in local parameter directions? | Numerical Hessians can be unstable, and boundary constraints complicate interior interpretations. |
| Profile likelihood or likelihood-based checks | How does the objective behave across a wider range of a chosen parameter? | These require careful treatment of boundaries and do not remove dependence on the model and data. |
The main alternative to these local diagnostics is to examine objective behavior more broadly, for example with a profile likelihood for a particular covariance parameter. This can reveal asymmetry or boundary behavior that a local Hessian misses, but it is still conditional on the fitted model. The earlier material on random-effects singularity and flat likelihood directions provides the relevant concepts; here the emphasis is on combining them to diagnose overall structural support.
A useful reading sequence is to inspect the covariance estimate and its eigen-directions, examine the gradient in the optimizer’s actual parameterization, and then assess Hessian curvature and its numerical reliability. The sequence is not an automatic simplification rule. It is a way to separate three possible findings: a converged but weakly identified fit, an optimizer that has not settled, or numerical diagnostics too unstable to interpret. A failed gradient check and a singular covariance can occur together, but neither implies the other.
Diagnostics, Limitations, and Common Mistakes
Before interpreting a diagnostic, verify that the model and computation make it meaningful. In particular, check the optimizer’s convergence tolerances, parameter scaling, covariance parameterization, and whether the reported Hessian is positive and stable under the software’s calculation method. Consider whether the grouping design actually contains information for the proposed slopes, including within-group predictor variation and the number and size of groups. A covariance estimate can look unstable because the design offers little information, because the optimization is poor, or because finite-precision calculations are problematic. These causes should not be conflated.
A common mistake is to read “the optimizer converged” as “the random-effects structure is supported.” A convergence message concerns the numerical stopping criteria, not whether all covariance directions are well identified. The correction is to inspect covariance eigenvalues, gradient behavior, and curvature as complementary evidence. The opposite mistake is to see one small eigenvalue and claim that the optimizer failed or that the true variance is exactly zero. Near-singularity can occur at a valid boundary solution, and a small eigenvalue is evidence of weak support in this data set—not proof of a particular population value.
- Do not use a gradient norm without its scale. Read it alongside the software’s parameterization and stated tolerance.
- Do not interpret a Hessian eigenvalue in isolation. Confirm its sign, numerical stability, and the parameter direction it represents.
- Do not rely only on pairwise correlations. In higher-dimensional random-effects models, a weak direction can involve several effects.
- Do not treat every singular fit as a numerical failure. A stable optimizer can produce a boundary estimate because the data do not support a full-rank covariance.
- Do not treat local diagnostics as a complete model assessment. They do not establish that the response distribution, fixed-effects structure, or scientific assumptions are correct.
These diagnostics have practical limits. Hessian estimates can be sensitive to finite precision and numerical differentiation, and approximate likelihood methods can affect reported curvature. A near-zero eigenvalue is also scale-dependent unless parameters have been placed on sensible comparable scales. Finally, evidence that a covariance structure is weakly supported does not by itself specify which modeling change is appropriate; that requires a separate, scientifically reasoned decision.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A fitted random-effects covariance matrix has a small eigenvalue but a small optimizer gradient. What can and cannot be concluded?
- Why might a covariance matrix have no extreme pairwise correlation yet still be nearly singular when it contains three or more random effects?
- For \(g=(0.001,0)^T\), \(H=\begin{pmatrix}0.04&0\\0&5\end{pmatrix}\), and \(\delta=(1,0)^T\), calculate the quadratic approximation to the objective change. What does it suggest?
- Why is a small Hessian eigenvalue not, by itself, proof that the optimizer failed?
Solutions. (1) The small eigenvalue suggests a weakly supported covariance direction, while the small gradient suggests local stationarity in the optimizer’s coordinates. Together they are consistent with a converged but weakly identified fit, subject to scaling and numerical checks; they do not prove an exact zero population variance. (2) A linear combination of three or more effects can have near-zero variance even when no single pairwise correlation is extreme. The corresponding eigenvector identifies that combination. (3) The gradient contribution is \(g^T\delta=0.001\); the curvature contribution is \(\frac12\delta^TH\delta=\frac12(0.04)=0.020\). The full approximation is \(0.001+0.020=0.021\). The first coordinate has relatively weak local curvature compared with the second, but interpretation depends on parameter scaling. (4) A small positive Hessian eigenvalue indicates a flat local direction, which can occur at a stationary point. It signals weak information about a parameter combination, not necessarily an optimization failure.