Setting Quadrature and Optimizer Controls for Maximum Likelihood

Why Numerical Controls Matter

Maximum-likelihood fitting of a joint model requires more than specifying the longitudinal and event submodels. Their shared random effects must be integrated out of each subject’s likelihood contribution, and that integration is generally numerical. An optimizer then searches for parameter values that maximize the resulting likelihood. Both stages introduce choices: integration accuracy affects the likelihood being optimized, while starting values and stopping controls affect how the optimizer searches it.

These controls are not cosmetic. Too few quadrature nodes can give a poor approximation to the likelihood; an overly strict tolerance can spend substantial time refining changes that are negligible for the analysis; and a limited iteration budget can stop an otherwise reasonable optimization before it has converged. Conversely, a program reporting convergence does not establish that the numerical integral was accurate or that the optimizer found the best relevant solution. A practical fit therefore checks sensitivity to key settings rather than relying on one successful software message.

The Likelihood Integral and Quadrature

Definition: For subject $i$, let $b_i$ denote the vector of shared random effects, with $b_i\sim N(0,D)$. Conditional on $b_i$ and covariates, suppose the longitudinal observations and event outcome contribute likelihood $L_i(b_i;\theta)$. The subject’s marginal likelihood contribution is $$ L_i(\theta)=\int L_i(b_i;\theta)\,\phi(b_i;0,D)\,db_i, $$ where $\theta$ collects the model parameters and $\phi(\cdot;0,D)$ is the multivariate normal density. This construction assumes independent subjects, the specified random-effects distribution, and conditional independence of the longitudinal and event outcomes given the shared random effects and covariates. The event contribution must also account for censoring under the model’s censoring assumptions.

The integral is often not available in closed form. Gaussian quadrature approximates an integral by a weighted sum of function evaluations. For a $d$-dimensional random-effects vector, write $D=LL^\mathsf{T}$, where $L$ is a Cholesky factor. Using $Q$ one-dimensional Gauss–Hermite nodes $x_q$ and weights $w_q$, the tensor-product approximation is

L_i(\theta)\approx \pi^{-d/2} \sum_{q_1=1}^{Q}\cdots\sum_{q_d=1}^{Q} \left(\prod_{j=1}^{d} w_{q_j}\right) L_i\!\left(\sqrt{2}Lx_q;\theta\right).

Here $x_q$ denotes the vector of selected nodes. The $\sqrt{2}$ and $\pi^{-d/2}$ factors correspond to the standard physicists’ Gauss–Hermite convention, which approximates integrals weighted by $e^{-x^2}$. Software may use a different internal convention, so the user should follow its documentation rather than manually combining nodes and weights from incompatible definitions.

The node count $Q$ controls the resolution of the approximation. Increasing $Q$ generally improves accuracy when the integrand is adequately represented by the quadrature rule, but accuracy is not guaranteed to improve monotonically for every parameter value. In a tensor-product rule, the number of node combinations is $Q^d$. Thus a modest increase in nodes can become expensive when the random-effects dimension $d$ is larger.

Worked Example: Checking Sensitivity to Quadrature Nodes

Suppose a hypothetical joint-model fit has two random effects. The analyst fits the same model using fixed Gauss–Hermite rules with 9, 15, 25, and 35 nodes per dimension. The software returns the following log-likelihoods and estimated association parameters:

Nodes per dimension $Q$	Node combinations $Q^2$	Log-likelihood	Association estimate
9	81	-312.846	0.412
15	225	-312.731	0.438
25	625	-312.728	0.440
35	1,225	-312.728	0.440

These values are hypothetical and illustrate a sensitivity check, not a universal pattern. From 15 to 25 nodes, the absolute log-likelihood change is $$ \left|-312.728-(-312.731)\right|=0.003. $$ From 25 to 35 nodes, the displayed change is $0.000$ at three decimal places, and the association estimate is unchanged at three decimal places. In this example, the 25-node fit appears adequate for these reported quantities, subject to checking other important estimates and standard errors.

The computational cost has increased substantially: $625/81\approx7.716$ times as many node combinations are evaluated per subject at 25 nodes as at 9. This ratio does not include implementation details such as vectorization or adaptive calculations, but it shows why choosing the largest available setting automatically is inefficient. The 9-node result, meanwhile, differs in both log-likelihood and association estimate, so treating it as adequate merely because optimization completed would be hard to justify.

Starting Values and Optimizer Controls

An optimizer begins from a starting parameter vector and iteratively updates it to improve the objective, usually the log-likelihood. Joint-model likelihoods can have complicated curvature and parameters on different scales. Starting values therefore matter: a poor initial point can slow optimization, produce unstable numerical derivatives, or lead to a different local solution. Reasonable initial values may come from simpler fits to the longitudinal and event data, provided they are transformed into the joint model’s parameterization and satisfy parameter constraints.

Controls commonly include a maximum number of iterations and one or more tolerances. A gradient tolerance assesses whether the objective is locally stationary; a parameter or step tolerance assesses whether successive updates are small; and an objective tolerance assesses whether the change in objective is small. These are not interchangeable. A small parameter update can occur because the optimizer has stalled, rather than because it has found a useful maximum. Tolerances also depend on software scaling and algorithm details, so their numerical values should not be treated as universal standards.

Practical principle: State the optimizer and its relevant controls, choose starting values that satisfy the model constraints, and allow enough iterations for the stopping criteria to be assessed. Then check whether the fitted estimates and likelihood are stable under reasonable changes to starting values and integration accuracy. Convergence is evidence about the numerical search, not proof that the model or approximation is correct.

Worked Example: Separating Convergence From an Iteration Limit

Consider a hypothetical fit using 25 quadrature nodes per dimension, an analytic gradient, a gradient tolerance of $10^{-5}$, and a maximum of 200 iterations. From an initial log-likelihood of $-338.4$, the optimizer reports a final log-likelihood of $-312.728$ after 84 iterations, with a scaled gradient norm of $3.2\times10^{-6}$.

The gradient criterion is met because $$ 3.2\times10^{-6}<10^{-5}. $$ The optimizer stopped before the iteration limit because $84<200$, and its reported gradient norm is below the chosen threshold. This supports the claim that the selected stopping rule was met; it does not by itself establish quadrature accuracy or a global maximum.

Suppose a second reasonable starting vector gives a final log-likelihood of $-312.727$ and association estimate 0.441, compared with $-312.728$ and 0.440 from the first start. The log-likelihood difference is $$ \left|-312.727-(-312.728)\right|=0.001. $$ The close results support numerical stability for these quantities in this example. They are not identical, and the analyst should consider whether differences at this scale matter for the intended interpretation.

Now suppose instead that the program stops at iteration 200 with a scaled gradient norm of $0.004$. Since $0.004>10^{-5}$, the gradient criterion has not been met. Reaching the iteration limit is not equivalent to convergence. The appropriate response is to investigate the starting values, scaling, parameterization, and numerical integration, and then refit with a justified control change—not simply to report the returned estimates as converged.

A Practical Tuning Sequence

Set a defensible baseline.
Use the intended model, a documented quadrature rule, sensible initial values, and tolerances and iteration limits appropriate to the optimizer. Keep a record of these settings.

Fit and inspect the stopping information.
Record the objective value, estimates, iteration count, stopping reason, and available gradient or step measures. Distinguish a tolerance-based stop from an iteration-limit stop.

Check the integration approximation.
Refit with more quadrature nodes, holding the model and other controls fixed as far as possible. Compare the log-likelihood, association estimates, variance estimates, and standard errors.

Check sensitivity to initialization.
Use more than one defensible starting vector when feasible. Materially different final values call for further investigation rather than selection of whichever fit is preferred.

Balance precision and cost.
Choose settings that give stable results for the quantities of interest at a practicable computational cost, and report the settings and sensitivity assessment.

Computational Trade-Offs and Alternatives

Fixed Gauss–Hermite quadrature uses nodes and weights fixed by the chosen rule. It is straightforward to compare across node counts, but may need many nodes when the integrand is concentrated away from the center of the random-effects distribution. Adaptive Gauss–Hermite quadrature shifts and scales the integration rule using information about the integrand, often its location and curvature near a mode. This can improve efficiency when that approximation is suitable, but adds computation and algorithm-specific choices. Fixed and adaptive results should not be assumed interchangeable without checking their agreement.

A Laplace approximation is another likelihood-integration approach: it approximates the integral using the integrand near its mode and local curvature. It can be computationally efficient, especially in higher dimensions, but its approximation may be inadequate when the relevant integrand is not well represented by a local quadratic form. These approaches trade computational cost against approximation behavior; none removes the need to assess whether the chosen calculation is adequate for the fitted model.

Dimensionality can dominate the decision. For three random effects, increasing a fixed tensor rule from 9 to 15 nodes changes the node-combination count from $$ 9^3=729 \quad\text{to}\quad 15^3=3{,}375, $$ a factor of $3{,}375/729\approx4.630$. With 200 subjects, that corresponds to 145,800 versus 675,000 node-level subject contributions for one full evaluation, before counting optimizer iterations and other computations. Real implementations may reuse calculations or use different integration schemes, but the arithmetic explains why adding random effects can make a previously routine setting costly.

Common Pitfalls and Limitations

Increasing the iteration limit instead of addressing the cause. More iterations may help if the optimizer is progressing but has not finished. They do not automatically fix poor scaling, unsuitable starting values, or an inaccurate integral.
Calling any stopping message convergence. A program can stop because it reached its iteration limit. Check the stopping reason and the relevant numerical criteria.
Changing several settings at once. If node count, starting values, and tolerances all change together, it becomes difficult to identify which change explains a different result. Vary one factor at a time where practical.
Assuming higher node counts always settle the question. Compare estimates and objective values across settings, and consider the software’s integration method and the model’s dimensionality. Rounded output can hide small differences, so retain adequate numerical precision in the comparison.
Treating numerical stability as model validation. Similar results across settings support computational robustness, but do not establish that the model assumptions are appropriate or that the estimated association is causal.

The appropriate controls depend on the optimizer, software implementation, model dimension, and precision required. A documented sensitivity assessment is more informative than claiming that one particular node count or tolerance is universally correct.

Worked Example: Choosing a Feasible Rule in Three Dimensions

A hypothetical analysis has three random effects and 200 subjects. The analyst is considering 9 or 15 fixed quadrature nodes per dimension. The number of node combinations per subject is $Q^d$, so the totals across subjects for one objective evaluation are:

For 9 nodes, $9^3=729$ combinations per subject, giving $200\times729=145{,}800$ subject-node contributions. For 15 nodes, $15^3=3{,}375$ combinations per subject, giving $200\times3{,}375=675{,}000$ contributions. The latter requires $675{,}000/145{,}800\approx4.630$ times as many contributions per evaluation.

Suppose the lower-cost fit gives an association estimate of 0.31 and log-likelihood $-842.16$, while the higher-cost fit gives 0.32 and $-842.12$. The likelihood difference is $\left|-842.12-(-842.16)\right|=0.04$. These figures alone do not prove that 15 nodes are necessary or that 9 are sufficient. The analyst should decide whether the estimate difference is material for the scientific question, consider a further node-count check if feasible, and weigh that evidence against the nearly 4.63-fold per-evaluation cost. The comparison is meaningful only if the same model and suitable optimization conditions were used for both fits.

Key takeaway: Quadrature settings determine how accurately the marginal likelihood is approximated; starting values and optimizer controls govern the search over that likelihood. Assess both parts, and report enough computational detail for the analysis to be understood and reproduced.

Check Your Understanding

Answer each question, then check the brief solutions below.

With four random effects and 11 quadrature nodes per dimension, how many tensor-product node combinations are evaluated per subject for one objective calculation?
A fit reaches its maximum of 150 iterations with a gradient norm of $0.002$, while the chosen gradient tolerance is $10^{-4}$. Has that gradient criterion been met?
Why can increasing the node count be especially costly when the number of random effects increases?
What does close agreement between two starting values establish, and what does it not establish?

Solutions. (1) $11^4=14{,}641$ combinations. (2) No: $0.002>10^{-4}$, and the run reached its iteration limit, so the gradient criterion was not met. (3) A tensor-product rule uses $Q^d$ combinations, so cost grows exponentially in the random-effects dimension $d$ for fixed $Q$. (4) Agreement supports stability of the numerical search for the compared fits; it does not establish quadrature accuracy, model adequacy, or that a global maximum was found.