Why a Forecast Needs More Than a Parameter Estimate
A concentration forecast depends on more than the estimated values of clearance, volume, or an elimination rate. Even after using the observed data to update an individual’s PK parameters, those parameters remain uncertain. The future concentration may also vary around the model prediction because of residual error. If the dosing record is uncertain, different plausible dose histories can produce different predictions as well.
These sources of uncertainty matter especially when the PK prediction is nonlinear in its parameters. For example, in a one-compartment model with first-order elimination, concentration contains the term \(\exp(-kt)\), where \(k\) is the elimination rate and \(t\) is time since a dose. Evaluating that function at a posterior mean for \(k\) generally does not give the posterior mean concentration. The appropriate forecast averages predictions over the uncertainty in the inputs.
The posterior predictive concentration forecast combines these sources rather than treating the most likely parameter values and dose history as certain. Here the focus is how to decompose and propagate the resulting uncertainty, conditional on the observed data \(D\).
Definition and Variance Decomposition
The first term is residual uncertainty: after fixing the parameters and history, the future observation may still differ from its conditional mean. The second is parameter uncertainty: predictions vary across plausible parameter values, averaged over the plausible histories. The third is dosing-history uncertainty: the parameter-averaged predictions differ across histories.
This is an application of the law of total variance, applied first to \(H\) and then to \(\theta\) within each history. It does not require \(\theta\) and \(H\) to be independent, because the parameter distribution is explicitly conditional on \(H\). The allocation of any parameter-history interaction depends on the order of conditioning; the displayed decomposition assigns between-history variation to the third term. The sum, \(\operatorname{Var}(C\mid D)\), is unchanged.
A useful special case has a known dose history. Then the between-history term disappears. If the model has no residual variation in the predicted quantity—for example, if the target is the conditional mean concentration rather than a future observation—the first term is zero. A forecast of an actual future measured concentration usually includes residual uncertainty as well.
Nonlinear Propagation of Parameter Uncertainty
Consider a hypothetical one-compartment intravenous bolus model with concentration \(C(t)=A\exp(-kt)\), where \(A\) is dose divided by volume and \(k\) is the elimination rate. Suppose the posterior distribution for \(k\), conditional on the data, places probability one-half on each of two values. This discrete posterior makes it possible to integrate the nonlinear prediction exactly.
Worked Example: Averaging Exponential Predictions
Suppose \(A=10\) mg/L, the forecast time is \(t=6\) h, and the posterior values of \(k\) are \(0.0800\) h\(^{-1}\) and \(0.1200\) h\(^{-1}\), each with probability \(0.5\). Assume no residual error for this example and a known dose history. Calculate the posterior mean and variance of the model-predicted concentration.
At \(k=0.0800\) h\(^{-1}\), \(C=10\exp[-(0.0800)(6)]=10\exp(-0.48)=6.1878\) mg/L. At \(k=0.1200\) h\(^{-1}\), \(C=10\exp[-(0.1200)(6)]=10\exp(-0.72)=4.8675\) mg/L.
The posterior mean prediction is \(E(C\mid D)=0.5(6.1878)+0.5(4.8675)=5.5277\) mg/L. The variance of these two equally weighted predictions is $$ \operatorname{Var}(C\mid D) =0.5(6.1878-5.5277)^2+0.5(4.8675-5.5277)^2 \approx 0.4358\;(\mathrm{mg/L})^2. $$
The posterior mean of \(k\) is \(0.1000\) h\(^{-1}\). Substituting only that value gives \(10\exp[-(0.1000)(6)]=5.4881\) mg/L, not the integrated mean of \(5.5277\) mg/L. The difference arises because the exponential prediction is nonlinear in \(k\). In this case \(C(k)\) is convex, since its second derivative with respect to \(k\) is \(At^2\exp(-kt)>0\) for \(A,t>0\). Thus averaging the predictions and predicting at the average parameter need not agree.
The example illustrates a general point: parameter uncertainty should be propagated through the prediction function, rather than appended to a forecast after evaluating the function at a single estimate. The same principle applies when several parameters enter jointly or have correlated posterior distributions. Their joint uncertainty must be retained when generating predictions.
Uncertain Dosing History and the Total Variance
When the dosing history is uncertain, each plausible history implies its own concentration prediction. Replacing an uncertain dose with a probability-weighted “average dose” is generally not equivalent to averaging the predictions under the possible histories. The following example uses a linear superposition model for bolus doses and first-order elimination.
Worked Example: A Dose May or May Not Have Been Taken
Consider a hypothetical patient with \(V=30\) L and \(k=0.100\) h\(^{-1}\). A 300 mg bolus dose was taken at time 0, and concentration is forecast at time 6 h. An additional 300 mg dose may have been taken at time 4 h: the posterior probability it was taken is \(0.25\), and the probability it was not taken is \(0.75\). Treat parameters as known and residual error as zero.
Each dose contributes \((300/30)\exp(-k\times\text{time since dose})\) mg/L. If the additional dose was not taken, the prediction is \(C_0=10\exp[-(0.100)(6)]=10\exp(-0.6)=5.4881\) mg/L.
If the additional dose was taken, its contribution at 6 h is \(10\exp[-(0.100)(2)]=10\exp(-0.2)=8.1873\) mg/L. The total is \(C_1=10\exp(-0.6)+10\exp(-0.2)=5.4881+8.1873=13.6754\) mg/L.
The history-averaged forecast mean is \(E(C\mid D)=0.75(5.4881)+0.25(13.6754)=7.5349\) mg/L. For a two-history distribution, the variance is \(q(1-q)(C_1-C_0)^2\), where \(q=0.25\). Therefore, $$ \operatorname{Var}_{H\mid D}\{E(C\mid H,D)\} =0.25(0.75)(13.6754-5.4881)^2 \approx 12.5685\;(\mathrm{mg/L})^2. $$
That variance is due solely to uncertainty about whether the dose was taken. It is not residual variability and does not arise from uncertain \(k\) or \(V\), which were fixed here. The large difference between the two history-specific predictions shows why the history should be represented as uncertain when the evidence does not resolve it.
Combining the Three Contributions
The decomposition is most useful when all relevant uncertainties are represented in one calculation. In practice, posterior simulation often does this directly: draw a history, draw parameters from the posterior conditional on that history, calculate the model prediction, and then draw residual variation if the target is an observed future concentration. The next example works through the variance components rather than relying on simulation software.
Worked Example: Parameter, Residual, and History Uncertainty Together
Use the dose histories above, with probabilities \(0.75\) for no additional dose and \(0.25\) for a dose at time 4 h. Now suppose \(k\) is uncertain, with posterior values \(0.0800\) and \(0.1200\) h\(^{-1}\), each equally likely conditional on either history. Let \(V=30\) L, each dose be 300 mg, and let the future observation have additive residual error with conditional mean zero and standard deviation \(0.250\) mg/L. Assume that residual variance is constant and that the history and parameter distributions factor as specified.
For the history with no additional dose, the two parameter-specific predictions are \(6.1878\) and \(4.8675\) mg/L. Their mean is \(5.5277\) mg/L and their variance is \(0.4358\;(\mathrm{mg/L})^2\).
For the history with the additional dose, the predictions are \(10\exp[-(0.0800)(6)]+10\exp[-(0.0800)(2)]=6.1878+8.5214=14.7093\) mg/L, and \(10\exp[-(0.1200)(6)]+10\exp[-(0.1200)(2)]=4.8675+7.8663=12.7338\) mg/L. Their mean is \(13.7215\) mg/L and their variance is approximately \(0.9756\;(\mathrm{mg/L})^2\).
The average parameter contribution is \(0.75(0.4358)+0.25(0.9756)=0.5708\;(\mathrm{mg/L})^2\). The two history-specific parameter-averaged means differ by \(13.7215-5.5277=8.1938\) mg/L, so the history contribution is \(0.75(0.25)(8.1938)^2\approx12.5886\;(\mathrm{mg/L})^2\). The residual contribution is \(0.250^2=0.06250\;(\mathrm{mg/L})^2\).
Consequently, the total forecast mean is \(0.75(5.5277)+0.25(13.7215)=7.5761\) mg/L, and the variance is \(0.5708+12.5886+0.06250=13.2219\;(\mathrm{mg/L})^2\). The corresponding standard deviation is \(\sqrt{13.2219}=3.6362\) mg/L. In this hypothetical calculation, dosing-history uncertainty is the largest variance contribution. That conclusion is specific to these assumed probabilities and model values; it is not a general ranking of uncertainty sources.
Derivation and Practical Calculation
The variance identity follows from conditioning. For any random variable \(C\) with finite variance and any conditioning variable \(Z\), \(\operatorname{Var}(C)=E\{\operatorname{Var}(C\mid Z)\}+\operatorname{Var}\{E(C\mid Z)\}\). Take \(Z=(\theta,H)\). The first term is the average conditional residual variance. For the second term, set \(M=E(C\mid\theta,H,D)\), then condition on \(H\): $$ \operatorname{Var}(M\mid D) =E_{H\mid D}\{\operatorname{Var}_{\theta\mid H,D}(M)\} +\operatorname{Var}_{H\mid D}\{E_{\theta\mid H,D}(M)\}. $$ Combining these two identities gives the stated three-part decomposition. This derivation also shows why the components add exactly when defined with the same conditioning structure.
A practical calculation can be organized as follows.
Distinguish a conditional mean concentration from a future measured concentration. The latter generally includes residual variation.
Use the joint posterior for PK parameters and, when needed, a probability distribution over plausible dosing histories. Preserve dependence represented by that posterior.
For each parameter-history draw, evaluate the full nonlinear PK model. Add a residual draw only when predicting an observation rather than its conditional mean.
Estimate the predictive mean and variance, then calculate the conditional variance components using the same draws and conditioning order. Confirm that their sum agrees with the total predictive variance within simulation error.
Diagnostics, Alternatives, and Limitations
A basic diagnostic is to compare the integrated forecast with a plug-in forecast at posterior mean parameters and the most probable history. A difference is not automatically evidence of a problem; it can be expected under nonlinear predictions or multimodal history uncertainty. It does indicate that the plug-in value is not a substitute for the integrated forecast.
With Monte Carlo calculations, inspect the distribution of predicted concentrations, not only its mean and variance. A few extreme draws can dominate the variance, and a multimodal distribution can be poorly summarized by a mean and standard deviation. Check that simulated parameter-history pairs follow the intended posterior, that predictions use consistent units and times, and that component estimates stabilize as the number of draws increases. If the component sum and the directly calculated predictive variance disagree beyond Monte Carlo error, check the conditioning, weighting, and whether residual variation was counted once or more than once.
The main alternatives have different roles:
| Approach | What it does | Main limitation |
|---|---|---|
| Plug-in prediction | Evaluates the model at one parameter estimate and one history. | Does not represent parameter or history uncertainty and can misstate the mean under nonlinearity. |
| Local approximation | Uses a local expansion around a parameter estimate to approximate uncertainty propagation. | Can be inaccurate for strongly nonlinear predictions, broad posteriors, or multiple plausible histories. |
| Posterior integration or simulation | Averages model predictions over joint parameter and history uncertainty; residual draws can be included for observation forecasts. | Requires adequate posterior representation and enough computation or simulation draws. |
The decomposition quantifies uncertainty under the assumed model; it does not automatically account for model misspecification, unrepresented dose-record errors, or a residual model that fails to describe future measurements. An additive residual model may also assign probability to negative concentrations, so it may be unsuitable when that probability is material. A residual model should be chosen and checked for the concentration scale and forecast purpose. Finally, a variance decomposition depends on the chosen conditioning order when parameter and history uncertainties interact. Report the order and the component definitions, rather than treating the allocation as unique.
Common Pitfall
Correction: For a nonlinear model, calculate predictions over the posterior distribution and average those predictions. If history is uncertain, average over its plausible values as well. Add residual variability when forecasting an observation, and use the variance decomposition to distinguish the contributions.
Another frequent error is to include residual error both in the conditional variance term and as an extra variance added after the full predictive variance has already been calculated. Define the forecast target first and include each source exactly once.
Check Your Understanding
Answer each question, then compare with the concise solutions below.
- If the dosing history is known and the forecast is the model-based conditional mean, which terms in the three-part decomposition remain?
- For \(C=A\exp(-kt)\), why can \(C\) evaluated at the posterior mean of \(k\) differ from the posterior mean of \(C\)?
- In the uncertain-dose example, what is the between-history variance contribution, and what does it represent?
- Suppose two histories have the same parameter-averaged predicted concentration. What is their contribution to the between-history variance?
- Why should residual error not be added again if the directly simulated predictions already include residual draws?
Solutions. (1) The residual term is zero for the model-based conditional mean, and the history term is zero when the history is known; the parameter term remains unless parameters are also known. (2) The exponential function is nonlinear, so averaging its values over posterior \(k\) need not equal evaluating it at the average \(k\). (3) It is \(0.25(0.75)(13.6754-5.4881)^2\approx12.5685\;(\mathrm{mg/L})^2\), representing variation between the predictions under the two possible dose histories. (4) Zero, because the variance across identical history-specific means is zero. (5) Those draws already include residual variation, so adding the residual variance again would count that source twice.