Simulation Checks For Population PD Model Stability

Why Simulate Before Relying on a Population PD Model?

A candidate population PD model can fit observed data and still be difficult to estimate reliably under the intended study design. For example, sparse observations near the concentration range where a response changes may leave an \(EC_{50}\) weakly identified. A covariate effect may also be hard to distinguish from between-subject variability when the covariate has little variation or is unevenly represented. Simulation checks make these vulnerabilities visible by asking a direct question: if this model and design generated the data, would the planned estimation procedure recover the parameters?

A parameter-recovery simulation generates many hypothetical datasets from a specified population model, using the planned sampling design, and refits the candidate model to each dataset. The results describe the behavior of the full estimation process under those assumptions. They do not prove that the model is true, nor do they establish that its assumptions match the actual population.

Specify the Simulation Before Interpreting It

The simulation needs to represent the features of the planned study that affect estimation: numbers of subjects and observations, dose or concentration profiles, observation times, covariate distributions, missing data assumptions, between-subject variability, residual variability, and the response model. When concentration is treated as an input to the PD model, its simulated or supplied values should still reflect the planned concentration range and observation schedule. The analysis step should use the same model structure, parameter constraints, and estimation procedure intended for the study.

For illustration, suppose response \(Y_{ij}\) for subject \(i\) at observation \(j\) follows an Emax model with a covariate on the population Emax:

Here \(C_{ij}\) is concentration, \(z_i\) is a covariate centered and scaled so that one unit is one standard deviation, and \(\eta_i\) represents subject-level variation in Emax. A typical simulation might specify \(\eta_i\sim N(0,\omega^2)\) and \(\epsilon_{ij}\sim N(0,\sigma^2)\), with independent subject effects and residual errors. Those distributional and independence assumptions are part of the simulation, not facts established by the fitted model. Other error structures or random-effects assumptions require corresponding changes to the generator.

A useful design specification preserves the planned distribution of \(z_i\), including relevant subgroup sizes and correlations with other covariates. It also uses realistic concentration and time profiles. Simulating an evenly balanced covariate when the planned study has a small subgroup, for example, can make the covariate effect look more precisely estimable than it will be in practice.

Summarize Parameter Recovery

Let \(R\) be the number of simulated datasets for which the estimation procedure produced an estimate of a parameter \(\theta\), and write the estimate from replicate \(r\) as \(\hat{\theta}_r\). If all \(R\) estimates are available, empirical bias and relative bias are:

Relative bias is useful for nonzero parameters, but can be misleading or undefined when the true value is near zero. In that setting, report bias in the parameter’s units and consider a scale suited to the parameter. Empirical precision is commonly summarized by the sample standard deviation of the estimates across replicates:

This is the variation of estimates across simulated datasets, not the standard error reported by a single fit. Root mean squared error (RMSE) combines variability and bias. For the sample-variance convention above, it can be written:

The identity follows by writing \(\hat{\theta}_r-\theta=(\hat{\theta}_r-\bar{\hat{\theta}})+(\bar{\hat{\theta}}-\theta)\), squaring and summing across replicates. The cross-product sums to zero because \(\sum_r(\hat{\theta}_r-\bar{\hat{\theta}})=0\). Thus the average squared error is the centered sum of squares divided by \(R\), plus the squared bias. RMSE is therefore small only when both bias and variability are limited.

Convergence Is Part of the Result

A fitted estimate is informative only if the estimation procedure has reached a defensible solution. Define convergence criteria before running the simulation. They might include a successful optimizer termination, acceptable gradient or objective-function checks, finite estimates and standard errors, and estimates away from inappropriate parameter boundaries. A software convergence flag alone is not a guarantee that the fitted solution is reliable; inspect whether different starting values or optimizers lead to materially different results.

Report the convergence rate as the number of fits meeting the prespecified criteria divided by the total number of simulated datasets. Also tabulate failure types, boundary estimates, and other concerning outcomes. Bias and precision calculated only from converged fits describe the distribution conditional on convergence. If difficult datasets fail disproportionately, those summaries can give an overly reassuring picture of performance. Do not silently discard failures or treat them as zero-valued estimates.

Check Covariate-Effect Stability Under the Planned Design

For a covariate parameter such as \(\beta\), assess recovery on its model scale and in an interpretable clinical scale where appropriate. In the Emax example, increasing \(z_i\) by one standard deviation multiplies the typical Emax by \(\exp(\beta)\). Thus \(\beta=0.30\) corresponds to a multiplier of \(\exp(0.30)=1.35\), approximately a 35.0% increase. This interpretation depends on the stated log-linear covariate model.

Stability is not just whether the average estimate is close to the generating value. Also examine spread, sign changes, boundary estimates, convergence, and sensitivity to the covariate distribution and observation schedule. If covariates are correlated, simulate their planned joint distribution rather than varying each one in isolation. A targeted set of scenarios can reveal whether the effect remains estimable when subgroup sizes, missingness, or concentration coverage are less favorable.

Diagnostics, Alternatives, and Limitations

Review simulation output as a set of diagnostics, not a single pass-or-fail statistic. Plot estimates against replicate number and compare their distributions with the generating values. Look for skewness, multimodality, long tails, repeated boundary estimates, or patterns linked to particular covariate or concentration scenarios. Compare the intended estimator with its reported standard errors: consistently small empirical spread alongside much larger reported standard errors, or the reverse, deserves investigation. Examine convergence failures by scenario, since aggregate rates can hide a design feature that makes a particular subgroup or parameter difficult to estimate.

Simulation-based parameter recovery answers a different question from common alternatives. A nonparametric bootstrap resamples observed data units to assess variation relative to the observed sample; it does not automatically test recovery under a fully specified generating model and planned future design. A visual predictive check compares simulated and observed data summaries for model adequacy, but a good predictive comparison does not by itself show that individual parameters are precisely recoverable. Standard errors from a single fitted model are local, model-based uncertainty measures; simulation can assess their behavior across repeated datasets under chosen assumptions.

The principal limitation is that the conclusions depend on the generator. If it reproduces the candidate model exactly, reassuring recovery shows that the design and estimation procedure can recover parameters under that model; it does not rule out misspecification. Add plausible alternative scenarios when important assumptions are uncertain, and state which assumption each scenario changes. Simulation results also have Monte Carlo uncertainty: a small number of replicates can make estimated bias, precision, or convergence rates themselves unstable.