Why a Stepwise Diagnosis Matters
A GLMM warning is evidence to investigate, not a verdict that the model is either unusable or wrong. A convergence warning may arise because predictors have awkward scales, an optimizer has not located a stable optimum, the likelihood contains a weakly identified direction, or a variance estimate lies at its allowed boundary. A singular random-effects fit is related to boundary behavior but is a distinct diagnostic: it indicates that the estimated covariance structure has less than full rank, or is numerically close to having less than full rank.
These explanations call for different responses. Rescaling a predictor may improve numerical conditioning without changing the statistical model. Trying another optimizer can reveal an optimization failure, but cannot create information absent from the data. Simplifying random effects may be justified by a stable, unsupported variance direction; it is not a cure for every warning. The practical goal is therefore to find the simplest explanation supported by diagnostics, while preserving the intended model unless evidence warrants changing it.
Model and Assumptions Behind the Workflow
For clarity, consider observations \(Y_{ij}\) from unit \(i\) in cluster \(j\), with fixed-effect design vector \(x_{ij}\) and random-effect design vector \(z_{ij}\). A GLMM specifies a conditional distribution from an exponential-family model and a link \(g\):
Here \(\beta\) is the fixed-effect vector, \(b_j\) is the cluster-specific random-effect vector, and \(D\) is its covariance matrix. The usual likelihood-based workflow assumes that the conditional response distribution, link, fixed- and random-effect structures, and random-effect distribution are appropriate for the intended analysis. It also assumes clusters are independent and, conditional on the random effects and specified predictors, responses have the stipulated dependence structure. Troubleshooting numerical fit does not verify these substantive assumptions.
The likelihood optimized for estimation is the marginal likelihood, which integrates out the random effects. If \(f(y_j\mid b_j;\beta)\) is the conditional likelihood contribution for cluster \(j\), then:
A covariance matrix \(D\) must be positive semidefinite. Its variances cannot be negative; a variance of zero is an allowed boundary value. Numerical convergence means an algorithm has met its stopping criteria near a candidate optimum. It does not, by itself, establish that the optimum is well identified, that the covariance estimate is nonsingular, or that the model is appropriate.
A Decision Sequence for Fit Problems
Proceed from checks that leave the model unchanged toward decisions that may change its structure. Keep the data, model specification, optimizer settings, and diagnostic output for each fit so comparisons are interpretable.
Check the response, coding, units, missing-data handling, cluster identifiers, and model formula. Confirm that the intended observations and random-effects terms are actually used. A data or formula error is not an optimizer problem.
Inspect predictor ranges and the reported gradient and Hessian diagnostics. Centering and scaling continuous predictors can make optimization more stable. Refit the same substantive model after a transparent change of units, then compare fitted quantities on the original scale.
With the same data and model, use a reasonable alternative optimizer or starting values. Compare log likelihoods and estimates, not warning labels alone. Materially different solutions suggest the optimum has not been located reliably; similar solutions with small gradients and stable curvature are reassuring, but do not rule out a boundary or singular fit.
Look for flat likelihood directions, large uncertainty, unstable estimates under small perturbations, or highly correlated parameters. These indicate that the data may not distinguish some parameter combinations well. A more aggressive optimizer cannot resolve a genuine lack of information.
Check whether a variance is estimated near zero and whether the random-effects covariance has a very small eigenvalue or near-perfect random-effect correlation. These can indicate a boundary optimum or a redundant random-effects dimension. Use profile behavior and stability evidence to distinguish a supported zero variance from a poorly determined direction.
Retain the model if the fit is stable and the warning has a defensible numerical explanation. If a random-effects component is unsupported, consider a scientifically defensible simpler structure and refit; report the original issue, diagnostic evidence, change, and its consequences.
Why Rescaling Can Help Without Changing the Model
Suppose a continuous predictor \(x\) is replaced by \(x^*=(x-c)/s\), where \(s\ne0\). In a linear predictor with a fixed intercept and slope, \(\beta_0+\beta_1x\), set \(\beta_1^*=s\beta_1\) and \(\beta_0^*=\beta_0+c\beta_1\). Then:
Thus the same linear predictor, and therefore the same likelihood, can be represented using different units. In a random-slope model the corresponding random-effect coefficients and covariance must also be transformed consistently. Rescaling can improve the numerical coordinates presented to an optimizer, but it does not add information or repair a misspecified model. When interpreting results, transform coefficients back to the original units.
Worked Example: Rescaling a Predictor and Its Slope
Suppose a hypothetical model uses age in years, with an estimated fixed slope of \(0.020\) log-odds per year and standard error \(0.0067\). Age has mean 52 years and standard deviation 12 years. Define \(z=(\text{age}-52)/12\). The slope per one standard deviation is:
The arithmetic uses the fact that one unit of \(z\) is 12 years. Conversely, the slope per year is \(0.240/12=0.020\), and its standard error is \(0.0804/12=0.00670\). The estimated change in log-odds per year is unchanged. If the original fit had a large gradient and the scaled fit has a small gradient while giving equivalent likelihood and transformed estimates, scaling is a plausible numerical improvement. It would be incorrect to call the scaled coefficient a \(0.240\) change per year; its unit is one standard deviation of age.
Worked Troubleshooting Cases
Worked Example: Distinguishing an Optimizer Problem From a Stable Fit
A hypothetical Poisson GLMM for weekly event counts in 30 clinics includes treatment as a fixed effect and a random intercept for clinic. The same data and model are fit with two optimizers. Optimizer A reports log likelihood \(-428.31\), maximum absolute scaled gradient \(0.018\), and a convergence warning. Optimizer B reports log likelihood \(-428.30\), maximum absolute scaled gradient \(0.0012\), and no warning. Both estimate the treatment coefficient as \(-0.31\), to two decimal places.
The log-likelihood difference is \((-428.30)-(-428.31)=0.01\), and both estimates agree at the stated precision. Optimizer A's gradient is \(0.018/0.0012=15\) times as large as B's. This is evidence that A stopped less satisfactorily, while B found a solution with a smaller gradient. It is not proof that every parameter is well identified: the covariance diagnostics and sensitivity of estimates still need inspection. The immediate response is to retain the better-supported fit as a candidate and examine curvature and random-effect diagnostics, not to alter the model merely to silence A's warning.
Worked Example: Reading a Near-Singular Random-Slope Covariance
A hypothetical model uses a centered, standardized predictor and estimates the random-intercept and random-slope covariance matrix as:
For a matrix of this form, the eigenvalues are \(0.040+0.039=0.079\) and \(0.040-0.039=0.001\). Both are positive, but the smaller eigenvalue is only \(0.001/0.079=0.0127\), or about 1.27% of the larger. The estimated correlation is \(0.039/\sqrt{0.040(0.040)}=0.975\). Together, the very small eigenvalue and near-perfect correlation suggest that the fitted random effects occupy nearly a lower-dimensional space; this is evidence of near-singularity.
This calculation does not by itself prove that the random slope should be removed. Eigenvalues depend on the units used for random-effect predictors, so the predictor scaling must be understood; here it was standardized. Next compare the covariance across optimizer starts and inspect the likelihood as the slope variance approaches zero. If the small eigenvalue and fit are stable, the model may be estimating a boundary or nearly redundant direction. If they vary widely, weak identification is a stronger concern. A singularity warning is not the same as evidence that the fixed effects are uninterpretable, nor does it identify the best simpler model automatically.
Worked Example: Applying the Full Sequence to a Hypothetical Fit
Suppose a hypothetical study records counts for 24 clinics over four visits at days 0, 7, 14, and 21. The proposed Poisson GLMM has treatment and visit as fixed effects, plus clinic-specific random intercepts and random slopes for visit. In the original day units, a fit reports maximum scaled gradient \(0.014\), Hessian condition number \(2.4\times10^9\), and a warning.
Setup and scale. The records and clinic identifiers are checked, and visit is confirmed to be measured in days. Re-express visit in weeks and center it at the second visit. This is an invertible change of predictor units, so the same mean patterns remain available when coefficients and random-effect covariance are transformed consistently.
Optimizer check. In the re-expressed model, two optimizers give log likelihoods \(-312.846\) and \(-312.845\), and treatment estimates that differ by \(0.003\). The better fit has maximum scaled gradient \(0.0011\) and Hessian condition number \(3.1\times10^5\). The improvement in gradient and conditioning supports the conclusion that original units contributed to numerical difficulty. The close likelihoods and treatment estimates are reassuring about optimizer stability, but do not establish a full-rank random-effects fit.
Covariance check and decision. The standardized random-effects covariance has eigenvalues \(0.062\) and \(0.00002\). Their ratio is \(0.00002/0.062=0.000323\), so the smaller direction is extremely weak relative to the larger one. Alternative starts produce essentially the same likelihood and near-zero slope-variance direction; a profile check places its maximum at or very near zero. The evidence now supports a stable fit at a boundary with a near-singular random-effects structure, rather than an unresolved optimizer failure alone. A reasonable next analysis is to fit a scientifically defensible random-intercept model and compare estimates and uncertainty, while preserving the original fit and reporting the reason for the simpler model. The available evidence does not justify claiming that the random slope is exactly absent in the population.
How to Interpret Diagnostics and Choose an Alternative
The distinctions are practical. Scaling trouble is suggested when a harmless change of units substantially improves gradients or conditioning while back-transformed estimates and likelihood remain equivalent. Optimizer failure is suggested when starts or algorithms reach materially different likelihoods or estimates, or leave large gradients. Weak identifiability is suggested by flat likelihood directions and unstable estimates, even when an optimizer has converged. A boundary estimate is a parameter estimate at an allowed edge, such as a variance of zero. A singular fit means the estimated random-effects covariance is rank deficient or effectively so under the numerical diagnostic being used.
These categories can overlap. A weakly identified variance may end up at the boundary, and awkward scaling can make a near-singular covariance harder to diagnose. Conversely, a stable boundary estimate can be the actual likelihood maximum, not a numerical failure. Software-specific gradient tolerances and singularity thresholds are not universal statistical constants; report the diagnostics and software criterion used rather than presenting a threshold as a theorem.
The main alternatives have different purposes. Rescaling retains the model and changes numerical representation. Trying another optimizer tests whether the solution depends on the numerical algorithm. Profile likelihood and perturbation checks probe how strongly the data support a parameter direction. Simplifying random effects changes the model and may trade flexibility for estimability. These approaches are complementary, not interchangeable. Do not select a simpler model only because it produces a warning-free fit; connect any change to a diagnosed unsupported component and the scientific question.
Common Pitfalls and Reporting
- Treating a convergence message as the diagnosis. A warning identifies a condition reported by the fitting routine, not its cause. Examine gradients, curvature, alternative fits, and covariance structure.
- Assuming a successful optimizer proves identification. An algorithm can converge to a flat or boundary optimum. Convergence is necessary evidence about computation, not sufficient evidence about information.
- Dropping a random slope immediately. First check scaling, optimizer stability, and whether the boundary direction is stable. A simpler model may be justified, but it should answer the same scientific question as well as possible.
- Comparing coefficients without units. After scaling, state whether a coefficient is per day, per week, or per standard deviation, and transform it before comparing estimates on different scales.
- Calling a near-zero variance proof of no variation. A boundary estimate describes the fitted likelihood under the specified model and data. It does not establish exact population homogeneity.
A concise report should name the original warning, the rescaling or optimizer checks performed, the relevant likelihood and gradient comparisons, and any boundary or singularity evidence. If the random-effects structure changed, state what was removed and why, and indicate whether the substantive fixed-effect conclusions were sensitive to that decision.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- A predictor is recorded in millimeters and is converted to centimeters. What should remain invariant if the model is reparameterized correctly?
- Two optimizers give nearly equal likelihoods but a random-effects covariance has one eigenvalue near zero. Does the agreement prove the fit is nonsingular?
- A variance estimate is zero in two stable refits, and the profile likelihood is maximized at zero. What does this support, and what does it not prove?
- A GLMM warning disappears after a random slope is removed. Why is this alone insufficient justification for removing it?
Solutions. (1) The fitted linear predictors, likelihood, and predictions should agree; coefficients must be transformed to reflect the changed units. (2) No. Optimizer agreement supports numerical stability of the solution, but the near-zero eigenvalue is evidence of a singular or nearly singular covariance. (3) It supports a boundary maximum for that variance under the specified model and data; it does not prove the true population variance is exactly zero. (4) Removing a term changes the model and may hide rather than diagnose the problem. Check scaling, optimizer stability, identification, and covariance evidence, then justify simplification scientifically.