Why VPC Quantiles Need Uncertainty Bands
A visual predictive check (VPC) compares quantiles of observed concentrations with quantiles calculated from datasets simulated under a population PK model. But a sample quantile is not a fixed value: it changes when the subjects observed in a study change. Simulated quantiles also vary from one simulated dataset to another, and estimated simulation reference limits change when the number of simulated datasets is finite. These are related but distinct sources of variation.
The distinction matters for interpretation. Uncertainty around an observed quantile describes how imprecisely that study estimates its underlying quantile. The spread of simulated quantiles describes the variation expected across datasets generated under the model and the study design. Monte Carlo uncertainty describes how imprecisely a finite number of simulations estimates that simulated distribution. A clear VPC should not treat all three as one kind of interval.
What Each Band Represents
Suppose the bin contains \(n_b\) observed concentrations. Write the observed quantile as \(\hat Q_{\mathrm{obs},p,b}\), where \(b\) identifies the bin and \(p\) the target probability. A confidence interval for this value can be estimated by resampling the observed sampling units and recalculating the quantile. For longitudinal PK data, the sampling unit is usually the subject: resample whole subject profiles, not individual concentrations, so within-subject dependence is retained.
Now generate \(R\) replicate datasets using the same design and calculate the target quantile in each one, giving \(\hat Q_{1,p,b},\ldots,\hat Q_{R,p,b}\). The spread of these \(R\) values describes how much the quantile can vary across datasets under the model. For example, the 5th and 95th percentiles of the replicate quantiles can form a pointwise 90% reference interval for the dataset-level quantile. This is a predictive reference interval, not a confidence interval for the quantile from the single observed dataset.
The endpoints of that reference interval are themselves estimated from only \(R\) replicates. Their Monte Carlo uncertainty can be displayed as confidence intervals around the endpoints or assessed numerically. Increasing \(R\) reduces this computational uncertainty; it does not remove the underlying variation among simulated datasets. Conversely, collecting more subjects can improve the observed quantile estimate, but does not improve the precision of an already-fixed simulation calculation.
These are pointwise statements for a specified quantile and bin. Bands across many bins and quantiles do not automatically provide a simultaneous confidence statement about the whole VPC. Also, unless parameter uncertainty is deliberately included, simulations conditional on a fitted model treat its parameter values as fixed. The resulting bands do not account for uncertainty in those parameter estimates.
Calculating Uncertainty
A practical method for an observed quantile is the subject-level bootstrap. Resample subjects with replacement, keeping each selected subject’s relevant observations together; apply the same binning and quantile definition; and repeat this process many times. The resulting bootstrap quantiles approximate the sampling distribution of the observed quantile. A percentile confidence interval uses the appropriate lower and upper percentiles of this bootstrap distribution. The interpretation relies on the resampling scheme representing the study’s sampling process and on there being enough independent subjects for the bootstrap distribution to be informative.
For a continuous outcome with density \(f\) that is positive and reasonably smooth at the population \(p\)-quantile \(Q_p\), an approximate standard error for an independent-sample quantile is
This formula makes two sources of imprecision visible: smaller \(n\) increases uncertainty, and a lower density near the target quantile increases it. The density must be expressed in inverse outcome units, so the standard error has the same units as the concentration. The formula is not a substitute for accounting for clustered observations: with repeated measurements from the same subject, a subject-level bootstrap or an appropriate cluster-aware variance estimate is needed.
The same approximation helps quantify Monte Carlo error for a simulated reference limit. Let \(G_b\) be the distribution of replicate quantiles in bin \(b\), and let \(q_\gamma\) be its \(\gamma\)-quantile. If the replicate quantiles are independent draws from a continuous distribution with density \(g_b(q_\gamma)>0\), then
Here \(R\), not the number of observed subjects, controls Monte Carlo precision. The density is that of the replicate quantiles near the reference limit. It can be estimated from the simulated values, although density estimates can be unstable in the tails. A bootstrap that resamples whole replicate datasets is another practical way to assess endpoint precision. Neither approach turns Monte Carlo uncertainty into biological or sampling variability.
Worked Example: Bootstrap Interval for an Observed Median
Suppose a hypothetical VPC bin contains one concentration from each of 60 independent subjects. The observed median is 5.10 mg/L. Investigators resample subjects with replacement 2,000 times and calculate the median in each bootstrap sample. Under a prespecified percentile convention, the 2.5th and 97.5th percentiles of the 2,000 bootstrap medians are 4.40 and 5.80 mg/L.
The estimated 95% bootstrap confidence interval for the observed median is therefore \([4.40,\ 5.80]\) mg/L. Its width is \(5.80-4.40=1.40\) mg/L. Under repeated sampling and the stated bootstrap assumptions, the procedure is intended to cover the corresponding population median about 95% of the time. It does not mean there is a 95% probability that the fixed population median lies in this particular interval.
The calculation assumes that the 60 subjects are independent sampling units and that the observed values represent the bin’s sampling process. If each subject contributed several concentrations to the bin, the resampling unit should instead be the whole subject profile. Treating those concentrations as independent could make the interval too narrow.
Worked Example: Monte Carlo Uncertainty in a Simulated Reference Limit
For a hypothetical bin and target quantile, \(R=500\) simulated datasets produce 500 replicate quantiles. Their 5th percentile is estimated as 3.80 mg/L, and the estimated density of replicate quantiles near that limit is \(0.10\) per mg/L. The simulated 5th percentile describes a lower reference limit for quantiles across datasets; it is not an interval for the observed quantile.
Using the Monte Carlo standard-error approximation,
An approximate 95% Monte Carlo interval for the estimated reference limit is \(3.80 \pm 1.96(0.0975)\) mg/L. Since \(1.96(0.0975)=0.191\), the interval is approximately \([3.61,\ 3.99]\) mg/L. This interval expresses uncertainty from estimating the 5th-percentile limit with 500 replicates; it does not replace the predictive spread of replicate quantiles.
The density and normal approximation are assumptions here, not guaranteed facts. In particular, tail quantiles can be poorly approximated when \(R\) is small or the replicate-quantile distribution is irregular. In that situation, more simulations and a replicate-level bootstrap can help assess stability.
Worked Example: Why More Replicates Do Not Narrow the Predictive Spread
In a hypothetical VPC bin, suppose the estimated 5th percentile of replicate quantiles is 4.10 mg/L, the density of replicate quantiles there is estimated as \(0.08\) per mg/L, and the estimated 95th percentile is 7.20 mg/L. The predictive reference interval has width \(7.20-4.10=3.10\) mg/L. The Monte Carlo standard error for the lower endpoint depends on \(R\).
With \(R=200\), the approximation gives
With \(R=800\), it gives
Four times as many replicates halve the approximate Monte Carlo standard error because the standard error scales as \(1/\sqrt{R}\). They can make the estimated 4.10 and 7.20 mg/L reference limits more precise, but they do not make the 3.10 mg/L predictive interval narrower in the underlying model. The latter reflects variability among simulated datasets.
Diagnostics, Alternatives, and Limitations
Before interpreting bands, check that each observed and simulated quantile uses the same bin definition, target probability, concentration handling, and quantile convention. The simulated datasets should reproduce the relevant design, including the sampling schedule and the number of subjects. For longitudinal data, uncertainty calculations should preserve subject-level dependence. If observations are sparse within a bin, a bootstrap interval may be highly discrete or unstable; the small-sample resolution of the VPC is then a substantive limitation, not something that a smooth-looking band fixes.
Check the stability of simulated reference limits by increasing \(R\), examining how the endpoint estimates change, and estimating their Monte Carlo error. A large change after adding simulations suggests the original reference limits were imprecise. Also inspect the replicate quantile distribution: a highly skewed, discrete, or multimodal distribution can undermine a symmetric standard-error approximation. Bootstrap resampling of whole replicate datasets can avoid that approximation, though it still cannot create information absent from a small \(R\).
The main practical alternatives have different purposes. The subject-level bootstrap directly estimates observed-quantile sampling uncertainty and can accommodate dependence if resampling is correctly clustered. The density-based standard-error formula is fast and useful as a diagnostic, but depends on smoothness and a reliable density estimate. Replicate-level quantiles or bootstrap resampling estimate the precision of simulated reference limits. Simply increasing the number of simulated datasets improves Monte Carlo precision but does not estimate uncertainty in the observed data. These approaches are complementary rather than interchangeable.
A common mistake is to call the spread of simulated quantiles “Monte Carlo error.” That is incorrect: it is the variability across datasets generated under the model. Monte Carlo error is the additional uncertainty in estimating that spread or its reference limits from a finite \(R\). Another mistake is to infer model failure because the observed quantile falls just outside a simulated band without considering uncertainty in the observed quantile, bin resolution, and the pointwise nature of the comparison. Bands are diagnostic summaries; they are not, by themselves, a formal test of overall model adequacy.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A bin has repeated concentrations from each subject. What should generally be resampled to estimate uncertainty in its observed quantile, and why?
- What does the spread of replicate quantiles across simulated datasets represent?
- If \(R\) increases from 100 to 400, by what factor does the approximate Monte Carlo standard error change, assuming the reference quantile and density remain the same?
- Why can a density-based quantile standard error be unreliable for a sparse bin with many tied values?
Solutions. (1) Resample whole subjects, retaining their relevant observations together, because concentrations from the same subject are generally dependent. (2) It represents dataset-to-dataset sampling variability under the simulated model and design, not finite-\(R\) Monte Carlo error. (3) The standard error is proportional to \(1/\sqrt{R}\), so it is multiplied by \(\sqrt{100/400}=1/2\): it is approximately halved. (4) The approximation assumes a continuous distribution with a positive, reasonably smooth density at the quantile. Many ties and sparse observations make that density assumption unsuitable or difficult to estimate, so the approximation may not describe the quantile’s uncertainty well.