Why Nest the Simulation?
A linked drug-development analysis can carry uncertainty through several kinds of models: a PBPK model may inform population PK parameters; the population PK model may generate patient exposures; an exposure-response model may map exposure to response; and a trial simulation may predict a decision-relevant outcome. The final prediction depends not only on the equations in those models, but also on uncertain parameters, differences among patients, and the finite sample in a future trial.
Nested simulation represents these sources at the level where they arise. Parameter uncertainty is sampled at the outer level; patient characteristics and individual effects are sampled within each parameter draw; and trial outcomes are sampled from the simulated patients. The nesting is substantive, not just a programming convenience: it specifies which quantities are shared and which are independently resampled.
Write the uncertain parameter vector as \(\theta\), a patient-level random quantity as \(U_i\), and a simulated trial outcome as \(T\). A linked chain can be represented schematically as \(X_{\mathrm{PBPK}}=f_{\mathrm{PBPK}}(\theta)\), \(\phi=g(\theta,X_{\mathrm{PBPK}})\), \(C_i=h(\phi,U_i)\), and \(T=k(C_1,\ldots,C_n,\epsilon)\). Here \(\phi\) denotes parameters or covariate functions handed to the population PK or exposure-response component, \(C_i\) is the exposure for patient \(i\), and \(\epsilon\) collects trial-outcome randomness. The equations and interfaces should use harmonized quantities and populations, as established in the earlier material on drawing model interfaces and harmonizing quantities.
The parameter draw must remain coherent across the chain. If the same uncertain clearance parameter affects a PBPK prediction and a subsequent PK prediction, draw it once for that simulated world and propagate its value. When parameters are correlated, draw their joint vector from the specified joint distribution rather than drawing each marginal independently. Patient-level random effects are different: they are drawn for each patient conditional on the selected population parameters. Trial sampling is different again: it represents which patients and outcomes appear in a finite simulated trial.
Set Up the Conditional Structure
Before writing simulation code, state the conditioning structure and the target. For example, the target might be the probability that a future trial meets a prespecified criterion, averaged over parameter uncertainty and future trial outcomes. That is a different target from the success probability conditional on one fixed parameter vector. For the unconditional target, each outer parameter draw should ordinarily be paired with multiple independently simulated trials, each containing its own patient sample.
For outer draw \(m\), sample \(\theta^{(m)}\) from the distribution representing parameter uncertainty. Carry that draw through every component that depends on it.
Using \(\theta^{(m)}\), generate patient covariates and individual random effects from the specified target-population distributions. Calculate each patient's exposure and response through the linked models.
For each simulated trial, sample the required number of patients and trial outcomes, including the randomization, residual, or event processes needed by the design.
Record the trial-level endpoint, such as success or failure. Summarize across outer draws and use independent outer draws, not every nested patient record, to assess Monte Carlo precision.
The distributions and conditional independence assumptions must match the scientific model. Patients may be conditionally independent given shared parameters, but outcomes from patients in the same trial can also depend on shared trial-level factors. Such factors should be drawn once per trial, not independently for every patient. Likewise, a parameter representing uncertainty about the model is generally shared across all patients in a simulated trial; drawing a new value for each patient would describe a different model.
Worked Example: Carrying a Parameter Draw Through Exposure and Response
Consider a hypothetical, deliberately small linked model. A parameter \(\theta\) is an uncertain clearance multiplier, with outer draws of 0.8 and 1.2, each having probability 0.5. Exposure is \(C=10V/\theta\) exposure units, where \(V\) is a patient multiplier taking values 0.8 and 1.2 with equal probability. The exposure-response model gives response probability \(p(C)=C/(C+10)\). Calculate the mean response probability over the parameter and patient draws.
First hold the parameter draw fixed and calculate patient exposures and response probabilities. For \(\theta=0.8\), the two patient draws give \(C=10(0.8)/0.8=10\) and \(C=10(1.2)/0.8=15\) exposure units. Their response probabilities are \(10/(10+10)=0.50\) and \(15/(15+10)=0.60\), so the conditional mean is \((0.50+0.60)/2=0.55\).
For \(\theta=1.2\), the exposures are \(10(0.8)/1.2=6.667\) and \(10(1.2)/1.2=10\) exposure units. The response probabilities are \(6.667/(6.667+10)=0.400\) and \(10/(10+10)=0.500\), giving a conditional mean of \((0.400+0.500)/2=0.450\).
Averaging over the two equally likely parameter draws gives \((0.55+0.45)/2=0.500\). This is the expected response probability under this hypothetical model and these specified distributions. It is not the response probability evaluated at the average parameter: substituting \(\theta=1.0\) and averaging patient multipliers happens to give 0.500 here, but nonlinear models do not generally preserve that equality. In a trial simulation, these probabilities would then be used to generate individual outcomes and finite-sample trial statistics.
Why the Levels Affect Predictive Variation
A short derivation shows how parameter uncertainty and trial sampling contribute to the variation in a predicted trial result. Suppose a simulated trial has \(n\) patients and a binary response for each patient. Conditional on a shared parameter draw \(\theta\), assume each patient's response is an independent Bernoulli variable with the same probability \(p_\theta\). Let \(\hat p_{\mathrm{trial}}\) be the fraction responding. Then, conditional on \(\theta\), its mean is \(p_\theta\) and its variance is \(p_\theta(1-p_\theta)/n\).
The law of total variance states that the variance of a quantity is the variance of its conditional mean plus the expected value of its conditional variance. Applying it here gives:
The first term reflects variation across parameter draws in the conditional mean trial response. The second reflects finite-trial sampling conditional on those parameters. In a heterogeneous population, the conditional variance term must instead be calculated from the patient-level response distribution and any shared trial effects; the equal-probability, conditionally independent Bernoulli assumption above is not automatic. For a nonlinear decision endpoint, an analytic variance decomposition may not be convenient, but the nested simulation still generates the relevant distribution directly.
Worked Example: Separating Parameter and Trial Contributions
In a hypothetical trial, a shared parameter draw implies response probability \(p_\theta=0.45\) or \(0.55\), each with probability 0.5. A trial enrolls \(n=100\) patients, whose binary responses are conditionally independent given the parameter. Calculate the mean response fraction and its variance across parameter uncertainty and trial sampling.
The mean conditional response probability is \(E_\theta(p_\theta)=(0.45+0.55)/2=0.500\). The parameter contribution is \(\operatorname{Var}_\theta(p_\theta)=\{(0.45-0.50)^2+(0.55-0.50)^2\}/2=0.002500\).
For either parameter draw, the conditional sampling variance is \(p_\theta(1-p_\theta)/100=(0.45)(0.55)/100=0.002475\). It is the same at 0.55, so its expectation is also 0.002475. The total variance is \(0.002500+0.002475=0.004975\), and the standard deviation is \(\sqrt{0.004975}=0.07053\), or about 7.05 percentage points.
The calculation concerns the predictive spread of a trial response fraction under the stated assumptions. It does not say that every component of this spread is parameter uncertainty: approximately half of the variance here comes from the between-parameter term and half from trial sampling. Patient heterogeneity or other shared trial factors would require additional structure.
Estimate a Decision Probability and Its Monte Carlo Error
Suppose the target is the probability \(q\) that a future simulated trial meets a defined success criterion, integrating over parameter draws and trial outcomes. Draw \(M\) independent outer parameter vectors. For each vector, simulate \(R\) trial replicates and let \(\hat q_m\) be the fraction of those replicates that succeed. The nested estimate is \(\hat q=M^{-1}\sum_{m=1}^{M}\hat q_m\). The outer groups are independent if their parameter draws and simulation streams are independent; trial replicates within a group share the same outer draw and are therefore not generally independent after averaging over parameter uncertainty.
A practical estimate of the Monte Carlo standard error is the sample standard deviation of the \(M\) outer-group estimates divided by \(\sqrt{M}\):
This estimates simulation noise in \(\hat q\), not uncertainty about whether the model is scientifically correct. More outer draws improve the precision of an integrated estimate and better explore parameter uncertainty. More trial replicates per outer draw can improve the accuracy of each \(\hat q_m\). For a fixed computational budget, allocation between \(M\) and \(R\) depends on how much variation comes from outer draws versus within-draw trial randomness; compare estimates from increased allocations rather than assuming one fixed ratio is adequate.
Worked Example: Calculating a Cluster-Based Monte Carlo Standard Error
Suppose four independent outer parameter draws each have ten simulated trial replicates. The numbers of successful trials are 5, 6, 5, and 8, so the group success fractions are \(\hat q_m=0.50, 0.60, 0.50,\) and \(0.80\). This small calculation illustrates the method; four outer draws are far too few for a reliable final analysis.
The overall estimate is \(\hat q=(0.50+0.60+0.50+0.80)/4=0.600\). The squared deviations are \(0.0100, 0, 0.0100,\) and \(0.0400\). Thus \(s_{\hat q_m}=\sqrt{0.0600/3}=0.1414\), and the estimated Monte Carlo standard error is \(0.1414/\sqrt{4}=0.07071\).
Treating all 40 trial outcomes as independent Bernoulli draws would give \(\sqrt{(0.60)(0.40)/40}=0.07746\). That calculation ignores the shared parameter draw within each group, so it is not the appropriate general error calculation for this nested design. In this particular small example, its numerical value is larger; dependence does not guarantee that a naive calculation is smaller. The outer-group estimate is itself highly unstable with \(M=4\), which is why precision should be checked with substantially more independent outer draws.
Diagnostics, Alternatives, and Limitations
Nested simulation is only as credible as its conditional distributions, interfaces, and implementation. Check that simulated parameter vectors reproduce the intended marginal distributions and covariance; that each shared parameter is used consistently across linked components; and that patient covariates, random effects, and trial processes are drawn at the correct level. Compare selected simulated outputs with analytic expectations or simpler calculations when available. If the simulated exposure distribution or response rate is inconsistent with the component models, troubleshoot the handoff before interpreting the final trial result.
Assess Monte Carlo stability by repeating the calculation with more outer draws and, where useful, more inner replicates. Track the estimate and its Monte Carlo standard error, and inspect the quantities driving a decision, especially tail probabilities or rare successes. A stable overall mean does not guarantee a stable probability near a decision threshold. Use independent random streams across outer groups, and ensure that any stratification, shared trial effect, or patient sampling rule in the simulated design is represented explicitly.
| Approach | Useful when | Key limitation |
|---|---|---|
| Nested Monte Carlo | Several nonlinear linked models and multiple conditional sources of randomness must be represented. | Can be computationally demanding; convergence and implementation require checking. |
| Plug-in calculation | A quick reference prediction at one chosen parameter vector is useful. | Does not integrate parameter uncertainty or, by itself, finite-trial variability. |
| Analytic approximation | A low-dimensional model has tractable conditional expectations or variances. | May be inaccurate for nonlinear links, complex endpoints, or nonstandard sampling. |
| Bootstrap | Resampling observed data is appropriate for the inferential target and available data structure. | Does not automatically represent future patient variation or uncertainty in every linked model component. |
A common pitfall is to draw an uncertain population parameter independently for every patient, then describe the resulting spread as a future-trial prediction. That changes a shared model parameter into patient heterogeneity and removes the common parameter shift within each simulated trial. The correction is to draw the parameter vector once per outer draw, hold it fixed across the linked chain and trial conditional on that draw, and generate patient-level variability separately. Conversely, drawing one patient random effect and reusing it for an entire trial would incorrectly turn individual variation into a shared trial effect.
Nested simulation has limitations. It can be expensive when each outer draw requires a costly PBPK calculation or many trial replicates. Rare decision events can require very large simulation runs to estimate precisely. Monte Carlo precision also does not measure structural model error, unsupported extrapolation, or uncertainty omitted from the distributions. Those limitations call for scientific assessment of the models and their interfaces; a large simulation count cannot repair a misspecified chain.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Question: A clearance parameter is shared across patients in a simulated trial. At which level should it be drawn?
Solution: Draw it at the outer parameter level and hold the draw fixed across the linked models and patients in that trial. Draw patient-specific effects separately for each patient. - Question: With \(n=50\) and fixed conditional response probability \(p=0.40\), what is the variance of the trial response fraction under independent Bernoulli responses?
Solution: \(\operatorname{Var}(\hat p_{\mathrm{trial}}\mid p)=p(1-p)/n=(0.40)(0.60)/50=0.004800\). - Question: Two equally likely parameter draws imply trial success probabilities 0.30 and 0.50. What is the integrated success probability?
Solution: \(q=(0.30+0.50)/2=0.400\). This is the average over the specified parameter distribution, not the success probability conditional on either draw. - Question: Why is it generally inappropriate to calculate Monte Carlo error as if all trial replicates sharing one outer parameter draw were independent draws from the fully integrated distribution?
Solution: They share the same parameter draw, so their outcomes can be dependent marginally over parameter uncertainty. Estimate error using independent outer groups, for example with the sample standard deviation of group estimates divided by \(\sqrt{M}\). - Question: A trial success probability estimate changes little when inner trial replicates increase, but its outer-group estimates vary substantially. What should be increased first to improve precision?
Solution: Increase the number of independent outer parameter draws \(M\), then reassess the Monte Carlo standard error and stability. Additional inner replicates alone do not provide as many independent draws from parameter uncertainty.