Nonadaptive Gauss-Hermite Quadrature at Zero

Why Use Fixed Nodes at Zero?

A GLMM marginal likelihood integrates the conditional likelihood over the random-effects distribution. In one dimension, the integral is conceptually simple, but the conditional likelihood may have a shape that makes direct calculation difficult. Gauss–Hermite quadrature replaces the integral with a weighted sum of evaluations at selected points. In the nonadaptive version considered here, the nodes are fixed in advance and centered at zero, the center of a standard normal distribution.

The method is especially useful when many likelihood evaluations are needed during estimation: the same nodes can be reused at each parameter value. But the fixed center is also its key limitation. If the integrand’s important mass lies far from zero or is concentrated in a narrow region, a low-order rule may place too few nodes where they matter.

Definition: Let \(b\sim N(0,\sigma_b^2)\) be a scalar Gaussian random effect and let \(L(b)\) denote the conditional likelihood contribution for the observed data at a fixed value of \(b\). The marginal likelihood contribution is \[ \int_{-\infty}^{\infty} L(b)\,\frac{1}{\sigma_b}\phi\!\left(\frac{b}{\sigma_b}\right)\,db =E\{L(\sigma_b Z)\}, \qquad Z\sim N(0,1), \] where \(\phi\) is the standard normal density. A nonadaptive Gauss–Hermite approximation evaluates \(L(\sigma_b z_j)\) at fixed standard-normal nodes \(z_j\), then combines those values using fixed weights \(w_j\).

From a Gaussian Integral to a Weighted Sum

The transformation \(b=\sigma_b z\) converts the random effect to standard-normal coordinates: \(db=\sigma_b\,dz\), and the Gaussian density and Jacobian combine to give \(\phi(z)\,dz\). Thus the marginal likelihood is a standard-normal expectation. This is the relevant setup from the earlier material on transforming random-effects integrals to standard normal coordinates.

For standard-normal quadrature, write the nodes and normalized weights as \(z_1,\ldots,z_m\) and \(w_1,\ldots,w_m\), respectively. The weights sum to one, and the approximation is

E\{L(\sigma_b Z)\} =\int_{-\infty}^{\infty}L(\sigma_b z)\phi(z)\,dz \ \approx\ \sum_{j=1}^{m}w_jL(\sigma_b z_j).

These are the standard-normal versions of the Gauss–Hermite nodes and weights; conventions that use a weight function involving \(\exp(-x^2)\) instead require a change of scale and weight normalization. Here, use the \(z_j,w_j\) convention above. For the three-node rule, the nodes are \(0,\pm\sqrt{3}\), with weight \(2/3\) at zero and weight \(1/6\) at each nonzero node.

The approximation follows directly from applying Gaussian-weighted polynomial quadrature to the function \(L(\sigma_b z)\). The Gaussian density is represented by the quadrature weight, so it is not multiplied into each likelihood evaluation a second time. This is a frequent implementation error: the weighted sum already accounts for the random-effects density.

Justification: If \(q(z)\) is a polynomial of degree at most \(2m-1\), the \(m\)-node Gauss–Hermite rule integrates \(q(z)\phi(z)\) exactly. For a general likelihood function, the rule is an approximation; exactness for polynomials does not imply exactness for an arbitrary GLMM integrand. Increasing \(m\) often improves an adequate approximation, but a low-order rule has no universal accuracy guarantee for a sharply peaked or displaced function.

Worked Calculations

Worked Example: A Three-Node Approximation and Its Error

Consider the hypothetical integrand \(g(z)=\exp(2z-z^2/2)\) under \(Z\sim N(0,1)\). Use the three-node rule to approximate \(E\{g(Z)\}\), and compare it with the exact value.

At the nodes, the function values are \(g(0)=1\), \(g(\sqrt{3})=\exp(2\sqrt{3}-3/2)\approx 7.1285\), and \(g(-\sqrt{3})=\exp(-2\sqrt{3}-3/2)\approx 0.00698\). Therefore,

\begin{aligned} Q_3 &=\frac{2}{3}g(0)+\frac{1}{6}g(\sqrt{3})+\frac{1}{6}g(-\sqrt{3})\\ &\approx \frac{2}{3}+\frac{7.1285}{6}+\frac{0.00698}{6}\\ &\approx 1.8559. \end{aligned}

For the exact expectation, combine \(g(z)\) with the standard normal density and complete the square:

\begin{aligned} E\{g(Z)\} &=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} \exp(2z-z^2)\,dz\\ &=\frac{e}{\sqrt{2\pi}}\int_{-\infty}^{\infty} \exp\{-(z-1)^2\}\,dz =\frac{e}{\sqrt{2}} \approx 1.9221. \end{aligned}

The approximation is about \(1.9221-1.8559=0.0662\) below the exact value. Using the unrounded values, its relative error is approximately \(0.0662/1.9221=0.03444\), or \(3.444\%\). This calculation illustrates both the use of fixed nodes and the fact that a valid quadrature rule need not be exact for a nonpolynomial integrand.

Worked Example: A Binary Random-Intercept Cluster

Suppose a hypothetical cluster has three binary responses \(1,0,1\). Conditional on a random intercept \(b\), the responses are independent Bernoulli variables with common probability \(p(b)=\{1+\exp(-b)\}^{-1}\). Let \(b\sim N(0,\sigma_b^2)\), with \(\sigma_b=\log(4)/\sqrt{3}\approx0.8004\). Approximate the cluster’s marginal likelihood using three nodes.

The conditional likelihood for these data is \(L(b)=p(b)^2\{1-p(b)\}\). The standard-normal nodes \(z=0,\pm\sqrt{3}\) correspond to random-intercept values \(b=0,\pm\log(4)\). Thus \(p(0)=0.5\), \(p(\log(4))=0.8\), and \(p(-\log(4))=0.2\). The likelihood values are consequently \(L(0)=0.5^2(0.5)=0.125\), \(L(\log(4))=0.8^2(0.2)=0.128\), and \(L(-\log(4))=0.2^2(0.8)=0.032\).

Applying the normalized three-node weights gives

\begin{aligned} \widehat{L}_{\mathrm{marg}} &=\frac{2}{3}(0.125)+\frac{1}{6}(0.128)+\frac{1}{6}(0.032)\\ &=0.08333+0.02133+0.00533\\ &=0.1100. \end{aligned}

The resulting approximate log marginal likelihood is \(\log(0.1100)\approx-2.207\). This is an approximation to the probability of the observed response pattern after integrating over the random intercept. The weights represent the Gaussian random-effects distribution; the likelihood values are conditional on the node-specific intercepts.

When the Zero Center Is Inefficient

The rule samples a fixed set of locations around zero regardless of the observed cluster data or the shape of \(L(\sigma_b z)\). It can therefore be inefficient when the conditional likelihood is concentrated away from zero, when it is narrow relative to the spacing between nodes, or when its tail behavior makes a few evaluations dominate the sum. Sparse clusters and extreme response patterns can produce likelihood shapes whose important region is not well covered by low-order nodes. A centered rule can still be accurate in such cases at sufficiently high order, but that may require many evaluations.

Worked Example: A Narrow Peak Away from Zero

Consider the hypothetical nonnegative integrand \(g(z)=\exp\{-(z-2)^2/(2(0.25)^2)\}\). It is concentrated near \(z=2\), with width \(0.25\), while the three fixed nodes are \(0\) and \(\pm\sqrt{3}\). Compare the three-node estimate with the exact standard-normal expectation.

At zero the value is \(\exp(-32)\), which is negligible. At \(\sqrt{3}\), it is \(\exp\{-(2-\sqrt{3})^2/0.125\}\approx0.563\); at \(-\sqrt{3}\), it is negligible. The fixed-node approximation is therefore

Q_3\approx \frac{2}{3}\exp(-32)+\frac{1}{6}(0.563)+\frac{1}{6}(0) \approx 0.0938.

For the exact value, completing the square in the product \(g(z)\phi(z)\) gives, for a Gaussian-shaped function with center \(\mu\) and width \(s\),

E\!\left[\exp\!\left\{-\frac{(Z-\mu)^2}{2s^2}\right\}\right] =\frac{s}{\sqrt{1+s^2}}\exp\!\left\{-\frac{\mu^2}{2(1+s^2)}\right\}.

Substituting \(\mu=2\) and \(s=0.25\) gives \(\{0.25/\sqrt{1.0625}\}\exp\{-4/2.125\}\approx0.0369\). The three-node estimate is substantially larger: its relative error is approximately \((0.0938-0.0369)/0.0369\approx1.54\), or \(154\%\). The node at \(\sqrt{3}\) lies near the peak but does not resolve its narrow shape, while the other nodes contribute almost nothing. This is a concrete case where fixed centering is inefficient even though all calculations use the correct standard-normal weights.

Assumptions, Diagnostics, and Alternatives

The method applies directly when the random effect is one-dimensional Gaussian and the target is the corresponding marginal integral. The conditional likelihood must be evaluable at the transformed nodes, and the integral must exist. The quadrature formula does not require the conditional likelihood itself to be Gaussian. However, an accurate low-order approximation also requires the integrand to be sufficiently well represented at the chosen nodes; Gaussian random effects alone do not ensure that.

Several checks help identify numerical trouble:

Check normalization. The weights should sum to one. Applying the rule to the constant function \(L(b)=1\) should return exactly one, up to numerical rounding.
Check standardization. Evaluate the conditional likelihood at \(b=\sigma_b z_j\). Using \(z_j\) directly as \(b\), or multiplying by the Gaussian density again, changes the integral.
Check sensitivity to the chosen order. Compare estimates from more than one reasonable node count. A material change signals that the current approximation is not stable enough for the intended use. Agreement is useful evidence, not a proof of accuracy.
Inspect the integrand’s location and shape. A conditional likelihood concentrated far from zero, unusually narrow, strongly asymmetric, or with multiple peaks may be poorly represented by a low-order fixed rule.
Protect the arithmetic. Likelihood values can be extremely small. Compute weighted sums on the log scale where needed, using log-sum-exp, rather than allowing underflow to turn positive contributions into numerical zeros.

The principal alternatives differ in where their computational effort is placed. The Laplace approximation uses a local quadratic approximation at the conditional mode and can be economical, especially when the integrand is well approximated near that mode. Fixed zero-centered quadrature evaluates a global set of predetermined points and avoids mode finding at every integration step, but can waste evaluations when the mass is displaced. Adaptive Gauss–Hermite quadrature instead locates and scales nodes using local information around the posterior mode; its construction is treated in the following material. For more than one random effect, one-dimensional fixed nodes are not by themselves a complete integration strategy.

Approach	Where evaluations are placed	Main consideration
Laplace approximation	Local quadratic approximation near a mode	Fast, but sensitive to how well the local shape represents the integral
Nonadaptive Gauss–Hermite	Fixed nodes centered at zero in standard-normal coordinates	Reusable nodes, but potentially inefficient for displaced or narrow integrands
Adaptive Gauss–Hermite	Nodes adjusted using local posterior shape information	Can target displaced mass, at the cost of locating and scaling around the mode

Common Mistake and Practical Interpretation

A common mistake is to say, “The random effect is centered at zero, so the likelihood’s important region must also be centered at zero.” The Gaussian prior is centered at zero, but the product of that prior and the conditional likelihood need not be. The data in a particular cluster can shift the region contributing most to the marginal integral. The correction is to distinguish the prior center from the shape of the full integrand and to check whether the fixed nodes adequately cover the latter.

A second pitfall is to treat a quadrature result as exact merely because the weights are correct. Gauss–Hermite quadrature is exact for polynomials up to a specified degree under its Gaussian weighting, not for every GLMM likelihood. In practice, a reported marginal likelihood or log likelihood should be understood as a numerical approximation, with its adequacy checked for the model and data at hand.

Key takeaway: Nonadaptive Gauss–Hermite quadrature approximates a one-dimensional Gaussian random-effects integral by a weighted sum at fixed standard-normal nodes. The formula is straightforward; deciding whether a zero-centered node set adequately represents the likelihood-weighted region is the essential numerical judgment.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

For a standard-normal random effect and the three-node rule, what are the nodes and weights? What should the approximation return when the integrand is the constant one?
A random intercept has standard deviation \(\sigma_b\). At which random-effect values should the conditional likelihood be evaluated for standard-normal node \(z_j\)?
For the hypothetical cluster in the binary example, explain why the quadrature sum uses conditional likelihood values without multiplying each one by a Gaussian density value.
What feature of an integrand can make a zero-centered rule inefficient even when the prior random effect is correctly specified as Gaussian?
Why is agreement between two quadrature orders useful but not a mathematical guarantee of accuracy?

Solutions. (1) The nodes are \(0,\pm\sqrt{3}\), with weights \(2/3,1/6,1/6\). Their weights sum to one, so the constant integrand is returned as \(2/3+1/6+1/6=1\). (2) Evaluate at \(b=\sigma_b z_j\), because \(b=\sigma_b Z\) standardizes the Gaussian random effect. (3) The normalized weights already incorporate integration under the standard-normal density; multiplying by that density again would count the weighting twice. (4) A narrow peak or a region of substantial likelihood displaced from zero may be poorly covered by fixed nodes centered at zero. (5) Similar results suggest stability over those orders, but both rules may miss the same important region or share a similar approximation error.