Why Numerical Integration Is Part of the Model Fit
In a frailty model, observations from the same cluster are conditionally related through a latent random effect. For a lognormal frailty, integrating that effect out generally does not produce a closed-form likelihood. Estimation therefore requires a numerical approximation to the integrated likelihood, followed by numerical optimization. Both steps matter: a well-converged optimizer can maximize a poor approximation, while an accurate integral can still be paired with an optimizer that has not found the maximum.
The introductory material on likelihood integration for lognormal frailty motivates the integral. Here the focus is how adaptive Gaussian quadrature evaluates it and how to check that the resulting estimates are numerically trustworthy. The method is an estimation procedure, not a change to the underlying frailty model.
The Integrated Likelihood
Let \(\theta\) collect the fixed-effect parameters, baseline-hazard parameters if they are estimated, and the frailty standard deviation \(\sigma\). Write \(\ell_i(b;\theta)\) for the conditional log likelihood of all observations in cluster \(i\), given \(b_i=b\). Integrating over the Gaussian distribution gives the marginal cluster likelihood
where \(\phi\) is the standard normal density. The product of cluster likelihoods, or equivalently the sum of their log likelihoods, follows from independence between clusters. Dependence among observations within a cluster is retained through the shared \(b_i\). The estimate is obtained by maximizing \(\ell_{\mathrm{marg}}(\theta)\), not by maximizing each conditional likelihood over a separate \(b_i\) and treating those fitted values as ordinary parameters.
For quadrature, include the frailty density in the integrand by defining \(h_i(b;\theta)=\ell_i(b;\theta)+\log\{\phi(b/\sigma)/\sigma\}\). The integral is then \(\int \exp\{h_i(b;\theta)\}\,db\). Adaptive quadrature locates the mode \(m_i\) of \(h_i\) and scales the integration region according to its curvature there. Thus the quadrature is adapted to the cluster’s likelihood contribution and the frailty distribution together.
Adaptive Gaussian Quadrature
Let \(q_i=-h_i''(m_i;\theta)>0\), the negative second derivative at the mode. This positive-curvature condition means that the mode is locally a maximum with a non-flat quadratic approximation. Let \(x_j,w_j\) be nodes and weights for physicists’ Gauss–Hermite quadrature, which approximates integrals weighted by \(\exp(-x^2)\). The change of variable \(b=m_i+\sqrt{2/q_i}\,x\) gives
The factor \(\exp(x_j^2)\) in the sum compensates for the \(\exp(-x^2)\) weight built into the Gauss–Hermite rule. The mode-centred scale places nodes in the region where the integrand is substantial; the curvature factor adjusts their spread. \(K\) is the number of quadrature nodes. With \(K\) increasing, the approximation typically improves, but a finite-node result is not automatically exact.
For a vector of \(d\) Gaussian random effects, the same idea uses the mode vector and the negative Hessian matrix. A factorization of that matrix supplies a local scale and orientation for the quadrature nodes. The number of nodes needed can grow rapidly with \(d\), so the scalar case is substantially easier to compute than high-dimensional adaptive quadrature.
Worked Example: Evaluating One Adaptive Quadrature Approximation
Consider a hypothetical cluster whose contribution to the log-integrand has mode \(m_i=0.20\), height \(h_i(m_i)=-2.00\), and curvature \(q_i=4.00\). Suppose a two-node Gauss–Hermite rule uses \(x_1=-1/\sqrt{2}\), \(x_2=1/\sqrt{2}\), and \(w_1=w_2=\sqrt{\pi}/2\). At the transformed locations \(b=0.20\pm0.50\), suppose evaluation of the specified cluster likelihood and Gaussian density gives log-integrand differences of \(-0.25\) and \(-0.35\), respectively.
The transformed locations follow from \(\sqrt{2/q_i}=\sqrt{2/4}=0.7071\) and \(0.7071(1/\sqrt{2})=0.5000\). For either node, \(x_j^2=0.5000\), so the two residual exponents in the sum are \(-0.25+0.50=0.25\) and \(-0.35+0.50=0.15\). Substitution gives
This is a two-node approximation to this cluster’s integrated likelihood, not an exact value. The calculation also shows why the adjustment \(+x_j^2\) is required: it compensates for the weighting convention used by the quadrature rule.
Maximization and a Numerical Accuracy Check
At every trial value of \(\theta\), the fitting routine must evaluate the cluster integrals, sum their log contributions, and optimize that approximate marginal log likelihood. The cluster modes and curvatures generally change as \(\theta\) changes. A practical implementation therefore recomputes or updates them during optimization. The frailty standard deviation is constrained to be nonnegative; implementations often optimize a transformed parameter so that trial values respect this constraint.
A local maximizer reported by software is not, by itself, evidence of a reliable fit. Numerical accuracy is assessed by repeating the fit with more nodes and comparing the maximized log likelihood and parameter estimates. The optimizer should also be checked from different starting values, with its gradient and convergence status inspected. These checks address different problems: changing node count probes integration error, whereas starting values and gradient checks probe optimization.
Worked Example: Comparing Fits at Increasing Node Counts
Suppose a hypothetical recurrent-event data set has 40 independent clusters, with cluster event counts ranging from 0 to 5 and observed follow-up between 0.8 and 2.0 years. A lognormal shared-frailty model is fitted using adaptive quadrature. The following are hypothetical software results; the log likelihood is the maximized approximate marginal log likelihood.
| Nodes \(K\) | Maximized log likelihood | Estimated log rate ratio | Estimated frailty SD \(\hat{\sigma}\) |
|---|---|---|---|
| 5 | -82.417 | 0.271 | 0.710 |
| 9 | -82.391 | 0.263 | 0.760 |
| 15 | -82.390 | 0.262 | 0.760 |
| 25 | -82.390 | 0.262 | 0.760 |
From 9 to 15 nodes, the log likelihood changes by \((-82.390)-(-82.391)=0.001\), the log rate ratio changes by \(0.262-0.263=-0.001\), and the frailty SD changes by \(0.760-0.760=0.000\) to the displayed precision. The 15- and 25-node results also agree to the displayed precision. This is evidence that the reported fit is not materially sensitive to further increases over this range; it does not prove exact integration. The log rate ratio estimate of 0.262 corresponds to an estimated conditional rate ratio of \(\exp(0.262)=1.300\), holding the frailty fixed. These estimates describe this hypothetical model fit, not a causal effect.
Derivation and Interpretation of the Local Scaling
The adaptive transformation follows from a local quadratic expansion of the log-integrand. At a mode \(m_i\), the first derivative is zero, so near that mode
Set \(b=m_i+\sqrt{2/q_i}\,x\). Then \(db=\sqrt{2/q_i}\,dx\), and the quadratic term becomes \(-x^2\). To express the original integral in the form used by Gauss–Hermite quadrature, write
Approximating the remaining integral with the nodes and weights produces the adaptive formula above. If the log-integrand were exactly quadratic, this local scaling would align the integration with its Gaussian shape. In general it is only locally quadratic, so node-count checks remain necessary, especially when the integrand is skewed or has substantial tails.
Worked Example: Distinguishing Integration Error From Optimizer Failure
For a second hypothetical data set of 30 clusters, a fit using 11 nodes gives \(\hat{\beta}=0.18\), \(\hat{\sigma}=0.52\), and maximized log likelihood \(-61.204\). Increasing to 19 nodes gives \(0.18\), \(0.52\), and \(-61.203\), respectively. At 19 nodes, an optimizer started from two different initial values reaches the same estimates to three decimal places, and the maximum absolute component of the scaled gradient is \(0.0004\).
The log-likelihood change is \((-61.203)-(-61.204)=0.001\), and the displayed estimates are stable across node counts. Agreement from distinct starting points and a small scaled gradient further support convergence to the same local maximum. Together these checks suggest numerical stability for this fit, though they do not establish that the model itself is adequate. If, instead, two starting values produced materially different estimates or the gradient remained large, the optimizer might have stopped away from a stationary point even if the integration rule were accurate.
Diagnostics, Alternatives, and Limitations
Useful checks should be reported as part of the fitting workflow, not treated as optional software details. Warning signs include estimates that shift noticeably with node count, a log likelihood that changes erratically as nodes are added, a non-positive curvature at the proposed integration mode, failed mode-finding for some clusters, or optimizer warnings and large gradients. A boundary estimate \(\hat{\sigma}=0\) may be a legitimate result, but it can also expose numerical difficulty near the boundary; it should not be interpreted as proof that the true frailty variance is exactly zero.
Calculations should be performed on the log scale where possible. Directly multiplying many small cluster likelihoods can underflow; summing cluster log likelihoods avoids that problem. The quadrature sum itself can also be stabilized by subtracting a large exponent before exponentiating and then restoring the scale. The number of nodes should be increased until the quantities important to the analysis, including the parameter estimates and relevant likelihood comparisons, are adequately stable.
| Approach | How it differs | Main consideration |
|---|---|---|
| Fixed Gaussian quadrature | Uses a common node grid rather than centering and scaling it at each cluster’s integrand mode. | Can be less efficient when cluster integrands are concentrated away from the origin; may need many nodes. |
| Laplace approximation | Uses the local quadratic expansion at the mode directly. | Often economical, but can be inaccurate when the integrand is not close to Gaussian. |
| Simulation-based integration | Approximates the integral using random or quasi-random draws. | Can handle more complex integration, but introduces simulation error and may require many draws. |
| Analytic integration | Uses a closed-form expression when the frailty distribution and likelihood permit one. | Avoids quadrature error, but is available only for suitable model combinations. |
Adaptive quadrature is most attractive when the number of random-effect dimensions is modest and each cluster integral is smooth and concentrated around a clear mode. It is not a universal solution: high-dimensional random effects can make the tensor-product node count impractical, and highly irregular or heavy-tailed integrands may need many nodes or a different integration strategy. Numerical stability also cannot compensate for a misspecified conditional likelihood, incorrect risk sets, or an unsuitable frailty distribution.
Common Mistakes
- Treating a finite-node approximation as exact. A successful fit at one node count does not certify integration accuracy. Repeat the fit with more nodes and assess stability.
- Checking only the optimizer’s convergence flag. That flag concerns the optimization routine, not whether the cluster integrals were accurate. Check node sensitivity separately.
- Optimizing cluster frailties as if they were observed effects. The marginal likelihood integrates over their assumed distribution. Maximizing a conditional likelihood over each frailty is a different procedure.
- Interpreting a stable fit as proof of a correct model. Stability supports numerical reliability for the specified model; it does not validate its assumptions or establish causality.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- Why does the marginal likelihood multiply cluster contributions, and what induces dependence among observations within a cluster?
- In the adaptive transformation, what roles do the mode \(m_i\) and curvature \(q_i\) play?
- A fit’s estimates change substantially when the node count increases from 7 to 15, although the optimizer reports convergence both times. What should be concluded and checked next?
- Why is agreement across node counts not evidence that the frailty model is correctly specified?
Solutions. (1) Conditional independence across clusters gives the product; the shared latent frailty induces within-cluster dependence. (2) The mode centres the nodes where the integrand is largest, and the curvature sets their local scale. (3) Numerical integration may not yet be stable; increase the node count, inspect the likelihood and estimates, and separately check gradients and starting values. (4) Node stability concerns numerical approximation under the specified model, not whether its likelihood, covariates, or frailty distribution describe the data adequately.