What a One-Dimensional Calculation Requires
A Gaussian random effect is integrated out to obtain a cluster’s marginal likelihood contribution. In one dimension, Gauss–Hermite quadrature replaces that integral with a finite weighted sum. The calculation is compact, but its accuracy depends on matching the integration variable, nodes, weights, and scaling convention. A factor of \(\sqrt{2}\), or a missing normalization constant, changes the calculation rather than merely its presentation.
Consider a cluster with a scalar random intercept \(b\sim N(0,\tau^2)\). If its observations are conditionally independent given \(b\), their conditional likelihood is \(L(b)\), and the marginal likelihood contribution is \(I=\int L(b)\phi(b;0,\tau^2)\,db\), where \(\phi(\cdot;0,\tau^2)\) is the \(N(0,\tau^2)\) density. The calculation below applies to any integrable \(L(b)\); the examples use Gaussian conditional likelihoods so that an exact reference can also be found.
Some tables instead specify nodes \(x_j\) and weights \(a_j\) for the physicists’ Hermite rule, which approximates integrals weighted by \(e^{-x^2}\). Its conversion to the standard-normal convention is \(z_j=\sqrt{2}x_j\) and \(w_j=a_j/\sqrt{\pi}\). The three-node physicists’ rule has nodes \(-\sqrt{3/2},0,\sqrt{3/2}\) and weights \(\sqrt{\pi}/6,2\sqrt{\pi}/3,\sqrt{\pi}/6\). Converting gives the standard-normal nodes and weights above. Do not combine nodes from one convention with weights from the other.
Worked Examples
Worked Example: Transforming Nodes for a Gaussian Random Effect
Suppose \(b\sim N(0,0.64)\), so \(\tau=0.8\), and the likelihood kernel is \(L(b)=\exp(-b^2/2)\). Use the three-node rule to approximate \(I=E\{L(b)\}\).
The standard-normal nodes are \(z_j=0,-\sqrt{3},\sqrt{3}\), so the corresponding random-effect values are \(b_j=0,-0.8\sqrt{3},0.8\sqrt{3}\), or approximately \(0,-1.3856,1.3856\). At the center, \(L(0)=1\). At either outer node, \(L(b_j)=\exp\{- (1.3856)^2/2\}=\exp(-0.96)\approx0.3829\).
Substitution into the weighted sum gives \(I\approx(2/3)(1)+(1/6)(0.3829)+(1/6)(0.3829)=0.7943\). The exact expectation is \(1/\sqrt{1+\tau^2}=1/\sqrt{1.64}\approx0.7809\), obtained by combining the two quadratic terms in the normal density and the likelihood kernel. Thus this three-node estimate is about \(0.01343\), or \(1.72\%\), above the exact value. This example checks both scaling the nodes by \(\tau\) and applying the weights to the likelihood evaluated at those scaled nodes.
Worked Example: A Three-Node Estimate and an Exact Reference
Now let \(b\sim N(0,1)\), and suppose one Gaussian observation has observed value \(y=1\), conditional mean \(b\), and conditional variance \(1\). Omitting the conditional density’s constant \(1/\sqrt{2\pi}\), the likelihood kernel is \(L(b)=\exp\{-(1-b)^2/2\}\). We first calculate \(I=\int L(b)\phi(b;0,1)\,db\); the marginal density of \(y\) is \(I/\sqrt{2\pi}\).
At the three standard-normal nodes, the likelihood values are \(L(0)=e^{-1/2}\approx0.6065\), \(L(\sqrt{3})=\exp\{-(1-\sqrt{3})^2/2\}\approx0.7650\), and \(L(-\sqrt{3})=\exp\{-(1+\sqrt{3})^2/2\}\approx0.0239\). Therefore,
For an exact reference, complete the square in the exponent: \(\frac12\{b^2+(1-b)^2\}=(b-\frac12)^2+\frac14\). Integration then gives \(I=(1/\sqrt{2})e^{-1/4}\approx0.5507\). The three-node estimate is low by about \(0.0149\), or \(2.70\%\) of the exact integral. The corresponding marginal density is \(0.5507/\sqrt{2\pi}\approx0.2197\). Comparing estimates of \(I\) with a marginal density would be a scale mismatch: the conditional normalizing constant must be included on both sides, or on neither.
Worked Example: Adapting Nodes to the Integrand
Use the same integral as the preceding example, but place the nodes around the mode of the full log integrand. Define \(q(b)=\phi(b;0,1)L(b)\) and \(h(b)=\log q(b)\), including the prior-density constant. Thus \(h(b)=-\frac12\log(2\pi)-\frac12b^2-\frac12(1-b)^2\).
Differentiating gives \(h'(b)=1-2b\), so the mode is \(m=0.5\). Also \(h''(b)=-2\); the negative curvature is \(H=2\). Adaptive three-node points are \(b_j=m+z_j/\sqrt{H}\), giving \(0.5\) and \(0.5\pm\sqrt{3/2}\). The adaptive approximation is
where \(w_j\) are the standard-normal weights \(1/6,2/3,1/6\). Here \(h(b)=h(m)-(b-m)^2\) exactly, so every correction factor inside the sum is \(\exp\{-(b_j-m)^2+(b_j-m)^2\}=1\). The weights sum to one. Since \(h(m)=-\frac12\log(2\pi)-\frac14\), the result is \(\widehat I_{\mathrm{AGH}}=\exp\{h(m)\}\sqrt{\pi}=(1/\sqrt{2})e^{-1/4}\approx0.5507\), exactly the reference integral. This exactness follows from the quadratic log integrand, not from a general guarantee that three adaptive nodes are exact.
Why the Adaptive Formula Has This Form
Adaptive quadrature uses the local Gaussian shape of the integrand to choose its center and scale. Let \(m\) be a mode of \(h(b)=\log q(b)\), and let \(H=-h''(m)>0\). A second-order Taylor expansion at \(m\) gives \(h(b)\approx h(m)-H(b-m)^2/2\). To retain the difference between the true integrand and this local quadratic approximation, write
Set \(b=m+z/\sqrt{H}\). Then \(db=dz/\sqrt{H}\), and the Gaussian factor becomes \(\exp(-z^2/2)\). Applying standard-normal quadrature to the remaining weighted expectation gives the adaptive formula in the example, including \(\sqrt{2\pi/H}\). The correction factor is essential: it accounts for departures from the local Gaussian approximation at each quadrature point. Omitting it would integrate only the local quadratic approximation.
This argument requires a finite integral, an interior mode at which \(h\) is twice differentiable, and positive curvature \(H\). If the integrand is not well represented by one local Gaussian shape, the approximation may need more nodes or a different integration strategy. Multiple separated modes, for example, cannot be reliably represented merely by centering at one of them.
Checking a Hand Calculation
A reliable calculation keeps the following objects distinct: the random effect \(b\), its standard-normal coordinate \(z\), the quadrature nodes, the quadrature weights, and any likelihood constants. For fixed-node quadrature, transform the nodes using \(b=\tau z\). For adaptive quadrature, use \(b_j=m+z_j/\sqrt{H}\), not \(m+\sqrt{H}z_j\). In either case, check that the weights sum to one under the standard-normal convention.
The earlier material on quadrature order and numerical convergence motivates comparing results at increasing orders. For a hand calculation, also check the weighted contributions individually: a contribution is \(w_jL(b_j)\), not just \(L(b_j)\). A small outer-node contribution can still matter, but it must be multiplied by its actual weight. For an independent reference, direct one-dimensional numerical integration is useful; an exactly integrable Gaussian case is especially valuable for checking node scaling, constants, and arithmetic.
| Approach | What is evaluated | Main consideration |
|---|---|---|
| Fixed-node Gauss–Hermite | Likelihood at nodes centered at zero and scaled by the prior standard deviation | Simple to calculate; may use nodes inefficiently if the integrand is concentrated away from zero |
| Adaptive Gauss–Hermite | Corrected likelihood at nodes centered and scaled using a mode and local curvature | Can target a concentrated integrand; depends on locating a relevant mode and calculating curvature correctly |
| Laplace approximation | The local Gaussian approximation at the mode, without a quadrature sum over correction terms | Compact; accuracy depends on how well the local quadratic approximation represents the integral |
| Direct one-dimensional numerical integration | The integrand over a chosen numerical integration range or transformed domain | Useful as a check; numerical truncation and integration tolerances must be controlled |
The Gaussian example shows a special case where adaptive quadrature and Laplace integration are exact: the log integrand is quadratic, so its local Gaussian form is the entire integrand. In a nonquadratic GLMM integral, adaptive nodes do not make the integrand Gaussian; the correction terms and quadrature order still matter. Nor does a close match in one cluster guarantee equally accurate marginal likelihoods across all clusters or parameter values.
Common Mistakes and Limitations
- Mixing Hermite conventions. Physicists’ nodes and weights require a conversion before using the standard-normal weighted sum. Check that the resulting standard-normal weights sum to one.
- Scaling in the wrong direction. Fixed nodes use \(b=\tau z\); adaptive nodes use \(b=m+z/\sqrt{H}\). These transformations reflect different integration coordinates.
- Dropping or adding a constant inconsistently. If the target is the integral of a likelihood kernel, compare it with a reference for that same kernel. For the full marginal density, include the conditional density constants consistently.
- Adding adaptive curvature twice or not at all. The correction exponent is \(h(b_j)-h(m)+H(b_j-m)^2/2\). It corrects the true log integrand relative to the local quadratic form.
- Assuming a low-order estimate is exact. Exactness in the worked adaptive example comes from a quadratic log integrand. For a general integrand, compare orders or use a suitably accurate numerical reference.
Hand calculations are most useful as transparent checks on a one-dimensional integral. They do not remove the need to examine numerical convergence, validate the model’s conditional likelihood and random-effects distribution, or ensure that the selected mode and curvature describe the part of the integrand carrying substantial mass.
Check Your Understanding
Try each question, then compare your reasoning with the solutions.
- For a three-node standard-normal rule, what are the nodes and weights? If \(b\sim N(0,0.25)\), what random-effect values correspond to those nodes?
- For the integral in the second worked example, what is the three-node estimate, and is it above or below the exact integral?
- In adaptive quadrature, what condition on \(H\) is needed for the stated node transformation, and why?
- Why must the adaptive correction term include \(H(b_j-m)^2/2\)?
- In the Gaussian adaptive example, why does the three-node estimate equal the exact integral?
Solutions. (1) The nodes are \(-\sqrt{3},0,\sqrt{3}\), with weights \(1/6,2/3,1/6\). Since \(\tau=0.5\), the corresponding values are \(-0.8660,0,0.8660\). (2) The estimate is approximately \(0.5358\), below the exact integral \(0.5507\). (3) One needs \(H>0\), so the mode has positive local curvature on the negative-log scale and \(1/\sqrt{H}\) is a real scale factor. (4) The term restores the quadratic exponent that was factored out to create the standard-normal weight; together with \(h(b_j)-h(m)\), it gives the correction to the local Gaussian shape. (5) The log integrand is exactly quadratic, so each correction factor is one and the adaptive Gaussian integral is the entire integral, not merely a local approximation.