Tensor-Product Quadrature for Multiple Random Effects

Why a Tensor Product Is Needed

A GLMM with two random effects requires integrating a conditional likelihood over both effects. Even if each one-dimensional integral is familiar, the integrand generally depends jointly on the two effects, so evaluating two separate one-dimensional approximations is not enough. A tensor-product rule supplies a grid of paired evaluation points, with one coordinate from each one-dimensional rule.

The construction below uses independent standard-normal coordinates. This covers models in which the two random effects are independent after standardization. A correlated pair requires a suitable transformation to independent coordinates before applying this rule; the formulas here do not themselves perform that transformation.

From a One-Dimensional Rule to a Bivariate Rule

A common starting point is the physicists’ Gauss–Hermite rule, which approximates integrals weighted by \(e^{-x^2}\). If its nodes and weights are \(x_j\) and \(w_j\), then the standard-normal expectation rule has nodes \(z_j=\sqrt{2}x_j\) and normalized weights \(a_j=w_j/\sqrt{\pi}\). Thus \(\sum_j a_j g(z_j)\) approximates \(E\{g(U)\}\) for \(U\sim N(0,1)\), and the weights satisfy \(\sum_j a_j=1\).

For two independent standard normals, the joint density factors into the product of the two standard-normal densities. Consequently, the expectation is a two-dimensional integral with a product density. Replacing each coordinate’s integral by the same one-dimensional rule gives the tensor formula.

The product weight \(a_j a_k\) is essential: it is the probability weight assigned by the product rule to that node pair. The weights across the entire grid sum to one, since \(\sum_j\sum_k a_j a_k=(\sum_j a_j)(\sum_k a_k)=1\). If the original random effects are independent but not standard normal, write them as \(B_1=\mu_1+\sigma_1U\) and \(B_2=\mu_2+\sigma_2V\), and evaluate the conditional integrand at those transformed values. The quadrature weights remain the standard-normal weights.

Applying the Grid to a GLMM Marginal Likelihood

For a cluster with observed responses \(y\), let \(L(b_1,b_2\mid y)\) be its conditional likelihood given the random effects. If the two random effects are independent standard normals, the cluster’s marginal likelihood is the expectation of that conditional likelihood over their joint distribution. Tensor quadrature evaluates the conditional likelihood at every paired node and averages those values using the product weights. With nonstandard independent effects, substitute their means and scales at each node as described above.

For example, in a binary-logit model, a conditional Bernoulli likelihood contribution for response \(y_r\) and linear predictor \(\eta_r\) is \(p_r^{y_r}(1-p_r)^{1-y_r}\), where \(p_r=\operatorname{expit}(\eta_r)=1/(1+e^{-\eta_r})\). The cluster likelihood is the product of the contributions, conditional on the shared random effects. The random effects are then integrated out; they are not replaced by their mean before computing that product.

Unequal Weights and a Normalization Check

Tensor weights need not all be equal. For the three-node standard-normal rule, the nodes are \(0,\sqrt{3},-\sqrt{3}\), with respective normalized weights \(2/3,1/6,1/6\). A paired node’s weight is still the product of its coordinate weights. For example, the pair \((0,\sqrt{3})\) has weight \((2/3)(1/6)=1/9\), whereas \((\sqrt{3},-\sqrt{3})\) has weight \((1/6)(1/6)=1/36\).

Assumptions, Checks, and Alternatives

The tensor formula applies directly when the transformed coordinates being integrated are independent standard normals and the target expectation exists. In a GLMM marginal likelihood, the integrand is the conditional likelihood as a function of those coordinates. The likelihood must be evaluated consistently at each node pair, including all responses in the cluster that share those effects. If a covariance structure makes the random effects dependent, applying independent product weights to the original correlated coordinates is not justified; first use an appropriate transformation to independent coordinates.

A practical diagnostic is to recompute the same integral at successively higher quadrature orders and compare the marginal likelihood, or the log marginal likelihood, at the precision needed for the analysis. Agreement is evidence of numerical stability, not a proof of exactness. Large changes with order, sensitivity to small changes in parameters, or non-finite conditional likelihood calculations warrant investigation. In particular, check that both node coordinates are substituted into every linear predictor, that all \(m^2\) pairs are included, and that the normalized weights sum to one. When likelihoods are extremely small, summing their logarithms with a stable log-sum-exp calculation can help avoid numerical underflow.

Tensor-product quadrature is a direct extension of one-dimensional Gauss–Hermite quadrature: it is systematic and easy to verify, but its accuracy depends on the chosen order and the integrand. Adaptive quadrature instead chooses nodes and weights in relation to the integrand’s local shape; tensor products can also be formed from such one-dimensional rules. The Laplace approximation is a different alternative, using local curvature near a mode rather than a weighted grid of function evaluations. These methods are not interchangeable by definition: compare numerical results when accuracy matters, and check that the integration approximation is stable enough for the intended inference.

A common mistake is to use physicists’ Hermite weights \(w_j\) directly as if they were normalized standard-normal weights, or to multiply the bivariate sum by another normalizing constant after already using \(a_j=w_j/\sqrt{\pi}\). That mixes conventions and changes the integral’s scale. Another error is to use weight \(a_j+a_k\) for a node pair; independence requires the product \(a_ja_k\). Write down the one-dimensional convention first, form product weights, and verify that the full grid weights sum to one.