Numerical Integration Settings for Random Effects

Why Integration Settings Matter

In a nonlinear mixed-effects model, the likelihood for a subject is obtained by averaging the conditional likelihood over that subject’s random effects. That average is generally an integral, and the model’s likelihood is often not available in closed form. An integration setting therefore affects both the numerical approximation to the likelihood and the amount of computation required at each parameter update.

The key practical issue is not simply that one algorithm is “more accurate” than another. Integration dimension, the shape of the conditional likelihood across random effects, and how the integration points are placed all matter. A rule that works well for a smooth, nearly Gaussian one-dimensional integrand may be inefficient or inadequate when several random effects interact or the likelihood is sharply concentrated away from its prior center.

Definition: Let $y_i$ denote the observations for subject $i$, let $\eta_i$ be a vector of $q$ random effects, and let $\theta$ collect model parameters. If $\eta_i$ has density $p(\eta_i\mid\theta)$, the subject-level marginal likelihood contribution is $$ L_i(\theta)=\int p(y_i\mid\eta_i,\theta)\,p(\eta_i\mid\theta)\,d\eta_i. $$ The integration dimension is $q$, the number of components of $\eta_i$ integrated over. For example, if $\eta_i\sim N(0,\Omega)$, the integral is over the $q$-dimensional random-effects distribution.

This expression assumes that the conditional observation model and the random-effects distribution have been specified, and that the integral exists. In the common formulation where subjects are independent given the population parameters, the full likelihood is the product of the subject-level contributions. The numerical integration is performed for each subject and is repeated as parameter values change during estimation.

Fixed and Adaptive Quadrature

Quadrature approximates an integral by evaluating its integrand at selected points and combining those evaluations with weights. Fixed, nonadaptive quadrature places the points at predetermined locations, often according to a Gaussian quadrature rule matched to the assumed random-effects distribution. The same point pattern is used regardless of the observed data for a subject.

In $q$ dimensions, a common construction forms a tensor product of a one-dimensional rule. If that rule uses $m$ points per dimension, the full rule evaluates the integrand at $m^q$ combinations of points. This can work well when the integral is smooth and its important region is reasonably covered by points near the random-effects distribution’s center. Its computational cost, however, grows exponentially in $q$.

Adaptive quadrature changes the point placement to reflect the subject-specific integrand. A typical procedure locates a mode of the conditional kernel $p(y_i\mid\eta_i,\theta)p(\eta_i\mid\theta)$, estimates its local curvature, and centers and scales the quadrature rule around that region. This can place points where the likelihood contributes most to the integral, rather than spending many points in regions with negligible contribution. Adaptation requires additional calculations and may be unreliable if the kernel is strongly skewed, has multiple important modes, or has poorly estimated curvature.

Worked Example: Fixed Quadrature Misses a Concentrated Region

Suppose a hypothetical one-dimensional subject contribution, apart from any parameter-free multiplier, is $L=E[\exp(-\eta^2)]$, where $\eta\sim N(0,1)$. Compare the exact integral with a three-point Gaussian rule centered on the random-effects distribution.

For a standard normal random effect, the three points are $0$ and $\pm\sqrt{3}$, with weights $2/3$ and $1/6$ at each of the two nonzero points. Thus the quadrature approximation is $$ \widehat L_{\mathrm{fixed}} =\frac{2}{3}\exp(0)+\frac{1}{6}\exp(-3)+\frac{1}{6}\exp(-3) =\frac{2}{3}+\frac{1}{3}\exp(-3). $$

Since $\exp(-3)\approx0.04979$, this gives $\widehat L_{\mathrm{fixed}}\approx0.6833$. The integral can be evaluated exactly by combining the normal density and the exponential term: $$ L=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\exp\left(-\frac{3}{2}\eta^2\right)d\eta =\frac{1}{\sqrt{3}}\approx0.5774. $$

The approximation is high by about $0.1059$, or $0.1059/0.5774\approx18.3\%$ relative to the exact value. The integrand is concentrated more tightly around zero than the original standard normal density. This example shows that a fixed rule’s points can be poorly matched to the shape of the integrand, even in one dimension. It does not establish that fixed quadrature will always have this direction or magnitude of error.

How Dimension Changes Computational Effort

The number of random effects is consequential even if the model otherwise remains unchanged. Suppose each one-dimensional rule uses $m=7$ points. A tensor-product rule for $q=3$ random effects requires $7^3=343$ integrand evaluations per subject. For $q=6$, it requires $7^6=117{,}649$, or $117{,}649/343=343$ times as many. With 500 subjects, those counts become 171,500 and 58,824,500 evaluations per objective-function calculation, respectively.

These counts measure integrand evaluations, not total run time: evaluating the model at each point can itself be costly, and implementation details affect overhead. Still, the arithmetic explains why simply increasing the number of points in every dimension can become impractical. Correlation among random effects may also make the important region tilted relative to the original coordinate axes, so an apparently generous grid can still allocate points inefficiently.

Worked Example: Estimating the Cost of a Tensor Rule

Consider a hypothetical model with $q=4$ random effects and 300 subjects. A tensor-product rule using $m=5$ points in each dimension requires $$ m^q=5^4=625 $$ integrand evaluations per subject, or $300\times625=187{,}500$ evaluations for one objective-function calculation.

If the model instead has $q=7$ random effects and the same five points per dimension, the count is $5^7=78{,}125$ per subject. For 300 subjects that is $300\times78{,}125=23{,}437{,}500$ evaluations. The ratio is $23{,}437{,}500/187{,}500=125$, equal to $5^{7-4}$. Thus adding three dimensions multiplies the cost by 125 under this tensor construction. This is a planning calculation, not a prediction that run time will increase by exactly the same factor.

Local Gaussian Approximation and Other Strategies

A related method is the Laplace approximation. It is useful to understand its connection to adaptive integration, but it is an approximation based on local curvature rather than a quadrature rule that sums evaluations at a chosen set of nodes. Define the negative log kernel $h_i(\eta)=-\log\{p(y_i\mid\eta,\theta)p(\eta\mid\theta)\}$. If $h_i$ has a well-defined interior minimum $\hat\eta_i$, and its Hessian $H_i$ there is positive definite, a second-order Taylor expansion gives a local quadratic approximation.

Derivation: Near $\hat\eta_i$, the gradient of $h_i$ is zero. Retaining terms through second order gives $h_i(\eta)\approx h_i(\hat\eta_i)+\frac{1}{2}(\eta-\hat\eta_i)^T H_i(\eta-\hat\eta_i)$. Exponentiating the negative of this expression and integrating the resulting Gaussian kernel yields $$ L_i(\theta)\approx \exp\{-h_i(\hat\eta_i)\}(2\pi)^{q/2}|H_i|^{-1/2}. $$ The approximation follows from replacing the full log kernel by its local quadratic form; it is most credible when the important region is well represented by that form.

Laplace approximation can be computationally attractive in higher dimensions because it avoids a full $m^q$ grid, though finding the mode and computing the curvature still require work. It can be inaccurate when the kernel is highly skewed, has substantial non-Gaussian tails, or has multiple important modes. Adaptive quadrature also uses location and curvature information, but retains a weighted integration calculation around the adapted region. A local Gaussian approximation and adaptive quadrature should therefore not be treated as interchangeable.

Sparse-grid quadrature combines selected points from lower-dimensional rules rather than evaluating every combination in a full tensor product. It can reduce the point count for smooth integrands whose effective dimension is modest, but the gain depends on the integrand’s structure and the construction used. Simulation-based methods, including importance sampling, estimate the integral using draws from a proposal distribution and weights. They avoid a deterministic tensor grid, but introduce sampling variability and depend on an adequate proposal distribution. The earlier material on importance sampling covers that approach in more detail.

Worked Example: Why Adaptation Can Help After a Shift in the Kernel

Consider a hypothetical integrand formed from a standard normal random-effects density and a conditional likelihood proportional to $\exp[-2(\eta-2)^2]$: $$ L=E\left[\exp\{-2(\eta-2)^2\}\right] =\frac{1}{\sqrt{2\pi}}\int\exp\left\{-\frac{1}{2}\eta^2-2(\eta-2)^2\right\}d\eta. $$

The negative log kernel, ignoring a constant, is $h(\eta)=\frac{1}{2}\eta^2+2(\eta-2)^2$. Its derivative is $h'(\eta)=\eta+4(\eta-2)$, so the mode satisfies $5\eta-8=0$, giving $\hat\eta=1.6$. Its curvature is $h''(\eta)=5$, so the local standard deviation implied by the quadratic kernel is $1/\sqrt{5}\approx0.4472$. The important region is therefore shifted from the prior center at zero and is relatively narrow.

Completing the square gives $h(\eta)=\frac{5}{2}(\eta-1.6)^2+1.6$. Hence $$ L=\frac{\exp(-1.6)}{\sqrt{5}}\approx\frac{0.2019}{2.236}=0.09030. $$ An adaptive rule can center and scale its points around $1.6$ with local scale $1/\sqrt{5}$, matching the location and curvature of this integrand. A fixed rule centered at zero does not make that adjustment. In this particular example the kernel is exactly quadratic, so its local Gaussian representation is exact; real subject-level kernels need not be quadratic, and adaptation alone does not guarantee accuracy.

Choosing a Strategy and Diagnosing Problems

Integration settings should be chosen with the model’s random-effects dimension and the shape of its subject-level kernels in mind. Fixed quadrature is a natural candidate when dimension is low and the integrand is smooth and well covered by points centered on the assumed random-effects distribution. Adaptive quadrature is worth considering when subject-specific kernels are concentrated or shifted, provided reliable modes and curvature can be found. Sparse grids may help with smooth problems of moderate effective dimension. Laplace approximation may offer a computationally economical local approximation when the kernels are close to Gaussian. Simulation-based integration is an alternative when deterministic grids become unwieldy, with Monte Carlo precision and proposal adequacy as additional concerns.

Useful warning signs include large changes in parameter estimates or objective values when the integration strategy changes, unstable subject-level contributions, failure to find a mode, or computational effort that grows sharply as random effects are added. Such signs do not by themselves identify the cause: the model may be poorly conditioned, the integration may be inadequate, or both. Compare methods under the same model and data, and record the settings so that the comparison is reproducible.

A common mistake is to assume that a large number of integration points automatically solves every problem. In a tensor rule, adding points in each dimension can multiply the total work dramatically, while still placing them poorly for a shifted or irregular kernel. The correction is to consider both where the integrand is important and how many dimensions must be integrated. Changing to an adaptive or alternative strategy can address placement or cost, but its assumptions and failure modes also need checking.

No integration strategy removes the need to inspect the fitted model and its numerical behavior. A local approximation can miss non-Gaussian features; adaptive procedures can be misled by unreliable curvature; sparse-grid savings depend on smoothness and effective dimension; and simulation introduces sampling error. Integration is a computational component of estimation, not evidence by itself that the model is scientifically appropriate.

Check Your Understanding

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

A tensor rule uses six points per dimension. How many integrand evaluations are required for one subject with five random effects?
Why might adaptive quadrature be preferable to fixed quadrature when a subject’s conditional kernel is sharply concentrated away from the random-effects mean?
State one condition that supports the Laplace approximation and one kernel feature that can make it unreliable.
Give one advantage and one limitation of simulation-based integration compared with a deterministic tensor rule.

Solutions. (1) The count is $6^5=7{,}776$ evaluations per subject. (2) Adaptation can center and scale the integration points around the subject-specific region that contributes most to the integral; this avoids relying on a grid centered only on the random-effects mean. Its advantage depends on obtaining a reliable mode and curvature. (3) A well-defined interior mode with positive-definite curvature and a kernel close to quadratic in its important region supports the approximation. Strong skewness or multiple important modes can make a single local quadratic approximation unreliable. (4) Simulation avoids the exponential point count of a tensor grid as dimension grows, but introduces sampling variability and can perform poorly if its proposal distribution does not cover important regions of the integrand.