Monte Carlo Precision in Predictive Checks

Why Simulation Count Matters

A predictive check uses simulated data to estimate what summaries the fitted population PK model would produce under the observed design. The simulation process is random: if the same model and design are simulated again, the resulting data and summaries will usually differ. With a finite number of simulations, a displayed reference value is therefore an estimate, not an exact property of the model.

This variation is Monte Carlo error. It is distinct from uncertainty in the estimated model parameters, uncertainty due to the finite observed study, and model misspecification. In the calculations below, treat the fitted model, its parameter values, and the design being simulated as fixed. The question is how precisely a finite set of random simulations estimates a target under those fixed conditions.

Definition: Monte Carlo error is the random difference between a simulation-based estimate calculated from a finite number of random draws and the target quantity it estimates, conditional on the model and design used to generate those draws. Its standard deviation over repeated simulation runs is the Monte Carlo standard error (MCSE).

Increasing the number of independent simulations generally reduces MCSE, but it does not correct a model that predicts the wrong distribution. A narrow simulation-error interval means that the target has been estimated precisely under the simulation setup; it does not establish that the model is adequate.

Standard Errors and Confidence Intervals for Simulation Estimates

Suppose \(Z_1,\ldots,Z_M\) are independent draws of a quantity from the simulation target, and the aim is to estimate its expectation. The estimate is the average of the draws. If the target has finite variance, the variance of that average is the variance of one draw divided by \(M\). In practice, estimate the variance of one draw with the sample variance \(s_Z^2\).

\widehat{\theta}=\bar{Z}=\frac{1}{M}\sum_{m=1}^{M}Z_m, \qquad \widehat{SE}_{MC}(\bar{Z})=\frac{s_Z}{\sqrt{M}}.

For a sufficiently large number of independent draws, an approximate 95% Monte Carlo confidence interval is \(\bar{Z}\pm 1.96\widehat{SE}_{MC}\). With fewer independent draws and an approximately normal sampling distribution, a \(t\) interval using \(M-1\) degrees of freedom is often more suitable. These intervals describe repeated-simulation uncertainty for the simulation target. They are not intervals for an unknown model parameter.

A quantile requires a different calculation from a mean. For \(M\) independent draws from a continuous distribution with cumulative distribution function \(F\), let \(q_p\) be its \(p\)th quantile, and let \(f(q_p)\) be the density at that quantile. The empirical quantile has approximate variance \(p(1-p)/(M f(q_p)^2)\). The density matters: if observations are sparse around the quantile, small changes in the empirical cumulative distribution can move the quantile substantially.

\widehat{SE}_{MC}(\widehat{q}_p) \approx \frac{\sqrt{p(1-p)}}{\sqrt{M}\,f(q_p)}.

This approximation follows from the fact that the empirical cumulative distribution at \(q_p\) is an average of \(M\) indicator variables, each equal to 1 when a draw is at or below \(q_p\). Its variance is approximately \(p(1-p)/M\). Near \(q_p\), a small change \(\Delta F\) in cumulative probability corresponds to a change \(\Delta q\approx\Delta F/f(q_p)\) in the quantile. Applying that local conversion gives the variance above. In applications, \(f(q_p)\) is not known and must be estimated; the approximation can be unreliable when the density estimate is unstable or the distribution is discrete.

Worked Example: Approximate MCSE for a Simulated VPC Quantile

Suppose 1,000 independent simulated concentrations at a fixed design point are used to estimate the median concentration in a VPC. The estimated median is 10.0 mg/L, and a density estimate at the median is 0.20 per mg/L. Use the quantile approximation to assess simulation precision.

For the median, \(p=0.5\). Substituting \(M=1{,}000\) and \(f(q_p)=0.20\) per mg/L gives:

\widehat{SE}_{MC} \approx \frac{\sqrt{0.5(1-0.5)}}{\sqrt{1{,}000}(0.20)} = \frac{0.5}{31.62(0.20)} = 0.0791\ \text{mg/L}.

An approximate 95% simulation-error interval is \(10.0\pm1.96(0.0791)\), or 9.85 to 10.2 mg/L after rounding. Under the stated independent-draw and density-approximation assumptions, this interval describes the Monte Carlo precision of the estimated median. It does not account for uncertainty in the fitted model or the observed data.

VPC Quantiles From Replicate Datasets

In a VPC, simulated records from the same replicate dataset may be dependent, particularly when a dataset contains repeated measurements on the same individuals. It is therefore important not to assume automatically that every simulated record is an independent draw. One practical approach is to retain the replicate-dataset structure and calculate the summary of interest separately within each simulated dataset.

For example, if \(Q_m\) is the median concentration in a specified time bin for replicate dataset \(m\), the average of \(Q_1,\ldots,Q_M\) estimates the expected replicate-specific median under the simulation setup. The MCSE of this average is the sample standard deviation of the replicate-specific medians divided by \(\sqrt{M}\), provided the replicate datasets are independent and the replicate-specific medians have finite variance.

This target is the mean of replicate-specific medians. It is not necessarily the same as the median obtained by pooling all simulated records. Pooling can give more weight to datasets or individuals that contribute more records, and within-dataset dependence affects its precision. State which target is being estimated, and use a method that respects the independent replicate datasets when assessing its MCSE.

Worked Example: MCSE for Replicate-Specific VPC Medians

In a hypothetical VPC, each of 10 independent simulated datasets produces one median concentration in a time bin. The replicate-specific medians, in mg/L, are 9.4, 10.1, 9.8, 10.6, 10.0, 9.7, 10.4, 9.9, 10.3, and 9.8. Estimate their mean, MCSE, and a 95% \(t\) interval for the expected replicate-specific median.

The values sum to 100.0 mg/L, so the estimated mean is \(100.0/10=10.0\) mg/L. The sum of squared deviations from 10.0 is 1.16 \((\text{mg/L})^2\). Thus:

s_Q=\sqrt{\frac{1.16}{10-1}}=0.359\ \text{mg/L}, \qquad \widehat{SE}_{MC}=\frac{0.359}{\sqrt{10}}=0.114\ \text{mg/L}.

Using \(t_{0.975,9}=2.262\), the interval is \(10.0\pm2.262(0.114)\), which is 9.74 to 10.3 mg/L. This interval quantifies the simulation precision of the estimated mean replicate-specific median, subject to independent replicate datasets and the usual \(t\)-interval assumptions. Ten replicates give a fairly visible amount of Monte Carlo uncertainty; more independent replicates would generally narrow it.

Monte Carlo Precision for NPDE Summaries

An NPDE summary, such as the mean of the normalized prediction distribution errors, can also vary across simulation runs because the predictive distributions used to calculate the errors are estimated by simulation. To isolate this simulation component, repeat the calculation using \(B\) independent simulation batches under the same fitted model and observed design. Let \(A_b\) be the resulting summary from batch \(b\). The mean of the batch summaries estimates the target, and their sample standard deviation divided by \(\sqrt{B}\) estimates its MCSE.

This calculation needs a clear definition of a batch. If each batch produces a complete, independent estimate of the same NPDE summary, the batch estimates are the units used to calculate the MCSE. Do not use the standard deviation of individual NPDE values divided by the square root of the number of observations as a substitute: that quantity concerns variation across observations and does not, by itself, measure the error caused by finite simulation.

Worked Example: MCSE for the Mean NPDE

For a hypothetical dataset, eight independent simulation batches give the following estimates of the mean NPDE: 0.04, −0.03, 0.02, 0.01, −0.05, 0.03, 0.00, and −0.02. Estimate the average and its MCSE.

The batch estimates sum to zero, so their average is 0.000. The sum of squared deviations from the average is 0.0068. Therefore:

s_A=\sqrt{\frac{0.0068}{8-1}}=0.0312, \qquad \widehat{SE}_{MC}=\frac{0.0312}{\sqrt{8}}=0.0110.

A 95% \(t\) interval with \(7\) degrees of freedom is \(0.000\pm2.365(0.0110)\), or approximately −0.0261 to 0.0261. It describes the simulation precision of this mean-NPDE estimate under the batching procedure. It does not demonstrate that the model is correct, nor does it measure the full uncertainty in the NPDE summary from collecting a different study sample.

Choosing Replicate Counts and Checking Reliability

A useful simulation plan starts with a precision goal, such as an MCSE no larger than a specified amount in the units of the summary. For an estimated probability, the target itself gives a simple planning calculation. If \(M\) independent simulated outcomes are classified as either meeting or not meeting a condition, the estimated probability is a sample proportion \(\hat{p}\), with estimated MCSE \(\sqrt{\hat{p}(1-\hat{p})/M}\).

Worked Example: Replicate Count for a Tail Probability

In a hypothetical predictive check, 204 of 400 independent simulated datasets have a summary at or below a fixed observed value. Estimate the simulated probability and its MCSE, then estimate how many simulations would be needed for an MCSE near 0.01 if the probability remains about 0.51.

The estimated probability is \(\hat{p}=204/400=0.510\). Its MCSE is:

\widehat{SE}_{MC} = \sqrt{\frac{0.510(1-0.510)}{400}} = 0.0250.

For a planning value of \(p=0.51\), set \(\sqrt{p(1-p)/M}=0.01\) and solve for \(M\):

M=\frac{0.51(0.49)}{0.01^2}=2{,}499.

About 2,500 independent simulations would therefore be needed to target an MCSE near 0.01 under this approximation. This is a planning estimate, not a guarantee: the achieved MCSE depends on the probability actually being simulated and on the independence assumption.

Several checks help identify unreliable precision estimates:

Check independence at the right level. Repeated records within a simulated dataset are not automatically independent Monte Carlo units. Use independent replicate datasets or another variance calculation that accounts for clustering.
Check stability across batches. Large differences between independent batch estimates suggest that the total simulation count may be insufficient, the target may have high variability, or a batch may be constructed incorrectly.
Check quantile behaviour. Quantiles can be unstable when few values fall near the target quantile, when there are ties, or when the density estimate used in an asymptotic MCSE is poor. Compare independent replicate estimates or use a method that respects replicate-level dependence.
Check for rare events. A simulated probability near zero or one can be difficult to estimate precisely. If no simulated events occur, a plug-in MCSE of zero is misleading; it does not establish that the target probability is zero.
Separate numerical error from scientific conclusions. A small MCSE means repeated simulation would give similar estimates. It does not mean that the model adequately describes the observed study.

Independent replicate batches are a straightforward alternative to formulas that assume independent individual records. They can be computationally convenient and make the simulation unit explicit, but require enough independent batches to estimate their variability. Direct quantile formulas can be more efficient when their independent-draw and density assumptions are credible. Neither approach removes the need to define the target summary precisely.

As a practical rule, doubling the number of independent simulations reduces the MCSE of a mean or proportion by a factor of \(1/\sqrt{2}\), approximately 0.707, when the assumptions hold. Halving the MCSE requires about four times as many simulations. This diminishing return can guide resource decisions, but computational effort should be concentrated on summaries whose simulation uncertainty could affect interpretation.

Common Mistake

A common mistake is to regard the number of simulated records as the number of independent simulations. For example, 1,000 concentration records from a small number of simulated individuals do not necessarily provide the precision of 1,000 independent outcomes. If records from the same individual share simulated random effects, treating them as independent typically understates MCSE.

Key takeaway: Match the MCSE calculation to the simulation unit and the summary being estimated. Report the replicate count, the estimated summary, its MCSE, and—when useful—a simulation-error interval. Interpret these as precision conditional on the model and design, not as evidence that the model is correct.

Check Your Understanding

Work through each problem, then compare your reasoning with the solutions.

A simulation estimates a mean summary as 4.2 with a sample standard deviation of 1.5 across 100 independent replicates. What is its MCSE?
Why might a VPC quantile calculated from pooled repeated records have a different MCSE from one calculated under an independent-record assumption?
A probability estimate is 0.20 from 500 independent simulations. Calculate its approximate MCSE.
Does a very small MCSE show that a fitted population PK model is adequate? Explain.

Solutions. (1) The MCSE is \(1.5/\sqrt{100}=0.150\), in the units of the summary. (2) Repeated records within a dataset may be dependent, so the number of records can overstate the number of independent simulation units; a replicate-level or cluster-aware calculation is needed. (3) The MCSE is \(\sqrt{0.20(0.80)/500}=0.0179\). (4) No. A small MCSE indicates that the simulation target is estimated precisely under the specified model and design. Model adequacy is a separate question.