Why predict at more than one future time?
A patient-specific forecast is often needed at several points in a dosing schedule: for example, near a future peak, at a later concentration measurement, and just before a subsequent dose. These forecasts are connected, but their uncertainty need not be the same. A patient’s parameter uncertainty can affect concentrations differently as time passes or doses accumulate, while residual variability contributes additional uncertainty to each future observation.
A prediction interval concerns a future concentration or measurement, not just the model’s expected concentration. It must therefore account for uncertainty in the patient-specific parameters and, when predicting a measured concentration, the observation-level error represented by the model. A band of intervals at several times is usually a collection of pointwise intervals: each interval refers to its own time, rather than guaranteeing simultaneous coverage of the whole future concentration profile.
The conditional distribution \(p(C(t)\mid\theta,D)\) depends on what is being predicted. For a latent model-predicted concentration, it may be a point at the value \(f(t;\theta)\), where \(f\) is the structural PK model. For a future measured concentration, it also includes the observation or residual-error model. If the residual error depends on concentration, its conditional variance may differ across parameter draws and times.
Combining parameter and residual uncertainty
Suppose a future observation is represented as \(Y(t)=f(t;\theta)+\varepsilon(t)\), where \(f(t;\theta)\) is the model prediction and \(\varepsilon(t)\) is residual error. The posterior distribution \(p(\theta\mid D)\) describes remaining uncertainty about the patient’s parameters. Even if the residual error has mean zero conditional on \(\theta\), a future observation varies both because \(\theta\) is uncertain and because a new residual error occurs.
The variance can be separated using the law of total variance. First consider the variation of the future observation conditional on \(\theta\); then account for how its conditional mean varies over the parameter distribution.
Under an additive, mean-zero residual model with variance \(\sigma^2(t,\theta)\), the first term is the posterior average of residual variance and the second is the posterior variance of the structural prediction. The identity follows by writing \(Y-\operatorname{E}(Y\mid D)\) as the sum of \(Y-\operatorname{E}(Y\mid\theta,D)\) and \(\operatorname{E}(Y\mid\theta,D)-\operatorname{E}(Y\mid D)\). The cross-product has expectation zero because the first term has conditional mean zero given \(\theta,D\). Thus neither source of uncertainty can generally be omitted when predicting a future measurement.
Worked Example: Decomposing Uncertainty at One Future Time
Suppose a hypothetical patient’s posterior distribution is represented by three equally weighted parameter draws. At a specified future time, the structural model predicts 8, 10, and 12 mg/L from these draws. Assume independent additive residual error with mean zero and standard deviation 1.5 mg/L.
The predictive mean of the structural prediction is $$ \bar{C}_{\mathrm{model}}=\frac{8+10+12}{3}=10\ \mathrm{mg/L}. $$ The posterior variance of the model prediction, using the equally weighted draws as the posterior distribution, is $$ \operatorname{Var}(f\mid D) =\frac{(8-10)^2+(10-10)^2+(12-10)^2}{3} =\frac{8}{3}=2.667\ (\mathrm{mg/L})^2. $$
The residual variance is \(1.5^2=2.25\ (\mathrm{mg/L})^2\). The total predictive variance is therefore \(2.667+2.25=4.917\ (\mathrm{mg/L})^2\), and the predictive standard deviation is \(\sqrt{4.917}=2.217\) mg/L. If a normal approximation is appropriate, an approximate 95% prediction interval is $$ 10\pm1.96(2.217)=(5.655,\ 14.345)\ \mathrm{mg/L}. $$ This interval concerns a future observation. An interval for the structural prediction alone would omit residual variance and would be narrower.
Building intervals across future doses
A practical method is to draw parameter values from the patient’s posterior distribution, evaluate the structural PK model at every time of interest for each draw, then draw residual errors and add them when predicting future observations. The resulting simulated values approximate the posterior predictive distribution. Quantiles at each time give pointwise prediction intervals. This procedure retains nonlinearity and skewness in the model predictions, unlike a normal approximation based only on a mean and variance.
Specify future dose amounts and times, and the concentration times to predict. Use the same dosing history assumptions throughout the calculation.
Sample \(\theta^{(b)}\) from \(p(\theta\mid D)\) for draw \(b\). Preserve posterior dependence among parameters rather than drawing each parameter independently.
Calculate \(f(t;\theta^{(b)})\) at each future time. For an observation prediction, draw an error from the specified conditional residual model and add it on the appropriate scale.
At each time, use the simulated predictive values to estimate the desired quantiles. Report that the intervals are pointwise and state the prediction scale and units.
Worked Example: Interval Width as Doses Accumulate
Consider a hypothetical one-compartment bolus model with elimination fraction \(q=0.5\) over each dosing interval. Let \(A\) be the concentration increment from one dose, and suppose the patient’s three equally weighted posterior draws for \(A\) are 4, 6, and 8 mg/L. At the instant after dose \(m+1\), the concentration predicted from the accumulated doses is $$ f_m(A)=A(1+q+\cdots+q^m) =A\frac{1-q^{m+1}}{1-q}. $$ Assume additive residual error with standard deviation 0.5 mg/L at each prediction time.
After the first dose (\(m=0\)), the three structural predictions are 4, 6, and 8 mg/L. Their mean is 6 mg/L and their posterior variance is \(8/3=2.667\ (\mathrm{mg/L})^2\). Adding residual variance \(0.5^2=0.25\) gives total predictive variance 2.917 \((\mathrm{mg/L})^2\) and standard deviation 1.708 mg/L. The moment-matched normal interval is \(6\pm1.96(1.708)\), or approximately 2.654 to 9.346 mg/L.
After the second dose, the accumulation factor is \(1+0.5=1.5\), so the predictions are 6, 9, and 12 mg/L. The mean is 9 mg/L and the posterior variance is 6.000 \((\mathrm{mg/L})^2\). Including residual variance gives 6.250 \((\mathrm{mg/L})^2\), a standard deviation of 2.500 mg/L, and an approximate interval of 4.100 to 13.900 mg/L.
After the third dose, the factor is \(1+0.5+0.25=1.75\), giving predictions of 7, 10.5, and 14 mg/L. Their mean is 10.5 mg/L and their posterior variance is $$ \frac{(7-10.5)^2+(10.5-10.5)^2+(14-10.5)^2}{3} =\frac{24.5}{3}=8.167\ (\mathrm{mg/L})^2. $$ The total variance is \(8.167+0.25=8.417\ (\mathrm{mg/L})^2\), with standard deviation 2.901 mg/L. The corresponding approximate interval is 4.814 to 16.186 mg/L. In this example, uncertainty in \(A\) is amplified by accumulation, so the interval widens. This pattern is not universal: changes in elimination, sampling phase, or residual-error variance can alter how width changes over time.
Intervals on a transformed scale
Concentration predictions are positive, and their predictive distributions can be skewed. An additive normal approximation on the concentration scale may then be unsuitable, particularly if it gives appreciable probability to negative concentrations. If the model uses a lognormal observation error, intervals can instead be constructed on the log scale and transformed back, provided the predictive distribution on that scale is adequately represented.
Worked Example: Transforming a Log-Scale Predictive Interval
At one future time, suppose the hypothetical predictive distribution satisfies \(\log C\mid D\sim N(\log 8,0.3^2)\), where concentration is in mg/L and the logarithm is natural. Since a normal distribution’s central 95% interval is its mean plus or minus 1.96 standard deviations, the interval on the log scale is \(\log 8\pm1.96(0.3)\). Exponentiating both limits gives $$ \left(8e^{-0.588},\ 8e^{0.588}\right) \approx (4.443,\ 14.405)\ \mathrm{mg/L}. $$
At another future time, suppose \(\log C\mid D\sim N(\log 5,0.5^2)\). Its interval is \(5e^{-0.98}\) to \(5e^{0.98}\), or approximately 1.875 to 13.318 mg/L. The second interval is wider because predictive uncertainty on the log scale is greater. These intervals are asymmetric around their medians (8 and 5 mg/L). The assumed normal distributions on the log scale are essential to this calculation; if posterior predictions form a more complex mixture, simulate the predictive distribution and take its quantiles instead.
Diagnostics, alternatives, and limitations
Check that predictions at every time use the intended dose history, dose times, and units. Inspect the posterior draws of predicted concentrations across time: implausible jumps, unexpected negative concentrations, or rapidly increasing dispersion may indicate a coding error, an unsuitable residual model, or sensitivity to poorly informed parameters. Verify that residual variability is added on the scale specified by the model; adding a log-scale error directly to concentrations is not equivalent to adding it to log concentrations.
The principal alternative to draw-based prediction is a moment-based approximation. It can be efficient when the predictive distribution is close to normal on the chosen scale and its mean and variance are calculated appropriately. However, it may conceal skewness, multiple modes, or parameter-dependent residual variance. A plug-in forecast that uses only a single point estimate of \(\theta\) is simpler still, but it omits parameter uncertainty and is not a full patient-specific prediction interval. It should not be presented as though it represents all predictive uncertainty.
A pointwise 95% interval at each of several times does not imply a 95% probability that the entire future profile lies inside all intervals simultaneously. If the clinical question concerns a whole trajectory or whether a concentration ever crosses a threshold, the joint predictive distribution across times is needed. Prediction quality also depends on the structural model, the posterior distribution, and assumptions about future dosing and adherence. An interval can be computed precisely under a model that is nevertheless an inadequate description of the patient.
A common mistake is to calculate uncertainty in the model prediction, add a residual standard deviation rather than a residual variance, or ignore residual error altogether. The correct variance decomposition adds variances, not standard deviations, when the assumptions for the decomposition hold. Another mistake is to call a 95% prediction interval a 95% interval for the unknown parameter: it describes future outcomes under the predictive model, not parameter uncertainty alone.
Check Your Understanding
Answer the questions, then compare with the solutions below.
- A future model prediction has posterior variance 4 \((\mathrm{mg/L})^2\), and residual variance is 1 \((\mathrm{mg/L})^2\). What is the predictive variance under the additive, mean-zero error assumptions?
- Why is a pointwise interval at each of six future times not automatically a 95% interval for the entire six-time profile?
- If an additive normal approximation gives a lower prediction limit below zero, what should be examined before reporting that interval?
- In the bolus example, what is the structural prediction after the fourth dose for a draw with \(A=6\) mg/L and \(q=0.5\)?
Solutions. (1) The predictive variance is \(4+1=5\ (\mathrm{mg/L})^2\), because the conditional residual variance and the variance of the conditional mean are added. (2) Each interval concerns one time separately; simultaneous coverage depends on the joint distribution across times and is not guaranteed by the individual coverages. (3) Examine whether the normal approximation and additive error model are appropriate on the concentration scale; a log-scale model or direct predictive simulation may better represent positive, skewed concentrations. (4) The accumulation factor is \(1+0.5+0.25+0.125=1.875\), so the prediction is \(6(1.875)=11.25\) mg/L.