Constructing Prediction Intervals for Clinical Outcomes

Start With the Predictive Target

A simulated distribution is not yet a prediction interval. First specify what outcome is being predicted, for whom, under which treatment and clinical scenario, and at what time or set of times. Then decide which sources of uncertainty belong in that prediction. The resulting interval may describe plausible outcomes for a future patient, uncertainty about a population mean, or joint uncertainty across several outcomes. Those are different targets, even if they are produced by the same QSP model.

Forward simulation maps parameter ensembles and virtual patients to clinical outcomes. For a future-patient prediction, the simulated values should represent outcomes for individual patients in the intended target population, with parameter uncertainty and any relevant patient-level or residual variability included according to the stated prediction. For a population-mean prediction, the target is instead the mean outcome in that population under the specified model and scenario. Variation among individual patients is not, by itself, uncertainty in that mean.

These are model-based predictive statements: the coverage is with respect to the predictive distribution used to generate the simulations. It is not a guarantee that the model is correctly specified or that observed clinical outcomes will have the stated coverage. The prediction scenario and target population must be defined before calculating quantiles.

Pointwise Intervals From Simulated Outcomes

Suppose \(M\) simulated outcomes \(y^{(1)},\ldots,y^{(M)}\) represent draws from the predictive distribution for one future patient at one prespecified endpoint and time. Sort them as \(y_{(1)}\le \cdots \le y_{(M)}\). An equal-tailed interval with nominal coverage \(1-\alpha\) uses the predictive quantiles at \(\alpha/2\) and \(1-\alpha/2\). One simple empirical quantile convention is the nearest-rank rule, \(q_p=y_{(\lceil Mp\rceil)}\), for \(0<p\le1\). Other software may interpolate between order statistics, so the convention should be reported when it matters.

The equal-tailed construction allocates approximately \(\alpha/2\) predictive probability to each tail. It does not require a normal predictive distribution. Its empirical accuracy depends on the simulated values adequately representing the intended predictive distribution; with finite \(M\), the reported endpoints are estimates of its quantiles.

Population Means Are Different Targets

To predict a population mean, each simulation replicate should produce a value of the mean estimand, not simply one randomly selected patient's outcome. For example, under parameter draw \(r\), calculate or approximate the mean across the specified target population, denoted \(\mu^{(r)}\). Quantiles across \(r=1,\ldots,M\) then summarize uncertainty in that population mean induced by the simulated uncertainty represented in those draws.

If the target is the mean for a finite future sample of patients, that is a different quantity: its distribution also includes the variability from which patients enter that sample. State whether the target is a population mean or a future sample mean. If a virtual population is used to approximate the population mean, its finite simulation size can add numerical error to each replicate; that error should not be mistaken for biological variation or parameter uncertainty.

Simultaneous Intervals Across Outcomes

Separate 95% pointwise intervals do not generally provide 95% coverage that all specified outcomes lie inside their intervals at once. If \(A_k\) is the event that outcome \(k\) falls inside its interval, pointwise coverage concerns \(P(A_k)\); simultaneous coverage concerns \(P(A_1\cap\cdots\cap A_K)\). The latter depends on the joint predictive distribution, including dependence among outcomes. A joint simulation must therefore retain the simulated outcomes together as a vector for each replicate.

One construction for a finite set of \(K\) outcomes uses a maximum standardized deviation. Choose a center \(c_k\) and a positive scale \(s_k\) for each outcome, in its own units. In each joint simulation \(r\), calculate the largest absolute standardized deviation across outcomes. If \(c_{1-\alpha}\) is the corresponding \(1-\alpha\) quantile of these maxima, the simultaneous region is the rectangle with limits \(c_k\pm c_{1-\alpha}s_k\). A common choice is to estimate centers and scales in a separate pilot simulation and calibrate the maximum using joint calibration simulations.

The approximation reflects finite simulation and empirical quantile estimation. The construction is calibrated to the simulated joint distribution; it is not distribution-free protection against model misspecification. Because all outcomes from a replicate contribute to the same maximum, their dependence affects the calibrated multiplier. The scales also matter: standardization allows outcomes with different units to contribute on comparable relative scales, while a zero or inappropriate scale makes this method unsuitable.

Why the Maximum Construction Works

For a proposed multiplier \(c\), all \(K\) outcomes fall inside their limits exactly when every absolute standardized deviation is at most \(c\). That condition is equivalent to the maximum deviation being at most \(c\). Choosing \(c\) as the \(1-\alpha\) quantile of that maximum therefore targets joint coverage of \(1-\alpha\) under the simulated distribution.

A simpler alternative is Bonferroni adjustment. For \(K\) outcomes, construct each pointwise interval with error probability \(\alpha/K\). Since the probability that at least one interval misses is no greater than the sum of the individual miss probabilities, the joint coverage is at least \(1-\alpha\). This needs no dependence estimate, but can be conservative. Maximum-deviation calibration uses the simulated joint distribution and can reflect dependence, but relies more directly on that distribution and the chosen standardization.

For a trajectory, a simultaneous band over a finite set of prespecified visit times can be constructed by treating every outcome-time coordinate as one of the \(K\) coordinates in the maximum. Coverage then applies to that grid. It does not automatically establish coverage at every time between grid points. A denser grid or a separately justified continuous-time method may be needed if the clinical claim concerns the entire trajectory.

Assumptions, Diagnostics, and Limitations

Before interpreting any interval, check that the simulation scheme matches the target. In particular, the simulated patient characteristics should represent the intended population; uncertain parameters and patient-level variation should be included or held fixed according to the stated prediction; and joint outcomes must remain paired by simulation replicate when constructing a simultaneous region. The simulated scenario, endpoint definitions, visit times, units, interval level, and quantile convention should be explicit.

Useful checks include plotting or tabulating the simulated outcome distribution, verifying that the reported order statistics match the stated ranks, and checking that interval limits are on a clinically possible scale. For a simultaneous construction, inspect whether one outcome dominates the maximum statistic because its scale is too small or its units were entered incorrectly. If joint outcomes were simulated separately and then paired arbitrarily, their dependence has been lost, so the resulting simultaneous calibration may not represent the intended joint predictive distribution.

Intervals based on empirical quantiles can be unstable when simulation draws are too few, especially in the tails. A bounded outcome may also produce limits outside its allowable range if a symmetric maximum-deviation construction is used. Such behavior calls for reviewing the target, transformations, and interval construction—not mechanically clipping limits without reassessing coverage. A finite-grid band has a further limitation: its simultaneous statement applies only to the included coordinates. Finally, all these intervals inherit model assumptions; simulation alone cannot establish that the QSP model predicts clinical outcomes accurately.

Common Pitfall

Wrong reasoning: “The simulated individual outcomes are broad, so their 95% prediction interval is also a 95% interval for the population mean.” This confuses patient-to-patient variability with uncertainty in an average. Correction: generate one value of the mean estimand per uncertainty replicate, averaging over the defined target population as appropriate, and take quantiles of those mean values. Use individual-outcome draws when the target is a future patient. Also, do not label separate pointwise intervals “95% simultaneous”: joint coverage is a separate probability that requires joint calibration or a stated adjustment.