Common Monte Carlo Reporting and Interpretation Errors

Why Monte Carlo Results Are Easy to Misstate

Monte Carlo propagation turns specified input distributions and a model into a distribution of simulated outputs. A mistake in describing that distribution can change the scientific meaning of a result even when the simulation code is correct. In PBPK applications, frequent errors include calling an input range a prediction interval, presenting a poorly resolved tail quantile as precise, sampling correlated inputs independently, or describing uncertainty in estimated parameters as variability between patients.

A useful starting point is to distinguish three things: the distributions or ranges assigned to model inputs, the distribution of outputs induced by those inputs and the model, and the population or individual variation the simulation is intended to represent. They are related, but they are not interchangeable. The introductory material on prediction intervals distinguishes predictive intervals from confidence and credible intervals; here the emphasis is on avoiding misleading Monte Carlo interpretation.

Definition: Let \(X\) denote a simulated model output, such as AUC in \(\text{mg}\cdot\text{h}/\text{L}\). A prediction distribution is the distribution of \(X\) under specified model structure, input distributions, dependencies, and simulation conditions. A prediction interval summarizes a specified portion of that output distribution. An input range describes possible or assumed values for an input; it is not itself an interval for \(X\).

A Monte Carlo prediction distribution is conditional on how it was constructed. Its interpretation depends on which inputs were uncertain, the distributions assigned to them, their dependence, whether between-subject variation was included, and whether numerical error was controlled. A result should not be given a broader interpretation than those assumptions support.

Input Ranges Are Not Output Intervals

An input range can be a plausible range, a measurement range, a parameter support, or a range used in a sensitivity analysis. None of these automatically gives the probability that a model output lies between particular values. The model maps inputs to outputs, often nonlinearly, and the output distribution also depends on the joint distribution of the inputs.

For example, a stated range of clearance values does not determine a concentration interval by itself. Dose, other model parameters, covariates, the time point or summary being reported, and the probability distributions assigned to the inputs all matter. A prediction interval is calculated from simulated outputs, not copied from an input range or obtained by relabeling its endpoints.

Worked Example: Separating a Clearance Range From an Output Interval

Suppose a hypothetical PBPK analysis assigns clearance values between 4 and 8 L/h to an uncertain input. The model output of interest is a concentration, \(X\), at a specified time. For illustration, suppose the simulated output distribution is approximated by \(X\sim N(10,2^2)\) mg/L. This normal approximation is an assumption for this calculation, not a consequence of the clearance range.

For a continuous normal distribution, the central 90% interval uses the 5th and 95th percentiles. The corresponding standard-normal quantiles are approximately \(-1.645\) and \(1.645\), so:

10 \pm 1.645(2) = 10 \pm 3.290

Thus the interval is approximately 6.71 to 13.29 mg/L. These endpoints are output quantiles under the stated output-distribution assumption. They are not the clearance values 4 and 8 L/h, and they cannot be derived from those two input endpoints alone. In an actual Monte Carlo analysis, the interval would be computed from the simulated concentration outputs, with uncertainty about the estimated quantiles assessed separately.

Quantile terminology also needs care when the output distribution has atoms, as can occur with discrete outputs, censoring rules, or model logic that assigns positive probability to a fixed value. Using the conventional quantile definition \(Q(p)=\inf\{x:F(x)\geq p\}\), the interval \([Q(0.05),Q(0.95)]\) need not contain exactly 90% of probability. For a continuous distribution, it has 5% probability below the lower endpoint and 5% above the upper endpoint. With atoms, the actual tail probabilities depend on the probability at the endpoints and should be calculated or reported using a clearly stated convention. If \(X=0\) with probability 1, both quantiles are 0 and the interval contains all the probability, not exactly 90%.

Tail Quantiles Can Be Much Less Stable Than They Look

A Monte Carlo estimate of a quantile is based on the simulated outputs near its rank in the sorted sample. For \(M\) independent draws, an estimated 1st percentile is informed by roughly the lowest \(0.01M\) outputs. If \(M=1{,}000\), that is about 10 observations. A decimal-rich endpoint does not mean that the tail is well resolved.

One way to see the sampling limitation is to consider the number of draws below the true \(p\)-quantile. Under independent sampling from a continuous distribution, this count has a binomial distribution with parameters \(M\) and \(p\). The estimated fraction below the quantile therefore has standard error \(\sqrt{p(1-p)/M}\). Quantile uncertainty in the output’s units also depends on the density near that quantile.

Justification: For independent draws from a continuous distribution with density \(f\) that is positive at the true quantile \(q_p\), the large-sample standard error of the sample \(p\)-quantile is approximately \(\sqrt{p(1-p)}/(\sqrt{M}f(q_p))\). The numerator reflects variation in the fraction of draws below \(q_p\); dividing by the local density converts a probability-scale error into an error on the output scale. This approximation is not reliable when the density is near zero, the distribution has an atom, or the tail sample is too small.

Worked Example: Interpreting a Noisy Tail Estimate

Consider \(M=1{,}000\) independent hypothetical simulations and a target lower-tail probability \(p=0.01\). The expected number of draws below the true 1st percentile is \(Mp=1{,}000(0.01)=10\). The standard error of the simulated fraction below that quantile is:

\sqrt{\frac{0.01(0.99)}{1{,}000}} = 0.00315

This is about 0.315 percentage points. It describes uncertainty in the estimated tail fraction at the fixed true quantile; it is not an output-scale confidence interval for the quantile itself.

Suppose, additionally, that the output density near the true 1st percentile is approximately \(0.002\) per mg/L. The approximate standard error of the estimated quantile is then:

\frac{\sqrt{0.01(0.99)}}{\sqrt{1{,}000}(0.002\ \text{per mg/L})} = \frac{0.0995}{0.0632\ \text{per mg/L}} \approx 1.57\ \text{mg/L}

Under the approximation’s assumptions, this indicates substantial uncertainty in the tail endpoint. Reporting that endpoint to three decimal places would imply unjustified precision. Increasing the number of independent simulations can improve Monte Carlo precision, but it does not fix an inappropriate input distribution, a model limitation, or numerical error in the model outputs.

Useful diagnostics include checking how the quantile changes as the simulation count increases, examining the number of draws in the relevant tail, and comparing independent simulation batches. A quantile that moves appreciably across batches is not made reliable by rounding it neatly. For rare-event probabilities, the same basic warning applies: if very few events occur, the estimated probability can be imprecise even when the total number of runs appears large.

Dependence Between Inputs Must Be Preserved

Sampling each uncertain input from its correct marginal distribution does not guarantee that the joint input distribution is correct. If two inputs are correlated, drawing them independently changes their joint distribution. That can change the output variance and tail behavior, particularly when the model responds to both inputs in the same part of the calculation.

For random variables \(X_1\) and \(X_2\), the variance of their sum is:

\operatorname{Var}(X_1+X_2) = \operatorname{Var}(X_1)+\operatorname{Var}(X_2) +2\operatorname{Cov}(X_1,X_2)

The covariance term is zero only when the variables are uncorrelated; independence is sufficient, but not necessary, for zero covariance. Omitting a real dependence structure can therefore widen or narrow the simulated prediction distribution. Whether it does so in a particular model depends on the signs of the dependence and on how the output responds to the inputs.

Worked Example: The Effect of Ignoring Positive Dependence

For a deliberately simple illustration, let an output be \(Y=X_1+X_2\), measured in arbitrary output units. Suppose \(X_1\) and \(X_2\) each have mean 5 units and standard deviation 1 unit, and suppose their joint distribution is normal with correlation \(\rho=0.8\). Then \(\operatorname{Cov}(X_1,X_2)=0.8(1)(1)=0.8\) square units, giving:

E(Y)=5+5=10,\qquad \operatorname{Var}(Y)=1+1+2(0.8)=3.6

The standard deviation is \(\sqrt{3.6}=1.897\) units. Under the assumed normal distribution, a central 90% interval is approximately:

10\pm1.645(1.897) = 10\pm3.120 = [6.880,\ 13.120]\ \text{units}

If the same marginal distributions were instead sampled independently, the covariance would be set to zero. The resulting variance would be \(1+1=2\), and the corresponding interval would be \(10\pm1.645\sqrt{2}=[7.674,\ 12.326]\) units. In this example, independently sampling the inputs produces a narrower interval because it removes positive dependence. The example assumes a linear output and a bivariate normal joint distribution; a PBPK model can be nonlinear, so the direction and size of the effect should be assessed by propagating the appropriate joint distribution through that model.

Parameter Uncertainty Is Not Patient-to-Patient Variability

Uncertainty about a model parameter concerns limited knowledge of a parameter’s value. Between-subject variability describes genuine differences among individuals under a population model. The two sources can both contribute to a predictive distribution, but they answer different questions. Calling one by the name of the other can make a model appear to predict a population spread that was never included, or can conceal how much of a prediction’s uncertainty reflects limited parameter knowledge.

A parameter draw is often shared across a simulated population because it represents uncertainty about a common model quantity. By contrast, individual variability is represented by individual-level draws from a population distribution. The appropriate structure depends on the intended prediction: uncertainty about a population mean is not the same target as the distribution of outcomes for a new patient.

Worked Example: Combining Parameter Uncertainty and Individual Variation

Suppose a hypothetical concentration for a new patient is represented by \(Y=\theta+\eta\), in mg/L. Here \(\theta\) is an uncertain population-level concentration with \(\theta\sim N(10,1^2)\) mg/L, and \(\eta\) is individual variation with \(\eta\sim N(0,2^2)\) mg/L. Assume the two components are independent. Then:

E(Y)=10+0=10\ \text{mg/L},\qquad \operatorname{Var}(Y)=1^2+2^2=5\ (\text{mg/L})^2

The predictive standard deviation for a new patient is \(\sqrt{5}=2.236\) mg/L. Since \(Y\) is normal under these assumptions, its central 95% prediction interval is approximately:

10\pm1.96(2.236) =10\pm4.383 =[5.617,\ 14.383]\ \text{mg/L}

If only uncertainty in \(\theta\) were propagated, the 95% interval would instead be \(10\pm1.96(1)=[8.04,\ 11.96]\) mg/L. That interval describes uncertainty in the population-level quantity in this simplified example; it does not describe the full distribution for a new patient. Conversely, including \(\eta\) does not turn parameter uncertainty into patient variability: the two components retain their distinct interpretations even when combined in one prediction.

Diagnostics, Alternatives, and Common Pitfalls

When reviewing a Monte Carlo result, first identify what was sampled and what the output summarizes. Then check whether the reported interval or probability was calculated from simulated outputs, whether the number of draws supports the claimed tail precision, whether known input dependence was represented, and whether individual variability was actually included. If outputs show spikes, repeated values, or a discrete component, do not assume continuous-distribution tail behavior. Calculate the probability below and above the reported endpoints under the chosen quantile convention.

Different summaries serve different purposes. A mean or median describes a central location but does not show the spread or the tails. A prediction interval summarizes a specified portion of an output distribution, while a one-sided bound focuses on one tail. A point estimate with a Monte Carlo error assessment addresses simulation precision for that estimate, not all uncertainty in the model or inputs. No summary can compensate for a mismatch between the simulated quantity and the question being asked.

Wrong reasoning: “The input parameter ranges from 4 to 8 L/h, so the prediction interval is 4 to 8 L/h.” Correction: Those are input values with units of clearance. Compute an interval from the simulated output values, which have the output’s units.
Wrong reasoning: “The 1st percentile is reported to three decimal places, so it is precise.” Correction: Precision depends on simulation and model conditions, not the number of displayed digits. Check tail counts and stability across batches.
Wrong reasoning: “Each input distribution is correct, so independently sampling all inputs is correct.” Correction: Correct marginals do not establish a correct joint distribution. Represent justified dependence.
Wrong reasoning: “The spread from parameter draws is the amount patients differ from one another.” Correction: Parameter uncertainty and between-subject variability are distinct sources; state which are included and what population the output represents.
Wrong reasoning: “The 5th-to-95th-percentile interval always has exactly 90% probability.” Correction: This exact-tail description applies to continuous distributions. With atoms, state the quantile convention and report actual tail probabilities; the interval may contain more than 90% probability.

Key takeaway: A Monte Carlo result is meaningful only together with its target quantity and simulation assumptions. Distinguish input ranges from output intervals, treat tail precision as an estimable property rather than a formatting choice, preserve relevant dependence, and label parameter uncertainty and patient variability separately.

Check Your Understanding

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

A parameter is sampled between 2 and 6 mg/L, and a model produces a concentration output. Does the parameter range define a prediction interval for concentration? Why?
In 500 independent simulations, approximately how many draws are expected below the true 2nd percentile? What does that imply about the information available in that tail?
Two inputs each have variance 4, with covariance \(-1\). What is the variance of their sum? How would independent sampling change that variance?
A simulation includes uncertainty in a population parameter but no individual-level random effect. Is its output spread a measure of patient-to-patient variability?
An output is always zero. Using \(Q(p)=\inf\{x:F(x)\geq p\}\), what are its 5th and 95th percentiles, and how much probability lies between them?

Solutions. (1) No. The range describes an input, while a prediction interval must be calculated from output draws and has output units. (2) \(500(0.02)=10\) draws are expected. The estimate depends on a small tail sample, so its precision should not be assumed from the total simulation count alone. (3) \(\operatorname{Var}(X_1+X_2)=4+4+2(-1)=6\). Independent sampling sets the covariance to zero and gives variance \(8\), changing the propagated spread. (4) No. It describes output variation induced by parameter uncertainty under the simulation setup; individual-level variability was not included. (5) Both quantiles are zero, and the interval \([0,0]\) contains probability 1. The tails outside the endpoints contain zero probability, not 5% each.