Fit a nonlinear mixed effects model

Description

Fits a NLME model using FOCE or FOCEI estimation with a Rust backend.

Usage

ferx_fit(
  model,
  data = NULL,
  method = NULL,
  covariance = NULL,
  verbose = NULL,
  bloq_method = NULL,
  threads = NULL,
  mu_referencing = NULL,
  sir = NULL,
  gradient = NULL,
  optimizer_trace = FALSE,
  scale_params = FALSE,
  inits_from_nca = FALSE,
  fd_hessian_step = NULL,
  settings = NULL,
  output = NULL,
  include_data = FALSE,
  ignore = NULL,
  accept = NULL,
  ignore_ids = NULL,
  ...
)

Arguments

• model: Path to a .ferx model file, or a [ferx_model](ferx_model.qmd) object.
• data: Path to a NONMEM-format CSV file. When omitted, it defaults to the dataset declared in the model file’s [data] block (path = ...); passing data here overrides that. Required columns: ID, TIME, DV, EVID, AMT, CMT. Optional columns recognised by the engine: RATE (infusion rate; RATE = -1 infuses AMT at a modeled rate given by a per-subject parameter R{n}, and RATE = -2 over a modeled duration given by D{n}, on dose compartment n), MDV (missing-DV flag), II (dosing interval, required when SS > 0), SS (steady-state flag: 1 = pre-dose at steady state, 2 = add SS concentration to current state), CENS (LOQ censoring flag for M3 method: 1 below LLOQ, -1 above ULOQ), OCC (occasion index for IOV), and any covariate columns referenced in the model. See the steady-state section below for SS/II details.
• method: Estimation method(s). NULL (the default) uses whatever the model file’s [fit_options] method specifies, falling back to "focei" if the model file sets none (the engine warns when it defaults). Pass a value here to override the model file: either a single string or a character vector of methods to run in sequence. Each stage is seeded with the previous stage’s converged parameters, and only the final stage produces the reported covariance/diagnostics. Supported methods: "foce", "focei", "laplace" (alias "laplacian"; the Laplace approximation with the exact Hessian - NONMEM $EST METHOD=1 LAPLACIAN INTER. This is not the same estimator as "focei", which builds its Gaussian from the Gauss-Newton Hessian and reports a different OFV. At the default n_agq = 1 it is a single node; settings = list(n_agq = N) with N > 1 turns it into adaptive Gauss-Hermite quadrature - the exact conditional likelihood evaluated on a Gauss-Hermite grid around each subject’s empirical-Bayes mode, which makes no Gaussian-residual assumption and so handles non-Gaussian endpoints (time-to-event, categorical). Cost is n_agq^n_eta per subject per iteration, so higher node counts suit models with few random effects; n_agq is capped at 21, and a fit whose tensor grid would exceed 100000 nodes is rejected at check time. Odd values are conventional. Supports IOV at any occasion count. There is no separate "agq" method - adaptive quadrature is this method with n_agq > 1; "focei" likewise accepts n_agq > 1 for the Gauss-Newton-anchored quadrature), "saem", "gn" (Gauss-Newton / BHHH), "gn_hybrid" (Gauss-Newton followed by a FOCEI polish step), "imp" (also accepted as "importance_sampling" or "importance-sampling"; the NONMEM METHOD=IMP importance-sampling Monte-Carlo EM estimator), or "impmap" (also accepted as "importance_sampling_map"; Importance Sampling assisted by Mode A Posteriori, the NONMEM METHOD=IMPMAP Monte-Carlo EM estimator). "imp" is an estimator by default (it updates parameters) and may run standalone (method = "imp"), lead, or sit mid-chain. Set settings = list(imp_eval_only = TRUE) (NONMEM EONLY=1) to make it instead evaluate the marginal -2 log L at the fixed input parameters; in that mode it must be the last entry of the chain, e.g. c("focei", "imp"). Plain "imp" re-centers its proposal from the previous iteration’s samples and so is fragile on rich data; prefer "impmap" or warm-start with c("focei", "imp") there. "impmap" is a full estimator: it may run standalone (method = "impmap") or as a chain stage (c("focei", "impmap")), requires a mu-referenced parameterization, and does not yet support IOV. Example chain: c("saem", "focei"). "bayes" (also accepted as "bayesian" or "mcmc") runs full MCMC Bayesian estimation (Gibbs-within-HMC, NONMEM METHOD=BAYES parity): it returns posterior means with credible intervals and convergence diagnostics on $bayes rather than a point estimate, and runs standalone. Supports BSV and zero-mean inter-occasion variability (per-occasion kappa); the IOV variance posterior appears as OMEGA_IOV(...) in $bayes. SAEM fully supports inter-occasion variability (IOV / kappa) models.
• covariance: Logical, or NULL (the default). Whether to compute the covariance step for standard errors. NULL uses the model file’s [fit_options] covariance (engine default TRUE when unset); a logical overrides the model file. Previously defaulted to TRUE and silently overrode a model file that set covariance = false (#558).
• verbose: Logical, or NULL (the default). Print progress during estimation. NULL uses the model file’s [fit_options] verbose (engine default TRUE when unset); a logical overrides it.
• bloq_method: Handling of observations outside quantification limits. NULL (default) keeps whatever the model file specified; "m3" enables Beal’s M3 likelihood (requires a CENS column in the data, with DV carrying the limit value: LLOQ on CENS=1 rows and ULOQ on CENS=-1 rows); "drop" disables M3 and treats censored rows as ordinary observations – it does not remove them. Each censored row is fitted at the limit value in DV as though it had been measured there, which biases the fit. To genuinely exclude them, use ignore = "CENS == 1".
• threads: Number of worker threads for the per-subject parallel loops in the Rust backend (inner EBE search, SAEM, SIR). NULL (default) uses the engine’s default: available cores minus one (floored at 1), capped at 8 – most fits gain little from spreading across every core, and not all cores are equal on asymmetric platforms (e.g. Apple Silicon E-cores). Pass a positive integer to pin the count. The setting is per-call, so successive fits in the same R session can use different values.
• mu_referencing: Logical, or NULL (the default). When TRUE, automatically detects mu-referencing from the model structure for faster and more accurate convergence. NULL uses the model file’s [fit_options] mu_referencing (engine default TRUE when unset). Applies to all estimation methods. Set to FALSE to disable for comparison purposes. Detection works automatically for standard parameterizations such as PARAM = THETA * exp(ETA); unusual parameterizations fall back silently to zero-centred ETA initialisation with no error. No changes to the .ferx model file are needed. Check fit$warnings to see which ETAs were detected.
• sir: Logical, or NULL (the default); run Sampling Importance Resampling after the fit to produce non-parametric parameter uncertainty intervals. Requires covariance = TRUE. NULL uses the model file’s [fit_options] sir (engine default FALSE when unset). Tuning knobs (sir_samples, sir_resamples, sir_seed) still flow through settings.
• gradient: Inner-loop (per-subject EBE) gradient method. One of "auto" or "fd", or NULL (the default) to use the model file’s [fit_options] gradient (engine default "auto" when unset). "ad" is no longer supported: the Enzyme automatic-differentiation path was retired in favour of the analytic Dual2 sensitivities, so a model file (or call) carrying gradient = ad now fails validation with E_AD_RETIRED rather than being tolerated.The inner optimizer is BFGS; what differs is how the gradient of the individual NLL w.r.t. ETA is computed:
• "fd": central finite differences, 2 * n_eta forward evaluations per gradient call. Always available.
• "auto" (default): let the engine choose the best available gradient path for the model. This is the right choice for almost all uses.When to set "fd" explicitly. Mainly for validation or for reproducibility against a known finite-difference baseline.Set FERX_TIME_GRADIENTS=1 in the environment to print per-gradient-call timings at the end of a fit, which is the easiest way to check which method is faster on your specific model/data.
• optimizer_trace: Logical. If TRUE, write a per-iteration CSV trace to a temporary file, store its path in fit$trace_path, and read it into fit$trace as a data frame. Pass the result to [ferx_trace](ferx_trace.qmd) or [plot.ferx_fit](plot.ferx_fit.qmd) to inspect optimizer progress. Default FALSE.
• scale_params: Logical. If TRUE, apply a per-coordinate scaling layer on top of the existing log/Cholesky parameterization, dividing each transformed coordinate by |x0[i]| (when |x0[i]| > 0.1, otherwise by 1.0) so the outer optimizer works in a near-unit-magnitude space. Default FALSE.Why it defaults to FALSE. The scaling is not trajectory-transparent. Although it leaves the OFV value unchanged at any fixed point, it rescales the gradient the optimizer sees, and that gradient feeds the SLSQP overshoot cap, the quasi-Newton Hessian estimate, and the xtol/ftol termination - all of which act in the scaled coordinate system. The scaling layer only ever runs on log/Cholesky-packed coordinates (it auto-disables when any identity-packed theta is present), and for those, dividing by |log value| is counterproductive: a coordinate like log(V) = log(20) ~ 3 gets scale 3, so the optimizer’s unit step becomes a 3-unit move in log space - an e^3 ~ 20x multiplicative jump in V. That large step both overshoots and, through the uniform gradient cap, starves the step in every other dimension (notably OMEGA), so the fit can halt well short of the minimum. See ferx-core issue #99. The FALSE default reproduces the well-tested pre-scaling-layer behaviour.When to set TRUE. Left as an opt-in for experimentation - e.g. to A/B the scaled vs unscaled trajectory on a specific model. When enabled it applies to all outer optimizers: NLopt FOCE/FOCEI (BOBYQA, SLSQP, L-BFGS, MMA), the hand-rolled BFGS, Gauss-Newton / BHHH, and the SAEM M-step.
• inits_from_nca: Derive NCA-based starting values from the data before the optimizer runs, overriding the model file’s defaults. Either a logical (FALSE, the default, disables it; TRUE is an alias for "nca_sweep") or one of "nca", "nca_sweep", "nca_ebe" to pick a strategy explicitly. Most useful with settings = list(optimizer = "trust_region") or method = "gn", where bad starting values can stall the optimizer. See [ferx_inits_from_nca](ferx_inits_from_nca.qmd) to inspect the values without fitting.
• fd_hessian_step: Positive finite numeric. Relative step size for the finite-difference Hessian used in the covariance step (default 0.01). The actual perturbation for parameter i is fd_hessian_step * (1 + |x_hat[i]|). Increase (e.g. 0.1) when the fit warns about ill-conditioned Hessian entries; decrease (e.g. 1e-3) on smooth OFV surfaces where FD noise is the primary concern. Has no effect when covariance = FALSE, nor when the exact analytic R-matrix serves the fit – which is now the default for in-scope models, so on those this argument is inert. See settings = list(analytic_cov_hessian = ...) below for the scope and for how to force finite differences back on.
• settings: Optional named list of fine-grained options forwarded to the Rust FitOptions. Use this to tune knobs that do not have a dedicated ferx_fit() argument. Keys are validated: values that duplicate a dedicated argument (method, covariance, verbose, bloq_method, bloq, threads, sir, gradient, gradient_method) are rejected – pass them via the dedicated argument. Unknown keys and malformed values also raise an error.Every key the engine accepts is listed below, with two deliberate exceptions: frem_predictions and frem_sigma are structural maps written into the generated model file by [ferx_model_to_frem](ferx_model_to_frem.qmd) and cannot sensibly be hand-authored.Precedence: dedicated ferx_fit() arguments win over settings, which in turn win over the model file’s [fit_options] block. A warning() is issued whenever a call-time value overrides a different value from [fit_options]. Inspect fit$model_file_settings and fit$call_settings to audit the full picture.Shared options (FOCE / FOCEI / Laplace / GN / GN-hybrid / SAEM / IMP / IMPMAP / Bayes)The engine accepts the two inner-loop keys for every method that runs an inner EBE solve, which includes "laplace", "imp", "impmap" and "bayes" as well as the five above.
• inner_maxiter: Per-subject EBE iteration cap (default 200).
• inner_restarts: (default 1) Guarded multi-start count for the inner (per-subject EBE) optimizer, to escape a multimodal individual objective – saturable protein binding is the classic case. A subject re-solves its EBE from this many Omega-scaled alternate seeds and keeps the lowest-objective mode; an alternate seed that reconverges to the same mode is not accepted, so unimodal subjects are bit-identical to inner_restarts = 0. The covariance step reconverges with the same setting. FOCE / FOCEI / Laplace / GN / GN-hybrid.
• inner_optimizer: "auto" (default), "bfgs", "lbfgs" or "nelder_mead". Inner-loop algorithm. "auto" uses dense BFGS, switching to limited-memory L-BFGS above 32 random effects; an explicit value pins one algorithm with no switching. All gradient-based variants converge to the same EBE, so this trades per-step cost against memory rather than accuracy.
• cov_inner_tol: Inner EBE-reconvergence tolerance used only by the covariance step, decoupled from inner_tol (default: whatever inner_tol is; for LTBS models min(inner_tol, 1e-8)). The covariance R-matrix is a second-difference of the reconverged objective and so is more sensitive to EBE precision than the fit itself – on a flat surface an EBE converged only to inner_tol can visibly perturb the standard errors. Worth reaching for on heavily-censored M3 + IOV models (try 1e-11). Note the engine may emit a spurious “not used by method … will be ignored” warning for this key; the value is applied.
• parameter_scaling: "auto" (default), "none", "abs" or "rescale2". Parameter-scaling strategy for the outer optimizer; supersedes scale_params when not "none". "auto" applies "rescale2" to the gradient-based optimizers that benefit (nlopt_lbfgs, slsqp) and leaves the derivative-free bobyqa unscaled, where it would distort the trust region. "rescale2" is the nlmixr2-style normalisation – each packed coordinate divided by its bound half-range – which markedly improves cold-start convergence. "abs" is the legacy |packed value| normalisation, equivalent to scale_params = TRUE.
• ebe_warm_start: (default FALSE) When an inner EBE solve fails its BFGS step and falls back to Nelder-Mead, warm-start the simplex from the BFGS partial eta-hat rather than cold-starting from eta = 0. Substantially fewer prediction walks on fallback-heavy fits (~1.7x on a 2-cpt unidentifiable-V2 benchmark). Opt-in because it moves the fallback subjects’ EBEs, which perturbs the outer optimizer’s trajectory – harmless for the derivative-free default, but it can derail a gradient-based outer optimizer into a worse basin. Validate OFV and estimates on your own model before enabling.
• checkpoint: (default TRUE) and checkpoint_interval_secs (default 300). Periodically save a resume point so an interrupted run restarts from where it stopped. Restart is coarse: the population estimates at the last save become the new starting point, so a resumed run may need a few extra iterations to re-converge. A change to the model or data invalidates the checkpoint via a hash check and the run starts fresh; a run shorter than the interval leaves nothing behind.
• iov_column: Name of the occasion column in the dataset (e.g. "OCC") for inter-occasion variability.
• iov_occasion: Derive the occasion partition from the model instead of a data column: "column" (default, use iov_column), "dose" to start a new occasion at each administration, or "time(24, 48)" for time-window breakpoints. "dose" and "time(...)" override iov_column (with a warning). Supported by FOCE, FOCEI, Laplace and SAEM – not IMP/IMPMAP.
• npde_nsim: (default 0, i.e. disabled) Monte-Carlo replicates per subject for the simulation-based NPDE/NPD diagnostics computed after the fit; with 0 no NPDE/NPD columns are produced. A typical value is 1000, and cost scales linearly. npde_seed makes the simulation reproducible. See also [ferx_calc_npde](ferx_calc_npde.qmd) to add them post hoc.
• inner_tol: Gradient-norm convergence tolerance for the inner (per-subject EBE) loop (default 1e-5). A looser tolerance leaves residual noise in each subject’s EBE solution, which propagates into the marginal objective the outer optimizer sees; the engine moved from 1e-4 to 1e-5 to reach NONMEM’s minimum, at roughly 1.5x the per-fit cost. Note this changes the converged point, so it is not a pure cost knob. Tighter is not uniformly better – at 1e-6 some ill-conditioned fits over-converge the inner Hessian and the outer optimizer can land in a worse basin.
• fd_hessian_step: Also available as the dedicated fd_hessian_step argument above.
• analytic_cov_hessian: TRUE (default) or FALSE. When TRUE and the model is in scope, the covariance step assembles the exact analytic R-matrix from third-order sensitivities of the closed-form prediction instead of differencing the objective. That removes both the fd_hessian_step tuning knob and the eps/h^2 differencing noise, and costs 2 * (n_theta + n_eta) + 1 sensitivity evaluations per subject rather than roughly 2 * n_free^2 objective evaluations that each re-solve every subject’s inner loop. Serves method = "focei" and method = "foce" alike, from two separate assemblies rather than one shared formula – the non-interaction case is built on the Sheiner-Beal gradient and carries no log|H~| term – so both are exact, and interaction does not silently drop you back onto finite differences. Scope is otherwise deliberately narrow – plain analytical (closed-form) Gaussian models only: no IOV, LTBS, [scaling], Form-C readout, [initial_conditions], iiv_on_ruv, block_sigma or custom-magnitude residual, M3 censoring, FREM, non-Gaussian (TTE / categorical / Markov) endpoint, and not method = "laplace" at any node count, method = "focei" with n_agq > 1, a covariate-Selected error spec, or gradient = "fd". Anything outside it – including a single out-of-scope subject – silently keeps the finite-difference stencil, all-or-nothing across the population. Set FALSE to force finite differences for an in-scope model, e.g. to reproduce standard errors from an earlier run.
• covariance_method: Covariance estimator, mirroring NONMEM $COV MATRIX=: "r" (inverse Hessian, the default), "s" (score cross-product), or "rsr" (the Huber-White sandwich, robust to model misspecification). "s" and "rsr" are supported for FOCEI, FOCE and IOV fits alike, all three anchored against NONMEM $COV MATRIX=S/RSR within ~10
  to "r" while NONMEM’s $COVARIANCE default is "rsr" – set "rsr" when reconciling against a NONMEM run that used the default $COV. No effect when covariance = FALSE.
• covariance_fallback: "none" (default) or "sir". When the finite-difference Hessian is not positive definite, "sir" runs SIR with an absolute-eigenvalue-rectified proposal instead of leaving the covariance step failed; covariance_status is then "sir_fallback" and SIR-based credible intervals are reported.ODE models: solver method and tolerance
• ode_method: Which stepper integrates the [odes] block: "rk45" (default), "vern7", "rosenbrock23", "rodas4" or "rodas5p". There are two independent reasons a fit’s step size is small, and they want opposite fixes, so diagnose before switching. If steps stay tiny whatever ode_reltol you ask for, the model is stability-limited (stiff – fast reversible binding / TMDD, Michaelis-Menten with KM far below observed concentrations, long transit chains, QSP cascades): use one of the linearly implicit Rosenbrock methods, "rosenbrock23" at crude tolerances, "rodas4" at a typical 1e-6–1e-9, or "rodas5p" at 1e-9 and tighter. If instead nearly every step is accepted and it is tightening ode_reltol that drives the step count up, the fit is accuracy-limited and a stiff method buys nothing – "vern7" (higher order) is the lever, worth about 2.3x at 1e-9 on ferx-core’s transit benchmark but about 1.4x slower at default tolerances. A stiff step is not free: it costs n + 1 extra right-hand-side evaluations for the finite-difference Jacobian plus an n x n factorization, where n is the size of the system actually integrated (for a continuous-time Markov endpoint that is the s^2 occupancy system, not the [odes] state count). Every method is a full peer – analytic sensitivities, time-to-event and categorical endpoints, simulation and adaptive dosing all work with all of them. Can also be set in the model file’s [fit_options] block.
• ode_reltol: Relative tolerance for ODE models (default 1e-4; ignored for analytical PK). The default reproduces analytical closed forms in PRED to about 1e-4, but the FOCE objective amplifies solver error, so an ODE-form model’s OFV can differ from its analytical equivalent by several units. Set tighter (e.g. 1e-10) when the ODE-form OFV must match an analytical reference; expect slower fits. Can also be set in the model file’s [fit_options] block.
• ode_abstol: Absolute tolerance for ODE models (default 1e-6).
• ode_max_steps: Maximum solver steps per integration segment (default 10000). Raise if a tight ode_reltol exhausts the step budget on stiff multi-compartment systems – or switch ode_method to a stiff stepper, which is the actual fix when stiffness rather than accuracy is capping the step.FOCE / FOCEI / Laplace / GN / GN-hybrid: iteration cap
• maxiter: Maximum outer-optimizer iterations (default 500). Not applicable to SAEM, which controls iterations via n_exploration and n_convergence.FOCE / FOCEI / Laplace / GN-hybrid: outer optimizer
• optimizer: Population-level optimizer. One of "auto" (default; picks "nlopt_lbfgs" when the exact analytic FOCE/FOCEI gradient is available and "bobyqa" when only finite differences are – the fit reports the resolved pick as "auto (<resolved>)"), "bobyqa" (derivative-free, robust), "slsqp" (sequential quadratic programming, gradient-based), "lbfgs" / "nlopt_lbfgs" (limited-memory BFGS via NLopt, gradient-based), "mma" (method of moving asymptotes, gradient-based), "bfgs" and "lbfgs" (deprecated aliases that now select "nlopt_lbfgs" – the hand-rolled built-in BFGS is no longer reachable and is slated for removal), or "trust_region" (trust-region Newton, gradient-based). Gradient-based optimizers use the inner gradient method set by the gradient argument; "bobyqa" does not. Not accepted by pure "gn".
• outer_xtol: Relative step tolerance for the derivative-free "bobyqa" outer optimizer (NLopt xtol_rel; default 1e-4). Gradient-based optimizers use their own fixed internal tolerances. Can also be set in the model file’s [fit_options].
• outer_ftol: Relative objective tolerance for the derivative-free "bobyqa" outer optimizer (NLopt ftol_rel). Unset = auto: 1e-8 for a pure time-to-event model (its hazard objective is evaluated exactly) and 1e-6 otherwise. The TTE tightening lands the frailty variance on the NONMEM/nlmixr2 minimum across a near-flat omega-squared ridge (ferx-core #469); the 1e-6 floor elsewhere avoids grinding on noisy ODE / FD-inner objectives, where 1e-8 is unreachable. Set an explicit value to pin it for every model.
• global_search: Logical. When TRUE, run a global search phase before local refinement (default FALSE). Not accepted by pure "gn".
• global_maxeval: Function evaluations budget for the global search phase (default 0, i.e. disabled when global_search = FALSE). Not accepted by pure "gn".
• stagnation_guard: Logical (default TRUE). Terminates the NLopt outer loop early when the OFV plateau is numerically flat. Set FALSE to let SLSQP / L-BFGS run to their own xtol/ftol or maxiter. Not accepted by pure "gn".
• reconverge_gradient_interval: Integer (default 0). How often to re-solve each subject’s inner EBE during the population gradient, instead of holding it fixed. 0 = never (cheap fixed-EBE gradient); 1 = every gradient evaluation; N = every N-th. Reconverging recovers the full gradient surface at roughly 5-6x the per-call cost and can help gradient-based optimizers (e.g. slsqp) past ill-conditioned non-IOV plateaus. IOV models always reconverge and ignore this setting. Not accepted by pure "gn".Trust-region optimizer (optimizer = "trust_region")
• steihaug_max_iters: Conjugate-gradient iteration budget for the trust-region subproblem (default chosen by the engine). Only consumed when optimizer = "trust_region" is set under FOCE / FOCEI / GN-hybrid. Increase if the subproblem solver exits too early on ill-conditioned problems.SAEM
• n_exploration: Stochastic exploration phase iterations (default 150).
• conddist: (default FALSE; alias saem_conddist) Run a post-fit conditional-distribution pass that estimates each subject’s p(eta_i | y_i) by MCMC, reporting per-subject conditional mean and SD plus distribution-based eta-shrinkage – the shrinkage-unbiased basis for eta diagnostics. Read it with [ferx_conddist](ferx_conddist.qmd). The three keys below apply only when this is TRUE.
• conddist_nsamp: Retained draws per subject (default 200; production diagnostics usually want thousands).
• conddist_burnin: Draws discarded before accumulation (default 20), to forget the EBE-mode warm start.
• conddist_keep_samples: (default FALSE) Retain the raw draws rather than just the summaries.
• n_convergence: Averaging / convergence phase iterations (default 250).
• n_mh_steps: (default 20) Metropolis-Hastings steps per subject per iteration. Also sizes the componentwise decorrelating kernel that prevents block-Omega collapse. When n_leapfrog > 0 it applies to the subjects that fall back to MH. Consumed by "saem" and "bayes".
• adapt_interval: How often (in iterations) the MH proposal covariance is adapted (default 50).
• omega_burnin: Initial iterations during which Omega is held fixed while the sampler warms up (default 20). Guards against Omega collapse on sparse data. Set to 0 to disable.
• seed / saem_seed: RNG seed for the SAEM Metropolis-Hastings sampler (default 12345). Independent of multi_start_seed.
• n_leapfrog / saem_n_leapfrog: Leapfrog steps for HMC proposals (default 0 = Metropolis-Hastings). A positive value replaces MH with one HMC proposal per subject per iteration.Bayes ("bayes")
• bayes_warmup: Warmup (burn-in + adaptation) sweeps per chain, discarded from the posterior (default 1000).
• bayes_iters: Retained sampling sweeps per chain, before thinning (default 1000).
• bayes_chains: Number of independent chains (default 4); used for split-R-hat.
• bayes_thin: Keep every bayes_thin-th sampling draw (default 1).
• bayes_seed: Base RNG seed for the Bayes sampler. Independent of seed / saem_seed.Gauss-Newton ("gn" / "gn_hybrid")
• gn_lambda: Levenberg-Marquardt damping factor (default 0.01). Larger values make steps more conservative. Accepted by both "gn" and "gn_hybrid".Importance Sampling ("imp")By default "imp" is a Monte-Carlo EM estimator (NONMEM METHOD=IMP); set imp_eval_only = TRUE to evaluate the marginal -2 log L at fixed parameters instead (NONMEM EONLY=1).
• imp_eval_only: Logical; TRUE evaluates -2 log L at the fixed input parameters without estimating (NONMEM EONLY=1; must be the terminal chain stage). FALSE (default) estimates.
• imp_iterations: Number of Monte-Carlo EM iterations, ignored when imp_eval_only (default 200).
• imp_averaging: Number of final iterations whose parameters are averaged into the reported estimate, ignored when imp_eval_only (default 50).
• imp_samples: Importance samples drawn per subject (default 1000). Halving the Monte-Carlo SE requires quadrupling this value.
• imp_proposal_df: Degrees of freedom for the Student-t proposal distribution (default 5); the string "normal" (or "mvn") selects a multivariate-normal proposal.
• imp_seed: RNG seed for the IS step (default chosen by the engine).
• imp_auto: (default TRUE) Adaptive sample count (NONMEM AUTO). With this on, imp_samples is the starting count and is ramped up (doubling per iteration, capped at 10000) while the objective’s Monte-Carlo SE exceeds 1.0. Strongly recommended for FREM and other high-dimensional models, where a fixed count leaves a sample-count-dependent bias in the typical-value and Omega estimates. Note this makes the “quadruple the samples to halve the MC SE” rule of thumb apply only when it is FALSE.
• imp_defensive_alpha: (default 0, i.e. off; alias impmap_defensive_alpha) Defensive-mixture weight in [0, 1). Each subject draws this fraction of its samples from the prior N(0, Omega) rather than the mode-centred proposal, and every sample is scored under the mixture density. That bounds the importance weights, so a weakly-identified subject cannot hijack the weighted M-step and walk theta to its bounds. Try 0.1. Enabling it disables Sobol QMC and raises the per-subject ESS floor.
• iscale_min: (default 0.1) and iscale_max (default 10). Bounds on the proposal scaling factor for adaptive importance sampling (NONMEM ISCALE_MIN / ISCALE_MAX): the proposal covariance is multiplied by s^2 with s chosen from this interval to maximise per-subject ESS. Set both to 1 to disable.
• frem_rao_blackwell: (default TRUE) FREM only: Rao-Blackwellise the covariate ETAs – integrate them analytically and importance-sample only the PK ETAs. Strongly recommended, since brute-force sampling of the near-singular covariate dimensions has very poor ESS. Set FALSE only to diagnose the Rao-Blackwell path against the full-dimensional sampler.
• imp_low_ess_threshold: ESS fraction below which a subject is flagged in fit$importance_sampling$low_ess_subject_ids (default 0.1, i.e. 10% of imp_samples).IMPMAP ("impmap" estimator)
• impmap_iterations: Number of Monte-Carlo EM iterations (default 200).
• impmap_samples: Importance samples drawn per subject per iteration (default 300).
• impmap_proposal_df: Proposal degrees of freedom (default 4). A finite value >= 1 gives a heavier-tailed Student-t, which is the ferx default; the string "normal" (or "mvn") gives a multivariate-normal proposal, which is NONMEM’s default. ferx diverges deliberately: a Gaussian’s lighter tails under-cover the posterior of weakly-identified parameters, so the importance weights blow up in the tail and bias the M-step moments.
• impmap_averaging: Number of final iterations whose parameters are averaged into the reported estimate (default 50).
• impmap_seed: RNG seed for the IMPMAP sampling (default chosen by the engine).
• impmap_auto: (default TRUE) Adaptive sample count, as imp_auto above – impmap_samples is the starting count.
• impmap_mceta: (default 0) Additional random starting points for the per-subject MAP optimization (NONMEM MCETA). Each start draws eta from N(0, Omega) and the lowest individual NLL wins; 0 means a single warm start. 3 is a good choice for high-dimensional models (e.g. FREM with 5 or more ETAs).
• impmap_sobol: (default FALSE) Use Sobol quasi-random sequences (with Cranley-Patterson randomization) for the importance draws instead of pseudo-random, giving more uniform posterior coverage per sample. Applies only to multivariate-normal proposals (impmap_proposal_df = "normal"); under the Student-t default it is inert.
• impmap_low_ess_threshold: ESS fraction below which a subject is flagged as poorly sampled (default 0.1).
• impmap_trace: Logical; when TRUE, collect per-iteration parameter values into fit$impmap_trace (analogous to NONMEM .ext output). Default FALSE.SIR (Sampling Importance Resampling)
• sir_samples: Candidate draws for SIR (default 1000).
• sir_resamples: Resampled draws retained for CI computation (default 250).
• sir_seed: RNG seed for the SIR step (default chosen by the engine). Independent of seed / saem_seed.
• sir_df: (default 5) Degrees of freedom for the SIR multivariate Student-t proposal; higher values approach a normal proposal.
• sir_keep_samples: (default FALSE) Retain the resampled parameter vectors, which [ferx_simulate_with_uncertainty](ferx_simulate_with_uncertainty.qmd) requires.Multi-start optimization
• n_starts: Number of optimizer starts (default 1, i.e. single start). When > 1, starts 2..n are initialized from log-space perturbations of the model’s initial values; the best OFV wins.
• start_sigma: Perturbation spread in log-space (default 0.3, approx. 30% CV). Log-packed thetas are multiplied by exp(N(0, start_sigma)); identity-packed thetas are shifted by start_sigma * N(0,1).
• multi_start_seed: RNG seed for the start-point perturbation (default 42). Independent of seed / saem_seed.
• output: Optional path to a .fitrx file. When non-NULL, [ferx_save_fit](ferx_save_fit.qmd) is invoked on the result so the fit is persisted to disk in a portable, cross-language bundle (zip of JSON + CSV) before the function returns. Equivalent to calling ferx_save_fit(fit, output) after the fit. See the format reference in the ferx-core docs for the on-disk schema.
• include_data: Logical. Only meaningful with output. When TRUE, embeds the input data CSV verbatim inside the .fitrx bundle so the file is self-contained. Default FALSE.
• ignore: Character vector of filter expressions that exclude a record when the expression evaluates to true (NONMEM $DATA IGNORE= equivalent). Expressions use standard comparison operators (==, !=, <, <=, >, >=) on column names. Use && within one expression for AND; for OR use multiple elements. These conditions are merged with any [data_selection] ignore rules in the model file. Example: ignore = c("DV < 0.001", "EVID != 0 && MDV == 1").
• accept: Character vector of filter expressions. A record is kept only when all conditions pass; excluded if any fail (NONMEM $DATA ACCEPT= equivalent). Merged with model-file accept rules.
• ignore_ids: Numeric or character vector of subject IDs to exclude entirely. Sugar for ignore = "ID == <id>" applied per-subject. Merged with [data_selection] ignore_subjects from the model file.
• ...: Reserved for future use. Unrecognised arguments raise an error.

Process noise (SDE / diffusion)

ODE-based PK/PD models occasionally produce autocorrelated IWRES when the structural model is misspecified (missing compartment, unmodelled feedback). Adding continuous within-subject process noise via an SDE framework - solved by an Extended Kalman Filter (EKF) - absorbs this drift and yields a better-calibrated likelihood. Add a [diffusion] block to the .ferx model file to enable SDE mode. Each line declares the diffusion variance for one ODE state:

[diffusion]
  central ~ 0.5   # initial estimate for DIFF_CENTRAL (variance units)

Key points:

• The declared value and the fitted estimate are variances, not standard deviations. A value of 0.5 means \(\sigma^2 = 0.5\), not \(\sigma = 0.5\).
• Each diffusion parameter appears in fit$theta as DIFF_<STATE> (e.g. DIFF_CENTRAL). Standard errors and fit$estimates treat them as regular thetas.
• SDE models use finite differences for the EKF covariance propagation.
• SAEM is not supported with SDE models; a hard error is raised.
• fit$uses_sde is TRUE whenever a [diffusion] block was present.

c(“## Unit conversion and scaling (an optional [scaling] block to the .ferx model file tomodel predictions before comparing them to the observations in DV.is useful when the model is parameterised in one unit (e.g. amounts in) but DV is recorded in another (e.g. concentrations in ng/mL).engine divides each prediction by the scale factor before computing, so the fitted sigma is on the DV scale.Form A – scalar divisor:\n[scaling]\n obs_scale = 1000 # divide every predicted value by 1000\nForm B – expression divisor:\n[scaling]\n obs_scale = V # divide by the individual volume (theta/eta expression)\nmay reference thetas, etas, and [individual_parameters]variables. Differentiated exactly on the defaultgradient = auto – both the outer and inner loops serve an expressionanalytically via the eta-quotient rule (ferx-core #486). Do not setgradient = fd for it.Per-CMT scaling (Form A or B per compartment):\n[scaling]\n obs_scale[CMT=1] = 1000 # CMT 1 in ng/mL, dose in ug\n obs_scale[CMT=2] = 1 # CMT 2 already in correct units\nForm C – ODE readout expression:For ODE models, y = <expr> defines the observation as a function ofstate variables, thetas, etas, and individual parameters. This replacesobs_cmt= argument in [structural_model] and is requiredthe observation is not a raw compartment amount:\n[structural_model]\n ode(states = [depot, central], ...) # no obs_cmt= here\n\n[scaling]\n y = central / V # observation = central compartment / volume\n-CMT Form C: y[CMT=1] = central / V1, y[CMT=2] = central2 / V2.Constraints:[scaling] is not yet supported on SDE ([diffusion]) models.A uniform* Form B or Form C scale is differentiated exactlythe default gradient = auto. The per-CMT variants(obs_scale[CMT=N], y[CMT=N]) fall back to finitesilently* – there is no error and no warning, sofit$gradient_used if it matters.”, “## list("[scaling]")an optional [scaling] block to the .ferx model file tomodel predictions before comparing them to the observations in DV.is useful when the model is parameterised in one unit (e.g. amounts in) but DV is recorded in another (e.g. concentrations in ng/mL).engine divides each prediction by the scale factor before computing, so the fitted sigma is on the DV scale.Form A – scalar divisor:\n[scaling]\n obs_scale = 1000 # divide every predicted value by 1000\nForm B – expression divisor:\n[scaling]\n obs_scale = V # divide by the individual volume (theta/eta expression)\nmay reference thetas, etas, and [individual_parameters]variables. Differentiated exactly on the defaultgradient = auto – both the outer and inner loops serve an expressionanalytically via the eta-quotient rule (ferx-core #486). Do not setgradient = fd for it.Per-CMT scaling (Form A or B per compartment):\n[scaling]\n obs_scale[CMT=1] = 1000 # CMT 1 in ng/mL, dose in ug\n obs_scale[CMT=2] = 1 # CMT 2 already in correct units\nForm C – ODE readout expression:For ODE models, y = <expr> defines the observation as a function ofstate variables, thetas, etas, and individual parameters. This replacesobs_cmt= argument in [structural_model] and is requiredthe observation is not a raw compartment amount:\n[structural_model]\n ode(states = [depot, central], ...) # no obs_cmt= here\n\n[scaling]\n y = central / V # observation = central compartment / volume\n-CMT Form C: y[CMT=1] = central / V1, y[CMT=2] = central2 / V2.Constraints:[scaling] is not yet supported on SDE ([diffusion]) models.A uniform* Form B or Form C scale is differentiated exactlythe default gradient = auto. The per-CMT variants(obs_scale[CMT=N], y[CMT=N]) fall back to finitesilently* – there is no error and no warning, sofit$gradient_used if it matters.”, “## block)an optional [scaling] block to the .ferx model file tomodel predictions before comparing them to the observations in DV.is useful when the model is parameterised in one unit (e.g. amounts in) but DV is recorded in another (e.g. concentrations in ng/mL).engine divides each prediction by the scale factor before computing, so the fitted sigma is on the DV scale.Form A – scalar divisor:\n[scaling]\n obs_scale = 1000 # divide every predicted value by 1000\nForm B – expression divisor:\n[scaling]\n obs_scale = V # divide by the individual volume (theta/eta expression)\nmay reference thetas, etas, and [individual_parameters]variables. Differentiated exactly on the defaultgradient = auto – both the outer and inner loops serve an expressionanalytically via the eta-quotient rule (ferx-core #486). Do not setgradient = fd for it.Per-CMT scaling (Form A or B per compartment):\n[scaling]\n obs_scale[CMT=1] = 1000 # CMT 1 in ng/mL, dose in ug\n obs_scale[CMT=2] = 1 # CMT 2 already in correct units\nForm C – ODE readout expression:For ODE models, y = <expr> defines the observation as a function ofstate variables, thetas, etas, and individual parameters. This replacesobs_cmt= argument in [structural_model] and is requiredthe observation is not a raw compartment amount:\n[structural_model]\n ode(states = [depot, central], ...) # no obs_cmt= here\n\n[scaling]\n y = central / V # observation = central compartment / volume\n-CMT Form C: y[CMT=1] = central / V1, y[CMT=2] = central2 / V2.Constraints:[scaling] is not yet supported on SDE ([diffusion]) models.A uniform* Form B or Form C scale is differentiated exactlythe default gradient = auto. The per-CMT variants(obs_scale[CMT=N], y[CMT=N]) fall back to finitesilently* – there is no error and no warning, sofit$gradient_used if it matters.” )c(“## Parameter declaration syntax (Theta (fixed effects):\ntheta CL(0.134, 0.001, 10.0) # (initial, lower, upper)\n theta CL(0.1, 0.001) # lower-bound only (no upper bound)\n theta CL(0.134, 0.001, 10.0) (FIX) # FIX at end (traditional)\n theta CL (FIX) (0.134, 0.001, 10.0) # FIX anywhere (flexible placement)\nOmega (inter-individual variability):\nomega ETA_CL ~ 0.07 # variance parameterisation (default)\n omega ETA_CL ~ 0.07 (FIX) # fixed omega, FIX at end\n omega ETA_CL (FIX) ~ 0.07 # fixed omega, FIX before the tilde\nsame flexible (FIX) placement applies to sigma andkappa (IOV) declarations.Unused-parameter warning:* a warning severity message with\"unused_parameter\" is emitted by the parser when ais declared in [parameters] but never referenced in[individual_parameters] or [error_model]. This usuallya commented-out expression or a typo in the parameter name.ferx_get_warnings(fit) or check ferx_model_validate()before fitting.”, “## list("[parameters]")Theta (fixed effects):\ntheta CL(0.134, 0.001, 10.0) # (initial, lower, upper)\n theta CL(0.1, 0.001) # lower-bound only (no upper bound)\n theta CL(0.134, 0.001, 10.0) (FIX) # FIX at end (traditional)\n theta CL (FIX) (0.134, 0.001, 10.0) # FIX anywhere (flexible placement)\nOmega (inter-individual variability):\nomega ETA_CL ~ 0.07 # variance parameterisation (default)\n omega ETA_CL ~ 0.07 (FIX) # fixed omega, FIX at end\n omega ETA_CL (FIX) ~ 0.07 # fixed omega, FIX before the tilde\nsame flexible (FIX) placement applies to sigma andkappa (IOV) declarations.Unused-parameter warning:* a warning severity message with\"unused_parameter\" is emitted by the parser when ais declared in [parameters] but never referenced in[individual_parameters] or [error_model]. This usuallya commented-out expression or a typo in the parameter name.ferx_get_warnings(fit) or check ferx_model_validate()before fitting.”, “## block)Theta (fixed effects):\ntheta CL(0.134, 0.001, 10.0) # (initial, lower, upper)\n theta CL(0.1, 0.001) # lower-bound only (no upper bound)\n theta CL(0.134, 0.001, 10.0) (FIX) # FIX at end (traditional)\n theta CL (FIX) (0.134, 0.001, 10.0) # FIX anywhere (flexible placement)\nOmega (inter-individual variability):\nomega ETA_CL ~ 0.07 # variance parameterisation (default)\n omega ETA_CL ~ 0.07 (FIX) # fixed omega, FIX at end\n omega ETA_CL (FIX) ~ 0.07 # fixed omega, FIX before the tilde\nsame flexible (FIX) placement applies to sigma andkappa (IOV) declarations.Unused-parameter warning:* a warning severity message with\"unused_parameter\" is emitted by the parser when ais declared in [parameters] but never referenced in[individual_parameters] or [error_model]. This usuallya commented-out expression or a typo in the parameter name.ferx_get_warnings(fit) or check ferx_model_validate()before fitting.” )c(“## Steady-state dosing (a dosing row as a steady-state dose by setting SS = 1 (orSS = 2) in the NONMEM CSV and providing the dosing intervalII (in the same time units as TIME).SS = 1: the subject is assumed to be at steadybefore this dose. The engine computes the steady-stateconditions analytically (for 1/2/3-cpt models) or viaexpansion (for ODE models) and uses them as the startingstate.SS = 2: the steady-state concentration is addedthe current compartment state (superposition).II: dosing interval (required when SS > 0).match the units of TIME. Rows with SS = 0 or missingSS ignore II.changes to the .ferx model file are needed. Steady-stateis purely data-driven and works for all PK model families(1-cpt, 2-cpt, 3-cpt oral/IV/infusion, ODE-based).”, “## list("SS")a dosing row as a steady-state dose by setting SS = 1 (orSS = 2) in the NONMEM CSV and providing the dosing intervalII (in the same time units as TIME).SS = 1: the subject is assumed to be at steadybefore this dose. The engine computes the steady-stateconditions analytically (for 1/2/3-cpt models) or viaexpansion (for ODE models) and uses them as the startingstate.SS = 2: the steady-state concentration is addedthe current compartment state (superposition).II: dosing interval (required when SS > 0).match the units of TIME. Rows with SS = 0 or missingSS ignore II.changes to the .ferx model file are needed. Steady-stateis purely data-driven and works for all PK model families(1-cpt, 2-cpt, 3-cpt oral/IV/infusion, ODE-based).”, “## anda dosing row as a steady-state dose by setting SS = 1 (orSS = 2) in the NONMEM CSV and providing the dosing intervalII (in the same time units as TIME).SS = 1: the subject is assumed to be at steadybefore this dose. The engine computes the steady-stateconditions analytically (for 1/2/3-cpt models) or viaexpansion (for ODE models) and uses them as the startingstate.SS = 2: the steady-state concentration is addedthe current compartment state (superposition).II: dosing interval (required when SS > 0).match the units of TIME. Rows with SS = 0 or missingSS ignore II.changes to the .ferx model file are needed. Steady-stateis purely data-driven and works for all PK model families(1-cpt, 2-cpt, 3-cpt oral/IV/infusion, ODE-based).”, “## list("II")a dosing row as a steady-state dose by setting SS = 1 (orSS = 2) in the NONMEM CSV and providing the dosing intervalII (in the same time units as TIME).SS = 1: the subject is assumed to be at steadybefore this dose. The engine computes the steady-stateconditions analytically (for 1/2/3-cpt models) or viaexpansion (for ODE models) and uses them as the startingstate.SS = 2: the steady-state concentration is addedthe current compartment state (superposition).II: dosing interval (required when SS > 0).match the units of TIME. Rows with SS = 0 or missingSS ignore II.changes to the .ferx model file are needed. Steady-stateis purely data-driven and works for all PK model families(1-cpt, 2-cpt, 3-cpt oral/IV/infusion, ODE-based).”, “## columns)a dosing row as a steady-state dose by setting SS = 1 (orSS = 2) in the NONMEM CSV and providing the dosing intervalII (in the same time units as TIME).SS = 1: the subject is assumed to be at steadybefore this dose. The engine computes the steady-stateconditions analytically (for 1/2/3-cpt models) or viaexpansion (for ODE models) and uses them as the startingstate.SS = 2: the steady-state concentration is addedthe current compartment state (superposition).II: dosing interval (required when SS > 0).match the units of TIME. Rows with SS = 0 or missingSS ignore II.changes to the .ferx model file are needed. Steady-stateis purely data-driven and works for all PK model families(1-cpt, 2-cpt, 3-cpt oral/IV/infusion, ODE-based).” )## Specifying model and data

There are two equivalent ways to supply the model and data: 1. Path strings (classic style):

ferx_fit("pk.ferx", "data.csv")
ferx_fit(model = "pk.ferx", data = "data.csv")

2. ferx_model object (pipe style):

"data.csv" |> ferx_model("pk.ferx") |> ferx_fit()

Both dispatch to the same Rust backend. The ferx_model form is convenient when combined with [ferx_model_set_section](ferx_model_set_section.qmd) to modify model options in the same chain (see Examples). The data path stored in the ferx_model object can always be overridden by supplying dataexplicitly to ferx_fit(). 3. Data declared in the model file: a .ferx file may carry a [data] block with path = <csv> (resolved relative to the model file’s directory). When data is omitted, ferx_fit() falls back to that path, so ferx_fit("pk.ferx") fits against the declared dataset. An explicit data argument always wins.

Controlling estimation

Estimation method - pass a single method or a vector to chain methods in sequence (each stage seeds the next with its converged parameters):

ferx_fit(m, d, method = "focei")
ferx_fit(m, d, method = c("saem", "focei"))  # SAEM warm-start, FOCEI polish

Standard errors - the covariance step is on by default:

ferx_fit(m, d, covariance = TRUE)   # default: produces SE / 
ferx_fit(m, d, covariance = FALSE)  # skip for speed during development

# Tune the FD step for the Hessian (default 0.01):
ferx_fit(m, d, fd_hessian_step = 0.1)    # increase when Hessian is ill-conditioned
ferx_fit(m, d, fd_hessian_step = 1e-3)   # decrease for very smooth surfaces

Parallelism - cap the Rust thread pool:

ferx_fit(m, d, threads = parallel::detectCores(logical = FALSE))

LOQ censoring:

ferx_fit(m, d, bloq_method = "m3")    # M3 likelihood (CENS: 1=LLOQ, -1=ULOQ)
ferx_fit(m, d, bloq_method = "drop")  # disable M3; use rows as observed

Gradient method for the inner EBE loop:

ferx_fit(m, d, gradient = "auto")  # default: engine-selected
ferx_fit(m, d, gradient = "fd")    # force finite differences

Optimizer trace - write a per-iteration CSV for convergence diagnostics:

fit <- ferx_fit(m, d, optimizer_trace = TRUE)
plot(fit)

c(“## Fine-tuning withsettings argument forwards a named list of low-level options toRust backend without requiring a new wrapper release. Arguments withparameters (e.g. method, covariance) cannot bein settings - pass them via the named argument.the complete settings reference - including which options apply to which, valid combinations, and a convergence troubleshooting guide - see(https://ferx-nlme.github.io/model-dsl/fit-options.html).Outer optimizer selection:\nferx_fit(m, d, settings = list(optimizer = \"auto\")) # default\nferx_fit(m, d, settings = list(optimizer = \"bobyqa\")) # derivative-free\nferx_fit(m, d, settings = list(optimizer = \"slsqp\")) # gradient-based\nferx_fit(m, d, settings = list(optimizer = \"lbfgs\"))\nferx_fit(m, d, settings = list(optimizer = \"bfgs\"))\nferx_fit(m, d, settings = list(optimizer = \"trust_region\"))\nferx_fit(m, d, settings = list(optimizer = \"mma\"))\nIteration cap and inner loop:\nferx_fit(m, d, settings = list(\n maxiter = 500L,\n inner_maxiter = 100L,\n inner_tol = 1e-6\n))\nSAEM tuning:\nferx_fit(m, d, method = \"saem\", settings = list(\n n_exploration = 200,\n n_convergence = 400,\n n_mh_steps = 3,\n omega_burnin = 20L,\n seed = 42L\n))\n\n# HMC proposals:\nferx_fit(m, d, method = \"saem\", settings = list(n_leapfrog = 3L))\nSIR uncertainty (requires sir = TRUE, covariance = TRUE):\nferx_fit(m, d, sir = TRUE, settings = list(\n sir_samples = 2000L,\n sir_resamples = 500L,\n sir_seed = 1L\n))\nGlobal search before local refinement:\nferx_fit(m, d, settings = list(\n global_search = TRUE,\n global_maxeval = 1000L\n))\nTrust-region CG budget:\nferx_fit(m, d, settings = list(\n optimizer = \"trust_region\",\n steihaug_max_iters = 100L\n))\nConvergence tolerance for difficult subjects:max_unconverged_frac (numeric in [0, 1], default 0.1) is theof subjects allowed to have unconverged EBEs before the outerrejects the trial step (scoring it as Inf). It is an-loop guard on the search, not a post-hoc relaxation of the converged flag,it never raises an error; set it to 1.0 to disable the guard. min_obs_for_convergence_check (non-negative integer) is thenumber of observations a subject must have before its inner-loopcounts toward that fraction; sparser subjects are excluded.\nferx_fit(m, d, settings = list(\n max_unconverged_frac = 0.1,\n min_obs_for_convergence_check = 2L\n))\nStagnation guard (NLopt-based optimizers):NLopt-based outer optimizers (BOBYQA, SLSQP, L-BFGS, MMA) short-circuit byonce recent evaluations show no OFV improvement above 1e-3 over awindow, letting them terminate quickly instead of burning throughmaxiter at full inner-loop cost. Set stagnation_guard = FALSEto disable this and let the optimizer run to its natural termination- useful when debugging or for problems with genuinely-but-real OFV improvements.\nferx_fit(m, d, settings = list(stagnation_guard = FALSE))\nMulti-start optimization:Runs n_starts independent fits from perturbed initial values in(via rayon) and returns the converged result with the lowest OFV.n_starts = 1 disables multi-start (single run, no overhead).for models prone to local minima: Michaelis-Menten elimination,-block omega, or many covariate parameters. On an 8-core machinen_starts = 8 costs approximately the same wall-clock time as a single run. always uses the exact initial values from the model file...n apply a log-space perturbation of size start_sigma(default 0.3) to each theta.\n# 8 parallel starts (approx. same wall time as one run on 8 cores)\nferx_fit(m, d, settings = list(n_starts = 8L))\n\n# Wider perturbation for ridge-shaped surfaces (e.g. Michaelis-Menten):\nferx_fit(m, d, settings = list(n_starts = 8L, start_sigma = 0.5))\n\n# Pin the perturbation RNG seed for reproducibility (independent of the SAEM seed):\nferx_fit(m, d, settings = list(n_starts = 8L, multi_start_seed = 123L))\n”, “## list("settings")settings argument forwards a named list of low-level options toRust backend without requiring a new wrapper release. Arguments withparameters (e.g. method, covariance) cannot bein settings - pass them via the named argument.the complete settings reference - including which options apply to which, valid combinations, and a convergence troubleshooting guide - see(https://ferx-nlme.github.io/model-dsl/fit-options.html).Outer optimizer selection:\nferx_fit(m, d, settings = list(optimizer = \"auto\")) # default\nferx_fit(m, d, settings = list(optimizer = \"bobyqa\")) # derivative-free\nferx_fit(m, d, settings = list(optimizer = \"slsqp\")) # gradient-based\nferx_fit(m, d, settings = list(optimizer = \"lbfgs\"))\nferx_fit(m, d, settings = list(optimizer = \"bfgs\"))\nferx_fit(m, d, settings = list(optimizer = \"trust_region\"))\nferx_fit(m, d, settings = list(optimizer = \"mma\"))\nIteration cap and inner loop:\nferx_fit(m, d, settings = list(\n maxiter = 500L,\n inner_maxiter = 100L,\n inner_tol = 1e-6\n))\nSAEM tuning:\nferx_fit(m, d, method = \"saem\", settings = list(\n n_exploration = 200,\n n_convergence = 400,\n n_mh_steps = 3,\n omega_burnin = 20L,\n seed = 42L\n))\n\n# HMC proposals:\nferx_fit(m, d, method = \"saem\", settings = list(n_leapfrog = 3L))\nSIR uncertainty (requires sir = TRUE, covariance = TRUE):\nferx_fit(m, d, sir = TRUE, settings = list(\n sir_samples = 2000L,\n sir_resamples = 500L,\n sir_seed = 1L\n))\nGlobal search before local refinement:\nferx_fit(m, d, settings = list(\n global_search = TRUE,\n global_maxeval = 1000L\n))\nTrust-region CG budget:\nferx_fit(m, d, settings = list(\n optimizer = \"trust_region\",\n steihaug_max_iters = 100L\n))\nConvergence tolerance for difficult subjects:max_unconverged_frac (numeric in [0, 1], default 0.1) is theof subjects allowed to have unconverged EBEs before the outerrejects the trial step (scoring it as Inf). It is an-loop guard on the search, not a post-hoc relaxation of the converged flag,it never raises an error; set it to 1.0 to disable the guard. min_obs_for_convergence_check (non-negative integer) is thenumber of observations a subject must have before its inner-loopcounts toward that fraction; sparser subjects are excluded.\nferx_fit(m, d, settings = list(\n max_unconverged_frac = 0.1,\n min_obs_for_convergence_check = 2L\n))\nStagnation guard (NLopt-based optimizers):NLopt-based outer optimizers (BOBYQA, SLSQP, L-BFGS, MMA) short-circuit byonce recent evaluations show no OFV improvement above 1e-3 over awindow, letting them terminate quickly instead of burning throughmaxiter at full inner-loop cost. Set stagnation_guard = FALSEto disable this and let the optimizer run to its natural termination- useful when debugging or for problems with genuinely-but-real OFV improvements.\nferx_fit(m, d, settings = list(stagnation_guard = FALSE))\nMulti-start optimization:Runs n_starts independent fits from perturbed initial values in(via rayon) and returns the converged result with the lowest OFV.n_starts = 1 disables multi-start (single run, no overhead).for models prone to local minima: Michaelis-Menten elimination,-block omega, or many covariate parameters. On an 8-core machinen_starts = 8 costs approximately the same wall-clock time as a single run. always uses the exact initial values from the model file...n apply a log-space perturbation of size start_sigma(default 0.3) to each theta.\n# 8 parallel starts (approx. same wall time as one run on 8 cores)\nferx_fit(m, d, settings = list(n_starts = 8L))\n\n# Wider perturbation for ridge-shaped surfaces (e.g. Michaelis-Menten):\nferx_fit(m, d, settings = list(n_starts = 8L, start_sigma = 0.5))\n\n# Pin the perturbation RNG seed for reproducibility (independent of the SAEM seed):\nferx_fit(m, d, settings = list(n_starts = 8L, multi_start_seed = 123L))\n” )## Post-fit outputs and pipe extensions

ferx_fit() returns a ferx_fit object. All of the following work both as standalone calls and as the tail of a |> pipe:

fit |> print()              # full parameter table
fit |> summary()            # compact diagnostic summary
fit |> ferx_model_inspect() # model structure auto-derived by the engine
fit |> plot()               # convergence trace (needs optimizer_trace = TRUE)

Diagnostics data frame (PRED, IPRED, CWRES, ETAs, …) lives in fit$sdtab; tidy estimates, correlations, and ETA-covariate correlations are plain fields (no accessor call needed):

fit$sdtab
fit$ebe_etas
fit$individual_estimates
fit$estimates    # tidy data frame: theta / omega / sigma + SE / 
fit$cor_matrix   # parameter correlation matrix (needs covariance = TRUE)
fit$eta_cov      # ETA-covariate correlation table

Seealso

• Fit options reference - full documentation of every [fit_options] key and settings knob, organised by method with a convergence-troubleshooting guide and a method x option compatibility table.
• [ferx_fit_async](ferx_fit_async.qmd) - non-blocking version for long runs.
• [ferx_check_init](ferx_check_init.qmd) - 5-iteration pilot to validate starting values before a full run.
• [ferx_inits_from_nca](ferx_inits_from_nca.qmd) - NCA-derived starting values.
• [ferx_get_warnings](ferx_get_warnings.qmd) - structured warnings from the fit.
• fit$estimates - tidy parameter table with SE /
  
• [plot.ferx_fit](plot.ferx_fit.qmd) - convergence trace plot.

Other fitting: [ferx_check_init](ferx_check_init.qmd), [ferx_collect](ferx_collect.qmd), [ferx_covariance](ferx_covariance.qmd), [ferx_inits_from_nca](ferx_inits_from_nca.qmd), [ferx_sir](ferx_sir.qmd), [ferx_stop](ferx_stop.qmd), [print.ferx_job](ferx_fit_async.qmd)

Concept

fitting

Value

A named list. Where to find commonly-needed quantities:

What you want	Accessor	Shape	When present
Residuals, PRED, IPRED, diagnostics	fit$sdtab	one row per observation	always
Covariates echoed per observation (via [output])	fit$sdtab	one row per observation, LOCF	declared in [output]
Raw covariate values for all dataset records	fit$covtab	one row per dataset record (doses + obs)	model has [covariates] block
ETA / EBE values per observation row	fit$sdtab (ETA_CL, ETA_V, …)	one row per observation, value repeated	always
ETA / EBE values per subject	fit$ebe_etas	one row per subject	model declares etas
Individual PK parameters per subject	fit$individual_estimates	one row per subject	always

Full slot descriptions:

Examples

ex <- ferx_example("warfarin")
fit <- ferx_fit(ex$model, ex$data, method = "gn", covariance = FALSE)
fit$theta
fit$ofv


ex <- ferx_example("warfarin")

# -- Classic style (path strings) ------------------------------------------

# Minimal call: FOCEI with default optimizer, covariance step on
fit <- ferx_fit(ex$model, ex$data)

# Named arguments (fully equivalent)
fit <- ferx_fit(model = ex$model, data = ex$data, method = "focei")

# All common options at once
fit <- ferx_fit(ex$model, ex$data,
  method      = "focei",
  covariance  = TRUE,
  threads     = 4L,
  gradient    = "auto",
  verbose     = TRUE,
  scale_params = FALSE
)

# SAEM warm-start, then FOCEI polish
fit <- ferx_fit(ex$model, ex$data, method = c("saem", "focei"))

# Fine-tune via settings
fit <- ferx_fit(ex$model, ex$data,
  method   = "focei",
  settings = list(
    optimizer     = "slsqp",
    maxiter       = 500L,
    inner_maxiter = 100L,
    inner_tol     = 1e-6
  )
)

# LOQ-censored observations (M3 method)
bloq <- ferx_example("warfarin_bloq")
fit  <- ferx_fit(bloq$model, bloq$data, method = "focei", bloq_method = "m3")

# SIR parameter uncertainty
fit <- ferx_fit(ex$model, ex$data,
  sir      = TRUE,
  settings = list(sir_samples = 2000L, sir_resamples = 500L)
)

# -- Pipe style (ferx_model object) -----------------------------------------

# Basic pipe: data flows into ferx_model() which bundles the model path
ex$data |>
  ferx_model(ex$model) |>
  ferx_fit(method = "focei", covariance = TRUE)

# Modify a model section, then fit
ex$data |>
  ferx_model(ex$model) |>
  ferx_model_set_section("fit_options", c(
    "  method     = focei",
    "  maxiter    = 500",
    "  covariance = true"
  )) |>
  ferx_fit()

# Full pipeline including post-fit outputs
ex$data |>
  ferx_model(ex$model) |>
  ferx_model_set_section("fit_options", c("  method = focei")) |>
  ferx_fit(covariance = TRUE, threads = 4L,
           settings = list(optimizer = "slsqp")) |>
  summary()

# Inspect a section, then fit
ferx_model_get_section(ex$model, "parameters")
ex$data |>
  ferx_model(ex$model) |>
  ferx_fit() |>
  (\(fit) fit$estimates)()

# Override data stored in ferx_model at fit time
# (substitute the path to your own dataset for "other_cohort.csv")
m <- ferx_model(ex$data, ex$model)
ferx_fit(m, data = "other_cohort.csv", method = "focei")

# -- Post-fit outputs -------------------------------------------------------

fit <- ferx_fit(ex$model, ex$data, covariance = TRUE, optimizer_trace = TRUE)

summary(fit)              # compact diagnostic table
ferx_model_inspect(fit)   # model structure auto-derived by the engine
plot(fit)                 # convergence plot (optimizer_trace = TRUE required)

fit$sdtab                 # per-observation diagnostics (PRED, IPRED, CWRES, etc.)
fit$ebe_etas              # per-subject empirical Bayes ETAs
fit$individual_estimates  # per-subject individual PK parameters
fit$estimates             # tidy data frame with SE and %RSE
fit$cor_matrix            # parameter correlation matrix
fit$eta_cov               # ETA-covariate correlation table
fit$eigenvalues           # sorted eigenvalues of parameter correlation matrix
fit$condition_number      # > 1000 flags potential ill-conditioning

# -- SDE model (Extended Kalman Filter) -------------------------------------

# The [diffusion] block in the .ferx file enables SDE mode.
# DIFF_CENTRAL is a within-subject process-noise variance on the central
# compartment.  SDE models use finite differences for this step.
sde <- ferx_example("warfarin_sde")
fit_sde <- ferx_fit(sde$model, sde$data)
fit_sde$uses_sde                  # TRUE
fit_sde$theta["DIFF_CENTRAL"]     # fitted diffusion variance
fit_sde$estimates                 # DIFF_CENTRAL appears in theta block

# -- Multi-start (avoid local minima) ----------------------------------------
ex <- ferx_example("warfarin")
fit_ms <- ferx_fit(ex$model, ex$data, settings = list(n_starts = 4L))
fit_ms$ofv

# -- NCA-based starting values (inits_from_nca) ------------------------------
# `inits_from_nca` is designed to mitigate stalling and local-minimum traps
# in **gradient-based estimation methods**, which navigate the likelihood
# surface using a local Jacobian / Hessian and so are sensitive to where the
# starting thetas sit. In ferx those are:
#   - method = "gn"         (FOCE Gauss-Newton)
#   - method = "gn_hybrid"  (FOCE-GN with SLSQP polish)
#   - method = "foce" / "focei" run with settings = list(optimizer =
#     "trust_region") or settings = list(optimizer = "slsqp") (the
#     default "auto" resolves to derivative-free "bobyqa" on these
#     poorly-started ODE/PD fits, which tolerates poor starts the best;
#     switch to slsqp / trust_region for gradient-based behaviour, then
#     NCA-derived starts help most)
# For these, NCA-derived starts often turn a non-converging or stagnating
# fit into a clean convergence.
#
# Stochastic / sampling-based methods ("saem", "imp") explore the space
# globally and are far less sensitive to starting values, so `inits_from_nca`
# gives a smaller benefit. The better tool for them is multi-start (run
# several fits from perturbed starts and keep the lowest OFV) - see the
# Multi-start section above (`settings = list(n_starts = 4L)`), or chain
# SAEM into FOCEI: `method = c("saem", "focei")`.

# TRUE = default strategy ("nca_sweep"); rescues a trust_region fit that
# would otherwise stall on poor defaults.
fit <- ferx_fit(ex$model, ex$data,
  settings       = list(optimizer = "trust_region"),
  inits_from_nca = TRUE
)

# Pick a strategy explicitly. "nca_ebe" handles large between-subject
# variability better than "nca_sweep" but falls back to it for ODE models.
fit <- ferx_fit(ex$model, ex$data, method = "gn",
                inits_from_nca = "nca_ebe")

# Fastest variant: NCA arithmetic only, no grid sweep.
fit <- ferx_fit(ex$model, ex$data, inits_from_nca = "nca")

# Same flag from the model file's [fit_options] (the call-time arg wins
# and warns on conflict). Inspect what NCA produces without fitting via
# `ferx_inits_from_nca()` -- see ?ferx_inits_from_nca.
#
# Combine with multi-start when even good inits aren't enough - good starts
# PLUS a small perturbed ensemble is often the most robust setup.
fit <- ferx_fit(ex$model, ex$data,
  method         = "gn",
  inits_from_nca = "nca_sweep",
  settings       = list(n_starts = 4L)
)

# -- Deep Compartment Model (covariate neural network) -----------------------
# Requires ferx-r built with the `nn` cargo feature.
# The [covariate_nn TYPICAL_PK] block replaces analytical covariate functions
# with a small MLP: WT + CRCL -> multiplicative modulator on TVCL/TVV1/etc.
# ferx_fit() call is identical to any other model.
dcm <- ferx_example("warfarin_dcm")
fit_dcm <- ferx_fit(dcm$model, dcm$data, method = "focei")

# NN metadata is in fit$neural_networks (one sub-list per [covariate_nn] block)
fit_dcm$neural_networks[[1]]$name          # "TYPICAL_PK"
fit_dcm$neural_networks[[1]]$shape         # e.g. c(2, 8, 8, 5)
fit_dcm$neural_networks[[1]]$n_weights     # total weight + bias count
fit_dcm$neural_networks[[1]]$input_names   # c("WT", "CRCL")
fit_dcm$neural_networks[[1]]$output_names  # c("CL", "V1", "Q", "V2", "KA")

# NN weight thetas are suppressed from the THETA table in print();
# a NEURAL NETWORKS summary block is shown instead.
print(fit_dcm)

# See inst/examples/ex_warfarin_dcm.R for a full worked example including
# interpretability heuristics and comparison against a no-covariate baseline.