Numerical Integration for Individual Forecast Distributions

Why Numerical Integration Matters

An individual forecast combines uncertainty about that individual’s pharmacokinetic parameters with uncertainty in a future concentration conditional on those parameters. Even when the model is simple, this combination can produce integrals without a convenient closed form. A calculation at a single “best” parameter value generally does not equal an average over the individual’s parameter uncertainty, and nonlinear concentration predictions make that distinction important.

This tutorial compares three ways to calculate individual predictive means and quantiles: posterior simulation, numerical quadrature, and local Gaussian approximation. The focus is computation, not choosing a point forecast or constructing prediction intervals. Those summaries begin with the same predictive distribution, but they are not interchangeable: each numerical method has different accuracy conditions and failure modes.

Definition: Let $H_i$ denote the information available for individual $i$, including the relevant population model, covariates, and that individual’s observed concentrations. Let $\theta_i$ be the individual parameters and $Y^*$ a future concentration at a specified time and dosing history. The individual predictive distribution is $$ p(y^*\mid H_i)=\int p(y^*\mid \theta_i,H_i)\,p(\theta_i\mid H_i)\,d\theta_i. $$ Its mean is the integral of the conditional mean, and its cumulative distribution function (CDF) is the integral of the conditional CDF. A predictive quantile is obtained by inverting that CDF.

For example, if the forecast concerns the latent concentration rather than a future measured value, then $Y^*=g(\theta_i)$, where $g$ is the pharmacokinetic prediction at the forecast time. In that case the predictive distribution is the distribution of $g(\theta_i)$ under $p(\theta_i\mid H_i)$. For a measured future concentration, the conditional distribution $p(y^*\mid\theta_i,H_i)$ also includes residual variability. One must decide which quantity is being forecast before calculating any summary.

Three Ways to Evaluate the Predictive Distribution

Posterior simulation draws parameter values from $p(\theta_i\mid H_i)$, then draws a future concentration from its conditional distribution. The resulting forecast draws approximate the full predictive distribution. Their sample mean estimates the predictive mean, and their empirical CDF estimates predictive quantiles.

Quadrature approximates integrals directly by evaluating the integrand at selected parameter values and combining those evaluations with weights. It can be highly efficient when the effective integration dimension is small and the integrand is smooth. To calculate a predictive CDF, for example, quadrature integrates the conditional CDF over the parameter distribution; the resulting CDF can then be numerically inverted.

Local Gaussian approximation replaces the posterior near its mode or another appropriate center with a multivariate normal approximation, and often linearizes the concentration function near that center. This can yield fast approximate means, variances, and quantiles. Its accuracy depends on whether the posterior and the transformed concentration distribution are sufficiently close to the local approximation.

E(Y^*\mid H_i)=\int E(Y^*\mid\theta_i,H_i)\,p(\theta_i\mid H_i)\,d\theta_i, \qquad F(y\mid H_i)=\int F(y\mid\theta_i,H_i)\,p(\theta_i\mid H_i)\,d\theta_i.

These equations clarify a useful computational choice. For a mean, integrate the conditional mean rather than repeatedly constructing the entire predictive CDF, if the conditional mean is available. For a quantile, integrate conditional CDF values and solve $F(y\mid H_i)=q$, or simulate predictive draws and calculate the empirical $q$ quantile. Averaging conditional quantiles is generally not a valid substitute: in general, the quantile of a mixture is not the weighted average of its component quantiles.

Worked Examples

Worked Example: Local Gaussian Approximation for a One-Compartment Forecast

Suppose a hypothetical individual receives an intravenous bolus of $D=500$ mg. The one-compartment concentration function is $C(t)=D/V\exp(-kt)$, with volume $V$ in L and elimination rate constant $k$ in h$^{-1}$. At $t=4$ h, suppose the posterior for $(\log V,\log k)$ is approximated by independent normal distributions centered at $(\log 50,\log 0.10)$, with standard deviations 0.10 and 0.20, respectively. Ignore residual error here and forecast the latent concentration.

At the center of this approximation, $$ C(4)=\frac{500}{50}\exp(-0.10\times4)=10\exp(-0.4)=6.7032\ {\rm mg/L}. $$ To obtain a first-order approximation, write $a=\log V$ and $b=\log k$. The concentration function is $C= D\exp(-a)\exp(-t\exp(b))$, so its derivatives at the center are $\partial C/\partial a=-C$ and $\partial C/\partial b=-ktC=-0.4C$.

Under the stated independent Gaussian approximation, the delta-method variance is the sum of the squared derivatives multiplied by their respective variances: $$ \operatorname{Var}(C\mid H_i)\mathrel{\approx} C^2\left[(-1)^2(0.10)^2+(-0.4)^2(0.20)^2\right] =(6.7032)^2(0.0164)\mathrel{\approx}0.7369\ ({\rm mg/L})^2. $$ The approximate standard deviation is $\sqrt{0.7369}=0.8584$ mg/L. A local normal approximation for concentration therefore has mean 6.7032 mg/L and standard deviation 0.8584 mg/L. Its approximate 90% quantiles are $6.7032\mathbin{\pm}1.645(0.8584)$, or about 5.291 and 8.115 mg/L.

This calculation assumes the posterior is reasonably represented by a Gaussian in the log parameters and that a first-order expansion is adequate across the posterior’s relevant range. The estimated quantiles are approximations, not exact quantiles of the nonlinear transformation. In particular, a normal approximation on the concentration scale allows negative values, whereas this one-compartment latent concentration is positive.

Worked Example: Weighted Quadrature for Concentration Quantiles

For a second hypothetical forecast, let $D=1000$ mg, $V=50$ L, and $t=4$ h. Suppose the individual posterior for $k$ is represented by a four-node weighted rule. The nodes are chosen so the one-compartment concentration predictions at those nodes are 4, 6, 8, and 12 mg/L, respectively, with weights 0.10, 0.40, 0.30, and 0.20. These are positive weights summing to one; this finite rule illustrates the weighted calculation. The corresponding rate constants $k=\log(20/C)/4$ are approximately 0.4024, 0.3010, 0.2291, and 0.1277 h$^{-1}$, so each node is consistent with $C(0)=D/V=20$ mg/L and the one-compartment model. Forecast the latent concentration.

The quadrature estimate of the predictive mean is the weighted sum of the node predictions: $$ E(C(4)\mid H_i)\mathrel{\approx}0.10(4)+0.40(6)+0.30(8)+0.20(12) =0.4+2.4+2.4+2.4=7.60\ {\rm mg/L}. $$

For this finite rule, the approximate CDF is the sum of weights at nodes no greater than the proposed concentration. It is 0.10 at 4 mg/L, 0.50 at 6 mg/L, and 0.80 at 8 mg/L. Using the quantile definition $Q(q)=\inf\{c:F(c)\ge q\}$, the median is 6 mg/L and the 80th percentile is 8 mg/L. This example shows why the weighted mean and quantiles answer different questions: the mean is 7.60 mg/L, while the median is 6 mg/L.

A practical quadrature rule for a continuous posterior would use more carefully selected nodes and weights than this illustrative finite rule. The weighted-sum logic is the same. If future residual error is included, the CDF calculation must also average the conditional CDF at each parameter node; the latent-concentration quantiles above would then no longer describe the measured future concentration.

Worked Example: Posterior Simulation and Monte Carlo Summaries

A third hypothetical one-compartment forecast uses posterior simulation to include both parameter uncertainty and future residual variability. Suppose 10 simulated future concentrations at a fixed forecast time, in mg/L, are $4.2,\ 4.8,\ 5.0,\ 5.4,\ 5.7,\ 6.1,\ 6.5,\ 7.2,\ 8.0,\ 9.1$. These deliberately small numbers of draws make the arithmetic transparent; they are not enough to support a precise tail estimate.

The simulated mean is $$ \bar{Y}^*=\frac{4.2+4.8+5.0+5.4+5.7+6.1+6.5+7.2+8.0+9.1}{10} =\frac{62.0}{10}=6.20\ {\rm mg/L}. $$ For the empirical CDF, each draw contributes weight $1/10$. The empirical median, using the midpoint of the two central order statistics for this even sample, is $(5.7+6.1)/2=5.90$ mg/L. The empirical 80th percentile under the inverse empirical-CDF convention is the eighth ordered value, 7.20 mg/L.

These are estimates based on simulated draws, not exact values of the individual predictive summaries. With only 10 draws, the mean and especially the quantiles may vary substantially across simulation runs. Increasing the number of effectively independent draws reduces simulation error, but does not correct a misspecified model or a simulation algorithm that fails to represent the posterior.

Why the Methods Work—and Where They Can Fail

For posterior simulation, the justification is the law of large numbers: if $Y_1^*,\ldots,Y_M^*$ are independent draws from the predictive distribution and have finite variance, their sample mean converges to the predictive mean as $M$ increases. The simulation standard error of that sample mean is estimated by the sample standard deviation of the draws divided by $\sqrt{M}$. Predictive quantiles are estimated from the empirical CDF, which converges to the predictive CDF. For correlated draws, as in a Markov chain, correlation can make the effective sample size differ from the raw draw count; it is often smaller for positively autocorrelated draws and should be considered when judging Monte Carlo precision.

Quadrature is justified by numerical approximation of the defining integrals. If the posterior density and integrand are smooth, a well-chosen rule can achieve high accuracy with relatively few evaluations. But accuracy can deteriorate when the posterior is concentrated in a region poorly covered by the nodes, the integrand changes rapidly, or the integration dimension is high. Comparing successively finer rules, changing the integration range, or checking against posterior simulation can reveal numerical sensitivity; agreement is evidence of computational stability, not proof that the model is appropriate.

A local Gaussian approximation is fastest when the posterior is unimodal and approximately symmetric near its center, and when the concentration function changes nearly linearly across the region carrying most posterior mass. It can fail with skewness, heavy tails, boundaries, multiple modes, strong parameter dependence not captured by the approximation, or a strongly nonlinear concentration transformation. The local Gaussian example also makes a specific diagnostic visible: compare the approximate distribution with positive-valued posterior simulations. If a substantial part of a normal approximation lies below zero, its shape is not a credible representation of concentration uncertainty.

Derivation: For a differentiable concentration function $g(\theta_i)$, expand it around a posterior center $\theta_0$: $$ g(\theta_i)\mathrel{\approx}g(\theta_0)+\nabla g(\theta_0)^\mathsf{T}(\theta_i-\theta_0). $$ If $\theta_i\mid H_i$ is approximately $N(\theta_0,\Sigma)$, the linear term has approximate variance $\nabla g(\theta_0)^\mathsf{T}\Sigma\nabla g(\theta_0)$. This yields the local Gaussian variance calculation used above. The derivation depends on the first-order expansion being adequate over the posterior’s spread; it does not establish that condition by itself.

The methods have different computational trade-offs. Simulation is flexible for complicated conditional distributions, nonlinear models, and joint parameter-residual uncertainty, but quantile precision can require many draws. Quadrature can be more efficient for low-dimensional, smooth integrals and can calculate a CDF without adding simulation noise, but its cost grows quickly with dimension and its node placement matters. A local Gaussian approximation is inexpensive and useful as a benchmark or rapid calculation, but can hide skewness and other non-Gaussian features. When feasible, comparing a local approximation with quadrature or simulation is a practical check, especially if the forecast is sensitive to the upper or lower tail.

No numerical method can repair an incorrect forecast target, an inappropriate individual posterior, or a misspecified concentration or residual model. Always state whether the result concerns latent concentration or a future measurement, retain units, and distinguish numerical approximation error from uncertainty represented by the predictive distribution. A stable algorithm can still calculate the wrong target accurately.

Common Mistakes

Plugging in a posterior center and calling it the predictive mean. In general, $E\{g(\theta_i)\mid H_i\}\ne g\{E(\theta_i\mid H_i)\}$ for nonlinear $g$. Integrate over the posterior or explicitly label a plug-in calculation as an approximation.
Averaging component quantiles. A mixture quantile is determined by the mixture CDF, not generally by the average of component quantiles. Integrate CDF values or use predictive draws.
Treating a large raw simulation count as proof of precision. Correlation can make the effective number of draws much smaller, and tail quantiles need adequate representation of the tails.
Using a local normal approximation without checking the transformed scale. A Gaussian parameter approximation does not imply a Gaussian concentration forecast. Inspect the nonlinear transformation and compare methods when possible.

Check Your Understanding

Answer each question, then compare with the short solutions below.

Why does a predictive mean generally require averaging conditional means over the individual posterior?
For a weighted quadrature rule with values 3 and 9 mg/L and weights 0.25 and 0.75, calculate the weighted mean and median using the inverse-CDF quantile definition.
What is one condition supporting a local Gaussian approximation, and one sign that it may be unsuitable for a concentration forecast?
How does adding future residual variability change the calculation for a predictive quantile?
What does increasing the number of independent posterior predictive simulations improve, and what does it not fix?

Solutions. (1) Nonlinearity means $E\{g(\theta_i)\mid H_i\}$ need not equal $g\{E(\theta_i\mid H_i)\}$, so uncertainty in the parameters must be integrated. (2) The mean is $0.25(3)+0.75(9)=7.50$ mg/L. The CDF reaches 0.25 at 3 and 1 at 9, so the median is 9 mg/L. (3) A unimodal, approximately symmetric posterior and a nearly linear transformation over its important region support the approximation; substantial negative mass in a normal concentration approximation is a warning. (4) Integrate the conditional CDF, including residual variability, over the parameter posterior before finding the quantile. (5) More independent draws generally reduce Monte Carlo error, but do not correct model misspecification or a faulty posterior simulation algorithm.