Quantifying Monte Carlo Error in Prediction Quantiles

Why Quantiles Need Their Own Monte Carlo Error

A PBPK simulation can produce a distribution of predictions by repeatedly drawing uncertain inputs and propagating them through the model. A reported median or prediction-distribution boundary is then estimated from a finite sample of simulated outcomes. Increasing the number of draws may change that estimate, even when the input distributions and PBPK model are unchanged. Monte Carlo error is this finite-simulation uncertainty in the estimated quantity; it is not the biological or predictive spread of the simulated outcomes themselves.

The error is especially important for quantiles. A mean averages information across the full sample, whereas a quantile depends on the outcomes near a particular rank. The uncertainty of a quantile estimate therefore depends both on sample size and on how densely outcomes lie near the target quantile. A tightly concentrated distribution near the median can yield a precise median, while a sparse tail can make a boundary estimate imprecise even with many simulations.

The empirical quantile convention above is one standard choice; software can use other interpolation rules, particularly when \(np\) is not an integer. The convention should be stated when reproducibility or small-sample comparisons matter. The large-sample approximations below apply to common quantile conventions under their regularity conditions.

Large-Sample Quantile Error

Suppose \(Y_1,\ldots,Y_n\) are independent draws from the same continuous distribution, and that the density \(f\) is positive and sufficiently smooth near \(q_p\). In particular, \(F(q_p)=p\) and \(f(q_p)>0\). For large \(n\), the quantile estimate is approximately normal around \(q_p\), with variance determined by the local density:

The units of \(f(q_p)\) are inverse units of \(Y\), so the MCSE has the same units as the prediction. For instance, if \(Y\) is a concentration in mg/L, its quantile MCSE is also in mg/L. This standard error quantifies simulation precision, not the width of a predictive interval.

This approximation requires more than a large total simulation count. The draws must be independent, or the dependence must be handled with an appropriate method; the target distribution must be adequately represented; and the density near the target quantile must be positive and reasonably smooth. A point mass, a sharp discontinuity, or a density close to zero can invalidate the approximation or make the quantile intrinsically difficult to estimate.

Estimating Error From Order-Statistic Spacings

The density \(f(q_p)\) is usually unknown. One practical approach estimates it from the span of nearby ordered outcomes. Let \(k=\lceil np\rceil\), choose a positive integer \(r\) such that \(k-r\) and \(k+r\) are valid indices, and estimate the local density by the probability mass between those ranks divided by their outcome-value span:

The spacing is an estimate of how much outcome value corresponds to a probability change of about \(2r/n\). A larger \(r\) gives a more stable spacing but averages density over a wider region; a smaller \(r\) is more local but can be noisy. The choice should be checked for sensitivity rather than treated as automatic.

Estimating Error With Independent Batches

A second approach divides independent simulation draws into \(B\) batches of equal size \(m\), so \(n=Bm\), and calculates the same quantile separately in each batch. If the batch estimates are \(\hat q_{p,1},\ldots,\hat q_{p,B}\), calculate their sample standard deviation \(s_q\). The batch-based estimate of the standard error of their average is:

This directly estimates the standard error of the average of batch quantiles. Under the large-sample approximation, each batch quantile has variance about \(p(1-p)/(m f(q_p)^2)\). Dividing this variance by \(B\) gives \(p(1-p)/(Bm f(q_p)^2)\), the same leading-order variance as the pooled-sample quantile. Thus the formula is a useful approximation to the pooled quantile MCSE when batches are independent and sufficiently large. The average of batch quantiles is not generally identical to the pooled quantile, especially for small batches or nonlinear quantile behavior.

Diagnostics, Alternatives, and Limitations

Order-statistic and batch methods address the same target but expose different weaknesses. The spacing method uses all outcomes in one sample and can be informative when a reasonable local density estimate is available. Its weakness is sensitivity to the spacing window. Batch estimates avoid explicit density estimation and provide a direct picture of between-batch variation, but require enough independent batches and enough draws within each batch to estimate the quantile meaningfully. Agreement between methods is reassuring, but neither method can compensate for a misspecified PBPK model or an unrepresentative input distribution.

Check the following before relying on an MCSE:

• Local spacing stability. Recalculate the spacing estimate using several reasonable values of \(r\). Large changes suggest that the density estimate is unstable or the distribution changes rapidly near the quantile.
• Batch stability. Inspect the batch quantiles for outliers, trends with simulation order, or much greater spread than expected. Such patterns can indicate dependence, nonstationarity, or an implementation problem.
• Effective information near the target. For a quantile with probability \(p\), there are approximately \(np\) outcomes below it and \(n(1-p)\) above it. If either count is small, the estimate may be sensitive to a few observations. This is a warning about quantile precision, not a rare-event probability analysis.
• Target and units. Verify that batches use the same model, input distributions, outcome definition, and quantile convention. Report the MCSE in the same units as the prediction.

A common mistake is to report the width of the simulated prediction distribution as the Monte Carlo error of its median or boundary. These answer different questions: the distribution width describes variation in the predicted outcomes, whereas the MCSE describes how much the estimated quantile would vary across repeated finite simulations. Another mistake is to assume that a very large \(n\) always ensures a precise tail boundary. The MCSE also depends on \(f(q_p)\), which may be small in a sparse tail.

Both methods have limitations. The spacing approximation can fail near a support boundary, at a discontinuity, or in a region with very low density; a spacing-based density estimate can also be biased by curvature. Batch estimates are unreliable if batches are too small, too few, or not independent. In either case, a reported MCSE is an estimated standard deviation, not a confidence interval for the true quantile and not a guarantee that the estimate lies within one MCSE of its target. For Monte Carlo workflows that use dependent or stratified sampling, the simple independent-draw formulas should not be applied without accounting for that design.