Simulate from a NLME model

Description

Simulates observations from a parsed model with between-subject variability and residual error. When a fit is supplied, the fitted theta, omega, and sigma replace the model file’s initial values - which is the usual flow after [ferx_fit](ferx_fit.qmd) (e.g. for posterior-predictive checks or VPCs).

Usage

ferx_simulate(
  model,
  data = NULL,
  n_sim = 1L,
  seed = 42L,
  fit = NULL,
  match = FALSE,
  horizon = NULL
)

Arguments

  • model: Path to a .ferx model file
  • data: Path to a NONMEM-format CSV (provides population structure: doses, obs times). When omitted, the model file’s [data] block (path = ...) is used. The DV column may be left empty (. / NA) on the sampling rows - the DV is what the simulation produces, so an empty cell means “simulate here” (a placeholder value is not needed). Rows marked MDV = 1 are excluded, as always. Kept empty-DV records are counted in the simulation_warnings attribute and re-emitted as an R warning: ferx_fit() skips those same records, so simulated rows at those times have no counterpart in a fit’s sdtab (do not overlay the two, e.g. in a VPC).
  • n_sim: Number of simulation replicates
  • seed: Random seed for reproducibility
  • fit: Optional ferx_fit result. When provided, simulation uses fit$theta, fit$omega, and fit$sigma instead of the model file’s initial values.
  • match: Propensity-score matching method (default FALSE). When enabled, each replicate’s drawn etas are reassigned to subjects by propensity-score matching against the subjects’ fitted (posthoc) etas - Mahalanobis matching under the model omega. This pairs each subject’s observed dosing/sampling design with a similar drawn eta, correcting VPC bias from treatment adaptation in real-world data (e.g. longer dosing intervals for high-clearance patients). Accepts:
  • FALSE / "none": No matching (default).
  • TRUE / "optimal": Global linear-assignment minimum Mahalanobis distance (MatchIt(method = "optimal")); best on average in simulation, the recommended method.
  • "nearest": Greedy nearest-neighbour (MatchIt(method = "nearest", distance = "mahalanobis")).
  • "rank": Pair by the rank of each eta’s Mahalanobis norm.Requires data to be real observed data (every subject must have observations, so its posthoc eta can be computed). The posthoc etas are computed at the fitted parameters when fit is supplied, otherwise at the model file’s initial values.
  • horizon: Optional administrative censoring time for time-to-event (TTE) endpoints. A finite, positive horizon is required to simulate a drug-driven (joint PK-TTE) model: the augmented hazard ODE is integrated until the cumulative hazard reaches -log U, censoring at horizon if no event fires (ferx-core #564). Purely-Gaussian models ignore it. NULL (default) leaves it unset.

Which predictive distribution you get

ferx_simulate() draws a fresh set of random effects for every ID in data, in every replicate, and adds residual error on top. At the default match = FALSE it does not read the observed DV at all (see the data argument), so each subject’s draw is unconditional: the prior predictive under the estimated parameters, not a posterior one. match is the exception, on the assignment rather than on the draw. The pool of etas per replicate is the same unconditional sample with the same marginal distribution, but which member of it lands on a given subject is chosen by matching against that subject’s posthoc eta - computed from that subject’s observed DV, which is why matching requires every subject to carry observations. So under match the eta paired with a particular observed design is informed by that design’s data; that is the point of it (it restores the design-eta association that adaptive dosing puts in real-world data), and it means the per-subject draw is no longer the unconditional one described above. What that distribution is of follows from what one ID means in your data and what level the etas were estimated at. The two are the same thing in individual-level PK, where an ID is a patient - but not in a model-based meta-analysis, where a row is a trial-arm summary, an ID is a study or an arm and the etas are between-study random effects. There, each replicate is a set of new studies, and DV_SIM is the predictive distribution of the next trial’s readout - not of the next patient. A patient-level interval needs a patient-level random effect in the model and rows at patient level in the data; ferx will not manufacture one, and quoting a between-study interval as a between-patient one is the error this note exists to prevent. Within a run, pick the column that matches the question:

  • IPRED: the mean response of the newly drawn study/subject - random effects, no residual error.
  • DV_SIM: that study’s/subject’s observed readout - random effects plus residual error. This is the column a VPC and a predictive interval are built from.

Two neighbouring functions cover the other levels: [ferx_predict](ferx_predict.qmd)gives the typical-value curve (all etas at zero, no residual error), and [ferx_simulate_with_uncertainty](ferx_simulate_with_uncertainty.qmd) adds parameteruncertainty (a theta/omega/sigma draw per replicate) on top of the random effects, which is what a decision interval for a future trial usually wants. For the individual etas conditioned on observed data, see fit$ebe_etas and [ferx_conddist](ferx_conddist.qmd).

Seealso

Other simulation: [ferx_calc_npde](ferx_calc_npde.qmd), [ferx_predict](ferx_predict.qmd), [ferx_predict_survival](ferx_predict_survival.qmd), [ferx_simulate_adaptive](ferx_simulate_adaptive.qmd), [ferx_simulate_with_uncertainty](ferx_simulate_with_uncertainty.qmd)

Concept

simulation

Value

A data.frame. Gaussian rows carry DRAW, SIM, ID, TIME, CMT, IPRED, DV_SIM (with OBSERVED = NA). For a joint PK-TTE model each subject also yields a TTE row on the event CMT, where TIME is the sampled event/censor time and OBSERVED is 1 (event before horizon) or 0 (right-censored at it); its IPRED and DV_SIM are NA. Use is.na(OBSERVED) to separate continuous rows from event rows. For a [binary_model] endpoint the categorical row on the binary CMT carries the simulated 0/1 outcome in DV_SIM (coded as the input CSV codes DV), with IPRED and OBSERVED both NA; select it by its CMT. The returned frame carries a simulation_warnings attribute (a character vector, empty for a clean run) listing any per-subject simulation diagnostics from ferx-core (e.g. a degenerate or pathological hazard that censored a subject with no event); these are also raised as an R warning.

Examples

ex <- ferx_example("warfarin")
fit <- ferx_fit(ex$model, ex$data, method = "gn", covariance = FALSE)
sim <- ferx_simulate(ex$model, ex$data, n_sim = 10L, seed = 1L, fit = fit)
head(sim)

# Propensity-score-matched simulation for a real-world-data VPC:
sim_pm <- ferx_simulate(ex$model, ex$data, n_sim = 10L, seed = 1L,
                        fit = fit, match = "optimal")