Decomposing Forecast Uncertainty Across Sources

Why Decompose Forecast Uncertainty?

A future concentration forecast can be uncertain for several different reasons. The individual’s pharmacokinetic parameters may not be precisely known from their observed concentrations. The population parameters used to describe those individual parameters may also be estimated with uncertainty. Even if the parameters were known, residual variability can make a future observation differ from its model-predicted value. Finally, a future input—such as adherence, dose timing, or a covariate value—may itself be uncertain.

These sources have different meanings and may call for different action. More individual concentration measurements can reduce uncertainty about an individual’s parameters, while additional population data can improve population parameter estimates. Residual variability is not removed merely by learning the individual’s parameters, and an uncertain future input may require scenarios or better information about what will happen. Decomposition makes clear which uncertainty is being represented, without assuming that all sources can be reduced in the same way.

Definition: Let \(Y\) be a future observation for an individual, and let \(D\) denote the information already available, including the population data and that individual’s observations. Forecast uncertainty is the conditional variability of \(Y\) given \(D\). Its sources may include uncertainty in population parameters \(\theta\), uncertainty in the individual’s latent parameters \(\eta_i\), residual error \(\epsilon\), and uncertain future inputs \(U\). Individual parameter uncertainty means uncertainty about this individual’s latent parameters; it is not the same as the variation among parameters of different individuals.

A Nested Variance Decomposition

Write the future observation as \(Y=f(\theta,\eta_i,U)+\epsilon\). Here \(f\) is the model-predicted value before residual error, and \(\epsilon\) is the residual term. The relevant distributions are conditional on the available information: population parameter uncertainty is represented by a distribution for \(\theta\), and the individual parameter distribution is conditional on \(\theta\) and that individual’s data. The future input \(U\) can also have a conditional distribution. These distributions might come from a Bayesian analysis or, with an appropriate interpretation, a sampling or bootstrap approximation.

Assume that the residual error has conditional mean zero given \(\theta,\eta_i,U,D\), and has conditional variance \(\sigma_\epsilon^2(\theta,\eta_i,U)\). If the residual mean is not zero, include that mean in \(f\) before applying the decomposition. Define the input-averaged model prediction

\mu(\theta,\eta_i) = E\!\left[f(\theta,\eta_i,U)\mid \theta,\eta_i,D\right].

The law of total variance first separates residual error from variation in the model-predicted value. Applying it again, in the order population parameters, individual parameters, then future inputs, gives:

\begin{aligned} \operatorname{Var}(Y\mid D) ={}& \underbrace{E\!\left[\sigma_\epsilon^2(\theta,\eta_i,U)\mid D\right]}_{V_{\mathrm{res}}}\\ &+ \underbrace{\operatorname{Var}_{\theta\mid D} \left(E_{\eta_i\mid\theta,D}[\mu(\theta,\eta_i)]\right)}_{V_{\mathrm{pop}}}\\ &+ \underbrace{E_{\theta\mid D} \left[\operatorname{Var}_{\eta_i\mid\theta,D} \left(\mu(\theta,\eta_i)\right)\right]}_{V_{\mathrm{ind}}}\\ &+ \underbrace{E_{\theta,\eta_i\mid D} \left[\operatorname{Var}_{U\mid\theta,\eta_i,D} \left(f(\theta,\eta_i,U)\right)\right]}_{V_{\mathrm{input}}}. \end{aligned}

This identity does not require a linear model or independent sources. Each conditional distribution must, however, represent the intended forecasting information and dependencies. The ordering specifies how variation is allocated: for a nonlinear model with interactions, a different ordering can assign a different share to each named component, even though the total variance remains the same. Report the ordering and conditional structure when component values are important.

Key distinction: Between-individual variability describes how latent parameters differ across individuals in the population. Individual parameter uncertainty describes how uncertain we are about one particular individual’s latent parameters after considering their data. A forecast for a new individual and a forecast for an observed individual therefore use different conditional distributions for \(\eta_i\).

Worked Calculations

Worked Example: Combining Four Variance Components

Suppose a hypothetical NLME forecast for one individual has a mean of \(8.20\) mg/L. A nested calculation, using the order above, gives residual, individual-parameter, population-parameter, and future-input variances of \(0.64\), \(0.49\), \(0.16\), and \(0.25\) \((\mathrm{mg/L})^2\), respectively. Assume the resulting predictive distribution is adequately approximated by a normal distribution.

The total predictive variance is the sum of the four components:

\operatorname{Var}(Y\mid D)=0.64+0.49+0.16+0.25 =1.54\;(\mathrm{mg/L})^2.

Thus the predictive standard deviation is \(\sqrt{1.54}=1.241\) mg/L. Under the normal approximation, a central 95% prediction interval is the forecast mean plus or minus \(1.96\) predictive standard deviations:

8.20\pm1.96\sqrt{1.54} =8.20\pm2.4323 =[5.768,\;10.632]\ \mathrm{mg/L}.

The interval describes uncertainty about a future observation under the stated model and information. It includes all four sources. It is not an interval for the population mean alone, and it is not a claim that there is a 95% probability that a fixed parameter lies in this particular interval.

Worked Example: Nonlinearity and the Allocation of Uncertainty

Consider a simplified hypothetical model prediction \(f(\theta,\eta,U)=\theta\eta U\), with no residual error. Let the uncertain population multiplier \(\theta\) be equally likely to be \(2\) or \(4\), the individual multiplier \(\eta\) equally likely to be \(1\) or \(3\) conditional on \(\theta\), and the future input \(U\) equally likely to be \(1\) or \(2\) conditional on \(\theta,\eta\). These are illustrative distributions, not data from a real study. We decompose in the order \(\theta,\eta,U\).

First, \(E(U)=1.5\), so \(\mu(\theta,\eta)=1.5\theta\eta\). For \(\theta=2\), averaging over \(\eta=1,3\) gives \(E_{\eta}[\mu\mid\theta=2]=6\); for \(\theta=4\), it gives \(12\). These two values are equally likely, with overall mean \(9\), so \(V_{\mathrm{pop}}=[(6-9)^2+(12-9)^2]/2=9\).

For fixed \(\theta\), the two values of \(\mu\) across \(\eta\) are \(1.5\theta\) and \(4.5\theta\). Their variance is \((1.5\theta)^2\), since their difference is \(3\theta\) and they are equally likely. Averaging over \(\theta\), whose squared values have mean \((2^2+4^2)/2=10\), gives \(V_{\mathrm{ind}}=2.25(10)=22.5\).

For fixed \(\theta,\eta\), the two predictions across \(U\) are \(\theta\eta\) and \(2\theta\eta\). Their variance is \((\theta\eta)^2/4\). Since \(E(\theta^2)=10\) and \(E(\eta^2)=(1^2+3^2)/2=5\), the input component is \(V_{\mathrm{input}}=(10)(5)/4=12.5\). The total variance is therefore \(9+22.5+12.5=44\), in squared units of the model prediction. Directly, \(E(Y^2)=E(\theta^2)E(\eta^2)E(U^2)=10(5)(2.5)=125\), and \(E(Y)=9\), so \(\operatorname{Var}(Y)=125-9^2=44\), confirming the decomposition.

The calculation illustrates why one cannot generally evaluate a nonlinear model only at average inputs: \(E(\theta\eta U)\) and the variability around it depend on the distributions of all three factors. It also illustrates why the stated order matters when interpreting component shares, even though the total is unchanged.

Worked Example: The Effect of Omitting Population Parameter Uncertainty

Use the four-component forecast from the first example, but suppose an analyst fixes the population parameters at their estimates and omits the \(0.16\;(\mathrm{mg/L})^2\) population parameter component. The other components and the mean are left unchanged. Under this simplified calculation, the variance becomes \(0.64+0.49+0.25=1.38\;(\mathrm{mg/L})^2\), and the standard deviation is \(\sqrt{1.38}=1.175\) mg/L.

The resulting normal-approximation interval is:

8.20\pm1.96\sqrt{1.38} =8.20\pm2.302478\ldots =[5.898,\;10.502]\ \mathrm{mg/L}.

This interval is narrower than the interval including population parameter uncertainty, but it answers a conditional question: what would the forecast variability be if the estimated population parameters were treated as known? It does not represent the full forecast uncertainty specified in the first example. Fixing population estimates may be a useful sensitivity analysis, but it should not be presented as though parameter-estimation uncertainty had been included.

Implementation, Diagnostics, and Alternatives

In a simulation-based forecast, the decomposition can be implemented with nested draws. Draw \(\theta\) from its uncertainty distribution, then draw \(\eta_i\) from its conditional distribution, then draw \(U\) from its conditional distribution, and finally draw residual error from its conditional model. To estimate the separate terms, retain the nesting: conditional averages and variances at each level are needed, not merely an unlabelled collection of simulated \(Y\) values.

Specify the target and conditioning information.
State what future observation is being forecast, what is already known, and which future inputs remain uncertain.

Define conditional distributions.
Specify uncertainty for \(\theta\), the individual’s \(\eta_i\) conditional on \(\theta\) and their data, future inputs \(U\), and residual error.

Compute nested predictions and variances.
Use conditional means and variances in the chosen order, preserving dependencies rather than drawing each source independently by convenience.

Check the total and the forecast distribution.
Verify that components sum to the total variance within simulation error, and inspect whether a normal approximation is adequate before using it for interval endpoints.

Several checks are important. Repeat the simulation with more draws and compare component estimates and interval endpoints; large changes suggest Monte Carlo error. Verify that conditional draws reflect the data actually available for the individual, and that population parameter uncertainty is not inadvertently represented twice. For example, combining a population-parameter draw with an extra variance term intended to represent the same parameter uncertainty can overstate the forecast spread. Conversely, using only a point estimate of \(\theta\) omits that source.

A sum of variances does not guarantee a normal predictive distribution. Nonlinear transformations, skewed random effects, mixture distributions, or uncertain inputs with a few distinct scenarios can produce asymmetry or multiple modes. In such cases, a normal interval based on the summed variance may have poor coverage or even extend into impossible values. Inspect predictive simulations and use empirical predictive quantiles when appropriate. Decomposition explains variance; it does not by itself determine the shape or calibration of an interval.

A full joint simulation is the direct alternative when the main goal is a predictive distribution or interval: simulate all relevant uncertainty sources together and summarize the resulting \(Y\) draws. It naturally preserves nonlinearities if the joint distribution is specified correctly, but does not automatically explain which source contributes what. A delta-method approximation can be faster for smooth models with modest uncertainty, but may be inaccurate for strong curvature or non-normal inputs. Fixing uncertain quantities at point estimates is simpler still, but produces a conditional forecast and can understate total uncertainty.

The decomposition is limited by the model and uncertainty distributions supplied to it. It does not capture an omitted structural mechanism, incorrect residual model, or unrecognized dependence merely by assigning variance components. Component estimates can also be sensitive to the ordering when sources interact. The decomposition is therefore both a calculation and an accounting convention: state what each component means, which conditional distributions were used, and whether the reported interval includes every intended source.

Common Pitfall

A common error is to call the spread of random effects in the population “individual parameter uncertainty” and then add it to a posterior spread for the same individual. These quantities are not interchangeable. Population heterogeneity describes variation across individuals; the conditional distribution of \(\eta_i\) after observing this individual’s data describes uncertainty about their latent parameters. Adding both without a coherent joint model can count variation twice. The correction is to define the forecast target first and use the appropriate conditional distribution for the individual’s latent parameters.

Key takeaway: A forecast interval should include the uncertainty sources relevant to its target and conditioning information. Decompose with conditional distributions, preserve dependencies, and distinguish uncertainty about an individual from variability across the population.

Check Your Understanding

Answer each question, then compare with the brief solutions.

A future concentration model includes random effects, residual error, and an uncertain dose time. Which source describes uncertainty about an observed individual’s latent clearance?
Why can the components assigned to individual parameters and future inputs depend on decomposition order, while their sum remains the same?
Suppose a forecast has total variance \(2.00\;(\mathrm{mg/L})^2\), including residual variance \(0.36\;(\mathrm{mg/L})^2\). What variance and standard deviation result if residual error is omitted, with all other components unchanged?
Why is a normal prediction interval not automatically justified by having computed a total variance?

Solutions. (1) The conditional uncertainty in the individual’s latent clearance, represented within the distribution for \(\eta_i\) given the population parameters and that individual’s data. (2) Nonlinear interactions can change conditional means and variances depending on which source is averaged over first; the law of total variance still partitions the same overall variance, so the total remains invariant. (3) The variance is \(2.00-0.36=1.64\;(\mathrm{mg/L})^2\), and the standard deviation is \(\sqrt{1.64}=1.281\) mg/L. This narrower forecast omits residual error. (4) Variance specifies spread but not distributional shape. Skewness, mixtures, or strong nonlinearities may make a normal approximation inaccurate, so the predictive distribution or simulated quantiles should be checked.