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 reported 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). Every other data-reader diagnostic the engine raises for data - a dose that never landed, a covariate with no value for some subjects, and so on - travels the same channel.
  • 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).

Errors raised by the engine

A model or dataset the engine refuses is an R error, never a NULLresult. Where the engine’s validation pass names the finding, the error has class ferx_engine_error and carries code, block, line and suggestion - the same condition a refused [ferx_fit](ferx_fit.qmd) raises, so one tryCatch(..., ferx_engine_error = ) handler covers every entry point. A failure the validation pass cannot attribute to one diagnostic is raised as an ordinary error carrying the engine’s message.

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, one row per simulated record, with these columns:

  • DRAW: Parameter draw. Always 1 here; it is the grouping variable in [ferx_simulate_with_uncertainty](ferx_simulate_with_uncertainty.qmd).
  • SIM: Replicate index, 1 to n_sim.
  • ID: Subject identifier, as a character string, exactly as the ID column of data spelled it.
  • TIME: Sampling time from data - or, on a TTE row, the sampled event or censoring time (see OBSERVED).
  • CMT: Observation compartment: the data file’s CMT for a continuous row, the [event_model] / [binary_model] compartment for an event or categorical row.
  • IPRED: Individual prediction - the drawn subject’s mean response, random effects but no residual error.
  • DV_SIM: The simulated observation: IPRED plus residual error, and the column a VPC or a predictive interval is built from.
  • OBSERVED: Event indicator for a time-to-event row only - 1 = the event fired before horizon, 0 = administratively right-censored at it - and NA on every other row. It is not the input DV: nothing in this frame echoes the data file’s observations, since simulation is what produces that column. Compare against the observed data by joining data on ID/TIME yourself.

Gaussian rows carry DRAW, SIM, ID, TIME, CMT, IPRED and DV_SIM, with OBSERVED = NA. For a joint PK-TTE model each subject also yields a TTE row on the event CMT, whose 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) together with the engine’s data-reader diagnostics for data - the same ones ferx_fit() returns in fit$warnings; 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")