Worked Strategy for Selecting an Integration Method

Why Integration Method Choice Matters

A GLMM likelihood averages the conditional likelihood over unobserved random effects. That average is an integral, and in many models it has no closed-form expression. The numerical method used to approximate it can affect the fitted likelihood, parameter estimates, and uncertainty calculations. A method that is accurate but too costly to run is not useful; a fast method that gives unstable or materially inaccurate estimates is not a sound substitute.

There is no universally best integration method. A useful choice depends on the dimension of the random-effects vector, the shape of the integrand, the accuracy required for the analysis, and the available computational budget. The earlier material on Laplace approximation, adaptive quadrature, order refinement, and integration diagnostics supplies the main tools; this tutorial organizes them into a selection strategy.

Definition: For cluster \(i\), let \(y_i\) denote its response data and \(b_i\) its \(q\)-dimensional random-effects vector. Write \(h_i(b_i)=\log p(y_i\mid b_i;\beta)+\log p(b_i;\theta)\), where \(\beta\) denotes fixed-effect parameters and \(\theta\) denotes random-effect distribution parameters. The cluster marginal likelihood contribution is \(L_i(\beta,\theta)=\int_{\mathbb R^q}\exp\{h_i(b_i)\}\,db_i\).

This expression presumes a specified conditional response model and random-effects distribution; in the usual GLMM setting the random effects are multivariate normal, with their covariance structure included in \(p(b_i;\theta)\). To write the full likelihood as \(\prod_i L_i\), two kinds of factorization are needed: responses from different clusters are conditionally independent given the random effects, and the random-effects distribution factorizes across clusters (for example, the \(b_i\) are mutually independent). Conditional independence of responses alone does not imply that the marginal likelihood factors: dependent random effects across clusters can couple their integrals.

The dimension \(q\) is the dimension of the random-effects vector for one cluster, not the number of observations in that cluster. A model with a random intercept and random slope has \(q=2\), even if each cluster contains many observations. Computational costs are often multiplied by the number of clusters, so both \(q\) and the number of clusters matter.

What the Main Choices Trade Off

Laplace approximation uses the height and curvature of \(h_i\) at its conditional mode. If \(\hat b_i\) is an interior mode and \(A_i=-\nabla^2 h_i(\hat b_i)\) is positive definite, then

L_i \approx \exp\{h_i(\hat b_i)\}(2\pi)^{q/2}|A_i|^{-1/2}.

This is usually computationally economical, especially when \(q\) is moderate or large. Its central limitation is that it represents the integrand locally by a Gaussian shape. Skewness, heavy tails, multiple modes, or appreciable non-Gaussian behavior away from the mode can make that local approximation inaccurate.

Adaptive Gauss–Hermite quadrature (AGHQ) places quadrature nodes around the mode and scales them using local curvature. In one dimension, if \(A_i=-h_i''(\hat b_i)>0\), the change of variable \(b_i=\hat b_i+z/\sqrt{A_i}\) gives

L_i=A_i^{-1/2}\int_{\mathbb R} \exp\left\{h_i\left(\hat b_i+\frac{z}{\sqrt{A_i}}\right)+z^2\right\} e^{-z^2}\,dz.

Gauss–Hermite nodes and weights approximate the integral with the weight \(e^{-z^2}\). With enough nodes, AGHQ can capture departures from the local Gaussian shape better than Laplace, but its cost rises with both the number of nodes and the random-effects dimension.

A tensor-product rule uses combinations of one-dimensional nodes in all \(q\) dimensions. With \(m\) nodes per dimension, it requires \(m^q\) node combinations per cluster. This can be effective for small \(q\), but the exponential growth makes it quickly impractical. Sparse-grid quadrature combines lower-dimensional rules to reduce the number of nodes. It can be attractive for smooth integrands whose important variation is concentrated in relatively few directions. It is not a general escape from high-dimensional integration: strong interactions among many directions or irregular integrands can undermine its efficiency or accuracy.

A Short Derivation: Why Laplace Is a Local Choice

Suppose \(h_i\) is smooth, has a unique interior maximum at \(\hat b_i\), and has negative-definite Hessian there. A second-order Taylor expansion around the mode is

h_i(b_i)\approx h_i(\hat b_i) -\frac{1}{2}(b_i-\hat b_i)^\mathsf{T}A_i(b_i-\hat b_i), \qquad A_i=-\nabla^2 h_i(\hat b_i).

Substituting this approximation into the integral and using the multivariate Gaussian integral gives

\begin{aligned} L_i &\approx \exp\{h_i(\hat b_i)\} \int_{\mathbb R^q} \exp\left\{-\frac{1}{2}(b_i-\hat b_i)^\mathsf{T} A_i(b_i-\hat b_i)\right\}\,db_i\\ &=\exp\{h_i(\hat b_i)\}(2\pi)^{q/2}|A_i|^{-1/2}. \end{aligned}

The approximation is justified when the region contributing substantially to the integral is well represented by this quadratic expansion. A well-defined mode and positive curvature matrix are necessary checks, but they do not prove that the quadratic approximation is accurate across the important integration region. That is why comparison with higher-order or alternative quadrature, when feasible, is valuable.

A Stepwise Selection Strategy

Confirm the integral and its factorization.
Specify the random-effects vector and its dimension. Verify conditional response independence and independence of cluster-specific random effects before treating the likelihood as a product of cluster contributions.

Check the mode and local curvature.
Use mode and Hessian diagnostics. A poorly located mode, a nearly singular curvature matrix, or materially different modes from different starts is a warning against trusting a simple local approximation.

Start with the least costly plausible method.
Laplace is a natural candidate when the integrand is locally close to Gaussian or \(q\) makes direct quadrature expensive. For a low-dimensional integral where accuracy is important, use AGHQ or tensor-product quadrature as a comparison.

Match quadrature structure to dimension.
For \(q=1\), order-refined AGHQ is often straightforward. For small \(q\), tensor-product rules can be feasible. Consider sparse grids at larger \(q\) when the integrand is smooth and relatively low-dimensional in its effective variation.

Refine and compare at the analysis target.
Increase quadrature order or sparse-grid level and compare log likelihoods, parameter estimates, and standard errors. Decide in advance what numerical change would be negligible for the scientific conclusions.

Worked Decisions

Worked Example: A One-Dimensional Binary Random-Intercept Model

Suppose a hypothetical study has 80 independent clusters, each with a binary response and a random intercept, so \(q=1\). A fitted model has a unique conditional mode in every cluster, and the negative second derivative is positive at each mode. The main goal is stable fixed-effect estimates, and the integrand is suspected to be somewhat skewed because several clusters have few events.

Laplace is computationally attractive, but the skewness concern makes it prudent to check an adaptive rule. Suppose an order-refinement run reports total log likelihoods of \(-128.436\), \(-128.421\), and \(-128.420\) at orders 7, 11, and 15, respectively. The successive changes are \( |-128.421-(-128.436)|=0.015\) and \( |-128.420-(-128.421)|=0.001\).

The last change is much smaller than the preceding one, suggesting numerical stabilization at the likelihood level. Suppose further that the largest fixed-effect estimate change from order 11 to 15 is \(0.002\) on the log-odds scale, against a prespecified tolerance of \(0.01\). These results support using order 15 AGHQ for this fit, subject to checking parameter standard errors and convergence as well. They do not establish that the model itself is adequate, nor do they prove that every possible quadrature error is below \(0.001\).

The choice is reasonable because \(q=1\) makes AGHQ inexpensive, while the shape concern gives a reason to validate Laplace rather than assume it is adequate. If the order-15 and order-11 estimates had materially differed, the next step would be further refinement or investigation of the modes and integrand, not automatic acceptance of the cheaper answer.

Worked Example: Four Random-Effect Dimensions

Consider a hypothetical model with random intercepts and slopes giving \(q=4\). A tensor-product rule with 9 nodes per dimension would require \(9^4=9\times9\times9\times9=6{,}561\) node combinations per cluster. With 120 clusters, that is \(6{,}561\times120=787{,}320\) cluster-node evaluations for one likelihood evaluation, before accounting for repeated evaluations during parameter optimization.

Suppose an implementation of a nested sparse-grid rule uses 321 distinct nodes at one selected level and reports total log likelihoods \(-93.782\), \(-93.771\), and \(-93.770\) over three successive levels. The final change is \(0.001\), and the previous change is \(0.011\). If fixed-effect estimates and standard errors also stabilize within their prespecified tolerances, the sparse-grid calculation is a plausible choice: it uses fewer evaluations than the 9-node tensor rule while showing level refinement.

The count 321 belongs to this particular grid construction and level; sparse-grid node counts depend on the rule and convention. A sound decision also checks smoothness and potential interactions among the four effects. If estimates shift substantially at the next level, the lower node count alone is not evidence of sufficient accuracy.

Worked Example: Ten Random-Effect Dimensions and a Tight Budget

A hypothetical model has \(q=10\), and the analyst considers a 7-node tensor rule. Its node count per cluster would be \(7^{10}=282{,}475{,}249\). For 50 clusters, that is \(282{,}475{,}249\times50=14{,}123{,}762{,}450\) cluster-node evaluations for one likelihood evaluation. This is generally prohibitive when the likelihood must be evaluated repeatedly during fitting.

Laplace is a sensible initial candidate if each conditional mode is stable and the Hessian is well conditioned. Suppose a sparse-grid calculation on a representative subset of clusters is feasible and suggests that likelihood and parameter estimates agree closely with Laplace, while increasing the grid level changes the result little. That comparison supports using Laplace for the full fit under a limited budget, with the subset and refinement checks reported as evidence about numerical accuracy.

The argument would be weaker if the subset showed substantial disagreement, or if the integrand had unstable modes or pronounced non-Gaussian features. In that case, one should seek a more accurate feasible calculation, simplify only with scientific justification, or report that numerical accuracy could not be established. High dimension does not make an inaccurate approximation accurate.

Diagnostics, Alternatives, and Limits

Integration diagnostics should address more than the final log likelihood. Compare refined fits using the quantities that matter to the analysis: fixed-effect estimates, variance parameters, standard errors, and any reported contrasts. A small change in the total log likelihood can coexist with a consequential change in a variance estimate, especially near a boundary. Also examine whether conditional modes are reproducible from different starting values and whether the negative Hessian is positive definite and numerically well conditioned.

For a low-dimensional problem, direct tensor-product quadrature is a transparent reference when its cost is manageable. AGHQ is often a better practical choice when local adaptation is useful and the dimension is small. Sparse grids can reduce node counts for smooth integrands, but performance depends on the specific grid and the way variation is distributed across dimensions. Laplace is typically cheaper, but it relies more strongly on a local quadratic representation. Simulation-based integration is another broad alternative when deterministic quadrature is impractical, although its Monte Carlo error and computational demands require their own assessment.

These methods address numerical integration, not misspecification. Agreement among methods does not show that the chosen link, conditional response distribution, or random-effects structure is scientifically appropriate. Conversely, disagreement is a warning to investigate numerical behavior and integrand shape; it does not by itself identify which method is correct. Numerical tolerances should be chosen in relation to the inferential purpose rather than by a universal cutoff.

Common Mistakes

Assuming conditional response independence guarantees a product likelihood. It does not if random effects are dependent across clusters. Check that the random-effects distribution also factorizes across clusters.
Choosing a method from \(q\) alone. Dimension affects cost, but integrand shape and required accuracy also matter.
Treating a stable log likelihood as the only accuracy check. Parameter estimates and standard errors can change even when the likelihood changes little.
Assuming sparse grids always solve high-dimensional problems. Their efficiency depends on smoothness and the distribution of variation across dimensions.
Equating numerical agreement with model validity. Integration methods can agree while the statistical model remains inadequate.

Key takeaway: Choose the least costly method that is credible for the integral’s dimension and shape, then demonstrate stability at the likelihood and inferential quantities that matter. Verify the assumptions behind cluster factorization before applying cluster-by-cluster integration.

Check Your Understanding

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

A model has a random intercept and random slope in each cluster. What is \(q\), and how many nodes per cluster does a 6-node tensor-product rule require?
Responses are conditionally independent across clusters given all random effects, but the random effects are correlated across clusters. Can the marginal likelihood generally be written as a product of cluster-specific integrals? Explain.
A scalar AGHQ fit changes its log likelihood by \(0.0005\) after order refinement, but a variance estimate changes materially. Is the integration result established as adequate?
When might a sparse grid be considered instead of a tensor-product rule, and what must still be checked?

Solutions. (1) There are two random effects, so \(q=2\); the rule uses \(6^2=36\) nodes per cluster. (2) No. Conditional response independence is not sufficient: dependence in the joint random-effects distribution generally couples the cluster integrals. (3) No. The likelihood stabilization alone is insufficient when an inferentially important variance estimate remains unstable; refine further or investigate the source of the change. (4) A sparse grid may be considered when the dimension makes a full tensor product costly and the integrand is smooth with variation concentrated in relatively few directions. Check level refinement, parameter stability, and whether interactions or irregular shape undermine that efficiency.