Posterior Predictive Checks for QSP Calibration

Why Check Predictions After Calibration?

Calibration produces a posterior distribution for model parameters given the measurements. But parameter summaries alone do not show whether the calibrated model can generate data resembling the observed time course. A posterior predictive check addresses that question by simulating replicated observations under the fitted model and comparing them with the measurements.

For a QSP model, it is essential to distinguish a predicted biological trajectory from a replicated measurement. An ODE solution may describe the model’s latent biomarker concentration, while the laboratory observations also reflect measurement error, transformations, correlations, or censoring specified by the observation model. A check against measured data should therefore reproduce the observation process, not merely plot the latent ODE solution.

Definition: Let $y$ denote the observed data, $\theta$ the mechanistic parameters, and $\psi$ any parameters in the observation model. A posterior predictive replicate $y^{\mathrm{rep}}$ is generated by first drawing $(\theta,\psi)$ from the posterior $p(\theta,\psi\mid y)$, then drawing new data from the observation model $p(y^{\mathrm{rep}}\mid\theta,\psi)$ at the same design points as the observations. The posterior predictive distribution is $$ p(y^{\mathrm{rep}}\mid y)=\int p(y^{\mathrm{rep}}\mid\theta,\psi)\,p(\theta,\psi\mid y)\,d\theta\,d\psi. $$

In practice, posterior draws may come from MCMC. For each retained draw, evaluate the mechanistic model at the measured times, then simulate observations according to the fitted observation model. If the model includes process noise, that noise must also be simulated. For hierarchical models, decide whether the check is conditional on the individuals already observed or concerns new individuals: the latter requires drawing new subject-level effects from their population distribution. These are different predictive questions.

The basic assumptions are that the posterior draws adequately represent the calibrated posterior, the simulator and observation model are implemented as intended, and the simulated data correspond to the same measurement schedule and data-processing rules as the observations. The check assesses implications of this specified model; it does not establish that the model is biologically correct.

Generating and Reading Predictive Time-Course Plots

A common first check overlays the observed trajectory with trajectories generated from posterior draws. Alternatively, show pointwise predictive intervals at each measurement time. These intervals include uncertainty in parameters and, when simulated, observation noise. They are not intervals for the latent mean trajectory alone.

Pointwise intervals are calculated separately at each time. For instance, a 90% pointwise interval contains the central 90% of replicated observations at a particular time. It does not promise 90% coverage of an entire trajectory. With many time points, occasional observations outside their pointwise intervals can occur even when the model is adequate. Conversely, lying inside every interval does not rule out a systematic pattern, such as consistently underpredicting late measurements.

Worked Example: Plotting Replicated Biomarker Measurements

Suppose a hypothetical calibration has four measurement times. The observed biomarker concentrations are 10.2, 8.1, 6.8, and 5.2 mg/L. Across 1,000 posterior predictive replicates, suppose the predictive medians are 10.0, 8.4, 6.8, and 5.3 mg/L, with 90% pointwise intervals of [9.0, 11.1], [7.1, 9.7], [5.4, 8.2], and [3.9, 6.6] mg/L.

Each observation falls within its interval: for example, $8.1$ mg/L lies between $7.1$ and $9.7$ mg/L at the second time. The measured course is therefore compatible with the pointwise predictive ranges in this particular check. This does not show that the joint trajectory is well reproduced; the pattern across times still needs inspection.

To see how a replicate is generated, consider one hypothetical posterior draw with predicted means $[9.8, 8.2, 6.5, 5.0]$ mg/L and an additive Gaussian observation error with standard deviation $0.5$ mg/L. If the simulated standard-normal errors are $[0.4,-1.2,0.6,0.8]$, then the replicate is

y^{\mathrm{rep}}_j=f(t_j;\theta)+\epsilon_j,\qquad y^{\mathrm{rep}}=[9.8+0.5(0.4),\;8.2+0.5(-1.2),\;6.5+0.5(0.6),\;5.0+0.5(0.8)] =[10.0,\;7.6,\;6.8,\;5.4]\ \mathrm{mg/L}.

This is one possible replicate, not a new estimate of the fitted curve. Repeating the process over posterior draws produces the predictive distribution used for the plot.

Prespecified Discrepancy Statistics

A plot is useful, but some forms of mismatch are easier to detect with a scalar discrepancy statistic. Choose the statistic before inspecting the check’s results, and choose it to express a scientifically relevant feature: for example, late-time bias, peak height, time to peak, or the number of threshold crossings. For each posterior draw, calculate the statistic both for the observed data and for a replicate generated using that same draw. Pairing the calculations this way accounts for parameter uncertainty in the comparison.

Definition: For a discrepancy statistic $T(y,\theta,\psi)$, the posterior predictive tail area for an upper-tail check is estimated from $M$ draws as $$ \hat{p}_{\mathrm{ppc}}=\frac{1}{M}\sum_{m=1}^{M} I\!\left\{T(y^{\mathrm{rep},m},\theta^{(m)},\psi^{(m)}) \geq T(y,\theta^{(m)},\psi^{(m)})\right\}, $$ where $I\{\cdot\}$ equals 1 when its condition is true and 0 otherwise. For a lower-tail check, reverse the comparison; for a two-sided check, define the extremeness rule explicitly.

This tail area describes how often the replicated statistic is at least as extreme as the observed statistic under the posterior predictive procedure. It is not the probability that the model is true or that a null hypothesis is true. Nor is it generally a classical p-value with a uniform reference distribution under a null model: the observed data inform the posterior used to generate the replicates, so the check reuses the data. Treat the result as a diagnostic signal, interpreted alongside plots and scientific context.

Worked Example: Checking for Late-Time Underprediction

Suppose the prespecified statistic is the mean residual at the final two measurement times, with residual defined as observed concentration minus predicted mean concentration. For one illustrative posterior draw, the observed late measurements are 6.4 and 5.8 mg/L, while the predicted means are 5.4 and 4.6 mg/L.

The residuals are $6.4-5.4=1.0$ mg/L and $5.8-4.6=1.2$ mg/L. Thus the observed statistic is

T_{\mathrm{obs}}=\frac{1.0+1.2}{2}=1.10\ \mathrm{mg/L}.

Suppose 7 of 100 paired posterior predictive replicates have a late-time mean residual at least 1.10 mg/L. Then $\hat{p}_{\mathrm{ppc}}=7/100=0.070$. The observed late measurements are higher than the fitted means, and relatively few replicates show as much late-time excess. This is a reason to investigate possible late-time underprediction, not a standalone verdict. The estimate’s Monte Carlo standard error would be about $\sqrt{0.070(1-0.070)/100}=0.0255$ if the 100 comparisons were independent; correlated MCMC draws can make the effective simulation size smaller.

Check Patterns, Dependence, and Failure Modes

The discrepancy must match the feature being investigated. Pointwise intervals can identify times with unusual measurements but often conceal trajectory-level structure. A statistic based on a prespecified feature can detect structure that pointwise coverage misses. If measurement errors are temporally correlated, replicates must preserve that correlation; simulating independent errors when the fitted model specifies correlated errors would produce the wrong reference distribution.

Worked Example: Checking Residual Alternation Over Time

At four equally spaced times, suppose a hypothetical fitted mean trajectory is $[10,8,6,4]$ mg/L and the observations are $[11,7,7,3]$ mg/L. The residuals are $[1,-1,1,-1]$ mg/L. Define the prespecified statistic as the number of sign changes between adjacent residuals. The observed statistic is 3, since each of the three adjacent pairs changes sign.

Suppose the observation model includes temporal error correlation and 4 of 200 paired predictive replicates have at least 3 sign changes. The estimated upper-tail area is $4/200=0.020$. This suggests the observed alternation is unusual under the fitted error model and warrants checking the residual structure, measurement process, and model specification. It does not by itself identify the cause. In particular, changing the model to independent errors would not be justified merely because the observed residuals alternate; the error model and scientific rationale need examination.

Several practical diagnostics help prevent misleading conclusions:

Inspect complete trajectories. Look for systematic deviations across time, not just isolated interval crossings. Include units and show the observation schedule.
Use the correct observation scale. If the likelihood is defined on log concentrations, simulate and assess on that scale where appropriate; a log-scale error model does not imply additive Gaussian errors on the original scale.
Check model and computation separately. Poorly mixing chains, insufficient effective sample size, or ODE solver failures can make predictive summaries unreliable. A predictive plot is not a substitute for convergence and numerical diagnostics.
Account for the sampling design. Replicate at the actual observed times and apply the same censoring or reporting process. For censored measurements, compare replicated censored data with the recorded-data rule rather than treating censored values as exact measurements.
Do not select a statistic after seeing its result and present it as prespecified. Exploratory checks are useful, but label them exploratory and seek confirmation or additional data when the finding will guide a model change.

What Posterior Predictive Checks Can and Cannot Do

Posterior predictive plots and statistics complement, rather than replace, posterior summaries for mechanistic quantities. Posterior summaries describe uncertainty in parameters or model-derived quantities; predictive checks ask whether simulated measured data reproduce selected features of the observations. Neither establishes that a mechanism is uniquely identified or that the calibrated model will predict well in a new setting.

Approach	Main question	Important limitation
Pointwise predictive intervals	Are observations unusual at individual times?	Can miss a coherent pattern across times; intervals are not simultaneous.
Prespecified discrepancy statistic	Is a targeted feature of the observed data unusual in replicated data?	Findings depend on the statistic and are not a general model-adequacy proof.
Residual plots	Do residuals show visible structure, bias, or changing spread?	Visual assessment can be subjective and should reflect the observation model.
Predictive assessment on held-out data	How well does the model predict data not used in that fit?	Requires a suitable validation design and answers a different question.

A posterior predictive check is especially useful for exposing a mismatch between the chosen likelihood and the observed trajectories, such as implausible variability, systematic timing errors, or residual dependence. Its limits are equally important: many different models can pass a small set of checks, and a check has little power to reveal a feature it was not designed to examine. Use a small set of scientifically motivated checks, state what each targets, and report the observed pattern as well as the predictive comparison.

Key takeaway: Generate replicated measurements by propagating posterior uncertainty through both the mechanistic model and the observation model. Compare those replicates with the measured trajectories using plots and prespecified statistics, and interpret discrepancies as diagnostic evidence rather than a pass-or-fail test.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A posterior draw predicts 12 mg/L at a measurement time, and the additive Gaussian observation error has standard deviation 2 mg/L. If the simulated standard-normal error is $-0.5$, what replicated measurement is generated?
At the final two times, observed concentrations are 5.4 and 4.8 mg/L, and predicted means for a given draw are 4.9 and 4.3 mg/L. Calculate the mean residual. If 3 of 40 paired replicates have a statistic at least that large, calculate and interpret the estimated upper-tail area.
All 19 of 20 measurements lie inside their 90% pointwise predictive intervals. Does that establish that the complete trajectory is well reproduced? Explain.
Why should each replicated observation be generated from a posterior parameter draw and the observation model, rather than from a single posterior mean parameter value alone?

Solutions. 1. The replicate is $12+2(-0.5)=11$ mg/L. 2. The residuals are $5.4-4.9=0.5$ mg/L and $4.8-4.3=0.5$ mg/L, so their mean is $0.5$ mg/L. The estimated upper-tail area is $3/40=0.075$. Thus 7.5% of these paired replicates had a mean residual at least as large; this is a diagnostic comparison, not the probability the model is true. 3. No. Pointwise intervals assess individual times, not simultaneous coverage or all trajectory patterns. Inspect the time-course shape and targeted discrepancies as well. 4. Posterior draws propagate parameter uncertainty, while the observation model adds measurement variability and any modeled dependence. A single posterior mean can suppress both forms of variation and may not even correspond to a representative parameter combination in a nonlinear model.