Sequential and Nested Quadrature Rules

Why Integrate in Stages?

In a nested GLMM, observations may share more than one level of latent variation. For example, measurements are grouped within patients, and patients are grouped within clinics. A clinic-level random effect is shared by all patients in that clinic, while each patient has a separate patient-level random effect. The marginal likelihood must average over both levels without discarding their dependence structure.

When the model has a conditional structure, the integrations need not be treated as one large, undifferentiated integral. Conditional on the clinic effect, patient contributions can often be integrated separately over their patient effects. The resulting patient-level integrals are then combined and integrated over the clinic effect. This is sequential quadrature: numerical integration follows the factorization of the model.

Definition: Suppose cluster $j$ has an outer random effect $b_j$, and nested unit $k$ has an inner random effect $u_{jk}$. If the $u_{jk}$ are conditionally independent given $b_j$, and observations are conditionally independent given all random effects, the cluster likelihood contribution is $$ L_j=\int \left\{\prod_{k=1}^{K_j} \left[\int L_{jk}(b_j,u_{jk})f(u_{jk}\mid b_j)\,du_{jk}\right]\right\}f(b_j)\,db_j, $$ where $L_{jk}(b_j,u_{jk})$ is the conditional likelihood for the observations in nested unit $k$. The inner integral averages over one unit’s random effect, conditional on the outer effect; the outer integral then averages the combined contribution over the outer effect.

This formula is not merely a computational rearrangement. It encodes which latent effects are shared and which are conditionally independent. If the inner effects are independent of $b_j$, then $f(u_{jk}\mid b_j)=f(u_{jk})$. If their conditional distribution depends on $b_j$, that dependence must remain in the inner integration. Either way, the product over $k$ is justified only when the inner effects and their observations separate conditionally on $b_j$.

From the Integral to a Sequential Quadrature Sum

For a standard-normal random vector $Z$, quadrature approximates an expectation $E\{g(Z)\}$ by weighted function evaluations. The nodes and weights may be fixed or adapted to the integrand. In sequential quadrature, first approximate each inner expectation at a given outer node; then use those inner approximations to evaluate the outer expectation.

\begin{aligned} I_{jk}(b)&=\int L_{jk}(b,u)f(u\mid b)\,du,\\ L_j&=\int \left\{\prod_{k=1}^{K_j} I_{jk}(b)\right\} f(b)\,db,\\ \widehat L_j&=\sum_{r=1}^{Q_b} w_r \prod_{k=1}^{K_j}\left\{\sum_{s=1}^{Q_u}v_{s\mid r} L_{jk}(b_r,u_{s\mid r})\right\}. \end{aligned}

Here $b_r$ and $w_r$ are outer quadrature nodes and weights, and $u_{s\mid r}$ and $v_{s\mid r}$ are the inner nodes and weights used at outer node $b_r$. When the inner distribution is conditional on $b$, its nodes, weights, or both may need to change with $r$. The weights in this display are understood to approximate expectations under the relevant random-effects distributions; a density should not be multiplied in again if it is already incorporated into the quadrature rule.

The order of operations matters: at each outer node, calculate each nested unit’s inner integral, multiply those unit contributions, and only then average across outer nodes. A full tensor rule over every random effect would instead enumerate joint combinations of outer and all inner nodes. Sequential integration exploits conditional factorization: with $K_j$ scalar inner effects, it uses roughly $Q_b K_j Q_u$ evaluations rather than $Q_b Q_u^{K_j}$ combinations.

Worked Example: Two Successes in Nested Units

Consider a hypothetical clinic with two patients, each contributing one binary response $Y_k$. Let $b\sim N(0,1)$ be the shared clinic intercept and let $u_k\sim N(0,1)$ be independent patient intercepts, also independent of $b$. Given the random effects, suppose $Y_k$ is Bernoulli with success probability $\operatorname{expit}(b+u_k)$, where $\operatorname{expit}(x)=1/(1+e^{-x})$. Both observed responses are successes. Use two-node Gauss-Hermite quadrature for each standard-normal expectation: nodes $-1,1$, each with weight $0.5$.

The assumptions needed here are built into the hypothetical model: the responses are conditionally independent given the effects, and the two patient effects are independent given the clinic effect. For a fixed outer value $b$, the inner success probability is approximated by $$ \widehat m(b)=0.5\{\operatorname{expit}(b-1)+\operatorname{expit}(b+1)\}. $$ At the outer nodes, $\operatorname{expit}(-2)=0.1192$, $\operatorname{expit}(0)=0.5000$, and $\operatorname{expit}(2)=0.8808$. Thus $$ \widehat m(-1)=0.5(0.1192+0.5000)=0.3096,\qquad \widehat m(1)=0.5(0.5000+0.8808)=0.6904. $$

Conditional independence makes the clinic likelihood at fixed $b$ the product of the two identical inner contributions, $\widehat m(b)^2$. The outer quadrature therefore gives $$ \widehat L=0.5\{0.3096^2+0.6904^2\} =0.5(0.09585+0.4767)\approx 0.2863. $$

This is a two-node approximation to the marginal probability of observing two successes in this clinic. It is not the square of the marginal success probability: the shared $b$ induces dependence between the two responses after averaging over random effects.

A Short Justification of the Factorization

The sequential expression follows by conditioning on the outer effect and using the model’s conditional independence assumptions. For clarity, write the joint random-effects density as $f(b,\mathbf{u})$, where $\mathbf{u}=(u_1,\ldots,u_K)$, and suppose $f(b,\mathbf{u})=f(b)\prod_k f(u_k\mid b)$. Conditional independence gives a likelihood $\prod_k L_k(b,u_k)$. Substitution into the marginal likelihood yields

\begin{aligned} L &=\int\int \left\{\prod_{k=1}^{K}L_k(b,u_k)\right\} \left\{\prod_{k=1}^{K}f(u_k\mid b)\right\}f(b)\,d\mathbf{u}\,db\\ &=\int\left\{\prod_{k=1}^{K} \int L_k(b,u_k)f(u_k\mid b)\,du_k\right\}f(b)\,db. \end{aligned}

The second equality separates the integrations over $u_k$ because, at fixed $b$, each factor depends on only its own $u_k$. This is an application of iterated integration to a factorized integrand. If the conditional density of the inner effects does not factor, or if a likelihood term jointly involves multiple inner effects, this separation is not valid as written; the relevant joint block must be integrated together.

Worked Example: One Success and One Failure

Keep the model and quadrature rule from the preceding example, but suppose the first patient has a success and the second has a failure. For fixed $b$, let $m_1(b)$ denote the inner probability of success and $m_0(b)$ the inner probability of failure. Since the inner distribution is symmetric and the two-node rule is symmetric, $m_0(b)=1-m_1(b)$ for this calculation.

The values from the previous calculation give $m_1(-1)=0.3096$ and $m_1(1)=0.6904$, so $m_0(-1)=0.6904$ and $m_0(1)=0.3096$. Integrate the product over the shared clinic effect: $$ \widehat L=0.5\{0.3096(0.6904)+0.6904(0.3096)\} =0.3096(0.6904)\approx 0.2137. $$

A tempting but incorrect shortcut is to marginalize each patient separately and multiply the two results. Each marginal success or failure probability is $0.5$ under this symmetric approximation, so that shortcut gives $0.5(0.5)=0.2500$, not $0.2137$. It omits the shared clinic effect: the patient outcomes are conditionally independent given all effects, but not generally independent after integrating out a common effect.

Conditional Dependence and Numerical Choices

Sequential integration is especially useful when many inner units share a lower-dimensional outer effect. It does not eliminate the need to integrate over the full model; rather, it avoids evaluating every combination of inner nodes when conditional factorization makes that unnecessary. When the inner effects are correlated with $b$, their conditional distribution must be represented correctly. For example, if $u\mid b$ is multivariate normal with mean $\mu(b)$ and covariance $\Sigma(b)$, one can express $u=\mu(b)+A(b)z$, where $z$ is standard normal and $A(b)A(b)^\mathsf{T}=\Sigma(b)$, then apply inner quadrature in the $z$ coordinates for each outer node. For non-Gaussian conditional distributions, use a transformation or quadrature rule appropriate to the actual conditional distribution.

Adaptive quadrature can be used at both levels. The inner nodes may be centered and scaled for the conditional integrand at each $b_r$; the outer nodes may be adapted to the outer integrand, which includes the product of inner integrals. Adaptation does not change the target likelihood. It changes where the finite set of nodes is placed to approximate that likelihood more accurately. A rule adapted for one outer value should not silently be reused as though it were an exact rule for a different conditional distribution.

Worked Example: A Shared Random Intercept Versus Separate Effects

Suppose, hypothetically, that two binary observations are conditionally independent successes given a single shared random intercept $b\sim N(0,1)$, with success probability $\operatorname{expit}(b)$. There is no inner effect in this example; it isolates the outer integration that also appears in a nested calculation. Under the two-node rule, the probability of two successes is approximated by $$ \widehat L=0.5\{\operatorname{expit}(-1)^2+\operatorname{expit}(1)^2\}. $$

The node probabilities are $\operatorname{expit}(-1)=0.2689$ and $\operatorname{expit}(1)=0.7311$. Substitution gives $$ \widehat L=0.5(0.2689^2+0.7311^2) =0.5(0.07231+0.5345)\approx 0.3034. $$

If the shared intercept were incorrectly integrated separately for the two observations, the approximation would instead be $\{0.5(0.2689+0.7311)\}^2=0.2500$. The difference, $0.3034$ versus $0.2500$, reflects the shared latent effect, not a change in the conditional Bernoulli model.

Diagnostics, Alternatives, and Limitations

A numerical answer should be checked at both integration levels. Increase the inner order while holding the outer order fixed, then increase the outer order while holding the inner order fixed. If either refinement materially changes the marginal likelihood or the estimates obtained from it, the original rule was not yet stable. Also check that inner weights correspond to the conditional distribution at each outer node, that expectation weights are normalized as intended, and that the integrand is evaluated on the log scale when products of many small contributions risk underflow.

Factorization failure: If inner effects remain dependent given $b$, separate inner integrals cannot be multiplied. Integrate the dependent block jointly or use a valid conditional factorization.
Dependence erased: Multiplying marginal unit likelihoods can be wrong when units share an outer effect. Integrate the product conditional on that effect, then average over it.
Insufficient inner order: A good outer approximation cannot repair inaccurate inner integrals. Refine the inner rule separately.
Misapplied conditional nodes: If $f(u\mid b)$ changes with $b$, treating it as a fixed unconditional density can target the wrong integral.

Approach	Strength for nested models	Main caution
Sequential quadrature	Uses conditional factorization to avoid enumerating every combination of inner nodes.	Requires the factorization and conditional distributions to be specified correctly; inner errors propagate outward.
Full tensor-product quadrature	Directly evaluates combinations of nodes across all random effects.	Node counts grow rapidly with the number of effects and nested units.
Sparse-grid quadrature	Can reduce node growth for suitable moderate-dimensional integrands.	It does not automatically exploit a nested conditional factorization; accuracy still needs assessment.
Laplace approximation	Can be substantially less costly when the relevant integrand is well represented locally by a Gaussian shape.	It is a different approximation and may be less accurate when the integrand is skewed or otherwise poorly captured locally.

Sequential quadrature is not automatically exact or efficient in every model. Its cost can still be substantial when each inner effect is high-dimensional, when the conditional integrals are sharply peaked, or when many outer nodes require distinct inner calculations. It is most compelling when the model’s conditional structure yields a genuine reduction in dimension or node combinations. The central practical rule is to preserve the model’s dependence structure first, then assess numerical convergence at each integration level.

Common mistake: “The nested units are conditionally independent, so I can integrate each one out separately and multiply their marginal likelihoods.” The correction is that conditional independence permits multiplying the conditional inner integrals at a fixed outer effect. The product must still be integrated over the shared outer effect: $$ \int \prod_k I_k(b)f(b)\,db \ne \prod_k\int I_k(b)f(b)\,db $$ in general.

Check Your Understanding

Answer each question, then compare with the solutions below.

What assumptions permit the product of inner integrals in the nested likelihood?
Why must the product over nested units be formed before integrating over a shared outer effect?
If $u\mid b$ changes with $b$, what must be reflected in the inner quadrature?
In the two-success worked example, which refinement should be checked if the inner success probabilities are inaccurate?

Solutions. (1) The inner random effects must be conditionally independent given the outer effect, and the observations must have a likelihood that factors by nested unit conditional on the effects. (2) The shared effect induces dependence after marginalization, so in general the integral of a product differs from the product of integrals. (3) The conditional density, or an equivalent transformation that represents it, must be used at each outer node; fixed unconditional inner nodes and weights are not justified unless they correctly represent that conditional distribution. (4) Increase the inner quadrature order while initially holding the outer order fixed, then assess outer-order sensitivity separately.