Why Compare Integration Methods?
In a nonlinear mixed-effects model, a subject’s random effect is unobserved. The contribution of that subject to the population likelihood is therefore an average of the conditional data likelihood over the random-effects distribution. For a two-compartment pharmacokinetic (PK) model, the conditional likelihood can change sharply or asymmetrically with a random effect such as clearance. The numerical integration method matters: a convenient calculation is not necessarily an accurate one.
This example compares importance sampling with Gauss–Hermite quadrature for the same hypothetical subject. The numerical outputs below are illustrative rounded evaluations of the stated model. The purpose is to show how to calculate and compare the estimates, not to claim that either method is universally preferable.
Model and Marginal Likelihood
Suppose a single hypothetical subject receives an intravenous bolus dose \(D\). Let \(V_c\) and \(V_p\) be central and peripheral volumes, \(CL\) clearance, and \(Q\) intercompartmental clearance. With \(k_{10}=CL/V_c\), \(k_{12}=Q/V_c\), and \(k_{21}=Q/V_p\), the two-compartment concentration in the central compartment is
where \(k_{10}(\eta)=k_{10,\mathrm{pop}}e^\eta\), \(k_{12}\) and \(k_{21}\) are fixed, and
The individual clearance effect is \(\eta\sim N(0,\omega^2)\). Assume conditionally independent additive normal residual errors: \(Y_j\mid\eta\sim N(C(t_j;\eta),\sigma^2)\), where \(Y_j\) is the concentration at time \(t_j\). Under these assumptions the individual marginal likelihood is
where \(\phi\) is the standard normal density and \(m\) is the number of observations. Equivalently, set \(z=\eta/\omega\), so that \(z\sim N(0,1)\) and \(L=E\{L_c(\omega z)\}\). This standardized form is useful for both methods.
Worked Example: A Coarse Gauss–Hermite Estimate
Worked Example: A Coarse Gauss–Hermite Estimate
Consider a hypothetical subject with three concentration measurements at 0.5, 2, and 6 hours after a 100 mg dose. Set \(V_c=10\) L, \(V_p=20\) L, \(CL_{\mathrm{pop}}=3\) L/h, \(Q=2\) L/h, \(\omega=0.35\), and \(\sigma=0.39\) mg/L. Suppose the rounded conditional likelihood evaluations from this model are \(L_c(0)=1.0648\), \(L_c(-\sqrt{3}\omega)=0.3449\), and \(L_c(\sqrt{3}\omega)=0.0159\), in likelihood-density units of \((\mathrm{mg/L})^{-3}\).
For a standard normal expectation, the three-node Gauss–Hermite rule has nodes \(0,-\sqrt{3},\sqrt{3}\) and weights \(2/3,1/6,1/6\). The quadrature estimate is therefore
Suppose a much more refined integration is used as a numerical reference, giving \(L_{\mathrm{ref}}=0.7365\) in the same units. The three-node estimate is higher by \(0.7700-0.7365=0.0335\), about \(0.0335/0.7365=4.55\%\) of the reference. The arithmetic is straightforward, but the discrepancy shows why a small node count should not be trusted without a refinement check.
Importance Sampling: Weighted Likelihood Evaluations
In standardized coordinates, let the target density be \(\phi(z)\), the standard normal density, and choose a proposal density \(q(z)\) that covers the important region of the integrand. The identity underlying importance sampling follows by multiplying and dividing by \(q(z)\):
For independent draws \(z_1,\ldots,z_N\) from \(q\), the estimator is \(\widehat L_{\mathrm{IS}}=N^{-1}\sum_{i=1}^{N}W(z_i)\). It is unbiased for \(L\) when the proposal covers the target integrand and the expectation exists. In particular, \(q(z)\) must be positive wherever \(L_c(\omega z)\phi(z)\) is nonzero. The proposal need not equal the random-effects distribution: it is designed to sample efficiently where the likelihood-weighted density is substantial.
Worked Example: Calculating an Importance-Sampling Estimate
For illustration, suppose a five-draw importance-sampling batch produces the following weighted contributions, calculated as \(W(z_i)=L_c(\omega z_i)\phi(z_i)/q(z_i)\). The first row includes the component values so the weighting calculation is explicit.
| Draw | \(L_c(\omega z_i)\) | \(\phi(z_i)/q(z_i)\) | \(W(z_i)\) |
|---|---|---|---|
| 1 | 1.020 | 0.500 | 0.510 |
| 2 | Weighted contribution from the model evaluation | 0.660 | |
| 3 | Weighted contribution from the model evaluation | 0.740 | |
| 4 | Weighted contribution from the model evaluation | 0.890 | |
| 5 | Weighted contribution from the model evaluation | 0.882 | |
For draw 1, the likelihood contribution is \(1.020\times0.500=0.510\). The batch estimate is
This happens to be close to the reference \(0.7365\), but five draws do not establish accuracy. A small batch can agree by chance, particularly when the weights are variable. Importance sampling must be assessed using its sampling variability and the distribution of its weights, not by agreement of a single estimate alone.
For independent draws with finite weight variance, an estimated Monte Carlo standard error is \(SE(\widehat L_{\mathrm{IS}})=s_W/\sqrt{N}\), where \(s_W\) is the sample standard deviation of the weights. The effective sample size (ESS), often summarized as \((\sum_i W_i)^2/\sum_i W_i^2\), indicates how concentrated the weights are: a low ESS relative to \(N\) signals that a few draws dominate. ESS is a diagnostic, not a proof of accuracy; it can fail to reveal a region that the proposal almost never visits.
Adaptive Quadrature and a Practical Comparison
Ordinary Gauss–Hermite quadrature uses nodes and weights aligned to a Gaussian density. Adaptive quadrature first locates the mode of the subject-specific log integrand, \(\log L_c(\eta)+\log f_\eta(\eta)\), and uses its local curvature to set a location and scale for the nodes. It then evaluates the integrand at transformed nodes and applies the corresponding weights and transformation factors. This can place nodes near the dominant region even when the subject’s likelihood shifts the integrand away from the population random-effect mean.
For the hypothetical example, suppose adaptive quadrature with 15 nodes gives \(0.7368\), and increasing to 31 nodes gives \(0.7365\). The change is \(0.0003\), or about \(0.041\%\) of \(0.7365\). Against the same refined reference, the three-node standard rule gave \(0.7700\). Thus the adaptive calculations are more consistent with the reference in this illustration. The comparison does not establish that 15 nodes will suffice for another subject or parameter setting; node refinement is part of the evidence, not a guarantee.
Worked Example: Interpreting Accuracy and Computation Together
Suppose a computational run for this subject reports the following additional results: an importance-sampling estimate of \(0.7372\) from 20,000 independent draws, an estimated Monte Carlo \(SE\) of \(0.0021\), and ESS \(=8{,}460\). The estimate is \(0.0007\) above the reference \(0.7365\), a relative difference of about \(0.095\%\). Its difference from the reference is about \(0.33\) reported standard errors, since \(0.0007/0.0021\approx0.33\). That is compatible with ordinary Monte Carlo variation; it is not evidence that the reference is exact or that the estimate must be correct.
The ESS is \(8{,}460/20{,}000=42.3\%\) of the draw count. This indicates that weighting has reduced the information compared with equally contributing draws, but the estimate may still be useful. The adaptive 15- and 31-node results differ by only \(0.0003\), while the stochastic estimate has an estimated \(SE\) of \(0.0021\). Quadrature here is more stable at the displayed precision, whereas importance sampling provides a direct Monte Carlo uncertainty estimate and can be attractive when deterministic quadrature becomes costly in higher dimensions. Actual run time depends on the model evaluator, proposal construction, and implementation.
A short justification clarifies why refinement is informative but not conclusive. For a quadrature sequence \(\widehat L_K\) with increasing node count \(K\), agreement between \(\widehat L_K\) and \(\widehat L_{2K}\) means the reported integral is insensitive to that particular refinement. It does not bound the error unless additional error conditions are established. Similarly, an importance-sampling standard error estimates variability under the proposal’s sampled distribution; it cannot detect important regions that the proposal fails to explore. Both methods therefore require diagnostics suited to their distinct sources of error.
| Feature | Importance sampling | Adaptive quadrature |
|---|---|---|
| Evaluation pattern | Random draws and likelihood weights | Deterministic transformed nodes and weights |
| Main error signal | Monte Carlo SE, weight spread, and ESS | Change under node refinement and stability of the mode and curvature |
| Potential advantage | Can be used in dimensions where tensor-product node counts are prohibitive | Can be highly accurate with relatively few nodes for a smooth, concentrated low-dimensional integrand |
| Common vulnerability | A poorly matched proposal can yield a few dominant weights or miss important regions | A poor local approximation or inadequate nodes can miss tails, skewness, or multiple regions |
Diagnostics, Pitfalls, and Limitations
For importance sampling, inspect the weight distribution, ESS, and estimates across independent batches or larger draw counts. A few exceptionally large weights, a very low ESS, or estimates that shift materially as draws are added indicate instability. Check that the proposal is sufficiently broad and that it covers plausible tails as well as the mode. Heavy-tailed or defensive proposals can reduce the risk of severe weight instability, although they do not remove the need for diagnostics.
For adaptive quadrature, increase the node count and compare estimates at the same parameter values. Also check whether the estimated mode and curvature are stable and whether the integrand appears skewed or has more than one important region. Sharp curvature, a boundary-like tail, or a poorly located mode can make a local Gaussian adaptation inadequate. An apparent plateau across a few node counts should be treated cautiously if the likelihood evaluations themselves are unstable.
Common mistake: “The 15-node and 31-node answers are close, so the integral must be correct.” The correction is that agreement is evidence of numerical stability under that refinement, not a universal error bound; inspect the integrand and compare with an independent method when the result matters. Conversely, a large importance-sampling sample is not automatically reliable if a few weights dominate or the proposal misses a relevant region.
Neither method removes model uncertainty or fixes a misspecified PK or residual-error model. Their numerical reliability also depends on the smoothness and dimensionality of the integral and the quality of the model evaluations. In multiple random-effect dimensions, tensor-product quadrature can become expensive rapidly, while importance sampling can suffer weight degeneracy. Numerical integration should therefore be judged in the context of both the likelihood structure and the inferential task.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- Why is the subject likelihood an integral rather than simply \(L_c(0)\)?
- Using the three-node values in the example, what is the quadrature estimate and its absolute difference from the reference?
- What condition on \(q(z)\) is needed for the importance-sampling identity to cover the target integrand?
- What does an ESS of 8,460 from 20,000 draws suggest, and what does it not prove?
- Why does agreement between two node counts not establish a guaranteed quadrature error bound?
Solutions. 1. The random effect is unobserved, so the marginal contribution averages the conditional likelihood over its population distribution. 2. The estimate is \(0.7700\), and the absolute difference from \(0.7365\) is \(0.0335\). 3. The proposal density must be positive wherever \(L_c(\omega z)\phi(z)\) is nonzero, with the required weighted expectation finite. 4. The ESS is 42.3% of the draws, indicating that weights are unequal; it does not prove that the proposal explored every important region or that the estimate is accurate. 5. Refinement agreement shows stability for the tested node counts, but without further error conditions it does not bound error in regions the nodes may still inadequately represent.