Quadrature Order and Likelihood Accuracy

Why Quadrature Order Matters

A marginal likelihood for a model with random effects is obtained by integrating the conditional likelihood over the random-effects distribution. When this integral has no convenient closed form, numerical integration replaces it. The quadrature order determines how many weighted evaluations contribute to that approximation. Too few evaluations can distort the likelihood surface, shifting its maximum even when the optimizer has converged accurately for the approximation it was given.

The central practical question is therefore not just whether estimation converges, but whether the likelihood and parameter estimates are stable as quadrature order increases. This tutorial develops that check in a one-dimensional example, where the exact marginal likelihood is also available for comparison, and then shows why a satisfactory order in one dimension may not be satisfactory in several dimensions.

Definition: For a one-dimensional random effect $\eta\sim N(0,\omega^2)$ and conditional likelihood contribution $L(\theta\mid\eta)$, the marginal likelihood contribution is $$ L(\theta)=\int_{-\infty}^{\infty} L(\theta\mid\eta)\, \frac{1}{\omega\sqrt{2\pi}}\exp\left(-\frac{\eta^2}{2\omega^2}\right)\,d\eta. $$ A quadrature approximation replaces the integral by a weighted finite sum. Its order is the number of nodes used in that sum.

Gauss–Hermite Quadrature for a Normal Random Effect

For a standard normal random variable $Z$, an order-$m$ Gauss–Hermite rule can be written in normal-weight form as

E[g(Z)]\approx \sum_{j=1}^{m} w_j g(z_j), \qquad \sum_{j=1}^{m} w_j=1.

Here $z_j$ are the nodes and $w_j$ are the corresponding normalized weights. For $\eta\sim N(0,\omega^2)$, substitute $\eta=\omega Z$, giving

L(\theta)\approx L_m(\theta) =\sum_{j=1}^{m} w_j L(\theta\mid \omega z_j).

The nodes and weights are chosen to integrate polynomials accurately under the normal weight; an order-$m$ rule is exact for polynomials up to degree $2m-1$. That property does not imply that a general, non-polynomial likelihood contribution is exact. Increasing order usually improves the approximation for well-behaved integrands, but the sequence need not improve monotonically at every order or for every parameter value.

The approximation presumes that the random-effects distribution and conditional likelihood have been specified correctly, the integral is finite, and the quadrature rule is applied with the correct scale and weights. A computation can also be affected by numerical underflow or by a poorly optimized approximate likelihood. Thus a change in results across orders is evidence of integration sensitivity, but not by itself proof that quadrature is the only source of error.

A One-Dimensional Example With an Exact Benchmark

Worked Example: A Poisson Observation With a Normal Random Effect

Suppose a hypothetical subject has observed count $Y=1$. Conditional on a random effect $\eta\sim N(0,1)$, let $Y$ have a Poisson distribution with mean $\theta\eta^2$, where $\theta>0$. Compare the exact marginal likelihood with two low-order quadrature approximations. This is a constructed example for illustrating numerical integration, not a real study.

For a Poisson observation of one, the conditional likelihood is

L(\theta\mid\eta) =P(Y=1\mid\eta,\theta) =\theta\eta^2\exp(-\theta\eta^2).

Because $\eta\sim N(0,1)$, the marginal likelihood is

L(\theta)=\theta E\left[Z^2\exp(-\theta Z^2)\right].

To evaluate it, first use the normal integral $E[\exp(-\theta Z^2)]=(1+2\theta)^{-1/2}$. Differentiating with respect to $\theta$ gives

E[Z^2\exp(-\theta Z^2)] =-\frac{d}{d\theta}(1+2\theta)^{-1/2} =(1+2\theta)^{-3/2}.

Thus $L(\theta)=\theta(1+2\theta)^{-3/2}$. The log likelihood and its derivative are

\ell(\theta)=\log\theta-\frac{3}{2}\log(1+2\theta), \qquad \ell'(\theta)=\frac{1}{\theta}-\frac{3}{1+2\theta}.

Setting the derivative to zero yields $1+2\theta=3\theta$, so the exact maximum occurs at $\hat{\theta}=1$. The derivative is positive for $0<\theta<1$ and negative for $\theta>1$, confirming this is the maximum. The likelihood there is

L(1)=3^{-3/2}=0.1925.

For order 2, the standard-normal nodes are $-1,1$, each with weight $1/2$. Therefore

L_2(\theta) =\frac{1}{2}\theta(1)^2e^{-\theta} +\frac{1}{2}\theta(1)^2e^{-\theta} =\theta e^{-\theta}.

Its log derivative is $1/\theta-1$, so the order-2 approximation also has its maximum at $\hat{\theta}_2=1$. But its likelihood value there is $L_2(1)=e^{-1}=0.3679$, which differs substantially from the exact value $0.1925$. Agreement in the maximizing parameter has not made the likelihood value accurate.

For order 3, the nodes are $0,-\sqrt{3},\sqrt{3}$, with weights $2/3,1/6,1/6$. The zero node contributes zero because of the factor $\eta^2$, and the other two terms give

L_3(\theta) =\frac{1}{6}(3\theta e^{-3\theta}) +\frac{1}{6}(3\theta e^{-3\theta}) =\theta e^{-3\theta}.

The log derivative $1/\theta-3$ is zero at $\hat{\theta}_3=1/3$, with the derivative changing from positive to negative there. The order-3 estimate is therefore $0.3333$, rather than the exact $1$. This example shows that a change in order can affect the estimated parameter as well as the likelihood.

Interpreting Order Comparisons

Quadrature order controls the numerical approximation, not the statistical model. In the example, all three calculations used the same conditional Poisson model and the same random-effects distribution. They differed only in how the integral was approximated. The order-2 result happened to give the exact maximizing $\theta$, but missed the likelihood value; order 3 changed both. One low-order result cannot establish accuracy.

A useful assessment is to refit at a sequence of increasing orders while holding the model, data, parameterization, and other estimation settings fixed. Compare the estimated parameters and the marginal log likelihood. If these quantities move materially as order increases, the lower-order result is integration-sensitive. Stability across several larger orders is reassuring, although it is not a formal proof of accuracy. When feasible, compare with an exact calculation or a more accurate independent integration method.

The likelihood values above also illustrate an important interpretation: a quadrature likelihood is an approximation to the same marginal likelihood, not a new statistical likelihood with a different target. In the example, the order-2 maximum gives $0.3679$, exceeding the exact maximum $0.1925$, while the order-3 maximum gives $\frac{1}{3}e^{-1}=0.1226$. Quadrature approximations are not generally upper or lower bounds. Comparing their numerical values is meaningful as a sensitivity check, but a discrepancy should not be interpreted as evidence that one model fits better.

Worked Example: Comparing Likelihood Values and Maximizers

Using the preceding example, calculate the log likelihood at $\theta=1$ for the exact integral and for orders 2 and 3, then compare the order-specific maxima.

The exact log likelihood at $\theta=1$ is

\ell(1)=\log(3^{-3/2}) =-\frac{3}{2}\log(3) =-1.648.

At order 2, $L_2(1)=e^{-1}$, so $\ell_2(1)=-1.000$. At order 3, $L_3(1)=e^{-3}$, so $\ell_3(1)=-3.000$. Thus the approximate log likelihoods at the same parameter value range from $-3.000$ to $-1.000$, whereas the exact value is $-1.648$.

At their own maxima, order 2 gives $\ell_2(1)=-1.000$, while order 3 gives

\ell_3(1/3)=\log(1/3)-1=-2.099.

The exact maximum is $\ell(1)=-1.648$. The order-2 maximizing parameter agrees with the exact one in this particular example, but its maximum log likelihood does not; the order-3 approximation shifts the estimate and also changes the maximum log likelihood. A numerical optimization algorithm could report convergence for either approximation while neither approximation is sufficiently accurate for the intended use.

Why Dimension Changes the Calculation

In more than one dimension, a common construction applies a one-dimensional rule separately to each random-effect coordinate and forms a tensor product. If there are $d$ independent coordinates and $m$ nodes per coordinate, this requires $m^d$ joint node evaluations per likelihood contribution. For example, order 5 requires $5^2=25$ evaluations in two dimensions but $5^{10}=9{,}765{,}625$ in ten dimensions. This exponential growth is the main computational limitation of tensor-product quadrature.

Moreover, using the same order for each coordinate does not guarantee comparable accuracy as dimension grows. The joint integrand may have substantial curvature, skewness, or dependence across coordinates; accuracy then depends on the shape of the full integrand, not just the one-dimensional order chosen for each marginal coordinate. A rule that adequately resolves one coordinate in a simple model may poorly represent important regions of a higher-dimensional integral.

Worked Example: Two Random Effects and a Tensor-Product Rule

Consider two hypothetical independent observations, each with the conditional model from the earlier example but with its own independent $N(0,1)$ random effect. Both observations equal one, and they share the parameter $\theta$. Independence makes the exact joint marginal likelihood the product of the two one-observation likelihoods:

L(\theta)=\left[\theta(1+2\theta)^{-3/2}\right]^2 =\theta^2(1+2\theta)^{-3}.

Its log derivative is $2/\theta-6/(1+2\theta)$. Equating this to zero gives $2+4\theta=6\theta$, hence the exact maximum is $\hat{\theta}=1$.

Applying order 2 to each random effect yields the product of two order-2 approximations:

L_{2,2}(\theta)=(\theta e^{-\theta})^2 =\theta^2e^{-2\theta}.

The log derivative is $2/\theta-2$, so its maximum is also at $\theta=1$. With order 3 in each dimension, the tensor-product approximation is

L_{3,3}(\theta)=(\theta e^{-3\theta})^2 =\theta^2e^{-6\theta}.

Its log derivative $2/\theta-6$ is zero at $\theta=1/3$, so this approximation again shifts the maximizing parameter. The two-dimensional example shows that tensor-product quadrature can reproduce the one-dimensional sensitivity in a joint likelihood, while the number of evaluations grows multiplicatively with dimension. In a model with dependent random effects or a nonseparable likelihood, the joint integral need not factor into one-dimensional contributions, so each joint node must represent the combined integrand.

Diagnostics, Alternatives, and Limitations

When checking quadrature accuracy, vary order without changing other settings, and confirm that each run reaches a comparable optimization endpoint. Otherwise, apparent order sensitivity may be confounded with differences in starting values or optimizer behavior. Review the estimated parameters and the objective values at the fitted parameters; consider also evaluating the higher-order approximation at a common parameter vector to separate changes in the integral from changes in the maximizer.

Warning signs include material shifts in fixed-effect or variance estimates across orders, large changes in the marginal log likelihood, or estimates that move irregularly as order increases. A smooth progression toward stability is helpful but not required: quadrature error can vary non-monotonically. Extreme parameter estimates, very small or boundary-adjacent variance estimates, and sharply concentrated conditional likelihoods can also make a fixed quadrature grid poorly placed. These are reasons to investigate, not automatic proof that the integration is wrong.

Quadrature is deterministic for fixed nodes, weights, and inputs. Its error is therefore distinct from Monte Carlo error, which arises from random draws; Monte Carlo precision settings address the latter. Importance sampling is another numerical integration approach covered in the earlier material on Importance Sampling for Likelihood Evaluation. These approaches have different computational and error-control characteristics, but none eliminates the need to check whether the likelihood calculation is sufficiently accurate for the model and inferential task.

A frequent mistake is to treat optimizer convergence as evidence that the marginal likelihood is accurate. The optimizer can locate a maximum of $L_m(\theta)$ to high numerical precision even when $L_m(\theta)$ is a poor approximation to $L(\theta)$. The correction is to distinguish optimization accuracy from integration accuracy: check convergence of the optimization, then separately examine sensitivity to quadrature order or compare against an appropriate benchmark. Another mistake is to assume an order validated in one dimension automatically works in higher dimensions. Tensor-product cost and the shape of the joint integrand make that assumption unjustified.

Key takeaway: Quadrature order is an accuracy setting, not merely a computational preference. Assess it by checking whether both the marginal likelihood and its maximizing estimates are stable at higher orders, and reassess the choice when random-effects dimension or integrand shape changes.

Check Your Understanding

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

In the Poisson example, why does the order-3 rule produce a different maximizing $\theta$ from the exact marginal likelihood?
At $\theta=1$, calculate the order-2 likelihood and state how it compares with the exact likelihood.
For a tensor-product rule with 4 nodes per coordinate and $d=3$ random effects, how many joint node evaluations are required per contribution?
Why does optimizer convergence not establish quadrature accuracy?

Solutions.

1. The order-3 rule approximates the integral with its three nodes, producing $L_3(\theta)=\theta e^{-3\theta}$, whose log derivative is $1/\theta-3$. Its maximum is at $\theta=1/3$. The exact integral is $\theta(1+2\theta)^{-3/2}$, whose log derivative is $1/\theta-3/(1+2\theta)$, giving the maximum at $\theta=1$. The approximated integrands, and therefore their likelihood shapes, differ.

2. $L_2(1)=1\cdot e^{-1}=0.3679$. The exact value is $L(1)=3^{-3/2}=0.1925$, so the order-2 value is larger by about $0.1754$. This comparison is a quadrature discrepancy, not evidence that the approximation is an upper bound.

3. The number of evaluations is $4^3=64$ per contribution.

4. Convergence means the optimizer has adequately located a maximum of the objective it was given. If that objective is an inaccurate quadrature approximation, the converged maximum can still differ from the maximum of the true marginal likelihood. Optimization and integration accuracy require separate checks.