Convergence and Boundary Diagnostics

Why Convergence Needs More Than a Software Message

A shared-parameter model combines an outcome process and a dropout process through latent random effects. Its observed-data likelihood requires integration over those effects. Consequently, a numerical routine can stop because it has reached a stationary point, because it has encountered a difficult boundary, or simply because an iteration limit or numerical tolerance was reached. A message saying “converged” is useful, but it does not establish that the solution is stable, well identified, or insensitive to numerical integration.

This tutorial focuses on those diagnostics. It does not replace substantive assessment of the model’s fit. Instead, it asks whether the fitting procedure has located a credible solution to the specified likelihood, whether that solution lies near a problematic covariance boundary, whether the association estimate is numerically or inferentially unstable, and whether the integration approximation is adequate.

Setup: Write $\theta$ for outcome-process parameters, $\alpha$ for dropout-process and shared-association parameters, and $D$ for the random-effect covariance matrix. For a positive-definite $q\times q$ matrix $D$, the observed-data log-likelihood can be written $$ \ell(\theta,\alpha,D)=\sum_{i=1}^{n}\log\left\{\int L_i(\theta,\alpha\mid b_i)\,\phi_q(b_i;0,D)\,db_i\right\}, $$ where $L_i$ is the individual’s likelihood conditional on random effects $b_i$, and $\phi_q(\cdot;0,D)$ is the $q$-variate normal density. This density form assumes $D$ is positive definite. If a singular covariance is admitted as a boundary limit, the integral must instead be expressed as an expectation under the possibly degenerate Gaussian probability measure, or as an integral over its lower-dimensional support.

The numerical objective is usually an approximation to this likelihood, often based on adaptive quadrature. Let $\psi$ collect the free model parameters, and let $\ell_K(\psi)$ denote the approximation using a specified quadrature setting $K$. Diagnostics should therefore consider both the optimizer for $\ell_K$ and the adequacy of the approximation to $\ell$.

Interior Convergence and the Local Likelihood Surface

At an interior maximum of a differentiable log-likelihood, the score $\nabla\ell(\hat\psi)$ is zero. The Hessian $H(\hat\psi)=\nabla^2\ell(\hat\psi)$ must be negative semidefinite there. A negative-definite Hessian, with negative eigenvalues in every relevant direction, is a useful regular diagnostic and under usual smoothness conditions implies a strict local maximum. It is not a necessary condition for every strict maximum: a flat direction can have zero second derivative even at a strict maximum.

In practice, the computed score will rarely be exactly zero. Inspect a scaled score norm, not just the raw gradient, because parameters may have very different units. Also examine whether the Hessian is numerically stable and whether its eigenvalues suggest nearly flat directions. Near-singular curvature can make standard errors very large or unstable even if the optimizer reports success. At a boundary, the interior score condition need not hold; the relevant constrained optimum conditions must account for the allowed parameter space.

Worked Example: Assessing an Interior Convergence Check

Suppose a hypothetical fit has five parameters, transformed to comparable scales for diagnostics. At the reported solution, the maximum absolute scaled score component is $0.0008$, and the Hessian eigenvalues are $-12.0,-3.1,-0.40,-0.08,-0.02$. The estimate is interior to the parameter space.

The small scaled score is consistent with a stationary point: every component is close to zero on the chosen diagnostic scale. All five Hessian eigenvalues are negative, so the Hessian is negative definite at this estimate. This supports a regular local maximum, although it does not prove that the global maximum was found.

The eigenvalue closest to zero, $-0.02$, indicates a relatively flat direction compared with $-12.0$. Thus the fit can be locally maximizing while still having a weakly determined parameter combination. A responsible check would repeat optimization from dispersed starting values and compare solutions, then assess whether the weak curvature materially affects standard errors or key estimates. The score and eigenvalues alone do not certify stability.

A practical convergence assessment combines complementary evidence. Check the optimizer’s termination reason, the score or gradient norm, parameter changes across the final iterations, and the likelihood change. Repeat the fit from several reasonable starting values, including starts that vary the association parameters and random-effect scales. Similar parameter estimates and log-likelihoods across starts are reassuring; distinct solutions with similar likelihoods can indicate a flat or multimodal surface. When parameters are constrained, verify that reported values satisfy the constraints rather than relying on an unconstrained approximation near a boundary.

Random-Effect Covariance Near a Boundary

A random-effect covariance matrix must be positive semidefinite. A positive-definite estimate has strictly positive eigenvalues; a singular boundary estimate has at least one zero eigenvalue. In finite-precision calculations, a very small positive eigenvalue may signal that the fitted covariance is close to singularity. For a two-dimensional random effect, a correlation near $1$ or $-1$ is another warning: the two effects are then nearly linearly dependent under the fitted distribution.

These signs are diagnostic, not automatic proof of an error. A covariance boundary can arise because the data support little variation in one random-effect direction, because the model is weakly identified, or because the optimizer and integration are struggling. Standard errors based on an ordinary inverse Hessian can be misleading at or near a boundary; symmetric Wald intervals may extend into impossible covariance values. Report the boundary behavior and use an inference method appropriate to the constrained problem rather than treating a boundary estimate as a routine interior estimate.

Worked Example: Detecting a Nearly Singular Covariance

In a hypothetical fit, the estimated covariance matrix for a random intercept and random slope is $$ \hat D= \begin{pmatrix} 0.25 & 0.495\\ 0.495 & 1.00 \end{pmatrix}. $$ The intercept variance is $0.25$, the slope variance is $1.00$, and the covariance is $0.495$.

The implied correlation is $$ \hat\rho=\frac{0.495}{\sqrt{0.25(1.00)}}=\frac{0.495}{0.5}=0.990. $$ The determinant is $$ |\hat D|=(0.25)(1.00)-(0.495)^2=0.25-0.245025=0.004975. $$ The eigenvalues are approximately $1.246$ and $0.003993$. Both are positive, so this displayed matrix is positive definite, but the smaller eigenvalue is tiny relative to the larger one and the correlation is close to one. The estimate is therefore near-singular, not exactly singular.

This warrants checking sensitivity to starting values, integration settings, and reasonable covariance parameterizations. It also calls for caution with Hessian-based uncertainty summaries. It does not, by itself, justify deleting the random slope: that is a model decision requiring substantive consideration and a comparison that respects the boundary issue.

Extreme Association Estimates

In a shared-parameter model, an association coefficient describes how a specified latent effect enters the dropout process, conditional on the model’s other terms. Its numerical size depends on the latent-effect scale and link function. For a logistic dropout component, an association coefficient $\alpha$ multiplying a random effect measured in standard-deviation units corresponds to a conditional odds ratio $\exp(\alpha)$ for a one-unit increase in that effect, holding other terms fixed.

An extreme estimate may be plausible on a particular scale, but it can also flag sparse dropout information, near-separation in the dropout component, weak identification between the outcome and dropout parameters, or a poorly scaled latent effect. A large point estimate paired with a very large standard error is especially unstable. Wald intervals can be uninformative on the coefficient scale and asymmetric after exponentiation; inspect a profile likelihood or another appropriate uncertainty analysis rather than relying only on a Wald summary. The sign and size should always be interpreted using the model’s exact coding and scale.

Worked Example: Interpreting a Large Dropout Association

Suppose a hypothetical fitted dropout logit includes $\alpha b_i$, where $b_i$ is a standardized shared random effect, and gives $\hat\alpha=7.80$ with estimated standard error $2.40$. The conditional odds ratio for a one-standard-deviation increase in $b_i$ is $$ \exp(\hat\alpha)=\exp(7.80)\approx 2{,}441. $$

Thus the fitted model associates a one-standard-deviation increase in the latent effect with approximately 2,441 times the dropout odds, conditional on the other terms. This is a very large association, not a causal effect, and its interpretation depends on the model and latent-effect scale.

A rough Wald interval on the coefficient scale is $$ 7.80\mathbin{\pm}1.96(2.40)=7.80\mathbin{\pm}4.704=(3.096,12.504). $$ Exponentiating the endpoints gives an approximate odds-ratio interval of $(\exp(3.096),\exp(12.504))$, or about $(22.1,269{,}000)$. The enormous width shows that the point estimate alone gives a misleading impression of precision. Because the estimate is extreme and uncertainty is large, check convergence from alternative starts and examine a profile or other suitable interval before interpreting it substantively.

Quadrature Sensitivity

Adaptive quadrature approximates integration over the random effects. A fit can converge accurately to the maximum of $\ell_K$ while that approximation is inadequate for the desired inference. Refit with increased quadrature accuracy and compare the parameters of interest, covariance estimates, log-likelihood values, and conclusions. Stability is more informative than a single node count. The sequence need not change monotonically, so agreement between successive settings is a practical diagnostic rather than a mathematical guarantee.

The reason parameter estimates can shift is seen by a local expansion. Let $\hat\psi_K$ maximize $\ell_K$, and let $\hat\psi_{K'}$ be the fit under a changed quadrature setting. For an interior solution, expand the changed score around $\hat\psi_K$:

0\approx\nabla\ell_{K'}(\hat\psi_K) +H_{K'}(\hat\psi_K)(\hat\psi_{K'}-\hat\psi_K).

If $H_{K'}(\hat\psi_K)$ is invertible, rearrangement gives

\hat\psi_{K'}-\hat\psi_K \approx -H_{K'}(\hat\psi_K)^{-1}\nabla\ell_{K'}(\hat\psi_K).

This is a local approximation, not an exact identity. It shows why a small change in the objective’s score can move estimates substantially when curvature is weak. If the Hessians under the two settings are close, $H_K$ may approximate $H_{K'}$, but the changed score must still be evaluated under $\ell_{K'}$.

Worked Example: Comparing Quadrature Settings

A hypothetical analysis is refitted with increasing quadrature accuracy. The association estimate, treatment coefficient, and approximate log-likelihood are as follows.

Quadrature setting	Association estimate	Treatment coefficient	Log-likelihood
7 nodes	1.42	-0.31	-812.64
15 nodes	1.58	-0.30	-812.51
31 nodes	1.60	-0.30	-812.50

From 7 to 31 nodes, the association estimate changes by $1.60-1.42=0.18$; from 15 to 31 it changes by $1.60-1.58=0.02$. The treatment coefficient changes by only $0.01$ from 7 to 15 nodes and does not change to the reported precision from 15 to 31. The log-likelihood changes by $0.13$ from 7 to 15 and by $0.01$ from 15 to 31.

These comparisons suggest that the 7-node result is less stable for the association estimate, while the 15- and 31-node fits are close at the reported precision. They do not prove that 31 nodes is exact. The analyst should also check whether the key uncertainty summaries and substantive conclusions are stable at the higher settings.

A Diagnostic Sequence and Its Limits

Verify termination.
Record the optimizer’s stopping reason, final score or gradient, and changes in parameters and objective. A hard iteration limit is not evidence of convergence.

Check repeatability.
Refit from dispersed, reasonable starting values. Compare log-likelihoods and estimates, paying special attention to association and covariance parameters.

Inspect curvature and boundaries.
Review Hessian eigenvalues and covariance eigenvalues or correlations. Treat nearly flat directions and near-singular covariance estimates as warnings about uncertainty and identification.

Increase integration accuracy.
Repeat the fit with more accurate quadrature and compare important estimates, uncertainty summaries, and conclusions.

Report unresolved instability.
If solutions remain sensitive, state which quantities are unstable and avoid presenting a single numerical estimate as definitive.

These checks answer different questions. Multiple starts test optimizer repeatability; they do not test whether the latent-effect distribution or link is correctly specified. Quadrature comparisons assess numerical integration sensitivity; they do not establish model adequacy. Those broader questions require separate model assessment. Likewise, EM and direct likelihood optimization can behave differently computationally, and Bayesian estimation offers a different inferential framework, but neither automatically removes weak identification or sensitivity to modeling assumptions. The diagnostic priority is to identify what is unstable before choosing a remedy.

Common mistake: “The software converged, so the covariance and association estimates are reliable.” This reasoning conflates an algorithm’s stopping rule with the quality of the solution. A fit can satisfy a numerical tolerance while lying near a singular covariance boundary, having nearly flat likelihood curvature, or changing materially with quadrature accuracy. The correction is to report convergence evidence alongside boundary, curvature, repeatability, and integration checks.

Key takeaway: Convergence is a computational statement about a fitting procedure, not a guarantee of precise or robust inference. Check score behavior and repeatability, inspect covariance and likelihood curvature, scrutinize extreme association estimates, and compare results across quadrature settings.

Check Your Understanding

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

At an interior maximum, what second-order condition must the Hessian satisfy? Is negative definiteness necessary?
A two-random-effect covariance estimate has correlation $0.999$ and a small positive determinant. What does this indicate, and what does it not establish?
For a logistic dropout model with a standardized latent effect and association estimate $\alpha=2$, what is the conditional odds ratio for a one-unit increase in that effect?
Why should quadrature sensitivity be assessed by refitting under the changed quadrature setting rather than only checking that the original fit converged?
Give two distinct causes of a large association estimate that should prompt further diagnosis.

Solutions. (1) The Hessian must be negative semidefinite at an interior maximum. Negative definiteness is a useful regular diagnostic and sufficient for a strict local maximum under usual smoothness conditions, but it is not necessary. (2) The estimate is near-singular and its covariance structure may be weakly determined; this alone does not prove that the model is wrong or that the optimizer failed. (3) The odds ratio is $\exp(2)\approx7.39$, conditional on the other model terms. (4) Changing quadrature changes the approximate likelihood and its score, so a previous optimizer’s convergence says nothing by itself about the new approximation’s fitted solution or resulting parameter stability. (5) Examples include weak identification, sparse dropout information or near-separation, and poor scaling of the latent effect.