Why Adapt the Quadrature Grid?
A subject’s contribution to a joint-model likelihood is obtained by integrating its conditional longitudinal and event-time contribution over that subject’s random effects. In notation, let \(b_i\) be a \(q\)-dimensional random-effects vector, let \(g_i(b_i;\theta)\) be the subject’s conditional contribution given those effects, and let its density be \(f_b(b_i;\theta)\). The marginal contribution is
Fixed-node Gauss–Hermite quadrature evaluates this integral at nodes placed according to the random-effects distribution. But the integrand is the product of that distribution and the subject-specific conditional contribution. For a particular subject, that product may be concentrated far from the random-effects distribution’s center, or be much narrower than the distribution itself. A grid that does not reach the important region can approximate the integral poorly even when its number of nodes seems adequate.
Adaptive Gauss–Hermite quadrature addresses this by finding where the integrand is largest for each subject, then scaling the quadrature grid to the integrand’s local curvature there. It is a way of choosing integration nodes, not a different likelihood or a different random-effects model.
Posterior Mode and Curvature
The posterior here means the conditional distribution of the latent random effects given the subject’s observed data, at fixed parameter value \(\theta\). The log kernel can omit constants that do not depend on \(b_i\) when finding its mode and curvature. However, those constants must not be dropped from the marginal likelihood being approximated: they affect the integral’s value.
For the derivation, use the standard Gauss–Hermite convention \(\int_{-\infty}^{\infty}e^{-x^2}F(x)\,dx\approx\sum_{j=1}^m w_jF(x_j)\), where \(x_j\) and \(w_j\) are the nodes and weights. In one dimension, write \(H_i\) for the positive scalar curvature and change variables to \(b_i=\hat b_i+\sqrt{2/H_i}\,x\). The factor \(\sqrt{2/H_i}\) is the Jacobian. Substitution gives the adaptive approximation
The \(+x_j^2\) term is essential: it compensates for the \(e^{-x^2}\) weighting built into Gauss–Hermite quadrature. The Jacobian is essential too. Omitting either generally changes the integral, rather than merely changing the accuracy of its approximation.
For \(q\) dimensions, use a symmetric positive-definite matrix square root \(H_i^{-1/2}\) and set \(b_i=\hat b_i+\sqrt{2}H_i^{-1/2}x\). The corresponding tensor-product formula is
where \(x=(x_{j_1},\ldots,x_{j_q})^\mathsf{T}\), and \(|H_i|\) denotes the determinant. The factor \(2^{q/2}/|H_i|^{1/2}\) is the absolute Jacobian. The multivariate expression explains the role of curvature without requiring a particular matrix factorization.
Worked Examples
Worked Example: A Quadratic Kernel as a Calibration Check
Consider the hypothetical one-dimensional subject kernel \(h(b)=-1.3-\frac{9}{2}b^2\). Its mode is \(\hat b=0\), because the squared term is nonpositive and is zero at \(b=0\). Its curvature is \(H=9\). Thus \(I=\int e^{h(b)}\,db\) can be evaluated exactly as a Gaussian integral:
With one Gauss–Hermite node, \(x_1=0\) and \(w_1=\sqrt{\pi}\). The adaptive node in the original scale is \(b=0+\sqrt{2/9}(0)=0\). At this node, \(h(b)-h(0)+x_1^2=0\), so the approximation is \(e^{-1.3}\sqrt{2/9}\sqrt{\pi}\approx 0.2277\). It is exact because the log kernel is quadratic and the transformed integrand is constant after the \(e^{-x^2}\) weight has been accounted for.
Worked Example: Retaining the Kernel’s Scale
For another hypothetical kernel, let \(h(b)=-5.8-\frac{1}{2}b^2\). The mode is \(0\) and the curvature is \(1\). With one adaptive node, the transformed node remains \(b=0\), and the exponent correction is zero. Therefore
This is also the exact integral for this quadratic kernel. In particular, the constant \(-5.8\) contributes the multiplicative factor \(e^{-5.8}\); dropping it would give the wrong marginal contribution. Rounded to four significant figures, the integral is \(0.007589\).
Worked Example: Curvature Scaling for a Nonquadratic Kernel
Suppose a hypothetical subject’s one-dimensional log kernel is \(h(b)=-1.1-2(b-0.4)^2-0.15(b-0.4)^4\). Its derivative is \(-4(b-0.4)-0.6(b-0.4)^3\), which is zero at \(b=0.4\). The second derivative there is \(-4\), so the mode is \(\hat b=0.4\) and \(H=4\). The adaptive transformation is \(b=0.4+x/\sqrt{2}\).
Use the three-point Gauss–Hermite nodes \(0,\pm\sqrt{3/2}\), with weights \(2\sqrt{\pi}/3,\sqrt{\pi}/6,\sqrt{\pi}/6\), respectively. At the central node, the exponential correction in the sum is 1. At either side node, \((b-0.4)^2=x^2/2\), so the quadratic part of \(h(b)-h(0.4)\) is \(-x^2\), which cancels the \(+x^2\) correction. The remaining quartic part is \(-0.15x^4/4=-0.0375x^4\). Since \(x^4=(3/2)^2=2.25\), each side correction is \(e^{-0.084375}\approx0.9191\). Hence
The mode determines the center and the curvature determines the scale; the nonquadratic term remains in the summand and is approximated by the finite quadrature rule. This approximation is not exact, as it was in the first two quadratic examples.
Assumptions, Implementation, and Diagnostics
The transformation and formula require a finite, well-defined integral and an interior mode at which \(h_i\) is sufficiently smooth to calculate its Hessian. For the stated curvature scaling, \(H_i\) must be positive definite. In one dimension this means \(H_i>0\). These conditions make the local quadratic shape a usable scale for the grid; they do not guarantee that a finite number of quadrature nodes accurately captures the entire integral.
At the current parameter value, combine the conditional longitudinal and event-time contribution with the random-effects density. Include any relevant likelihood conditioning in the subject contribution.
Maximize \(h_i\) with respect to the random effects, then evaluate the negative Hessian at the solution. Check that the optimization has converged and the curvature is positive definite.
Center at \(\hat b_i\), scale by \(H_i^{-1/2}\), and include both the Jacobian and the exponential correction in the quadrature sum.
Repeat with a larger number of nodes and compare subject contributions and the total log likelihood. Investigate substantial changes rather than assuming the initial grid was adequate.
Useful warnings include a mode that changes sharply with small parameter updates, a Hessian that is singular or not positive definite, and an approximation that changes noticeably when the node count increases. Multiple local maxima can also make a single-mode grid misleading: an optimizer may find one peak while another important region contributes appreciably to the integral. A highly skewed or heavy-tailed kernel may be poorly represented by nodes scaled only to local curvature. Numerical optimization and derivative accuracy matter because errors in the mode or Hessian directly affect where the nodes are placed.
A practical check is to compare the adaptive approximation at successively larger node counts, subject by subject where feasible, and then compare the resulting sample log likelihood. Also inspect optimization convergence and the eigenvalues of the curvature matrix. These checks assess numerical adequacy; agreement at two node counts is reassuring but is not a proof of accuracy.
What Adaptation Does—and Does Not Do
| Approach | Node placement | Main consideration |
|---|---|---|
| Fixed Gauss–Hermite quadrature | Uses a grid centered and scaled according to the random-effects distribution. | Can work well when the subject kernels are close to that reference scale; can miss a shifted or concentrated subject kernel. |
| Adaptive Gauss–Hermite quadrature | Centers at each subject’s posterior mode and scales using local curvature. | Can use nodes more efficiently near the important region, but requires mode and Hessian calculations and still needs node-count checks. |
| Laplace approximation | Uses a local quadratic approximation at a mode rather than a finite Gauss–Hermite sum. | Is a distinct approximation; adaptive quadrature retains a numerical sum that can represent nonquadratic shape around the mode. |
Adaptation is especially useful when subject-specific posterior kernels are concentrated away from the random-effects distribution’s center. It does not remove integration error, guarantee a unique mode, or make a high-dimensional tensor grid inexpensive: with \(m\) nodes in each of \(q\) dimensions, the sum has \(m^q\) terms. Nor does it cure model misspecification. Its purpose is to allocate a fixed quadrature grid more effectively for each subject.
A common mistake is to treat the posterior mode as the only important point and report the kernel value there as the integral. A mode is a location, not an area: the integral also depends on the kernel’s spread and shape. Adaptive quadrature accounts for local scale through the Jacobian and evaluates multiple transformed nodes to capture departures from a purely quadratic shape.
Check Your Understanding
Answer each question, then compare with the solutions below.
- What is included in \(h_i(b_i;\theta)\), and why is its mode subject-specific?
- In one dimension, if the curvature at the mode is \(H=16\), what is the adaptive transformation from a standard node \(x\) to \(b\)?
- Why does the one-node rule give an exact result for a quadratic log kernel?
- Name two numerical diagnostics for an adaptive quadrature calculation.
Solutions. (1) It is the log of the conditional subject contribution plus the log random-effects density. It is subject-specific because each subject has different observed longitudinal and event-time data, which change the conditional contribution. (2) It is \(b=\hat b_i+\sqrt{2/16}\,x=\hat b_i+x/\sqrt{8}\), where \(\hat b_i\) is that subject’s mode. (3) After the curvature-based change of variables, a quadratic log kernel is exactly canceled by the \(+x^2\) correction, leaving a constant integrand under the Gauss–Hermite weight; the one-node sum integrates that constant exactly. (4) Examples include increasing the node count and comparing the approximations, checking optimizer convergence, and checking that the curvature is positive definite.