Why Adapt the Nodes?
Gauss-Hermite quadrature approximates an integral by evaluating its integrand at a finite set of nodes and combining those evaluations with weights. In a nonlinear mixed-effects model, a subject’s likelihood contribution can be sharply concentrated in a small region of random-effect space. A fixed set of nodes may place few nodes in that region. Adaptive quadrature addresses this mismatch by centering and scaling the nodes using the subject’s conditional random-effect mode and local curvature.
The adaptation is subject-specific: different observations can provide different information about their random effects, so their integrands can have different centers and widths. The procedure does not change the target marginal likelihood. It changes the locations and weights used to approximate its integral.
The standard-normal coordinates are obtained by expressing the original Gaussian random effects as a linear transformation of \(x_i\); this makes the prior density \(\phi_q(x_i)\) part of the kernel being integrated. In particular, the mode is the mode of the likelihood times the prior density, not generally the maximizer of the conditional likelihood alone. The prior term matters because it can pull the mode toward zero.
Deriving the Adaptive Formula
For one dimension, Gauss-Hermite nodes \(u_j\) and weights \(w_j\) are defined to approximate integrals of the form \(\int_{-\infty}^{\infty} e^{-u^2}f(u)\,du\) by \(\sum_j w_j f(u_j)\). In \(q\) dimensions, the ordinary tensor-product rule uses products of these nodes and weights. To adapt it, approximate the local shape of \(\ell_i(x)\) near its mode and transform the integration variable so that the local quadratic curvature matches the Gauss-Hermite weight.
A second-order Taylor expansion at the mode gives \(\ell_i(x)\approx \ell_i(m_i)-\tfrac12(x-m_i)^\mathsf{T}H_i(x-m_i)\), since the gradient at an interior mode is zero. Let \(A_i\) satisfy \(A_iA_i^\mathsf{T}=H_i^{-1}\), and set \(x=m_i+\sqrt{2}A_i u\). Then \(\tfrac12(x-m_i)^\mathsf{T}H_i(x-m_i)=u^\mathsf{T}u\), while the absolute Jacobian determinant is \(2^{q/2}|H_i|^{-1/2}\). Substitution into the integral and application of the product Gauss-Hermite rule yield the adaptive estimate.
The factor \(e^{u^\mathsf{T}u}\) inside the sum is essential: the Gauss-Hermite weights already approximate integration against \(e^{-u^\mathsf{T}u}\). The determinant factor is also essential; it accounts for the change in integration volume when the nodes are rescaled. Omitting either factor generally gives a wrongly scaled likelihood.
Worked Calculations
Worked Example: A Shifted Mode With an Exact Integral
Consider a hypothetical one-dimensional subject contribution with standard-normal random effect \(X\) and conditional likelihood kernel \(g(x)=\exp\{-(x-1)^2\}\). The marginal contribution is \(L=E[g(X)]\) for \(X\sim N(0,1)\). Find the conditional mode, curvature, and one-node adaptive estimate.
The log joint kernel, including the standard-normal density, is \(\ell(x)=-\tfrac12\log(2\pi)-\tfrac12x^2-(x-1)^2\). Differentiating gives \(\ell'(x)=-x-2(x-1)\), so \(\ell'(x)=0\) at \(m=2/3\). The negative second derivative is \(H=-\ell''(x)=3\), which is positive, confirming a local maximum with positive curvature.
For one Gauss-Hermite node, \(u_1=0\) and \(w_1=\sqrt{\pi}\). Thus the estimate is \(L\approx \sqrt{2/3}\sqrt{\pi}\exp\{\ell(2/3)\}\). Here \(\ell(2/3)=-\tfrac12\log(2\pi)-\tfrac13\), so \(L\approx \sqrt{2\pi/3}(2\pi)^{-1/2}e^{-1/3} = e^{-1/3}/\sqrt{3}\approx 0.4137\).
Completing the square in the original integral gives the exact value \(L=e^{-1/3}/\sqrt{3}\), so one adaptive node is exact here. The mode is shifted from the prior mean zero because this subject’s likelihood favors \(x=1\); the curvature determines the corresponding node scale.
Worked Example: Evaluating a Nonquadratic Kernel
For a numerical illustration of the node transformation, take a hypothetical integrable one-dimensional joint kernel whose log is \(\ell(x)=-1.78-2(x-0.5)^2-0.1(x-0.5)^4\). This constructed kernel is used to show the calculation, not to represent a particular study. Its maximum is at \(m=0.5\), since the first derivative vanishes there and the negative second derivative is \(H=4\).
The three-point Gauss-Hermite rule has nodes \(0,\ \pm\sqrt{3/2}\) and weights \(2\sqrt{\pi}/3,\ \sqrt{\pi}/6,\ \sqrt{\pi}/6\). Under the adaptive transformation, the corresponding \(x\) values are \(0.5,\ 0.5\pm\sqrt{3}/2\), approximately \(0.5000,\ 1.3660,\ -0.3660\). This uses \(x=0.5+\sqrt{2}(1/2)u\).
At the center, the summand’s exponential factor is \(e^{-1.78}\). At either outer node, \(x-0.5=\pm\sqrt{3}/2\), so the exponent inside the quadrature sum, after including \(u^2\), changes by \(-2(3/4)-0.1(9/16)+3/2=-0.05625\). The three-point adaptive estimate is therefore \[ \frac{1}{\sqrt{2}}e^{-1.78} \left\{\frac{2\sqrt{\pi}}{3} +2\frac{\sqrt{\pi}}{6}e^{-0.05625}\right\} \approx 0.2075. \]
For comparison, one adaptive node gives \(\sqrt{\pi/2}e^{-1.78}\approx 0.2114\). The three-point rule changes the estimate because it samples away from the mode and reflects the kernel’s nonquadratic tails. More nodes do not guarantee correctness in every problem, but this calculation shows how the transformed locations and curvature adjustment enter the integral.
Worked Example: Why Each Subject Gets Its Own Center
Suppose two hypothetical subjects have conditional likelihood kernels \(g_i(x)=\exp\{-(x-a_i)^2\}\), with \(a_1=1\) and \(a_2=-2\), and both have the same standard-normal random-effect prior. For a general value \(a_i\), the log joint kernel is \(\ell_i(x)=-\tfrac12\log(2\pi)-\tfrac12x^2-(x-a_i)^2\).
Differentiating gives \(\ell_i'(x)=-x-2(x-a_i)\). Setting this to zero gives \(m_i=2a_i/3\), and \(H_i=3\) for both subjects. Thus the centers are \(m_1=2/3\) and \(m_2=-4/3\), although the scales are equal. Completing the square gives the exact marginal contribution \(L_i=e^{-a_i^2/3}/\sqrt{3}\).
For subject 1, \(L_1=e^{-1/3}/\sqrt{3}\approx0.4137\). For subject 2, \(L_2=e^{-4/3}/\sqrt{3}\approx0.1522\). The difference in centers reflects the different subject-specific likelihood information. Reusing one subject’s mode for the other would place the second subject’s nodes around the wrong region, even though their curvatures happen to match.
Assumptions, Diagnostics, and Alternatives
The formula requires a finite marginal integral and a well-defined interior mode with positive-definite \(H_i\). The log kernel should be sufficiently smooth near that mode for its curvature to describe a meaningful local shape. These conditions justify the local centering and scaling; they do not imply that a low-order quadrature rule is accurate over the entire integration domain. In more than one dimension, the curvature matrix must be invertible, and a factor \(A_i\) with \(A_iA_i^\mathsf{T}=H_i^{-1}\) must be computed reliably.
Useful warning signs include a non-positive or nearly singular curvature matrix, a mode that changes substantially under small numerical perturbations, and quadrature estimates that change materially when the node count is increased. A single mode can also be misleading if the kernel is asymmetric, has multiple separated peaks, or has important tail mass far from the mode. In those cases the local curvature describes only one region. Optimization and integration should be checked separately: a well-converged mode search does not establish that the integral is accurate.
A common mistake is to say, “The likelihood is concentrated, so I can use fewer nodes without checking.” Concentration near one mode can help adaptive quadrature, but it does not rule out skewness, tails, or additional modes. Compare estimates at increasing node counts and inspect whether the mode and curvature are stable. When a normal approximation is poor, importance sampling can place draws around a proposal distribution and average weighted contributions instead; unlike deterministic adaptive quadrature, it has Monte Carlo error. Fixed Gauss-Hermite quadrature is simpler and reuses the same nodes for all subjects, but can be inefficient when subject kernels are far from the origin or differ greatly in width. Adaptive quadrature improves placement, not the fundamental suitability of a local Gaussian shape.
Include both the subject’s conditional likelihood and the random-effect prior density in \(\ell_i(x;\theta)\).
Optimize \(\ell_i\), verify the gradient is near zero, and confirm that \(H_i=-\nabla^2\ell_i(m_i;\theta)\) is positive definite.
Use \(x=m_i+\sqrt{2}A_i u\), retain the Jacobian and the \(e^{u^\mathsf{T}u}\) correction, then sum the weighted evaluations.
Repeat with a larger node count and compare the resulting subject contributions and total log likelihood.
Check Your Understanding
Try each problem, then compare with the brief solutions below.
- In one dimension, what does the curvature \(H_i\) measure, and how does a larger \(H_i\) affect adaptive node spacing?
- Why is the conditional mode generally defined using the likelihood multiplied by the prior density?
- For \(\ell(x)=-\tfrac12x^2-(x-1)^2+C\), find the mode and curvature.
- Name two conditions or diagnostics that should prompt caution when using adaptive quadrature.
Solutions. (1) \(H_i\) is the negative second derivative of the log joint kernel at its mode; larger curvature means a smaller local scale \(H_i^{-1/2}\), so nodes lie closer together. (2) The marginal integral contains the prior as well as the conditional likelihood, and the prior can shift the mode away from the likelihood maximizer. (3) The derivative is \(-x-2(x-1)\), giving \(m=2/3\); the negative second derivative is \(H=3\). (4) Examples include a singular or non-positive-definite curvature matrix, unstable mode estimates, substantial changes as node count increases, or evidence of multiple peaks and important tails.