Simulate with parameter-uncertainty propagation

Description

For each parameter set drawn from the uncertainty distribution, the per-subject random-effect / residual-error simulator runs n_sim_per_draw times. The result includes both individual variability (etas, epsilons) and parameter uncertainty - useful for uncertainty-aware VPCs, dose-recommendation intervals, and any analysis where treating the ML estimates as fixed would understate variability.

Usage

ferx_simulate_with_uncertainty(
  model,
  data,
  fit,
  n_uncertainty_draws = 100L,
  n_sim_per_draw = 1L,
  method = c("asymptotic", "sir"),
  seed = 42L
)

Arguments

model: Path to a .ferx model file
data: Path to a NONMEM-format CSV (provides population structure). 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.
fit: A ferx_fit result. Must carry either cov_matrix (asymptotic) or sir_resamples (SIR) depending on method.
n_uncertainty_draws: Number of parameter sets to draw from the uncertainty distribution
n_sim_per_draw: Number of eta/eps replicates per parameter draw
method: Either "asymptotic" (default) or "sir"
seed: Random seed for reproducibility

Details

Two uncertainty sources are supported:

method = "asymptotic": Multivariate normal around the ML estimate in the engine’s packed (log-theta, Cholesky-omega, log-sigma) parameter space, using fit$cov_matrix. Requires fit to come from a ferx_fit() call with covariance = TRUE.
method = "sir": Sample with replacement from fit$sir_resamples. Requires the fit to have been run with sir = TRUE and sir_keep_samples = TRUE (passed via settings).

c(“## Bounded thetas andtheta used as inv_logit(logit(THETA) + ETA) (alogit_probability theta, such as THETA_F in the bundledbioavailability example) is drawn on the logit scale, not the packedscale. Its centre is logit(THETA) and its variance comes fromfit$cov_matrix by the delta method, the same one behind the-scale confidence interval [ferx_estimates](ferx_estimates) reports fortheta. Each draw is then mapped back through inv_logit(), sodraw lies strictly inside (0, 1) and follows the logit-normalthe fit implies. SIR proposals are built the same way. Before-core #1548 these draws were log-normal: draws above a declared upperbelow 1 were rejected, which truncated the distribution, and draws were clamped to a probability of 1.FIX’d logit_probability theta is not drawn. A free one whoseis not strictly inside (0, 1) has no logit, so the engine refusesdraw it. Declare its upper bound below 1 (e.g. 0.999) and refit.”, “## list("logit_probability")theta used as inv_logit(logit(THETA) + ETA) (alogit_probability theta, such as THETA_F in the bundledbioavailability example) is drawn on the logit scale, not the packedscale. Its centre is logit(THETA) and its variance comes fromfit$cov_matrix by the delta method, the same one behind the-scale confidence interval [ferx_estimates](ferx_estimates) reports fortheta. Each draw is then mapped back through inv_logit(), sodraw lies strictly inside (0, 1) and follows the logit-normalthe fit implies. SIR proposals are built the same way. Before-core #1548 these draws were log-normal: draws above a declared upperbelow 1 were rejected, which truncated the distribution, and draws were clamped to a probability of 1.FIX’d logit_probability theta is not drawn. A free one whoseis not strictly inside (0, 1) has no logit, so the engine refusesdraw it. Declare its upper bound below 1 (e.g. 0.999) and refit.”, “## drawstheta used as inv_logit(logit(THETA) + ETA) (alogit_probability theta, such as THETA_F in the bundledbioavailability example) is drawn on the logit scale, not the packedscale. Its centre is logit(THETA) and its variance comes fromfit$cov_matrix by the delta method, the same one behind the-scale confidence interval [ferx_estimates](ferx_estimates) reports fortheta. Each draw is then mapped back through inv_logit(), sodraw lies strictly inside (0, 1) and follows the logit-normalthe fit implies. SIR proposals are built the same way. Before-core #1548 these draws were log-normal: draws above a declared upperbelow 1 were rejected, which truncated the distribution, and draws were clamped to a probability of 1.FIX’d logit_probability theta is not drawn. A free one whoseis not strictly inside (0, 1) has no logit, so the engine refusesdraw it. Declare its upper bound below 1 (e.g. 0.999) and refit.” )## 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.

Concept

simulation

Value

A data.frame with columns: DRAW, SIM, ID, TIME, CMT, IPRED, DV_SIM, OBSERVED. CMT is the observation compartment; OBSERVED is NA throughout (this path is Gaussian-only - a drug-driven ODE-TTE endpoint is not supported here, and OBSERVED never echoes the input DV; see [ferx_simulate](ferx_simulate.qmd) for what each column holds). Row count: n_uncertainty_draws * n_sim_per_draw * n_obs.

Examples

ex  <- ferx_example("warfarin")
fit <- ferx_fit(ex$model, ex$data, covariance = TRUE)

# Asymptotic (default): fast, MVN draws around the ML estimate
sims <- ferx_simulate_with_uncertainty(
  ex$model, ex$data, fit,
  n_uncertainty_draws = 200, n_sim_per_draw = 10
)
head(sims)         # SIM, ID, TIME, IPRED, DV_SIM
range(sims$DV_SIM) # sanity check on simulated range

# SIR method: requires sir = TRUE + sir_keep_samples = TRUE at fit time
fit_sir <- ferx_fit(ex$model, ex$data,
                    covariance = TRUE, sir = TRUE,
                    settings = list(sir_keep_samples = TRUE))
sims_sir <- ferx_simulate_with_uncertainty(
  ex$model, ex$data, fit_sir,
  n_uncertainty_draws = 200, n_sim_per_draw = 10,
  method = "sir"
)

# Summarise: 90% prediction interval per time point
pi90 <- aggregate(DV_SIM ~ TIME, data = sims,
                  FUN = function(x) quantile(x, c(0.05, 0.5, 0.95)))