Simulation

`simulate()`

Generate simulated observations from a model with random effects and residual error.

pub fn simulate(
    model: &CompiledModel,
    population: &Population,
    params: &ModelParameters,
    n_sim: usize,
) -> Result<Vec<SimulationResult>, String>

Parameters: - model: Compiled model - population: Template population (dose events and observation times are used; DV values are ignored) - params: True parameter values for simulation - n_sim: Number of simulation replicates

Returns: Ok with a vector of SimulationResult, one per observation per subject per replicate, or Err when the model and data fail a precondition — the same message fit() gives for that input (see Entry-point errors).

Example:

use ferx_core::api::read_population_for;

let model = parse_model_file(Path::new("model.ferx"))?;
let (population, _) = read_population_for(&model, &None, "data.csv", None, None, None, &[])?;

// Simulate 1000 replicates
let sims = simulate(&model, &population, &model.default_params, 1000)?;

for sim in &sims[..5] {
    println!("Sim {}, ID {}, TIME {}, IPRED {:.3}, DV {:.3}",
             sim.sim, sim.id, sim.time, sim.ipred, sim.dv_sim);
}

`read_population_for_simulation()`

Read a NONMEM-format CSV that will be simulated from rather than fitted.

pub fn read_population_for_simulation(
    model: &CompiledModel,
    covariate_decls: &Option<Vec<CovariateDecl>>,
    data_path: &str,
    fallback_columns: Option<&[&str]>,
    iov_column: Option<&str>,
    filter: Option<&SelectionFilter>,
    column_map: &[(String, String)],
) -> Result<(Population, Option<CovariateTable>), String>

Identical to read_population_for() — the reader every fitting entry point uses — except in how it reads a missing DV on a scored row (EVID=0, MDV=0). For a fit, a missing DV means “no observation here”, so the row is skipped and counted for the W_MISSING_DV warning. For a simulation the DV is the output, so the same row is a design point: a sampling time whose value has not been generated yet.

That makes the natural design template simulate as written:

ID,TIME,DV,EVID,AMT,CMT,MDV
1,0,.,1,100,1,1
1,0.25,.,0,.,1,0
1,1,.,0,.,1,0
1,4,.,0,.,1,0

let (population, _covariates) =
    read_population_for_simulation(&model, &None, "design.csv", None, None, None, &[])?;
let sims = simulate(&model, &population, &model.default_params, 1000)?;

Read with read_population_for() the same file yields zero simulated rows. Rows marked MDV=1 are excluded on both paths. Retained rows are reported once as a W_DESIGN_DV warning on the returned population, so a simulation run off an observed dataset makes visible that it kept rows a fit of the same data would have skipped.

The retained rows carry a placeholder DV (NaN, or the endpoint’s first declared state code for an integer-coded endpoint) that the simulated value replaces, so the returned population is for simulation only — do not pass it to fit(). That is enforced, not merely advised: a population carrying non-finite observations is rejected by fit() (and reported by ferx check) as E_NONFINITE_DV, and propensity-score-matched simulation (SimulateOptions::match_method) rejects it too, since a design template has nothing to compute posthoc etas from.

`simulate_with_seed()`

Same as simulate() but with a fixed random seed for reproducibility.

pub fn simulate_with_seed(
    model: &CompiledModel,
    population: &Population,
    params: &ModelParameters,
    n_sim: usize,
    seed: u64,
) -> Result<Vec<SimulationResult>, String>

`simulate_with_options()` — propensity-score matching

Same outputs as simulate(), with options for a seed and for propensity-score matching of the simulated random effects.

pub fn simulate_with_options(
    model: &CompiledModel,
    population: &Population,
    params: &ModelParameters,
    n_sim: usize,
    opts: &SimulateOptions,
) -> Result<Vec<SimulationResult>, String>

pub struct SimulateOptions {
    pub seed: Option<u64>,                  // None draws from entropy
    pub match_method: Option<MatchMethod>,  // None disables matching
    pub horizon: Option<f64>,               // TTE administrative censoring time; None = per-record window
}

pub enum MatchMethod {
    Optimal,   // global minimum total distance (linear assignment); recommended default
    Nearest,   // greedy nearest-neighbour in subject order
    Rank,      // pair by Mahalanobis-norm rank (k-th order statistic to k-th)
}

With match_method = None and horizon = None this is identical to simulate_with_seed() (or simulate() when seed is None).

horizon = Some(t) sets an administrative censoring time for time-to-event endpoints, decoupled from the observed event times: it overrides each TTE record’s per-record observation window so a re-simulated event-bearing subject (a competing-risks VPC) censors at the planned study end t rather than drawing unbounded. It has no effect on Gaussian endpoints. See TTE simulation.

Why match?

In real-world data, therapy is often adapted in response to a patient’s own PK — e.g. high-clearance patients are given longer dosing intervals or more frequent sampling. A standard VPC draws each subject’s eta independently of its dosing/sampling design, so this design↔︎eta association is lost and the VPC can show spurious model misspecification even when the model is correct.

Propensity-score matching restores the association. For each replicate:

Draw a pool of N etas from \(N(0, \Omega)\) (N = number of subjects).
Match the N drawn etas 1:1 (without replacement) to the N subjects’ fitted (posthoc) etas under the Mahalanobis distance over the model \(\Omega\): \(d^2(a,b) = (a-b)^\top \Omega^{-1} (a-b)\) using the chosen MatchMethod:
- Optimal — global minimum of the total matched distance (the linear assignment problem); mirrors MatchIt(method = "optimal", distance = "mahalanobis"). Best on average in simulation, the recommended default.
- Nearest — greedy nearest-neighbour in subject order; mirrors MatchIt(method = "nearest", distance = "mahalanobis").
- Rank — pair subjects and draws by the rank of their Mahalanobis norm \(\lVert \eta \rVert = \sqrt{\eta^\top \Omega^{-1} \eta}\) (k-th order statistic to k-th).
Each subject keeps its own observed design but is simulated with the drawn eta that matched its fitted eta — so a high-clearance draw lands on a subject whose adaptive design reflects high clearance.

The fitted etas are computed once via a posthoc (MAP) inner-loop pass over the observed population at params; only the drawn pool and the matching change across replicates.

Requires observed data: every subject must carry observations (so its posthoc eta is defined). The function returns Err if the population is empty or any subject has no observations. The matching is the only correction applied; binning and plotting the VPC is left to the caller (e.g. the vpc R package).

Validation against R

The matching numerics were cross-checked against R on a fixed 8-subject, 2-eta scenario:

the Mahalanobis cost matrix matches R’s mahalanobis() under Ω to 2e-11;
ferx-core’s optimal assignment is identical to two independent optimal solvers — clue::solve_LSAP and optmatch::pairmatch (the engine MatchIt(method = "optimal") uses) — with the same total cost to 10 dp.

A direct MatchIt(method = "optimal", distance = "mahalanobis") call (as in the PAGE-poster workflow) returns a slightly different pairing, because MatchIt scales the Mahalanobis metric by the empirical covariance of the combined etas rather than the model Ω. The optimal-matching algorithm is the same (optmatch); only the metric differs, which is the documented modelling choice here (Ω is exact for draws from N(0, Ω) and needs no extra estimation).

Example:

let pop = read_nonmem_csv(Path::new("rwd.csv"), None)?;   // observed data
let fit = fit(&model, &pop, &model.default_params, &opts)?;

// Recover the fitted ModelParameters (theta/Omega/sigma) from the result.
let params = fitted_params_from_result(&fit, &model);
let sims = simulate_with_options(
    &model, &pop, &params, 200,
    &SimulateOptions { seed: Some(42), match_method: Some(MatchMethod::Optimal) },
)?;
// → 200 replicates; build the pmVPC from `sims` as usual.

`predict()`

Population predictions without random effects (eta = 0). No simulation noise is added.

pub fn predict(
    model: &CompiledModel,
    population: &Population,
    params: &ModelParameters,
) -> Result<Vec<PredictionResult>, String>

Returns: Ok with a vector of PredictionResult holding population-level predictions, or Err on a failed precondition, as for simulate().

Example:

let preds = predict(&model, &population, &model.default_params)?;
for p in &preds {
    println!("ID {}, TIME {}, PRED {:.3}", p.id, p.time, p.pred);
}

`simulate_with_uncertainty()`

Simulate observations while propagating population parameter uncertainty. For each parameter draw from the uncertainty distribution, etas (from the drawn Omega) and epsilons (from the drawn Sigma) are sampled — giving simulation bands that include both individual variability and uncertainty in the population estimates.

pub fn simulate_with_uncertainty(
    model: &CompiledModel,
    population: &Population,
    fit_result: &FitResult,
    opts: &SimulateUncertaintyOptions,
) -> Result<Vec<SimulationResult>, String>

pub struct SimulateUncertaintyOptions {
    pub n_uncertainty_draws: usize,      // parameter sets drawn
    pub n_sim_per_draw: usize,           // eta/eps replicates per draw
    pub method: UncertaintyMethod,       // Asymptotic | Sir
    pub seed: Option<u64>,
}

pub enum UncertaintyMethod {
    Asymptotic,   // MVN around ML estimate using FitResult.covariance_matrix
    Sir,          // resamples from FitResult.sir_resamples_packed
}

Parameters: - fit_result: Result of a previous fit() call. Must include covariance_matrix (Asymptotic) or sir_resamples_packed (SIR). - opts.n_uncertainty_draws: number of parameter sets drawn. - opts.n_sim_per_draw: number of subject-level replicates per draw. - opts.method: how to draw the parameter sets.

Returns: Vec<SimulationResult> of length n_uncertainty_draws * n_sim_per_draw * n_subjects * n_obs. Each row carries a draw index (1..=n_uncertainty_draws) and a sim index (1..=n_sim_per_draw).

Prerequisites: - For UncertaintyMethod::Asymptotic: run fit() with run_covariance_step = true (or covariance = true in [fit_options]). - For UncertaintyMethod::Sir: run with sir = true and sir_keep_samples = true so the resampled vectors are kept on the FitResult.

Example:

let fit_opts = FitOptions {
    run_covariance_step: true,
    ..FitOptions::default()
};
let fit = fit(&model, &population, &model.default_params, &fit_opts)?;

let sims = simulate_with_uncertainty(
    &model,
    &population,
    &fit,
    &SimulateUncertaintyOptions {
        n_uncertainty_draws: 200,
        n_sim_per_draw: 10,
        method: UncertaintyMethod::Asymptotic,
        seed: Some(42),
    },
)?;

Result Types

pub struct SimulationResult {
    pub draw: usize,    // Uncertainty draw index (1-indexed; always 1 for simulate())
    pub sim: usize,     // Replicate number within a draw (1-indexed)
    pub id: String,     // Subject ID
    pub time: f64,      // Observation time
    pub ipred: f64,     // Individual prediction (no residual error)
    pub dv_sim: f64,    // Simulated observation (with residual error)
}

pub struct PredictionResult {
    pub id: String,
    pub time: f64,
    pub pred: f64,      // Population prediction (eta = 0)
}

Simulation Process

For each replicate and each subject:

Sample random effects: \(\eta_i \sim N(0, \Omega)\) using the Cholesky factor \(L\): \(\eta = L \cdot z\), where \(z \sim N(0, I)\)
Generate predictions through the same time-varying-covariate-aware dispatcher used by estimation and predict() (compute_predictions_with_tv). PK parameters are recomputed per event from each row’s covariate snapshot, so a covariate that changes over a subject’s record (e.g. body weight) drives the prediction at each time — not just the baseline value. Subjects with no time-varying covariates stay on the single-evaluation fast path (pk_param_fn(theta, eta, covariates)), so their output is unchanged.
Add residual error: \(DV = IPRED + \sqrt{V} \cdot \epsilon\), where \(\epsilon \sim N(0, 1)\) and \(V\) is the residual variance from the error model. FREM covariate pseudo-observations (FREMTYPE > 0) instead use the additive covariate error \(V = \sigma_\text{EPSCOV}^2\), matching the FREM likelihood.

The same per-event covariate snapshots drive the NPDE/NPD diagnostics and the NCA-based start-value sweep, so all simulation-backed surfaces agree.

Simulation with Uncertainty Process

simulate_with_uncertainty() wraps the steps above in an outer loop over parameter draws:

Outer loop — for each of n_uncertainty_draws, draw a population parameter set \((\theta_k, \Omega_k, \Sigma_k)\) from the uncertainty distribution.
- Asymptotic: \(x_k = \hat{x} + L_\text{cov} z\) in packed log-space (theta, Cholesky-omega, sigma share one packed vector), then unpack — theta, Omega, and Sigma are perturbed coherently from a single MVN draw.
- SIR: pick a parameter vector at random from FitResult.sir_resamples_packed (the resampled pool retained from the SIR step). Draws that fall outside the parameter bounds or yield non-positive theta/omega/sigma are rejected and resampled (up to 10 × n_uncertainty_draws attempts).
Inner loop — for the drawn \((\theta_k, \Omega_k, \Sigma_k)\), run n_sim_per_draw replicates of the standard simulation process above. Etas are drawn from \(N(0, \Omega_k)\) and epsilons from the drawn \(\Sigma_k\).

The resulting SimulationResult rows are tagged with both draw and sim so downstream code can compute either marginal bands (over all draws and sims) or hierarchical bands (e.g. median across sims within each draw, then percentiles across draws).

What the uncertainty distribution is

There is no inverse-Wishart draw here. Ω uncertainty is not sampled from a conjugate matrix distribution; the inverse-Wishart in ferx belongs to the Bayesian MCMC estimator, which is a separate code path (estimation/bayes.rs) and is not used by simulate_with_uncertainty().

What actually happens, for UncertaintyMethod::Asymptotic:

The proposal is a plain multivariate normal, and it lives in the packed parameter space — theta (log), the Omega Cholesky factor (log-diagonal), and sigma (log) in one vector. FitResult.covariance_matrix is the inverse covariance Hessian of the OFV in that same space, so it is the natural — and the only dimensionally consistent — place to draw. Its off-diagonal blocks are what makes each draw’s θ, Ω and Σ move together instead of independently.
Because the Omega segment is \(\log L_{ii}\), the implied marginal on an Omega diagonal is roughly log-normal on the standard-deviation scale, and every draw yields a positive-definite Ω by construction — no rejection of non-PD matrices needed. Off-diagonal Cholesky entries are drawn on the identity scale.
A probability-scale theta is drawn on the logit scale. A theta used as inv_logit(logit(THETA) + ETA) (theta_transform = LogitProbability, e.g. a bioavailability declared on (0, 1)) is packed as \(\log\theta\) like any other non-negative theta, and a normal draw there has no ceiling at 1. So that coordinate is moved to \(y = \operatorname{logit}\hat\theta\) first, its row and column of the covariance scaled by \(dy/dx = 1/(1-\hat\theta)\) (or \(1/(\hat\theta(1-\hat\theta))\) when the theta is identity-packed because its lower bound is negative), and each draw is mapped back with \(\theta = \operatorname{inv\_logit}(y)\). The draws are logit-normal with the delta-method SD, and never reach 1 whatever the declared upper bound (#1548). A FIX’d theta is pinned as usual. A free one whose estimate is not strictly inside (0, 1) (possible only with an upper bound above 1) has no logit, so the sampler, and SIR, returns an error asking for an upper bound below 1 rather than drawing values that logit would clamp. A logit draw above 34.5, the clamp logit itself applies, is rejected and redrawn like any other out-of-bounds draw: from about 36.7 on, \(\operatorname{inv\_logit}\) rounds to exactly 1.
If covariance_matrix is not positive-definite, it is symmetrised and eigenvalue-floored before the Cholesky (regularised_cholesky). A draw from the regularised factor is slightly wider than the raw covariance in the corrected directions.
FIX’d parameters carry no uncertainty: any packed index flagged FIX is reset to its fitted value on every draw, so a FIX’d theta, omega, sigma or kappa is identical across all draws. (compute_bounds() pins fixed indices with lower == upper, so an unpinned continuous draw would otherwise be rejected for every model with a FIX.)
A draw that lands outside the packed bounds, or that unpacks to a non-positive theta / sigma / Omega diagonal, is discarded and redrawn, up to 10 × n_uncertainty_draws attempts in total; exhausting that budget is an Err, and usually means an ill-conditioned covariance matrix or an estimate sitting on a bound.
If a valid draw then fails during simulation (e.g. the drawn parameters make the ODE solver fail), that single draw is skipped with a uncertainty draw k skipped — … warning rather than aborting the whole run, so the returned set may contain fewer than n_uncertainty_draws distinct draws.

UncertaintyMethod::Sir makes none of these distributional assumptions: it resamples whole parameter vectors from the SIR pool, so its shape is whatever the SIR step found. Prefer it when the asymptotic normal approximation is suspect (parameters near a bound, a banana-shaped likelihood) — at the cost of running SIR with sir_keep_samples = true at fit time.

Bootstrap (subject resampling + refit) is not implemented; the UncertaintyMethod enum is deliberately open for it.

Simulation

simulate()

read_population_for_simulation()

simulate_with_seed()

simulate_with_options() — propensity-score matching

Why match?

Validation against R

predict()

simulate_with_uncertainty()

Result Types

Simulation Process

Simulation with Uncertainty Process

What the uncertainty distribution is

`simulate()`

`read_population_for_simulation()`

`simulate_with_seed()`

`simulate_with_options()` — propensity-score matching

`predict()`

`simulate_with_uncertainty()`