Blockwise Integration for Nested Random-Effects Models

Why the Random-Effects Structure Matters

A GLMM marginal likelihood integrates over all random effects. In a nested design, however, not every observation shares every random effect: observations within a lower-level unit share that unit’s effect, while different lower-level units may share an upper-level effect. Preserving this structure can replace one very high-dimensional integral with a set of smaller integrals.

The factorization is an algebraic consequence of the model, not a numerical trick. It is valid only when the likelihood and random-effects density have the corresponding nested factorization. In particular, lower-level contributions can be integrated separately conditional on the shared upper-level effect; that does not usually mean they can be marginalized independently after the shared effect is removed.

Definition: Consider top-level clusters indexed by \(j=1,\ldots,J\), lower-level units \(k=1,\ldots,K_j\) within cluster \(j\), and observations \(i\) within unit \(jk\). A nested random-effects model has an upper-level random effect \(b_j\) shared by the cluster and a lower-level random effect \(u_{jk}\) shared by observations in unit \(jk\). Conditional on these effects and fixed parameters, the observation likelihood factors into the contributions of the observations.

Model and Factorization

Write the conditional density of observation \(Y_{jki}\) as \(f(y_{jki}\mid b_j,u_{jk};\theta)\), where \(\theta\) collects the fixed-effect, covariance, and any other model parameters. For example, a GLMM may specify a link-scale predictor \(x_{jki}^{\mathsf T}\beta+z_{jki}^{\mathsf T}b_j+w_{jki}^{\mathsf T}u_{jk}\), with the conditional response distribution determined by the chosen GLM family.

For the basic factorization, assume that the \(b_j\) are independent across top-level clusters, the \(u_{jk}\) are independent across lower-level units and independent of the \(b_j\), and that their densities are \(p_b(b_j;\theta)\) and \(p_u(u_{jk};\theta)\). These densities may be multivariate; covariance among components within a given effect is allowed. Also assume conditional independence of observations given the random effects and parameters. Define the conditional likelihood contribution of lower-level unit \(jk\) by

L_{jk}(b_j,u_{jk};\theta) = \prod_i f(y_{jki}\mid b_j,u_{jk};\theta).

The marginal likelihood contribution from cluster \(j\) is initially an integral over its upper-level effect and all its lower-level effects. Because the integrand is nonnegative and the lower-level densities and likelihood contributions factor across \(k\), the integrations over the \(u_{jk}\) can be performed separately for any fixed \(b_j\).

\begin{aligned} L_j(\theta) &= \int p_b(b_j;\theta) \prod_{k=1}^{K_j} \left\{ \int L_{jk}(b_j,u_{jk};\theta)\,p_u(u_{jk};\theta)\,du_{jk} \right\}db_j\\ &= \int p_b(b_j;\theta)\prod_{k=1}^{K_j}h_{jk}(b_j;\theta)\,db_j,\\ h_{jk}(b_j;\theta) &= \int L_{jk}(b_j,u_{jk};\theta)\,p_u(u_{jk};\theta)\,du_{jk}. \end{aligned}

The full marginal likelihood is \(\prod_{j=1}^J L_j(\theta)\). Thus the inner integral for a lower-level unit depends on the upper-level effect at which it is evaluated. The product of those inner results is then integrated against the upper-level density. The shared \(b_j\) is what links the lower-level blocks within a cluster.

Key condition: Blockwise integration applies when the conditional likelihood and the random-effects density separate into the stated lower-level blocks once the shared upper-level effect is fixed. If lower-level effects are dependent across units even conditional on \(b_j\), the product of separate \(h_{jk}\) integrals does not follow without a further valid factorization of their joint density.

Worked Factorization

Worked Example: Separating Two Nested Contributions

Suppose a hypothetical study has one cluster with two lower-level units. Each unit has observations conditionally independent given a cluster random intercept \(b\) and its own unit random intercept \(u_k\). Let \(b\sim N(0,\tau^2)\), let \(u_1,u_2\) be independent \(N(0,\sigma^2)\), and assume the three effects are mutually independent. For a binary response, suppose \(P(Y_{ki}=1\mid b,u_k)=\operatorname{expit}(x_{ki}^{\mathsf T}\beta+b+u_k)\).

For unit 1, with responses \(y_{1i}\), its conditional likelihood is

L_1(b,u_1;\theta) = \prod_i \left\{\operatorname{expit}(x_{1i}^{\mathsf T}\beta+b+u_1)\right\}^{y_{1i}} \left\{1-\operatorname{expit}(x_{1i}^{\mathsf T}\beta+b+u_1)\right\}^{1-y_{1i}}.

Define \(h_1(b;\theta)=\int L_1(b,u_1;\theta)\phi(u_1;0,\sigma^2)\,du_1\), and define \(h_2\) in the same way using unit 2’s data. Independence of \(u_1\) and \(u_2\) gives the cluster likelihood

L(\theta)=\int \phi(b;0,\tau^2)\,h_1(b;\theta)h_2(b;\theta)\,db.

The integral is not generally \(\{\int\phi(b;0,\tau^2)h_1(b;\theta)\,db\}\{\int\phi(b;0,\tau^2)h_2(b;\theta)\,db\}\). That expression would use two independently integrated cluster effects, although the model assigns the same \(b\) to both units. The distinction is preserved by multiplying \(h_1\) and \(h_2\) at the same value of \(b\), then integrating once.

A Numerical Blockwise Quadrature Calculation

The following deliberately small calculation illustrates the factorization, not a recommended quadrature order for fitting a model. Use a hypothetical cluster with two lower-level units, one binary observation in each, and both responses equal to 1. Let \(b\) and \(u_k\) be independent standard normal random intercepts, with conditional success probability \(\operatorname{expit}(b+u_k)\). Approximate each standard-normal integral using two-node Gauss–Hermite quadrature: nodes \(-1,1\), each with weight \(1/2\).

Worked Example: Two-Node Approximation for Two Successful Units

For a fixed outer node \(b\), the inner likelihood for either unit is its quadrature-averaged success probability:

h(b)\approx \frac12\left\{\operatorname{expit}(b-1)+\operatorname{expit}(b+1)\right\}.

At \(b=-1\), the two linear predictors are \(-2\) and \(0\). Since \(\operatorname{expit}(-2)\approx0.1192029\) and \(\operatorname{expit}(0)=0.5\),

h(-1)\approx \frac{0.1192029+0.5000000}{2}=0.3096015.

At \(b=1\), the predictors are \(0\) and \(2\), so

h(1)\approx\frac{0.5000000+0.8807971}{2}=0.6903985.

Both units contribute \(h(b)\). The outer two-node approximation therefore gives

\begin{aligned} L(\theta) &\approx \frac12\{h(-1)^2+h(1)^2\}\\ &=\frac12\{(0.3096015)^2+(0.6903985)^2\}\\ &\approx\frac12(0.0958531+0.4766501)=0.2862516. \end{aligned}

If each unit were incorrectly marginalized separately over its own cluster effect, each unit’s approximate marginal success probability would be \(\{h(-1)+h(1)\}/2=0.5\), and multiplying these would give \(0.25\). That is not the blockwise approximation for the stated shared-\(b\) model. The difference is a consequence of retaining the common cluster effect, not evidence that either calculation is an exact marginal likelihood: both use a finite quadrature rule.

Why the Computational Cost Can Fall

Suppose the upper-level effect has dimension \(d_b\), each lower-level effect has dimension \(d_u\), and a tensor rule uses \(q\) nodes per dimension. A direct tensor rule over all random effects in one cluster uses \(q^{d_b+K_jd_u}\) joint nodes. Under the block factorization, an outer rule uses \(q^{d_b}\) nodes, and at each outer node each of the \(K_j\) inner integrals uses \(q^{d_u}\) nodes. The approximate number of inner block evaluations is therefore \(q^{d_b}K_jq^{d_u}\), rather than \(q^{d_b+K_jd_u}\). The actual runtime also depends on the cost of evaluating each likelihood contribution and on implementation details.

Worked Example: Comparing Node Counts

Consider a hypothetical cluster with \(K_j=8\) lower-level units, a two-dimensional upper-level effect, a one-dimensional lower-level effect, and \(q=7\) nodes per dimension. A full tensor rule has

7^{2+8}=7^{10}=282{,}475{,}249

joint nodes. With blockwise integration, the outer rule has \(7^2=49\) nodes. At each outer node, eight inner integrals require \(8\times7=56\) lower-level node evaluations, for a total of

7^2\times 8\times 7=49\times56=2{,}744

inner block evaluations. This comparison assumes the same tensor order for every component and counts the evaluations implied by each strategy; it is not a guarantee of a particular runtime or accuracy. The substantial reduction is possible because the inner integrations do not require a joint grid over all eight lower-level effects.

Justification and Practical Assumptions

The mathematical justification can be stated directly. For fixed \(b_j\), the cluster integrand has the form \(p_b(b_j)\prod_k [L_{jk}(b_j,u_{jk})p_u(u_{jk})]\). Integrating with respect to the vector \((u_{j1},\ldots,u_{jK_j})\) gives a product of integrals because each factor depends on only its own \(u_{jk}\). Nonnegative likelihood contributions permit the order of integration to be changed; under the probability model, these integrals are also well-defined. The remaining integral over \(b_j\) is then the outer integral in the factorization.

For numerical blockwise integration, the model assumptions are not enough by themselves. The chosen rule must also resolve each inner and outer integrand adequately. A narrow peak, substantial skewness, or multiple important modes can require more nodes than a coarse rule supplies. Adaptive quadrature may help when its local approximation is appropriate, but it does not remove the need to check numerical convergence.

Check the nesting. Each observation should belong to the lower-level unit and upper-level cluster represented in the model. A lower-level effect must not be shared across clusters if the factorization treats those effects as independent.
Check conditional independence. The observation likelihood must factor as specified given the effects. Residual dependence left after conditioning can invalidate the likelihood product.
Refine the quadrature. Compare the cluster log likelihood at increasing quadrature orders. Material changes indicate that the approximation has not stabilized.
Check the integration target. Inspect modes and curvature in the relevant effects. A rule centered at one mode can miss meaningful mass elsewhere, a concern especially important for multimodal integrands.
Use stable arithmetic. Products of many small contributions can underflow. Evaluate sums of log likelihood contributions and use stable log-sum-exp calculations when combining weighted terms.

Alternatives, Limits, and Common Mistakes

A full tensor-product rule directly integrates over the complete random-effects vector. It is conceptually straightforward but becomes infeasible when many lower-level effects make the dimension large. Blockwise integration exploits the nested conditional factorization to avoid that full joint grid. Sparse-grid methods can reduce node growth in some moderate-dimensional settings, while a Laplace approximation replaces direct numerical integration with a local approximation around a mode. These alternatives have different accuracy and computational trade-offs; none makes an invalid dependence factorization valid.

Blockwise integration is most useful when many lower-level blocks are conditionally independent given a relatively small shared effect. If the shared effect itself is high-dimensional, its outer integral may still be costly. If lower-level effects are dependent across units conditional on the upper-level effect, separate inner integrals are generally not available in this form. A model with effects shared across otherwise separate clusters also requires a different joint integration structure. Factorization reduces computational dimension only to the extent that the model’s dependence structure permits it.

Common mistake: “The lower-level random effects are independent, so I can integrate each unit separately and multiply the marginal likelihoods.” The correction is that they can be integrated separately at a fixed shared upper-level effect. The resulting inner contributions must be multiplied before the single outer integration. Integrating each contribution over a separate copy of the shared effect changes the model.

Check Your Understanding

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

A cluster has three independent lower-level random effects conditional on its shared \(b_j\). What is the correct order of operations for its marginal likelihood?
Suppose \(d_b=1\), \(d_u=2\), \(K_j=4\), and each dimension uses \(q=5\) nodes. Compare the full tensor node count with the blockwise inner evaluation count.
What model feature would prevent the displayed product of separate inner integrals from following directly?
Why is multiplying separately marginalized lower-level likelihoods generally incorrect when the units share an upper-level effect?

Solutions. (1) For each fixed \(b_j\), integrate each lower-level likelihood contribution over its own random effect, multiply the resulting three values, and then integrate that product over \(b_j\). (2) The full tensor rule uses \(5^{1+4(2)}=5^9=1{,}953{,}125\) joint nodes. Blockwise integration uses \(5^1\times4\times5^2=5\times100=500\) inner block evaluations. (3) Conditional dependence among the lower-level random effects, or a likelihood that does not factor across lower-level units given the effects, prevents this direct factorization. (4) It effectively assigns separate upper-level effects to the units, whereas the model specifies one shared effect. The correct calculation multiplies the conditional inner contributions at the same value of that shared effect before integrating it.