Why Diagnose Convergence?
A PBPK simulation workflow can propagate uncertain inputs through a model and produce a distribution of predicted exposure. But a finite simulation estimates features of that distribution with finite precision. For example, a sample mean concentration-time area under the curve (AUC) may change as more parameter sets are simulated, and a sample variance may remain unstable even after the mean looks steady. Convergence diagnostics help assess whether the simulated estimates are settling down, rather than changing substantially as additional draws are added.
Here, convergence means practical stability of the Monte Carlo estimates as the number of simulations increases. It does not mean that a finite estimate is exactly equal to its population target, nor does it establish that the PBPK model or the input distributions are scientifically correct. The introductory material on reproducible Monte Carlo workflows addresses recording the simulation setup; this tutorial focuses on diagnosing the numerical stability of estimates from that setup.
The independent, identically distributed assumption describes ordinary Monte Carlo sampling: each input set is drawn from the specified joint uncertainty distribution, independently of the other sets, and the same deterministic model is applied to each. The target moments must also exist. A finite variance is needed for the usual standard error of the sample mean; diagnosing the sample variance’s sampling variability generally requires a finite fourth central moment. If the design instead uses stratification, Latin hypercube sampling, or nested draws, observations may not be independent in the ordinary sense. Use diagnostics that respect the design, such as independent replicated designs, rather than treating all rows as independent replicates.
Trace Plots and Running Summaries
A trace plot places the simulated output \(Y_i\) on the vertical axis and simulation index \(i\) on the horizontal axis. For independent Monte Carlo draws, the trace should look like an unstructured sequence around a broadly consistent range. Look for long stretches at a different level, gradual drift, abrupt changes, or isolated extreme outputs. A trace is particularly useful for finding features that a single overall mean or variance hides.
Running summaries show what happens to the estimates as simulations accumulate. At each \(k\), calculate \(\bar{Y}_k\) and \(s_{Y,k}^2\) using only the first \(k\) outputs. A stable running mean will generally fluctuate less as \(k\) grows, but it need not move monotonically. A running variance can be more erratic because squared deviations give extreme outputs considerable influence. Assess whether these summaries settle into a range narrow enough for the intended reporting precision; a visually flat segment alone is not proof of convergence.
Worked Example: Reading a Short Running Summary
Suppose eight hypothetical Monte Carlo simulations produce AUC values, in mg·h/L, of 8, 12, 10, 9, 11, 10, 8, and 12. The running mean after four simulations is
After six simulations, the cumulative total is 60 mg·h/L, so \(\bar{Y}_6=60/6=10.0\) mg·h/L. After all eight, the total is 80 mg·h/L and \(\bar{Y}_8=80/8=10.0\) mg·h/L. The running mean has returned to 10.0, but this tiny trace cannot establish convergence: the sample is short, and a later large output could still change the estimate.
For all eight outputs, the squared deviations from 10 are \(4,4,0,1,1,0,4,4\), which sum to 18. Thus,
This calculation gives the sample variance, not evidence by itself that the variance estimate is stable. The trace and running variance should be inspected over a much longer simulation.
Repeated Batches and Sampling Variability
A complementary check is to divide a simulation into batches and compare the estimates across them. If batches are independently generated, substantially different batch means or variances can indicate that the overall simulation is too short, or that outputs are highly variable. Batches should be large enough that each estimate is informative; many tiny batches can create noisy comparisons. Batch agreement is reassuring but not conclusive, because all batches can miss the same rare part of the output distribution.
For the mean, the usual estimated standard error under independent, identically distributed sampling is \(SE(\bar{Y})=s_Y/\sqrt{n}\). It describes the sampling variability of the estimated mean across repeated simulations of size \(n\), under the stated assumptions. Comparing changes in the running mean with this scale can help judge whether movement is consequential. It is a diagnostic aid, not a rule that a particular ratio proves convergence.
Worked Example: Comparing Four Independent Batches
Consider four hypothetical, independent batches of 1,000 PBPK simulations each. Their AUC means are 10.1, 9.8, 10.2, and 9.9 mg·h/L. The average of the batch means is
The deviations from 10.0 are \(0.1,-0.2,0.2,-0.1\) mg·h/L. The sample standard deviation across the four batch means is
The standard error of their average is \(0.1826/\sqrt{4}=0.0913\) mg·h/L. The individual batch means span 9.8 to 10.2 mg·h/L, which is reasonably consistent for a report rounded to about 1 mg·h/L, but would not support claiming stability to 0.01 mg·h/L. For comparison, if the within-batch sample standard deviation is about 3.0 mg·h/L, its estimated standard error for a batch mean is \(3.0/\sqrt{1000}=0.0949\) mg·h/L. This is similar to the between-batch standard error here, as expected for a plausible amount of sampling fluctuation.
Suppose the four batch variance estimates are 8.8, 9.4, 9.1, and 8.7 \((\text{mg·h/L})^2\). Their average is \(36.0/4=9.0\ (\text{mg·h/L})^2\), and their range is 8.7 to 9.4 \((\text{mg·h/L})^2\). This is a useful stability check, but four batches provide limited evidence. Also, the average of batch variances is not generally identical to the variance computed by pooling all observations: the pooled variance also reflects differences among batch means.
Why Variances Can Converge More Slowly
The sample variance depends on squared deviations, so rare high or low outputs can affect it more strongly than they affect the mean. This is one reason to inspect both trace plots and running summaries. A short simulation may have a stable-looking mean but a running variance that changes sharply when an extreme output appears.
A brief moment calculation explains the additional conditions involved. Write \(\mu_4=E[(Y-\mu_Y)^4]\), the fourth central moment. For independent draws with finite fourth moment, the sample variance is asymptotically based on the average of squared deviations. Its leading sampling variance is \((\mu_4-\sigma_Y^4)/n\), so its approximate standard error is
To see the source of this expression, use the identity \(s_Y^2=\frac{n}{n-1}(\overline{Y^2}-\bar{Y}^2)\), where \(\overline{Y^2}=\frac{1}{n}\sum_{i=1}^n Y_i^2\). For large \(n\), the factor \(n/(n-1)\) is close to 1. A first-order expansion around the population moments gives
The terms in this average have mean zero and variance \(E[(Y-\mu_Y)^4]-\sigma_Y^4=\mu_4-\sigma_Y^4\). Independence therefore gives variance approximately \((\mu_4-\sigma_Y^4)/n\), yielding the stated standard error. This justification requires a finite fourth moment; strongly heavy-tailed outputs can make variance estimates very unstable or make this approximation unreliable. In practice, compare running and batch variances and interpret any moment-based approximation cautiously.
Worked Example: An Extreme Output Changes Both Summaries
Imagine a hypothetical trace of eight AUC values, all 2 mg·h/L, followed by a ninth value of 20 mg·h/L. Before the ninth simulation, the mean is 2 mg·h/L and the sample variance is 0. After it is added, the mean becomes
The deviations from 4 are eight values of \(-2\) and one value of 16. The squared deviations sum to \(8(4)+16^2=288\), so
This simple sequence illustrates how an extreme draw can change the estimated mean and especially the variance. It does not establish whether 20 mg·h/L is an error or a legitimate tail value. Check the input draw, model output, and trace; do not delete an extreme output merely because it disrupts apparent stability. If it is a valid draw under the specified model and input distributions, it is part of the target distribution.
Practical Workflow, Alternatives, and Limitations
Specify the PBPK output, its units, the input distribution, and whether the target is its mean, variance, or both. A diagnosis is meaningful only relative to a defined target and reporting precision.
Plot outputs against simulation order. Investigate trends, long runs, abrupt shifts, extreme values, and model failures rather than relying on a single aggregate summary.
Plot \(\bar{Y}_k\) and \(s_{Y,k}^2\) against \(k\). Check whether their remaining movement is acceptably small for the intended use, not just whether the final segment looks flat.
Calculate each batch’s mean and variance and assess their spread. For structured sampling, use independent replicated designs or another comparison appropriate to that design.
If estimates remain unstable, increase simulation effort, then repeat the checks. Record the number of simulations and the resulting estimates and diagnostics.
These checks have different strengths. A trace can reveal when an unusual draw occurred; running summaries show how its influence accumulates; repeated batches show whether results vary across independent runs. The standard error of the mean provides a scale for expected mean-estimation fluctuation under its assumptions. No single check replaces the others. A trace that looks patternless can still miss an event that is extremely rare, and agreement among a few batches can be coincidental.
Do not interpret convergence as validation of the scientific model, the chosen input distributions, or the numerical solution of the PBPK equations. Those are distinct questions. Likewise, a stable mean and variance do not guarantee that every feature of a prediction distribution is well estimated. For this tutorial, the focus is on first and second moments; stability of other summaries requires diagnostics suited to those summaries.
Common Mistakes
- Calling a flat running mean proof of convergence. A finite trace can look flat before a rare influential output appears. Extend the simulation and compare independent batches.
- Judging the mean and ignoring the variance. Variances can be more sensitive to extreme outputs. Inspect both running summaries and batch estimates.
- Using independent-draw formulas for a dependent design without adjustment. Stratified or nested sampling may not produce independent rows. Match the diagnostic to the sampling design.
- Removing a large output to make estimates look stable. First determine whether it is a simulation error or a valid draw. Removing a valid draw changes the target being estimated.
- Equating Monte Carlo stability with model correctness. Stable estimates can still come from an unsuitable PBPK model or input distribution.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- A simulation’s running mean changes from 5.0 to 5.1 mg·h/L over its last 10,000 draws, while its running variance continues to rise. What should be checked before concluding the simulation has converged?
- Four independent batches of 500 simulations have means 7.8, 8.0, 8.1, and 8.1 mg·h/L. Calculate their average mean.
- Why can an extreme output affect a sample variance more than a sample mean?
- Can agreement across several batches establish that the PBPK model is scientifically valid? Explain.
Solutions. (1) Inspect the trace for extreme or changing outputs, check how the running variance responds, compare batch variances, and consider whether more simulations are needed. A nearly stable mean alone does not establish variance stability. (2) The average is \((7.8+8.0+8.1+8.1)/4=32.0/4=8.0\) mg·h/L. (3) The variance uses squared deviations from the mean, so a large deviation contributes a disproportionately large squared term; its influence can also shift the mean. (4) No. Batch agreement supports stability of the estimated summaries under the simulation setup, but does not validate the PBPK model or its input assumptions.