Simulation-Based Calibration for Bayesian Forecasting

Why Calibrate a Bayesian Forecasting Algorithm?

An individual forecasting algorithm may combine a population model, an individual’s concentration observations, and a posterior distribution for that individual’s parameters. It can then forecast a future concentration, AUC, Cmax, or the probability of meeting an exposure target. Even when the model is specified correctly, the computation that turns observations into a posterior can be wrong: a likelihood may be coded incorrectly, a parameter transformation may be mishandled, or posterior uncertainty may be understated.

Simulation-based calibration (SBC) is a way to test that computation under a known generative model. Repeatedly simulate an individual’s parameter values and observations, run the forecasting algorithm, and compare the simulated truth with its posterior distribution. If the model and computation meet the required conditions, the true parameter should occupy a uniform rank among posterior draws. Simulated interval coverage provides a complementary check of uncertainty summaries.

SBC is a test of the algorithm under the assumptions used to simulate. It is not, by itself, evidence that those assumptions describe patients in practice. That distinction matters in pharmacometrics: a forecasting method can pass SBC for its assumed individual PK model and still forecast poorly if real patients follow a different model.

What Must Be Simulated—and What Must Match?

For an individual forecasting algorithm, write \(\theta\) for the parameter or parameter vector being inferred for one individual, \(y\) for that individual’s observed data, and \(y^\ast\) for a future quantity to be forecast. A typical simulation replicate proceeds as follows:

1. Draw \(\theta_0\) from the prior or population distribution used by the inference model.
2. Generate the individual’s observed data \(y\) from the model conditional on \(\theta_0\).
3. Run the forecasting algorithm using \(y\), obtaining posterior draws \(\theta^{(1)},\ldots,\theta^{(L)}\).
4. Compare \(\theta_0\) with those draws, or evaluate posterior intervals and forecasts against their simulated truths.

The simulation and inference steps must correspond. If the inference model conditions on known population parameters, the simulation should use those same fixed values. If population parameters are estimated, simulate them according to the specified hierarchy and include them in the inference target as appropriate. Dose history, sampling times, residual error, and any assumed covariates also need to be generated consistently with the model being tested. This is a test of the implemented procedure under that model, not a test of every possible clinical data-generating process.

For the basic rank result, assume the simulated truth is drawn from the same prior used by the algorithm, the likelihood used to generate \(y\) is the likelihood used for inference, and the \(L\) posterior draws are independent draws from the exact posterior given \(y\). Assume also that the scalar quantity being ranked has a continuous posterior distribution, so ties occur with probability zero. These conditions make the truth and the posterior draws exchangeable conditional on \(y\): each has the same posterior distribution, and their order is equally likely.

Why the Rank Is Uniform

For a scalar parameter with continuous posterior distribution, define the rank of the simulated truth among \(L\) posterior draws by

Here \(I(\cdot)\) equals 1 when its condition is true and 0 otherwise. The rank is the number of posterior draws below the truth. It is not a p-value, nor does a particular rank establish that the algorithm is calibrated. The calibration statement concerns the distribution of ranks across many independent simulation replicates.

Conditional on the simulated observations \(y\), the truth \(\theta_0\) has the posterior distribution \(p(\theta\mid y)\). Each exact posterior draw also has distribution \(p(\theta\mid y)\). Under the independence assumption, the \(L+1\) values are exchangeable, so each of the \(L+1\) possible order positions for the truth has probability \(1/(L+1)\). Therefore

This conditional result also holds after averaging over simulated datasets. Across \(S\) independent replicates, each rank has expected count \(S/(L+1)\). A histogram that is roughly flat is consistent with calibration, while systematic shapes can indicate bias or other computational problems. Random fluctuations are expected, especially with few replicates.

For vector parameters, ranks can be computed separately for each component, and they can also be computed for a scalar derived quantity such as clearance or a forecasted AUC. A component-wise rank check does not fully assess joint dependence or every derived forecast. Check quantities that matter to the intended forecast, not only parameters that are convenient to extract.

Checking Rank Histograms in Practice

The rank calculation is simple, but a useful diagnostic requires enough replicates and a clear account of posterior sampling. As a hypothetical illustration, take \(L=99\) posterior draws per replicate and \(S=1{,}000\) independent replicates. Group the 100 possible ranks into ten consecutive bins of ten ranks each. If ranks are uniform, the expected count in each bin is \(1{,}000/10=100\).

A practical SBC run should track whether each replicate completed successfully, the parameterization in which ranks were calculated, and any numerical warnings. Examine ranks for important individual parameters and derived forecast targets separately. Repeat the exercise with enough replicates to distinguish a persistent pattern from Monte Carlo noise, and compare results across relevant designs, such as sparse and denser sampling schedules. A single pooled histogram can hide a failure confined to one design or subgroup.

The basic uniform-rank derivation assumes independent posterior draws. Markov chain Monte Carlo draws are usually autocorrelated, so treating a long correlated sequence as \(L\) independent draws can distort rank diagnostics even when the chain has the correct marginal distribution. Use an approach whose rank reference accounts for the sampling dependence, or use sufficiently independent posterior draws; thinning alone is not a guarantee of independence. For discrete parameters or quantities with posterior ties, use a prespecified randomized tie-breaking rule or a rank method designed for ties. Otherwise, apparent nonuniformity may arise from the rank construction rather than incorrect inference.

Checking Uncertainty Summaries

Ranks evaluate where the truth falls within posterior draws. Interval coverage asks a directly interpretable question: across repeated simulations, how often does an interval constructed by the algorithm contain the relevant simulated truth? If a procedure reports a nominal 95% posterior interval for a parameter, exact Bayesian computation under the prior-generative setup gives 95% prior-averaged coverage. This is a repeated-simulation property; it does not mean that a particular reported interval has a 95% probability of containing a fixed parameter.

For a forecast interval, the target must be clear. To assess a future concentration, simulate that future concentration from the model conditional on the simulated individual parameters and the planned dose and sampling time. Then check whether the algorithm’s predictive interval contains that simulated outcome. This tests the whole predictive computation, including the future observation model. An interval for a latent expected concentration and an interval for a future noisy measurement are different objects and need not have the same width or coverage.

What SBC Does Not Establish

SBC complements rather than replaces other forms of evaluation. The earlier material on calibration of individual predictive intervals addresses interval performance; SBC adds a controlled way to check whether the implemented Bayesian computation behaves as expected under its own generative assumptions. Posterior predictive checks instead compare features of observed data with data replicated from a fitted model. External validation evaluates forecasts against data not used to construct the model. These approaches answer different questions.

A model can pass SBC and still be misspecified for real patients, because the simulation uses the assumed model by construction. SBC can also miss a failure if the tested parameter or derived quantity is not relevant to it, if too few replicates are run, or if posterior sampling error obscures the pattern. Treat results as evidence about a specified computational pipeline under specified simulated conditions, not as a universal certification.

Common Mistakes and Failure Modes

• Interpreting one rank as a verdict. A single rank of 0 or \(L\) is allowed under a uniform rank distribution. Look for systematic patterns across replicates.
• Simulating from a different model than the one used for inference. That may be useful for robustness studies, but it is not the basic SBC uniform-rank test. For that test, simulation and inference must use matching assumptions.
• Calling prior-averaged coverage conditional coverage. Nominal interval coverage under prior-generated parameters averages over that prior. It does not guarantee nominal coverage for every fixed true parameter value or every clinical subgroup.
• Ignoring Monte Carlo error. A coverage estimate a few percentage points below nominal can occur by chance when the number of replicates is modest. Quantify its simulation uncertainty.
• Using correlated draws as if they were independent. The uniform-rank argument relies on exchangeability. Account for posterior-draw dependence when selecting or interpreting a rank procedure.