Why Move the Quadrature Nodes?
A random-effects integral can be difficult to approximate with nodes fixed at zero when most of its mass lies elsewhere. This can happen when the observed cluster data pull the conditional distribution of its random effect away from the population mean. Adaptive Gauss-Hermite quadrature responds by centering its nodes at the conditional mode and scaling their spread to reflect the local curvature there.
The procedure is useful when the integrand is concentrated around a location or width that a fixed rule does not represent well. It is still a numerical approximation: moving the nodes does not guarantee accuracy if the integrand is highly skewed, has important tails, or has multiple separated modes.
For a random-intercept GLMM, the integrand is commonly the conditional response likelihood multiplied by the random-effects density. For example, if \(b\sim N(0,\sigma_b^2)\), then \(I\) integrates \(p(y\mid b)\) against the \(N(0,\sigma_b^2)\) density. The conditional mode maximizes this full product, not the response likelihood alone. The mode-finding and curvature calculations rely on a sufficiently smooth log-integrand and a finite integral. A single-mode adaptive rule additionally requires care: the chosen mode should represent the region containing the important mass.
Adaptive Rule in Standard-Normal Coordinates
Let \(z_j\) and \(a_j\) be nodes and normalized weights for a Gauss-Hermite rule expressed as an approximation to expectation under \(Z\sim N(0,1)\). Thus the weights sum to one, and the rule approximates \(E[f(Z)]\) by \(\sum_j a_j f(z_j)\). Under the local scale transformation \(b=\hat b+\sqrt{V}z\), adaptive quadrature gives
The exponential factor inside the sum adjusts for the difference between the actual integrand and the local Gaussian shape used to choose the nodes. In particular, the \(z_j^2/2\) term compensates for the standard-normal weighting of the quadrature rule. Omitting it would integrate a different function.
A brief justification follows from the local quadratic approximation \(h(b)\approx h(\hat b)-(b-\hat b)^2/(2V)\). On setting \(b=\hat b+\sqrt Vz\), the integral is written as a standard-normal expectation times the factor \(\sqrt{2\pi V}\exp\{h(\hat b)\}\). The remaining expectation has integrand \(\exp\{h(\hat b+\sqrt Vz)-h(\hat b)+z^2/2\}\), which is approximated at the Gauss-Hermite nodes. If \(h\) is exactly quadratic, this remaining factor is identically one, so the rule integrates the Gaussian form exactly (apart from arithmetic error).
For the three-point standard-normal rule, the nodes are \(0,\sqrt 3,-\sqrt 3\), with respective weights \(2/3,1/6,1/6\). The scale is \(\sqrt V\), so the actual random-effect nodes are \(\hat b\), \(\hat b+\sqrt{3V}\), and \(\hat b-\sqrt{3V}\). For a general number of points, the same transformation applies to that rule’s nodes and weights.
Worked Example: A Non-Gaussian Integrand With a Shifted Mode
Consider \(I=\int_{-\infty}^{\infty}\exp\{h(b)\}\,db\), where \(h(b)=-(b-0.4)^2/2-0.05(b-0.4)^4\). This is a hypothetical integrand used to illustrate the calculation. Its derivative is \(h'(b)=-(b-0.4)-0.2(b-0.4)^3\), which is zero at \(\hat b=0.4\). Its second derivative at that point is \(-1\), so \(V=1\). The mode is unique because \(h'(b)\) is positive below 0.4 and negative above it.
Use the three-point rule. The transformed nodes are \(0.4\) and \(0.4\pm\sqrt 3\), approximately \(2.132\) and \(-1.332\). At the center the multiplier is 1. At either off-center node, \(z^2=3\), so \(\exp\{h(\hat b+z)-h(\hat b)+z^2/2\}=\exp(-0.05z^4)=\exp(-0.45)\approx0.6376\).
The weighted multiplier is \((2/3)(1)+(1/6)(0.6376)+(1/6)(0.6376)\approx0.8792\). Therefore \(I\approx\sqrt{2\pi(1)}(0.8792)\approx2.204\). This is a three-point numerical approximation, not a claim that the quartic integrand is exactly integrated. The example shows how centering and scaling are reflected in the nodes and how the correction multiplier accounts for non-Gaussian shape.
Application to a Binary Random-Intercept Cluster
For a binary cluster with a random intercept \(b\sim N(0,1)\), suppose there are four observations, three with response 1 and one with response 0. Let \(p(b)=1/(1+\exp(-b))\). Ignoring the constant normal-density factor for the moment, the log-integrand is \(h(b)=3b-4\log(1+\exp(b))-b^2/2\). The omitted factor is \(1/\sqrt{2\pi}\); it must be restored if the target is the full marginal likelihood.
Worked Example: Adaptive Nodes for a Binary Cluster
The score and curvature are \(h'(b)=3-4p(b)-b\) and \(h''(b)=-4p(b)\{1-p(b)\}-1\). The score decreases strictly because \(h''(b)<0\), so there is at most one mode. It changes from positive at sufficiently negative \(b\) to negative at sufficiently positive \(b\), giving a unique mode. Solving the score equation gives \(\hat b\approx0.5052\). At this value \(p(\hat b)\approx0.6237\), and \(3-4(0.6237)-0.5052\) is approximately zero (the small discrepancy is due to rounding).
The negative curvature is \(4(0.6237)(1-0.6237)+1\approx1.9388\), so \(V\approx1/1.9388\approx0.5158\). The three adaptive nodes are approximately \(0.5052\) and \(0.5052\pm\sqrt{3(0.5158)}\), or \(0.5052,\ 1.749,\ -0.7388\).
Substitution into the stated log-integrand gives \(h(\hat b)\approx-2.5213\), \(h(1.749)\approx-3.9203\), and \(h(-0.7388)\approx-4.0511\). The center multiplier is 1. The positive-node multiplier is \(\exp\{-3.9203-(-2.5213)+3/2\}\approx1.106\); the negative-node multiplier is \(\exp\{-4.0511-(-2.5213)+3/2\}\approx0.9706\).
Thus the weighted multiplier is \((2/3)+(1/6)(1.106+0.9706)\approx1.0128\). The integral with the normalizing constant omitted is approximately \(\exp(-2.5213)\sqrt{2\pi(0.5158)}(1.0128)\approx0.1465\). The full marginal likelihood contribution for this cluster is this value divided by \(\sqrt{2\pi}\), approximately \(0.0585\). Keeping track of the omitted density constant is essential when comparing likelihoods.
What the Curvature Does—and Does Not—Tell You
The value \(V\) is the inverse negative curvature of the log-integrand at its mode. A larger curvature magnitude means a narrower local peak and a smaller \(V\); the adaptive nodes then lie closer together. A smaller curvature magnitude means a wider local peak and a larger \(V\). This is a local scale choice, not an assertion that the true conditional distribution is Gaussian with variance \(V\).
The construction presumes an interior mode with finite, negative second derivative. If the curvature is zero or positive, \(V\) is not a valid positive scale. If the log-integrand is nearly flat at its mode, the scale can be very large and local curvature may poorly describe the relevant region. In vector problems, the corresponding curvature matrix must be positive definite after changing the sign of the Hessian; a singular or indefinite matrix does not define a valid local Gaussian scale.
Adaptive placement can be especially helpful when the mode is far from zero, since a fixed rule at zero may place too few nodes near the dominant part of the integrand. It does not, however, resolve every failure of a single local approximation. Strong skewness can leave one side of the mass poorly represented, and separated modes may require attention beyond nodes around one selected mode. Practical checks include inspecting the log-integrand over the node range, ensuring that meaningful mass is not concentrated beyond the outer nodes, and comparing estimates at more than one quadrature order.
Worked Example: Exact Integration of a Shifted Gaussian
Let \(h(b)=-(b-2)^2/(2(0.5)^2)\). The mode is \(\hat b=2\), and \(h''(b)=-1/(0.5)^2=-4\), so \(V=0.25\). The three adaptive nodes are \(2\) and \(2\pm\sqrt{3(0.25)}\), approximately \(2,\ 2.866,\ 1.134\).
Because \(h(b)=h(\hat b)-(b-\hat b)^2/(2V)\) exactly, every multiplier in the adaptive sum equals 1. The weights sum to one, and the approximation is \(\sqrt{2\pi(0.25)}=\sqrt{\pi/2}\approx1.253\). This is also the exact integral of the stated unnormalized Gaussian kernel. The example confirms that adapting the center and scale recovers a shifted Gaussian integral; it does not imply exactness for non-Gaussian integrands.
Choosing Among Integration Approaches
| Approach | Node placement or approximation | Main consideration |
|---|---|---|
| Nonadaptive Gauss-Hermite quadrature | Uses a fixed rule at zero in the chosen coordinates. | Simple to reuse, but may poorly cover an integrand displaced from zero. |
| Adaptive Gauss-Hermite quadrature | Centers and scales nodes using the conditional mode and local curvature. | Can better represent a concentrated, displaced integrand; still needs adequate node coverage and order. |
| Laplace approximation | Replaces the log-integrand by its quadratic expansion at the mode. | Uses only local mode and curvature information; it is exact for a Gaussian integrand but can be inaccurate when non-Gaussian features matter. |
Adaptive quadrature and the Laplace approximation use the same local mode and curvature, but they use them differently. Laplace replaces the integrand by its local Gaussian form. Adaptive quadrature uses that form to choose a coordinate system and then evaluates the actual integrand at several transformed nodes. It therefore can represent departures from the quadratic approximation, provided the selected nodes and weights adequately capture them. Its cost rises with the number of evaluations, and the mode and curvature must be recomputed when model parameters change during estimation.
Common Mistakes and Practical Checks
- Using the likelihood mode instead of the conditional mode. The random-effects density is part of the integrand, so the mode must maximize the likelihood times that density.
- Using curvature with the wrong sign. At a local maximum \(h''(\hat b)<0\), and the scale is \(V=-1/h''(\hat b)\), not \(1/h''(\hat b)\).
- Forgetting the correction term. The factor \(\exp(z_j^2/2)\) is required with standard-normal nodes and weights.
- Dropping normalization constants inconsistently. A constant may be omitted while finding the mode, but it must be included when the requested integral or likelihood requires it.
- Treating a good local fit as proof of global accuracy. A local mode and curvature do not rule out skewness, heavy tails, or other modes. Inspect coverage and compare numerical approximations.
Check Your Understanding
Answer each question, then compare with the solutions below.
- If \(h''(\hat b)=-2\), what are \(V\) and the scale \(\sqrt V\)?
- For three-point standard-normal quadrature, what are the adaptive nodes when \(\hat b=-0.3\) and \(V=0.64\)?
- Why does the adaptive multiplier include \(z_j^2/2\)?
- What happens to the adaptive rule if the negative curvature at the proposed mode is not positive?
Solutions. (1) \(V=-1/(-2)=0.5\), and \(\sqrt V=\sqrt{0.5}\approx0.7071\). (2) Since \(\sqrt{3V}=\sqrt{1.92}\approx1.386\), the nodes are \(-0.3\), \(-0.3+1.386=1.086\), and \(-0.3-1.386=-1.686\). (3) Standard-normal quadrature weights already incorporate the density factor proportional to \(\exp(-z^2/2)\); the multiplier’s \(\exp(z_j^2/2)\) cancels that weighting when expressing the original integral. (4) A nonpositive negative-curvature value does not yield a positive local variance scale. The proposed point may not be a suitable interior maximum, or the curvature calculation or mode-finding procedure may have failed, so the adaptive rule should not be applied with that scale.