Choosing Replicate Counts for NLME Simulation

Why Replicate Counts Depend on the Target

A simulation from a fitted nonlinear mixed-effects model can involve several layers of replication. One may draw population parameters, simulate new subjects conditional on each parameter draw, and generate multiple records for each subject. These counts are not interchangeable. More records from the same subjects may improve the detail of their trajectories, but they do not generally provide the same information as simulating more subjects. Likewise, many subjects simulated under only a few parameter draws may represent subject variability well while inadequately exploring parameter uncertainty.

The practical goal is not to find one universally sufficient number of replicates. It is to choose counts that make Monte Carlo error small enough for the summaries and decisions of interest. The introductory material on Monte Carlo error in simulated summaries establishes the central idea: a simulation estimate itself varies across repeated simulation runs. Here, the task is to locate which replication layer contributes most to that variation and allocate computation accordingly.

Counts should be selected for a specified target summary, such as a population mean concentration at a fixed time, a probability of exceeding a threshold, a quantile of an exposure measure, or an uncertainty-band endpoint. A simulation can be stable for one target and too noisy for another. For example, a mean may stabilize well before an upper-tail quantile.

Assumptions and a Variance Decomposition

Consider a target that can be calculated as a scalar \(Z_{ms}\) for subject \(s\) simulated under parameter draw \(m\). The scalar might be one subject's AUC, an indicator that the subject exceeds a threshold, or a concentration at a specified time. Write \(\theta_m\) for parameter draw \(m\), and let \(\mu(\theta_m)=E(Z_{ms}\mid\theta_m)\) be the target's conditional mean under that draw. The nested estimate of the overall mean is

For the variance calculation below, assume that parameter draws are independent draws from the intended parameter-uncertainty distribution; subjects are independently simulated conditional on their parameter draw; and the same \(S\) is used for each draw. Also assume the target has finite variance. Define \(\tau^2=\operatorname{Var}_{\theta}\{\mu(\theta)\}\), the variability in conditional target means across parameter draws, and \(\sigma^2=E_{\theta}\{\operatorname{Var}(Z\mid\theta)\}\), the average conditional subject-level variance.

This formula is an allocation guide, not a claim that the model's uncertainty is fully described by these two components in every application. In particular, it applies to a mean of subject-level summaries. If the target is a quantile or another nonlinear summary, use convergence checks for that target rather than treating the mean formula as its exact standard error.

Why the Nested Variance Formula Holds

Let \(\bar{Z}_m=S^{-1}\sum_s Z_{ms}\) be the subject average under draw \(m\). Conditional on \(\theta_m\), its mean is \(\mu(\theta_m)\) and its variance is \(\operatorname{Var}(Z\mid\theta_m)/S\), by conditional independence of the simulated subjects. The law of total variance gives

The \(M\) group averages are independent and equally weighted. The variance of their average is therefore their common variance divided by \(M\):

The decomposition explains an important allocation principle. Increasing \(S\) reduces the conditional subject-simulation term, but does not remove the parameter-draw term. If \(\tau^2/M\) dominates, simulating many more subjects under the same small set of parameter draws may yield little improvement in the overall estimate.

Worked Example: Selecting Counts for a Mean Exposure

Choosing Counts With a Precision Goal

A useful plan starts by stating an acceptable Monte Carlo error in units meaningful for the target. For a mean, this may be a maximum standard error or an approximate interval half-width. For a probability, it may be an absolute error in percentage points. For a quantile, it may be a tolerance in the outcome's units. The tolerance should reflect the intended use: a small numerical difference may be immaterial for a broad descriptive summary but consequential near a decision threshold.

Pilot simulations can estimate variance components for a mean. If a maximum standard error \(e\) is required, use the planning inequality

Here \(\hat{\tau}^2\) and \(\hat{\sigma}^2\) are pilot estimates, so the inequality is a planning approximation. If the parameter component dominates, consider increasing \(M\). If the subject component dominates, increasing \(S\) may be effective. If neither dominates, compare candidate allocations against computational cost. The minimum count that satisfies a variance estimate is not automatically adequate: check the resulting estimate empirically and allow for uncertainty in pilot variance estimates.

Observations Within Subjects Are Not Independent Subjects

The count \(R\) requires special care in longitudinal simulation. Records from one subject share that subject's random effects and often have correlated outcomes. Counting each row as an independent replicate exaggerates the amount of information. If the target is a subject-level summary calculated from the complete trajectory, the relevant independent units for subject variation are generally the simulated subjects, not the individual records used to construct each summary.

A simple special case makes the distinction clear. Suppose, conditional on a parameter draw, each record has a subject-specific random intercept with variance \(\omega^2\), plus independent record error with variance \(\epsilon^2\). For \(S\) independent subjects with \(R\) records each, the variance of the grand record mean, conditional on the parameter draw, is

More records reduce the independent record-error contribution, but they do not reduce the shared subject-intercept contribution. Real NLME trajectories may have more complex covariance, so this expression is illustrative rather than universal. The principle remains: within-subject dependence must be respected when deciding what counts as replication.

Practical Convergence Checks and Failure Modes

A planning calculation should be followed by a convergence check at the counts actually used. Examine estimates at increasing counts, such as successive cumulative checkpoints, and compare them with an explicit tolerance. For nested simulations, keep track of changes as \(M\) increases and as \(S\) increases; changing both at once makes it harder to identify the limiting source of noise. If computation permits, estimates from separate batches provide another view of run-to-run variability. Use the same target definition and simulated design at each checkpoint.

For each important summary, record its estimate, units, count allocation, and a measure of Monte Carlo error where one is available. Inspect more than the overall mean if the analysis will use quantiles, probabilities, or curves. Quantiles and interval endpoints can be unstable even when means look smooth, especially in tails. For a time-course summary, check representative times and relevant curve features rather than relying on one average over all records.

• Parameter estimates keep moving while the subject count is large. This suggests too few parameter draws or substantial variation across parameter draws; increasing subjects under the existing draws cannot remove that component.
• Substantial variation remains among subject averages under a draw. More subjects per parameter draw may be needed for the target mean.
• A curve changes when the time grid is refined. This indicates a resolution issue, not necessarily a shortage of simulated subjects.
• A smooth mean is treated as proof that every target has converged. Tail probabilities and quantiles require their own checks.
• Record count is reported as if every record were independent. Shared random effects and within-subject correlation make that interpretation generally invalid.

The main alternative to formula-based allocation is a purely empirical convergence strategy: run successively larger simulations and stop when target estimates change by less than a prespecified tolerance. This is easy to communicate, but small changes at one checkpoint can occur by chance and do not prove stability. Variance-based planning is more informative about where error comes from, but depends on pilot estimates and assumptions. Using both is usually more defensible than relying on either alone.

No finite replicate count corrects a misspecified model, an inappropriate target, or an incorrect simulation design. Nor does a stable estimate demonstrate that all aspects of a simulated distribution are well represented. Counts should be justified for the summaries that matter, and their limitations should be reported plainly.