Why Gaussian-Weighted Integration Calls for a Special Rule
A GLMM marginal likelihood requires integrating a conditional likelihood over random effects. The Gaussian random-effects density makes the integrand contain a Gaussian kernel, but the likelihood factor multiplying that kernel is generally not a polynomial. Gauss–Hermite quadrature is designed to integrate functions in precisely this weighted form: it handles the Gaussian kernel through its weights and evaluates the remaining factor at selected nodes.
The rule has two separate advantages. First, it is exact for a broad class of polynomial factors, with a precise degree guarantee. Second, it can approximate smooth nonpolynomial factors efficiently because the nodes are chosen for integration under a Gaussian weight, rather than distributed uniformly over an interval. Exactness for polynomials does not, however, guarantee a particular error for a GLMM likelihood; the approximation must still be assessed.
The Rule and Its Convention
There are two widely used Hermite conventions. This tutorial uses the physicists’ convention, with weight \(e^{-x^2}\). Its Hermite polynomial of degree \(m\) is \(H_m(x)=(-1)^m e^{x^2}\frac{d^m}{dx^m}e^{-x^2}\). The first two are \(H_0(x)=1\) and \(H_1(x)=2x\). Naming the weight is essential: nodes and weights from the alternative convention, which uses \(e^{-x^2/2}\), cannot be interchanged without adjustment.
The Gaussian weight is part of the integral, not part of \(f\) in the displayed approximation: the rule evaluates \(f\) at each node and multiplies by \(w_i\). The weights already account for the factor \(e^{-x^2}\). For instance, a constant factor \(f(x)=1\) integrates to \(\sqrt{\pi}\), so the weights sum to \(\sqrt{\pi}\).
Why the Rule Is Exact Through Degree \(2m-1\)
The key fact is that Hermite polynomials are orthogonal under the weight \(e^{-x^2}\): \(H_m\) has zero weighted integral against every polynomial of degree less than \(m\). The nodes make the quadrature sum vanish for any factor containing \(H_m\), since \(H_m(x_i)=0\). Together, these facts yield exactness beyond degree \(m-1\).
Let \(p\) have degree at most \(2m-1\). Polynomial division gives \(p=qH_m+r\), where \(q\) and \(r\) each have degree at most \(m-1\).
Orthogonality gives \(\int_{-\infty}^{\infty}e^{-x^2}q(x)H_m(x)\,dx=0\), because \(q\) has degree below \(m\).
At every node \(x_i\), \(q(x_i)H_m(x_i)=0\). Thus both the integral and quadrature sum of the product term are zero.
The \(m\) weights are chosen so that the rule integrates the Lagrange interpolation basis at the \(m\) nodes exactly. Since \(r\) has degree at most \(m-1\), its interpolation from those nodes is \(r\) itself. The rule therefore integrates \(r\) exactly, and hence integrates \(p=qH_m+r\) exactly.
This argument proves the degree guarantee. It does not say that a nonpolynomial function is integrated exactly, nor that the error necessarily decreases at a fixed rate as \(m\) increases. For a smooth nonpolynomial factor, the rule approximates the integral by sampling that factor at the nodes; practical accuracy depends on how well those evaluations represent the weighted integral.
Worked Example: Checking a Two-Node Rule on a Polynomial
For \(m=2\), \(H_2(x)=4x^2-2\). Its roots are \(x_1=-1/\sqrt{2}\) and \(x_2=1/\sqrt{2}\). The weight formula gives \(w_1=w_2=\sqrt{\pi}/2\). Use the rule to evaluate \(\int_{-\infty}^{\infty}x^2e^{-x^2}\,dx\).
The integrand factor is \(f(x)=x^2\), a polynomial of degree 2. Since a two-node rule is exact through degree \(2(2)-1=3\), exactness is guaranteed. Substitution gives $$ Q_2=\frac{\sqrt{\pi}}{2}\left[\left(-\frac{1}{\sqrt{2}}\right)^2+\left(\frac{1}{\sqrt{2}}\right)^2\right] =\frac{\sqrt{\pi}}{2}\left(\frac12+\frac12\right) =\frac{\sqrt{\pi}}{2} \approx 0.8862. $$ The exact integral is also \(\sqrt{\pi}/2\), confirming the result.
Using the Gaussian Weight for Random-Effects Integrals
For one random effect \(b\), a cluster’s marginal likelihood contribution has the form \(\int L(b)\phi(b;0,\sigma^2)\,db\), where \(L(b)\) is the cluster’s conditional likelihood and \(\phi(b;0,\sigma^2)\) is the Gaussian density with mean zero and variance \(\sigma^2\). This is a Gaussian-weighted integral: the likelihood factor is multiplied by a Gaussian density. Gauss–Hermite quadrature is well suited because its weights are designed for a Gaussian kernel.
The exact nodes and weights used for a particular random-effects integral depend on writing its Gaussian kernel in the same convention as the rule. In the physicists’ convention here, the weighted integral is \(\int e^{-x^2}f(x)\,dx\). The density normalization and the scale of the random effect must therefore be accounted for when expressing the desired integral in that form. The change of coordinates connecting a general Gaussian random effect to this standard weight is a separate step; using nodes from one convention with a weight from another is not valid.
Worked Example: A Gaussian-Weighted Random-Effects Factor
Suppose a hypothetical one-dimensional calculation has been expressed, with its Gaussian kernel separated, as \(\int_{-\infty}^{\infty}e^{-b^2}L(b)\,db\). At the two nodes \(b=\pm1/\sqrt{2}\), suppose the conditional-likelihood factors are \(L(-1/\sqrt{2})=0.20\) and \(L(1/\sqrt{2})=0.60\). Use the two-node rule to approximate the integral.
The two-node weights are both \(\sqrt{\pi}/2\). Thus $$ Q_2=\frac{\sqrt{\pi}}{2}(0.20)+\frac{\sqrt{\pi}}{2}(0.60) =\frac{\sqrt{\pi}}{2}(0.80) \approx 0.7090. $$ This is an approximation to the stated unnormalized weighted integral, not by itself a normalized probability or a complete marginal likelihood unless all density-normalization and scale factors have already been included in the expression. The example illustrates how the rule combines likelihood values with weights that account for the Gaussian kernel.
Approximation, Error, and Diagnostics
A nonpolynomial likelihood factor generally prevents exactness. A useful way to understand the rule’s limitation is to examine an integrand with a known exact integral. Consider \(f(x)=e^{-0.2x^2}\). Then \(\int e^{-x^2}f(x)\,dx=\int e^{-1.2x^2}\,dx=\sqrt{\pi/1.2}\). The two-node rule evaluates \(f\) at \(x=\pm1/\sqrt{2}\); these two values are both \(e^{-0.1}\).
Worked Example: Measuring Error for a Nonpolynomial Factor
Using the two nodes and weights above, the quadrature value is $$ Q_2=\frac{\sqrt{\pi}}{2}e^{-0.1}+\frac{\sqrt{\pi}}{2}e^{-0.1} =\sqrt{\pi}e^{-0.1} \approx 1.6038. $$
The exact value is $$ I=\sqrt{\frac{\pi}{1.2}}\approx 1.6180. $$ The approximation is below the exact value by about \(1.6180-1.6038=0.0142\). Using unrounded values, the relative error is $$ 100\frac{|I-Q_2|}{I} =100\left(1-\sqrt{1.2}e^{-0.1}\right) \approx 0.880\%. $$ The nonzero error is consistent with the theory: \(e^{-0.2x^2}\) is not a polynomial, so the degree-3 exactness guarantee does not apply.
In practice, compare quadrature estimates at successively higher orders and monitor the quantity that matters, often the log marginal likelihood or a difference in log likelihoods. Convergence should be judged against the required numerical precision, not merely by observing that two displayed values look similar after aggressive rounding. Also check that node evaluations are finite and that the weights, density normalization, and transformation are consistent. For likelihoods that vary sharply or have substantial mass far from the region sampled by the nodes, a seemingly stable low-order estimate can still be misleading; checking additional orders or an independent integration approach is prudent.
How It Relates to Other Integration Approaches
Gauss–Hermite quadrature and the Laplace approximation address the same broad problem—evaluating an integral over random effects—but use different information. Gauss–Hermite quadrature forms a weighted sum of function evaluations at fixed integration nodes and has exactness guarantees for Gaussian-weighted polynomials. The Laplace approximation instead uses a local quadratic approximation to the log integrand around a mode. It can be efficient when that local shape represents the important mass well, but its accuracy can suffer for skewed, sparse-cluster, or multimodal integrands.
A direct numerical integration method can be attractive in one dimension, but may require many evaluations over a sufficiently broad range. Monte Carlo integration uses random draws and is flexible in higher dimensions, but introduces simulation variability. Gauss–Hermite quadrature avoids simulation noise and can be highly effective when the Gaussian weight is appropriate and dimension is modest. Its limitations include sensitivity to the chosen scale and convention, approximation error for nonpolynomial factors, and rapidly increasing computational cost when multidimensional rules are formed by taking products of one-dimensional rules.
Common Pitfalls
- Mixing Hermite conventions. A rule for \(e^{-x^2}\) is not automatically a rule for \(e^{-x^2/2}\). Identify the weight in the integral before using any nodes and weights.
- Forgetting what the weights include. In the stated rule, \(e^{-x_i^2}\) is not multiplied into each term again. The weights are constructed for the Gaussian-weighted integral.
- Claiming exactness for a likelihood. The rule is exact when the factor \(f\) is a polynomial of degree at most \(2m-1\), not for every smooth likelihood factor.
- Dropping density constants. Approximating a kernel integral is not automatically the same as approximating a normalized expectation. Include the Gaussian density’s normalization and any scale factors required by the integral.
- Treating a stable low-order value as proof of accuracy. Numerical convergence checks help detect error, but the degree guarantee alone does not certify accuracy for a nonpolynomial integrand.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- Under the physicists’ convention, what weight is used, and what polynomial determines the nodes of an \(m\)-point rule?
- What is the maximum polynomial degree for which an \(m\)-point Gauss–Hermite rule is guaranteed to be exact?
- Why is it incorrect to multiply each quadrature term by \(e^{-x_i^2}\) again when applying the rule as defined here?
- A two-node rule estimates a Gaussian-weighted integral whose remaining factor is \(x^4\). Is exactness guaranteed? Explain.
- Name two checks that can help reveal a problematic approximation to a GLMM random-effects integral.
Solutions. (1) The weight is \(e^{-x^2}\), and the nodes are the roots of \(H_m(x)\). (2) Degree \(2m-1\). (3) The weights already account for the Gaussian weight; multiplying it again changes the target integral. (4) No. A two-node rule is guaranteed exact only through degree 3, whereas \(x^4\) has degree 4. (5) Examples include comparing estimates at increasing quadrature orders, checking node evaluations for numerical problems, and verifying density normalization and the Gaussian-weight convention.