Why a Convergence Report Needs More Than a Status Message
A software message saying “converged” is not enough to establish that a joint-model estimate is numerically reliable or well identified. A useful report lets another analyst understand what was optimized, what numerical approximations were used, which checks were passed, and whether plausible parameter values remain weakly distinguished by the data. These are related but different questions: an optimizer can stop successfully at a numerically stable point while the likelihood is nearly flat in one parameter direction.
For a generic shared-random-effects joint model, let \(y_i\) denote subject \(i\)'s longitudinal measurements, \(T_i\) the observed event or censoring time, and \(\delta_i\) the event indicator. Let \(b_i\) denote the subject-specific random effects and \(\theta\) the vector of model parameters. The observed-data likelihood combines the longitudinal and event contributions and integrates over the unobserved random effects. Numerical optimization and numerical integration can therefore both affect the reported fit.
What a Reproducible Report Specifies
Report the software and version, estimation method, optimizer and any important options, parameterization, starting-value strategy, and the criteria used to stop. For iterative methods, give the relevant iteration limit and tolerances, and say whether the algorithm stopped because it met a criterion or because it reached a limit. If several starts were tried, report how many, how they were generated, and whether the solutions agreed. Reporting only the best objective value hides possible dependence on initialization.
State how the likelihood was evaluated, including the integration method and settings, and identify any approximation used. In a shared-random-effects model, these details matter because the likelihood requires integration over \(b_i\). A fitted result should be associated with the settings that produced it; integration accuracy should not be inferred from an optimizer's status.
A concise report should also identify the parameter scale used for derivative diagnostics, the Hessian convention, and how uncertainty estimates were obtained. The Hessian of the log likelihood at a maximum is expected to be negative definite in an identifiable, regular interior case; equivalently, the observed information, the negative Hessian, is positive definite. These expectations require care with boundaries, redundant parameterizations, and numerical approximation.
Worked Example: Reporting Optimizer and Derivative Checks
Suppose a hypothetical joint model was fitted by maximizing the approximate log likelihood. The report states that software version 4.2 used a quasi-Newton optimizer, three prespecified starting values, and a maximum of 2,000 iterations. The stopping tolerances were \(10^{-8}\) for the relative objective change and \(10^{-6}\) for the scaled gradient. All three starts reached log likelihood \(-500.000\), with maximum absolute differences of \(0.0003\) across the resulting parameter estimates.
To make the gradient check interpretable across differently scaled parameters, define \(u_j=(\theta_j-\hat{\theta}_j)/s_j\), where \(s_j\) is a stated scale for parameter \(j\). Then the gradient in these coordinates is \(\partial\ell/\partial u_j=s_j\partial\ell/\partial\theta_j\). Suppose the largest absolute component of this scaled gradient at the reported estimate is \(2.1\times10^{-7}\), below the stated tolerance \(10^{-6}\). The report can say that the scaled-gradient criterion was met; it should not claim that the gradient is exactly zero.
Suppose the observed information, evaluated in the same scaled coordinates, has eigenvalues \(0.42, 1.7, 5.3,\) and \(12.0\). All are positive, supporting positive curvature of the negative log likelihood in these coordinates at the estimate. The condition number is \(12.0/0.42=28.6\). This is a description of local curvature, not a universal pass/fail threshold. The report should identify the scaling because eigenvalues and condition numbers depend on parameterization.
Checking Integration Sensitivity
A fit can satisfy its optimizer's stopping rule for a likelihood that is still materially affected by numerical integration error. Integration sensitivity is assessed by refitting the same model under increasingly accurate or otherwise varied integration settings, while holding the data, model, and other estimation choices fixed. Record the settings and compare estimates and uncertainty summaries that matter to the scientific interpretation, as well as the log likelihood when it is comparable across the settings.
For example, with Gaussian quadrature, one might increase the number of nodes. A stable estimate across those runs supports numerical robustness to that particular change in integration settings. It does not establish that the model is correctly specified, that all parameters are identified, or that the quadrature error is negligible in every direction. If estimates change appreciably, report the changes and investigate rather than selecting whichever run gives the most favorable result.
Worked Example: Comparing Integration Settings
In a hypothetical shared-random-effects fit, the association parameter \(\eta\) is a log hazard ratio per one standard deviation increase in the latent marker. The model is refitted with 9, 15, and 25 quadrature nodes per integration dimension, using the same optimizer settings and starting values. The results are:
| Nodes | Estimate of \(\eta\) | SE | Maximized log likelihood |
|---|---|---|---|
| 9 | 0.238 | 0.112 | -500.014 |
| 15 | 0.241 | 0.113 | -500.003 |
| 25 | 0.241 | 0.113 | -500.000 |
From 9 to 25 nodes, the estimate changes by \(0.241-0.238=0.003\), about \(0.003/0.241=1.24\%\) of the 25-node estimate. The standard error changes by \(0.001\), and the log likelihood increases by \(-500.000-(-500.014)=0.014\). These small changes support stability for these reported quantities across the tested settings. A suitable report gives the table and states this bounded conclusion; it does not call the quadrature error zero or claim stability for unreported quantities.
The hazard ratio corresponding to the 25-node estimate is \(\exp(0.241)=1.273\) per standard deviation increase in the latent marker. This transformation is included for interpretation, not as a substitute for the integration comparison. If the estimate or its uncertainty had shifted materially as nodes increased, the coarser fit would not provide adequate evidence of numerical stability.
Profile Evidence and Local Curvature
A Hessian summarizes local curvature near the estimate, so it can be useful for approximate standard errors but may not show the shape of the likelihood farther away. A profile likelihood gives a parameter-specific view. For a scalar parameter \(\eta\), its profile log likelihood is the maximum log likelihood over all other parameters while \(\eta\) is held fixed:
Comparing \(\ell_p(\eta)\) with its maximum shows how much fit is lost when the parameter is constrained. Under regularity conditions, an approximate likelihood-ratio confidence set for one parameter contains values satisfying \(2\{\ell_p(\hat{\eta})-\ell_p(\eta)\}\leq\chi^2_{1,0.95}\), where \(\chi^2_{1,0.95}\) is about 3.84. This approximation may be unreliable with weak identification, boundary parameters, or other violations of regularity; the report should state the method and its limitations.
Worked Example: Reporting a Profile Likelihood
Suppose a hypothetical fit has \(\hat{\eta}=0.24\) and maximized log likelihood \(-500.00\). When \(\eta\) is fixed at zero, re-optimizing the other parameters gives \(\ell_p(0)=-502.10\). The likelihood-ratio statistic for \(H_0:\eta=0\) is
Using the one-degree-of-freedom chi-square approximation, a statistic of 4.20 gives \(p\approx0.0404\). At the 0.05 level, this is evidence against \(H_0:\eta=0\) under the fitted model and the stated approximation. It is not the probability that the null hypothesis is true.
Suppose the profile log likelihood is \(-501.92\) at \(\eta=0.08\) and at \(\eta=0.40\), and the profile remains above that value between those points. The log-likelihood drop at either endpoint is \(1.92\), so the likelihood-ratio statistic is \(2(1.92)=3.84\), approximately the 95% cutoff. The resulting approximate profile interval is \([0.08,0.40]\). In hazard-ratio units, exponentiating the endpoints gives \([\exp(0.08),\exp(0.40)]=[1.083,1.492]\) per standard deviation increase in the latent marker. Under the model and approximation, this interval describes values not excluded by the 5% likelihood-ratio procedure; it does not mean there is a 95% probability that the fixed parameter lies in the interval.
What to Conclude—and What Not to Conclude
The evidence should be presented as a set of complementary checks. Optimizer termination and repeated starts address whether the reported solution was reached reliably under the chosen algorithm. Gradient and Hessian checks describe local first- and second-order behavior. Integration comparisons address sensitivity to a numerical approximation. A profile reveals parameter-specific likelihood shape. Agreement across these checks strengthens a numerical report, but none alone establishes global uniqueness or substantive adequacy.
A useful report distinguishes a parameter with a stable estimate and well-behaved profile from one whose profile is shallow, asymmetric, or interrupted by numerical failures. State the range examined, whether nuisance parameters were re-optimized at each fixed value, and how any interval was constructed. If profile calculations fail or multiple branches appear, disclose that. Do not present a symmetric standard-error interval as definitive when the profile indicates strong asymmetry or flatness.
There are alternatives, but they answer different questions. A Hessian-based standard error is efficient and compact near a regular maximum, while profiling can reveal non-quadratic shape for selected parameters at greater computational cost. Repeated starts help expose initialization sensitivity but do not prove a global maximum. Bootstrap variability can describe sampling instability, though it is not a replacement for checking numerical accuracy and may itself be unstable. Bayesian convergence diagnostics concern sampling from a posterior distribution, not the same optimization criteria; posterior summaries also depend on the prior.
Identifiability conclusions must remain conditional on the model, observed data, and estimation approach. A report can say that the association estimate is sensitive to integration settings, has a broad profile, or varies across starts. It should not convert one such symptom into an unsupported claim that a parameter is definitively unidentifiable. Likewise, a small gradient and positive-definite observed information do not rule out distant alternative maxima or weak information beyond the local neighborhood.
A Compact Reporting Workflow
Name the model, data analyzed, software and version, estimation method, and numerical integration method and settings.
Give optimizer options, starting values or start strategy, iteration limit, stopping criteria, and the actual termination status.
State the scaled-gradient result and Hessian or observed-information diagnostic, including parameter scaling and any irregularities.
Compare meaningful estimates, uncertainty measures, and likelihood values across stated integration settings and, where relevant, starting values.
Give profile results for key parameters and clearly separate observed evidence from conclusions the diagnostics cannot establish.
Common Reporting Pitfalls
- Treating the software flag as proof. “Converged” means the program met its specified stopping rule; it does not establish integration accuracy, global optimality, or identification.
- Reporting a gradient norm without its scale. The numerical size depends on parameter units. State how parameters were scaled and which norm or component summary was used.
- Calling a positive-definite Hessian conclusive. It supports a local curvature check at the estimate, but does not exclude distant modes or weak identification visible in a profile.
- Reporting only one integration setting. An optimizer can converge accurately for an inadequately approximated likelihood. Compare settings and report the comparison.
- Calling an interval a probability statement about a fixed parameter. Describe the procedure and its coverage interpretation under its assumptions, not a posterior probability unless a Bayesian analysis supplies one.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Why should a convergence report give the parameter scaling used for its gradient diagnostic?
- An estimate changes little across two integration settings. What conclusion is supported, and what conclusion is not?
- What does a positive-definite observed information matrix at the estimate support?
- A one-parameter profile gives a likelihood-ratio statistic of 2.7 for the null value. Is that value excluded by the approximate 5% test using a cutoff of 3.84?
- Name two details needed to make a reported profile-likelihood interval interpretable and reproducible.
Solutions. (1) Gradient magnitudes depend on parameter units and scaling, so the scaling is needed to interpret the diagnostic. (2) It supports stability of the quantities actually compared under those two settings; it does not prove zero integration error, correct model specification, or identifiability. (3) It supports positive local curvature of the negative log likelihood in the stated parameterization, but does not prove global uniqueness or rule out a flat profile farther away. (4) No. Since \(2.7<3.84\), the null value is not excluded by that approximate 5% likelihood-ratio test. (5) State how nuisance parameters were re-optimized at fixed parameter values and how the interval cutoff or approximation was chosen; also report the likelihood method and relevant regularity limitations.