Quantifying Monte Carlo Error in PTA

Why quantify Monte Carlo error?

A simulated probability of target attainment (PTA) is an estimate, not an exact value. For a fixed population PK model and a defined simulation protocol, the model implies a probability that a randomly generated virtual subject will meet the exposure target. A finite simulation estimates that probability using a finite number of simulated subjects. If the simulation is repeated with new random draws, the estimate will generally change, even though the model and target have not changed. Monte Carlo error describes this variation caused by using a finite number of simulations.

This error matters when a PTA is used to support a dose justification. For example, an estimated PTA of 90% from a few hundred simulated subjects may be too imprecise to distinguish a model-implied PTA near 88% from one near 92%. Reporting only the estimated proportion conceals that limitation. A Monte Carlo standard error and confidence interval show the precision attributable to simulation, but they do not capture every source of uncertainty in the analysis.

Definition: Suppose $M$ independent virtual subjects are simulated under a fixed population PK model, dose, and exposure-target definition. Let $I_i=1$ if subject $i$ meets the target and $I_i=0$ otherwise. If each subject has the same model-implied probability $p_{\mathrm{PTA}}$ of meeting the target, the estimated PTA is $\hat{p}_{\mathrm{PTA}}=K/M$, where $K=\sum_{i=1}^M I_i$. Under these assumptions, $K$ has a binomial distribution with parameters $M$ and $p_{\mathrm{PTA}}$.

The binary outcome must correspond to the target actually being evaluated. For a lower exposure threshold, for instance, $I_i=1$ when the simulated exposure for subject $i$ is at or above that threshold. Whether exposure is simulated with residual error, and which population or covariate distribution is represented, depend on the defined simulation protocol. Once that protocol is fixed, those choices define the probability being estimated; they do not change the calculation of binomial Monte Carlo error.

The binomial calculation assumes independent, identically distributed subject-level indicators: subjects are generated independently from the same target population under the same model and protocol, and each indicator is classified consistently. It also assumes that the simulations are not weighted or arranged using a sampling scheme that changes the variance of the proportion. These assumptions concern the simulation design and should be checked before applying the standard formula.

Monte Carlo standard error and interval estimation

For independent Bernoulli indicators, the variance of each $I_i$ is $p_{\mathrm{PTA}}(1-p_{\mathrm{PTA}})$. The variance of their average is that quantity divided by $M$. Since the unknown probability is not available, its estimate $\hat{p}_{\mathrm{PTA}}$ is substituted to estimate the standard error:

\widehat{SE}_{MC}(\hat{p}_{\mathrm{PTA}}) = \sqrt{\frac{\hat{p}_{\mathrm{PTA}}(1-\hat{p}_{\mathrm{PTA}})}{M}}

This is the estimated standard deviation of the PTA estimate over repeated simulations of size $M$, conditional on the fixed model and simulation protocol. It is measured in proportion units; multiply by 100 to express it in percentage points. It is not the between-subject standard deviation of exposure, nor the variability of individual target outcomes.

A common large-sample interval is the Wald interval, $\hat{p}_{\mathrm{PTA}}\pm 1.96\widehat{SE}_{MC}$ for a nominal 95% interval. It is easy to calculate but can perform poorly near zero or one, especially with small $M$ or few successes or failures. It can even extend below 0 or above 1. The Wilson score interval is generally preferable for reporting a binomial proportion, particularly when the observed proportion is near a boundary.

Definition: For a 95% Wilson score confidence interval, let $z=1.96$, the standard normal 97.5th percentile. The interval is $$ \frac{\hat{p}_{\mathrm{PTA}}+\frac{z^2}{2M}}{1+\frac{z^2}{M}} \ \pm\ \frac{z}{1+\frac{z^2}{M}} \sqrt{\frac{\hat{p}_{\mathrm{PTA}}(1-\hat{p}_{\mathrm{PTA}})}{M}+\frac{z^2}{4M^2}}. $$ It is an approximate interval for the fixed model-implied probability $p_{\mathrm{PTA}}$, based on binomial sampling.

The Wilson interval can be justified by inverting the score test for a proposed value $p_0$. Under the binomial model, the test statistic is $(\hat{p}_{\mathrm{PTA}}-p_0)/\sqrt{p_0(1-p_0)/M}$. Values of $p_0$ not rejected by the two-sided test satisfy

\frac{(\hat{p}_{\mathrm{PTA}}-p_0)^2}{p_0(1-p_0)/M}\le z^2.

Multiplying through by the positive denominator and collecting terms gives a quadratic inequality in $p_0$. Its two roots are the Wilson lower and upper limits. This inversion explains why the interval accounts for the fact that the standard error depends on the proposed probability, rather than simply adding and subtracting an estimated standard error as the Wald interval does.

Worked Example: PTA Near 72% With 1,000 Simulated Subjects

Suppose a hypothetical simulation uses $M=1{,}000$ independent virtual subjects. A prespecified exposure target is met by $K=720$ subjects. The estimated PTA is

$\hat{p}_{\mathrm{PTA}}=720/1{,}000=0.720$, or 72.0%.

The estimated Monte Carlo standard error is

$\widehat{SE}_{MC}=\sqrt{0.720(1-0.720)/1{,}000}=\sqrt{0.0002016}=0.01420$, or 1.420 percentage points.

For the Wilson interval, $z^2=1.96^2=3.8416$. Its center is \[ \frac{0.720+3.8416/(2{,}000)}{1+3.8416/1{,}000} = \frac{0.7219208}{1.0038416} \approx 0.7192. \] The half-width is \[ \frac{1.96}{1.0038416} \sqrt{\frac{0.720(0.280)}{1{,}000}+\frac{3.8416}{4(1{,}000)^2}} \approx 0.02779. \] Thus the approximate 95% Wilson interval is $0.7192\pm0.02779$, or approximately 0.6914 to 0.7470. Conditional on the fixed model and protocol, this method gives an interval with approximately 95% coverage over repeated simulations. It does not mean there is a 95% probability that the fixed probability lies in this particular interval.

Worked examples and practical interpretation

Worked Example: A Low Estimated PTA

In a second hypothetical simulation, $M=100$ independent subjects are generated and $K=8$ meet the target. Then $\hat{p}_{\mathrm{PTA}}=8/100=0.0800$, or 8.00%. The plug-in Monte Carlo standard error is \[ \sqrt{\frac{0.0800(0.9200)}{100}} =\sqrt{0.000736} =0.02713, \] or 2.713 percentage points.

For the 95% Wilson interval, the denominator is $1+3.8416/100=1.038416$. The center is \[ \frac{0.0800+3.8416/200}{1.038416} =\frac{0.099208}{1.038416} \approx 0.09554. \] The half-width is \[ \frac{1.96}{1.038416} \sqrt{\frac{0.0800(0.9200)}{100}+\frac{3.8416}{4(100)^2}} = \frac{1.96}{1.038416}\sqrt{0.00083204} \approx 0.05444. \] The interval is therefore approximately 0.0411 to 0.1500. The interval remains within the probability scale and is not centered at the observed proportion. With only eight simulated successes, the estimate is fairly imprecise relative to its size. A narrow-looking estimate should not be inferred from the point estimate alone.

Worked Example: Planning the Number of Simulations

Suppose the expected model-implied PTA is about 0.80, and the goal is an estimated Monte Carlo standard error no greater than 0.010, or one percentage point. Using the planning approximation \[ \sqrt{\frac{p(1-p)}{M}}\le 0.010, \] and substituting $p=0.80$, gives \[ M\ge\frac{0.80(0.20)}{0.010^2} =\frac{0.16}{0.0001} =1{,}600. \] Thus at least 1,600 independent simulations are indicated by this approximation. The achieved standard error will depend on the observed proportion, so this is a planning calculation, not a guarantee that every resulting estimate has an estimated standard error at or below 0.010.

If there is no reasonable prior estimate of the PTA and a conservative binomial planning value is wanted, $p(1-p)$ is largest at $p=0.50$, where it equals 0.25. For a one-percentage-point target, this gives $M\ge0.25/0.0001=2{,}500$. The conservative value is larger because proportions near 0.50 have the greatest binomial variance.

A related planning goal is to limit the approximate 95% interval half-width to $h$. Using the normal approximation gives $M\approx z^2p(1-p)/h^2$. For $p=0.80$ and $h=0.020$, this is $3.8416(0.16)/(0.020)^2=1{,}536.64$, so at least 1,537 simulations are suggested. This calculation is a useful approximation, not a substitute for checking the interval actually obtained, especially for rare target attainment.

What the interval does—and does not—measure

The PTA confidence interval quantifies simulation sampling error conditional on the specified model. It does not quantify uncertainty about the population PK parameter estimates, model structure, covariate distribution, target definition, or whether the model describes the intended population well. Those are distinct sources of uncertainty. A very large $M$ can make the Monte Carlo error small, but it cannot correct a misspecified model or make uncertain model parameters known.

It is also important to distinguish the simulated variation among subjects from uncertainty in the estimated proportion. Individual exposures may vary widely because of interindividual variability, covariates, and other components included in the simulation. That variation is part of the population distribution whose probability of target attainment is being estimated. The Monte Carlo standard error instead describes how much the estimated proportion varies because only $M$ subjects were simulated.

Assumptions, diagnostics, and alternatives

Before using the binomial standard error or Wilson interval, verify that each simulated subject contributes one binary indicator, that subjects are independently sampled from the same specified population, and that the denominator is the number of valid independent simulated subjects. Check for duplicated draws, unintended dependence, inconsistent target classification, and exclusions that alter the intended denominator. Reproducibility of the random-number seed makes results repeatable; it does not by itself establish independence or adequate precision.

The binomial formulas need modification when the estimator is not a simple unweighted proportion of independent, identically distributed indicators. For example, importance sampling, unequal weights, clustered simulation, or a fixed stratified allocation can change the variance. In such cases, use a variance estimator matched to the sampling design or assess simulation variability using appropriate independent batches. Repeated batches are useful only when they represent independent repetitions of the intended simulation procedure.

The Wald interval is a familiar alternative and can be adequate when both successes and failures are numerous and the proportion is not near a boundary, but the Wilson interval usually behaves better. An exact binomial interval is another option, particularly when counts are small; it avoids the normal approximation but can be conservative and wider. The choice among intervals does not change what is being estimated: under a fixed model, all are aimed at the same model-implied PTA, with different finite-sample properties.

Approach	Useful feature	Limitation
Plug-in Monte Carlo standard error	Simple measure of the estimated proportion's simulation precision	Does not itself give an interval; relies on the binomial design
Wald interval	Easy to calculate and explain	Can have poor coverage and limits outside 0 to 1 near boundaries
Wilson interval	Score-based interval with better boundary behavior in common settings	Approximate coverage; still assumes binomial sampling
Exact binomial interval	Uses the binomial distribution directly	Can be conservative and wider than necessary

Common mistake and key takeaway

A common mistake is to report a narrow Monte Carlo interval as evidence that the dose recommendation is certain. The interval only addresses finite simulation error under the fixed model and protocol. It does not include parameter uncertainty or model misspecification. Another mistake is to increase the number of virtual subjects to address every kind of uncertainty: increasing $M$ reduces simulation error, but does not resolve uncertainty about model inputs or assumptions.

Key takeaway: Report PTA together with the number of independent simulated subjects and a Monte Carlo precision measure. Interpret the interval as uncertainty from finite simulation, conditional on the specified population PK model and simulation protocol—not as a complete uncertainty interval for a clinical dose decision.

Check Your Understanding

Try each question, then compare with the solutions below.

A hypothetical simulation has $M=400$ independent subjects and 300 target successes. Calculate the PTA and its plug-in Monte Carlo standard error in percentage points.
Why does a 95% Monte Carlo confidence interval not represent a 95% probability that the fixed model-implied PTA lies inside the observed interval?
For an expected PTA of 0.50, estimate how many independent simulations are needed for an approximate 95% half-width of 0.02.
Name two situations in which the simple binomial standard error may not be appropriate.

Solutions. (1) $\hat{p}_{\mathrm{PTA}}=300/400=0.750$, or 75.0%. The standard error is $\sqrt{0.75(0.25)/400}=\sqrt{0.00046875}=0.02165$, or 2.165 percentage points. (2) The model-implied PTA is fixed under the stated model and protocol. The 95% coverage statement describes the long-run performance of the interval procedure over repeated simulations, not a probability assigned to the fixed PTA after observing one interval. (3) $M\approx1.96^2(0.50)(0.50)/(0.02)^2=0.9604/0.0004=2{,}401$ independent simulations. (4) Examples include unequal simulation weights and dependent or clustered subject draws; stratified allocation can also require a design-specific variance calculation.