Fit Options

The optional [fit_options] block configures the estimation method and optimizer settings.

Syntax

[fit_options]
  key = value

General Options

Key Values Default Description
method foce, focei, laplace, saem, gn, gn_hybrid, impmap, imp (only as final chain stage) focei Estimation method (or single-stage method in a chain). See chained fits. When neither this key nor a wrapper’s method argument is set, the engine falls back to focei and emits a warning so the implicit choice is visible — set method here to silence it.
threads positive integer, 0, auto auto Number of Rayon worker threads for per-subject parallelism. A positive integer pins the count. 0/auto (or omitting the key) uses the engine’s default: available cores − 1 (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). The CLI’s --threads flag mirrors this key (see ferx --help).
maxiter integer 500 Maximum outer loop iterations. maxiter = 0 is evaluation only (NONMEM MAXEVAL=0): the engine runs one inner EBE solve at the initial parameters and reports the objective there with no outer optimisation — useful for scoring a fixed parameter set or computing standard errors at known estimates (with covariance = true). The reported fit is marked not-converged, and the parameters are returned unchanged.
covariance true, false true Compute covariance matrix and standard errors
covariance_method r, s, rsr r Which covariance estimator to assemble, mirroring NONMEM $COVARIANCE MATRIX=. r = inverse Hessian R⁻¹ (model-based; the default). s = inverse score cross-product S⁻¹ (empirical information). rsr = the Huber–White sandwich R⁻¹SR⁻¹ (robust to model mis-specification; NONMEM’s default). ferx defaults to r, but NONMEM’s $COVARIANCE default is rsr — set covariance_method = rsr when reconciling SEs against a NONMEM run that used the default $COV. s/rsr are supported for FOCEI, FOCE, and IOV fits; all three anchored against NONMEM $COV MATRIX=S/RSR within ~10%. All methods require a positive-definite FD Hessian (a non-PD Hessian is not yet routed to the s cross-product). Requires covariance = true.
covariance_fallback none, sir none What to do when the FD Hessian is not positive definite (non-PD). none leaves the covariance step as failed. sir runs SIR with a fallback proposal covariance built from the \|eigenvalue\|-rectified Hessian, inflated 4×; the YAML reports SIR-based 95% CIs instead of Hessian SEs. Requires covariance = true. Ignored when the Hessian succeeds. Setting sir = true arms the same fallback on its own, so this key is only needed to get the fallback without requesting SIR generally.
analytic_cov_hessian true, false true Use the exact analytic R-matrix for the covariance step when the model is in scope, instead of the finite-difference-of-OFV Hessian (#436). The analytic route assembles the observed information from third-order sensitivities of the closed-form prediction, so it carries no fd_hessian_step to tune and no ε/h² differencing noise, and it costs 2·(n_theta + n_eta) + 1 sensitivity evaluations per subject rather than ~2·n_free² objective evaluations that each re-solve every subject’s inner loop. Scope is deliberately narrow — plain analytical (closed-form) Gaussian models, no IOV, no LTBS, no [scaling], no Form-C readout, no [initial_conditions], no iiv_on_ruv, no correlated (block_sigma) or custom-magnitude residual, no M3 censoring, no FREM, no covariate-Selected error spec, no non-Gaussian (TTE/categorical/Markov) endpoint, no method = laplace/agq (whose marginal this assembly does not describe), and not gradient = fd. Anything outside that scope, including a single out-of-scope subject, silently keeps the finite-difference stencil, which is correct for all of them — the choice is all-or-nothing across the population, because a Hessian half-assembled from two approximations would be neither. Set analytic_cov_hessian = false to force finite differences for an in-scope model (e.g. to reproduce standard errors from a pre-#436 run).
fd_hessian_step float 1e-2 Initial relative step size for the finite-difference Hessian used in the covariance step. Ignored when the analytic R-matrix (analytic_cov_hessian) serves the fit, which differences nothing at the population level. The actual perturbation for parameter i is fd_hessian_step * (1 + |x_hat[i]|). ferx automatically halves this value up to 8 times if any diagonal stencil evaluates to a non-finite value, so in most cases the default requires no tuning. Set explicitly only if the automatic reduction is insufficient (try 0.1) or if FD noise is the main concern on a very smooth surface (try 1e-3). Must be positive and finite.
optimizer auto, bobyqa, slsqp, nlopt_lbfgs, mma, trust_region (plus deprecated aliases bfgs/lbfgs → nlopt_lbfgs) auto Outer-loop optimisation algorithm. Applies to method = foce / focei and to the FOCEI polish phase of method = gn_hybrid. The default auto picks per model: nlopt_lbfgs when the exact analytic FOCE/FOCEI gradient is available, and bobyqa when only finite differences are (ODE/PD models, LTBS/SDE, or gradient = fd) — see the Outer Optimizers page for the applicability matrix and the benchmarking (#490) behind these defaults. The fit output records the resolved pick (e.g. auto (nlopt_lbfgs)). Ignored (with a runtime warning) under method = saem / imp / gn. bfgs and lbfgs are deprecated aliases that now select nlopt_lbfgs (the hand-rolled built-in BFGS/L-BFGS is slated for removal — see #483).
outer_xtol float 1e-4 Relative step tolerance: how far the optimizer’s step must shrink before it declares success. 1e-4 in scaled log-space is a ~0.01% move on the natural parameter scale. Affects bobyqa (the auto pick for FD-objective models — ODE/PD, LTBS/SDE, gradient = fd, and TTE), where it is NLopt’s xtol_rel, and trust_region, where it is the step limb of the no-progress stopping rule (#1000; see Trust-region convergence reporting). The remaining gradient-based NLopt optimizers use their own fixed internal tolerances. Must be positive and finite.
outer_ftol float auto Relative objective tolerance: the relative OFV change below which the optimizer stops. Applies to bobyqa (NLopt ftol_rel) and to trust_region, where it is the objective limb of the no-progress stopping rule (#1000; see Trust-region convergence reporting) — both resolve it the same way. Unset = auto: 1e-8 for a pure-TTE model and 1e-6 otherwise. The TTE tightening fixes #469 — on the near-flat frailty-ω² ridge of a nonlinear hazard parameter (the whole ridge spans <0.01 OFV) the historical 1e-6 stopped the optimizer short, leaving the variance above its own minimum (a Weibull shape-frailty read ω² 0.204 vs the NONMEM/nlmixr2 0.175 consensus; at 1e-8 it lands 0.176). The 1e-6 floor for non-TTE fits is deliberate: a TTE hazard objective is evaluated exactly, but on a noisy objective (ODE solver error, FD-inner FOCE such as LTBS) 1e-8 is unreachable, so BOBYQA grinds toward its maxeval budget instead of converging (≈3× the evaluations on an ODE fit). Set an explicit value to pin it for every model — e.g. outer_ftol = 1e-8 to tighten an LTBS flat-ridge fit, accepting the cost. Must be positive and finite.
inner_maxiter integer 200 Max iterations for the inner (per-subject EBE) optimizer
inner_optimizer auto, bfgs, lbfgs, nelder_mead auto Inner-loop (per-subject EBE) optimisation algorithm. auto uses dense BFGS, switching to limited-memory L-BFGS above 32 random effects. The explicit values pin a single algorithm with no automatic switching: bfgs (dense), lbfgs (limited-memory two-loop recursion — cheaper per step at high random-effect dimension), nelder_mead (derivative-free, for diagnosing gradient issues). All gradient-based variants converge to the same EBE; the choice trades per-step cost against memory. Applies to every method with an inner loop (foce/focei/gn/gn_hybrid/saem/impmap).
inner_tol float 1e-5 Gradient-norm convergence tolerance for the inner (per-subject EBE) optimizer. A looser tolerance leaves residual noise in each subject’s EBE solution, which propagates into the marginal objective the outer optimizer sees; on models with a noisy or flat marginal surface (FD-inner FOCE such as log-transform-both-sides) that noise can make the derivative-free BOBYQA outer optimizer false-converge above the true minimum. 1e-5 removes enough of that noise to match NONMEM’s minimum, at roughly 1.5× the per-fit cost of the older 1e-4 default. (For reference, NONMEM runs its inner conditional optimization far tighter than its outer convergence: the inner precision SIGL defaults to ~10 significant digits, versus the outer NSIG/SIGDIGITS default of 3.) 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. Loosen it (e.g. 1e-4) to speed up well-conditioned fits where it makes no difference; change it further only with validation against a reference fit.
cov_inner_tol float inner_tol (LTBS: min(inner_tol, 1e-8)) Inner EBE-reconvergence tolerance used only by the covariance step, decoupled from inner_tol. The covariance R-matrix is a second-difference of the reconverged OFV and is more sensitive to EBE precision than the fit itself, so on a flat or weakly-identified surface an EBE converged only to the fit’s inner_tol can leave enough residual noise to perturb the standard errors (e.g. the heavily-censored M3 + IOV case in #654 — set cov_inner_tol = 1e-11). This lets the covariance step reconverge tighter without slowing every outer iteration. Unset (default) uses inner_tol for ordinary models (SEs unchanged) and min(inner_tol, 1e-8) for closed-form LTBS models, whose g = ln(f) covariance Hessian needs the tighter reconvergence — including LTBS combined with [iov] since #486, when that combination’s inner loop moved onto the same analytic ln(f) gradient. The default is deliberately not tightened for every model — over-converging some ill-conditioned inner Hessians drives the covariance indefinite — so opt in per model with this key.
inner_restarts integer 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 high-volume/low-concentration fit and a low-volume/high-concentration fit can both explain a subject’s total-concentration data, so a single warm-started inner BFGS can settle in the wrong (shallower) basin, inflating that subject’s objective and biasing its EBE. With inner_restarts = N > 0 (default 1), a subject re-solves the EBE on a cold start from N Ω-scaled alternate seeds (±2·sd, ±3·sd, …) and keeps the lowest-objective mode; the outer loop’s warm start then carries the chosen basin through the rest of the fit, so the extra scan runs once per subject per fit (≈0 % overhead). Two situations are scanned. Subjects on the event-driven path (system resets EVID=3/4 or time-varying covariates — the saturable-binding case above) re-seed every random effect. Every other subject is checked, per random effect, for a weakly-identified coordinate — one whose individual objective is flat, i.e. the data adds less curvature than the prior already carries (high per-subject shrinkage). This is exactly the condition under which a distant, lower posterior mode can hide from a single BFGS descent (issue #891; the poorly-identified V1 in a saturable-clearance fluconazole model is the motivating case), and only those flagged coordinates are re-seeded. The flatness check is a two-point finite difference of the objective per coordinate, so a well-identified subject pays only that probe and its EBE is unchanged. An alternate seed that reconverges to the same mode is not accepted, so unimodal subjects are bit-identical to inner_restarts = 0; only a genuinely trapped subject changes (it moves to the deeper, correct basin). Applies to the gradient-based methods that share this inner EBE solve (focei, foce, laplace, gn, gn_hybrid, trust_region); the covariance step reconverges with the same setting so standard errors stay consistent with the fit. Set inner_restarts = 0 to disable; raise it above 1 only if a wider seed ladder is needed.
steihaug_max_iters integer adaptive Max CG iterations for the Steihaug subproblem (only used when optimizer = trust_region). Default (unset) uses ceil(sqrt(n_params)).clamp(5, n_params) — typically 5 for standard NLME models. Set explicitly (e.g. steihaug_max_iters = 50) to pin the budget.
ode_reltol float 1e-4 RK45 ODE solver relative tolerance. Applies to ODE models, and to the auto-generated ODE twin that closed-form transit / inverse-Gaussian absorption models (one_cpt_transit, two_cpt_transit, one_cpt_ig, two_cpt_ig) use to serve IOV, time-varying-covariate, and TIME-switch subjects (#719/#814) — so a call-time override reaches those rerouted fits too; otherwise ignored for analytical PK. The default reproduces analytical closed forms in PRED to ~1e-4, but the FOCE objective amplifies solver error, so the OFV of an ODE-form model can differ from its analytical equivalent by several units. Set a tighter value (e.g. 1e-10) when the ODE-form OFV must match an analytical reference; expect slower fits.
ode_abstol float 1e-6 RK45 ODE solver absolute tolerance (companion to ode_reltol). ODE models and the closed-form absorption ODE twin (as for ode_reltol).
ode_max_steps integer 10000 Max solver steps per integration segment. ODE models and the absorption ODE twin. Raise (e.g. 1000000) if a tight ode_reltol exhausts the step budget on stiff multi-compartment systems.
ode_method rk45, vern7, rosenbrock23, rodas4, rodas5p rk45 Which stepper integrates the [odes] block. Two axes: vern7 (Verner 7(6), 10 stages) buys order for fits that are accuracy-limited at tight tolerances — measured 2.8× fewer steps and ~2.3× faster than rk45 on the Savic transit NONMEM anchor at 1e-9, but ~1.4× slower at default tolerances, so it is not a blanket upgrade. The Rosenbrock methods buy stability. rk45 (default) is the explicit Dormand-Prince method and is the fastest choice on non-stiff models. The other three are linearly implicit Rosenbrock methods for stiff systems — fast reversible binding / TMDD, Michaelis-Menten with KM far below observed concentrations, long transit chains, QSP cascades — where an explicit method is stability-limited and grinds down to the minimum step (visible as a fit that crawls, or as ode_max_steps being exhausted). Pick by tolerance: rosenbrock23 (order 2, aliases ros23/ode23s) for crude tolerances or rough right-hand sides, rodas4 (order 4) at typical 1e-6–1e-9, rodas5p (order 5) at ode_reltol ≤ 1e-9. A stiff step costs n + 1 extra right-hand-side evaluations for the finite-difference Jacobian plus one n × n factorization, so this is not a free upgrade — switch only when stiffness, not accuracy, caps the step. Read n as the size of the system actually being integrated, which is not always the [odes] state count: a CTMM endpoint integrates the occupancy system dP/dt = P·Q(t), so an s-state chain is an s²-state integration and a Rosenbrock method there factors an s² × s² matrix per step attempt (a 5-state chain → a 25×25 LU). ode_method reaches that solve too, which is what you want when the chain has widely separated transition rates and is genuinely stiff, and is pure overhead otherwise. The stiff methods are full peers of rk45 — every feature works with every method, including analytic sensitivities, non-Gaussian endpoints (TTE / categorical / CTMM), time-to-event simulation and adaptive dosing — because each method carries its own continuous extension for in-step state readouts. See ODE Models.
global_search true, false false Run NLopt CRS2-LM (Controlled Random Search with Local Mutation) as a gradient-free global pre-search before the local optimizer. CRS2-LM samples within the parameter bounds; the local optimizer (e.g. bobyqa, slsqp) starts from the best point found. Useful for poorly-identified models — when the local optimizer can land in a degenerate basin (collapsed ETA, V/Q swap, parameters at bounds) from a far-from-truth start, the global pre-search usually escapes it. Adds the pre-search budget on top of the local optimisation, but typically more efficient than running multiple full fits from scratch. Requires a full NLopt build (e.g. brew install nlopt or apt install libnlopt-dev); a clear warning is emitted if CRS2-LM is unavailable.
global_maxeval integer 200 * (n_params + 1) Maximum evaluations of the FOCE objective during the global pre-search. Each eval is a full inner-loop pass over all subjects, so this is the dominant cost of global_search = true. The default (0 → auto) is empirically enough to escape bad basins on 10–20 parameter PK models without dominating the wall time of the subsequent local refine.
bloq_method drop, m3 drop How to handle rows with CENS=1 (below LLOQ) or CENS=-1 (above ULOQ). m3 enables Beal’s M3 likelihood (see BLOQ example).
npde_nsim integer ≥ 0 0 Number of Monte-Carlo replicates per subject used to compute the simulation-based NPDE/NPD diagnostics after the fit. 0 (default) disables the computation — no NPDE/NPD columns are emitted. A typical value is 1000. Cost scales linearly with the count.
npde_seed integer — RNG seed for the NPDE/NPD simulation, for reproducible diagnostics. Unset falls back to a fixed default. Only used when npde_nsim > 0.
mu_referencing true, false true Re-centre inner-loop ETA estimates on the current population mean (auto-detected from [individual_parameters]). See the FAQ entry for details. Set false to reproduce pre-automatic-mu behaviour. Under a [mixture] model this also switches off the class-aware mu-referencing θ shift for SAEM and estimating IMP/IMPMAP (see Mixture models).
iov_column string — Name of the occasion column in the dataset (e.g. OCC). Supplies occasions when the model uses kappa or block_kappa declarations. The column must contain integer occasion indices. Case-insensitive. Supported with foce, focei, and saem (not imp / impmap, which do not support IOV). See IOV documentation.
iov_occasion column, dose, time(t₁, …) column Derive the IOV occasion partition from the model instead of a data column. dose: each administration starts a new occasion. time(24, 48): time-window breakpoints. When set to dose/time(...), this overrides iov_column (with a warning). See Declaring occasions in the model.
optimizer_trace true, false false Write a per-iteration CSV to /tmp/ferx_trace_<pid>_<ts>.csv. The path is stored in FitResult::trace_path. Useful for diagnosing convergence problems or comparing optimizers. See Optimizer Trace.
checkpoint true, false true Periodically save a resume point so a run that is interrupted (e.g. Ctrl-C, a crashed job) can be restarted from where it stopped instead of from scratch (#755). The CLI writes {model}.tmp next to the other outputs and, on the next run of the same model + data, resumes from it automatically; the file is deleted on successful completion. The write is throttled to checkpoint_interval_secs, so a run that finishes before the first interval leaves nothing behind. Restart is coarse: the population estimates (θ/Ω/Σ) at the last save are used as the new starting point and the fit continues — not a bit-exact optimizer-state resume — so the resumed run may take a few extra iterations to re-converge. A change to the model or data invalidates the checkpoint (a hash check), and the run silently starts fresh. Pass the CLI flag --clean to force a fresh start (it deletes any existing {model}.tmp first). Set checkpoint = false to disable saving entirely.
checkpoint_interval_secs integer 300 Minimum wall-clock seconds between checkpoint writes. Larger values reduce I/O overhead but lose more progress on an interrupt; smaller values save more often. Only used when checkpoint = true.
reconverge_gradient_interval integer ≥ 0 0 How often to re-solve each subject’s inner EBE loop during the population gradient instead of holding the EBEs (η̂) and FOCE Hessian fixed. 0 (default) never reconverges — the cheap fixed-EBE gradient is used. 1 reconverges on every gradient evaluation; N reconverges on evaluations 0, N, 2N, … and uses the cheap gradient in between (amortizing the cost by ~N× while still periodically correcting the search direction). The fixed-EBE gradient omits the inner solution’s response to θ/Ω, so on ill-conditioned fits a gradient optimizer (e.g. slsqp) can stall well above the derivative-free (bobyqa) optimum; reconverging recovers the full surface at roughly 5–6× the per-gradient cost at interval 1. IOV models (kappa/block_kappa) always reconverge and ignore this setting. See the FOCE/FOCEI page.
inits_from_nca true, false, nca, nca_sweep, nca_ebe false Derive NCA-based starting values from the data before the optimizer loop. true (alias for nca_sweep) and false/off toggle the default strategy; the three named values pick a strategy explicitly (see NCA-based starting values). Fixed thetas are never overwritten; covariate-effect thetas (no mu-referencing link) keep the model default. Most useful with method = trust_region or method = gn / method = gn_hybrid, where bad starting values can cause early stalling. The same estimation is available without running a fit via the CLI flag --inits-from-nca[=METHOD] and (in ferx-r) the ferx_inits_from_nca() function.

NCA-based starting values

inits_from_nca estimates starting values directly from the data using non-compartmental analysis (NCA), then optionally refines parameters NCA cannot estimate. All three strategies are NCA-based; they differ only in how much refinement runs on top:

Value Strategy What it does Typical cost
nca NCA only Per-subject NCA (AUC, terminal slope, Wagner–Nelson Ka, biexponential peeling for 2/3-cpt) pooled to population geometric means. Leaves parameters NCA can’t reach (peripheral Q/V2, all ODE/PD thetas) at the model default. < 5 ms
nca_sweep NCA + sweep Runs nca, then sweeps every remaining non-fixed theta over a log-space rRMSE grid using population predictions (etas = 0). Model-agnostic — also covers ODE/PD models. This is what true selects. < 500 ms (analytical)
nca_ebe NCA + EBE sweep Like nca_sweep but evaluates the grid with empirical Bayes estimates (etas ≠ 0); more accurate under large IIV (omega > ~0.2). Falls back to nca_sweep for ODE models. < 500 ms (analytical)

The CL eta’s omega is also updated from inter-subject CV² when ≥ 3 subjects have a valid NCA estimate.

When nca_sweep is enabled but the fit fails to converge or the OFV looks suspiciously high, try nca_ebe.

Simulation-based diagnostics (NPDE / NPD)

Setting npde_nsim > 0 adds two simulation-based goodness-of-fit columns to the sdtab output: NPD (Normalized Prediction Discrepancies) and NPDE (Normalized Prediction Distribution Errors). Unlike CWRES — which linearises the model around the conditional mode and so inherits the bias of a first-order approximation — NPDE/NPD are built entirely from Monte-Carlo simulation under the fitted model, so they are robust to model nonlinearity and to non-Gaussian random effects (Brendel et al. 2006; Comets et al. 2008, the npde R package).

[fit_options]
method = focei
npde_nsim = 1000     # replicates per subject (off when 0, the default)
npde_seed = 12345    # optional, for reproducibility

The effective seed actually used — the explicit npde_seed when given, otherwise the built-in default — is recorded as model: npde_seed: in the {model}-fit.yaml output (and round-trips through .fitrx), so the NPDE/NPD columns can always be regenerated from the saved fit, even when no seed was set in the model file. The line is omitted when NPDE did not run (npde_nsim = 0).

For each observation y_ij, ferx simulates K = npde_nsim replicates under the fitted θ/Ω/Σ (sampling η ~ N(0, Ω) and the residual error), evaluated at the subject’s own observation design:

  • NPD (no decorrelation): pd_ij is the empirical CDF of the simulated values at observation ij, and npd_ij = Φ⁻¹(pd_ij).
  • NPDE (decorrelated): within each subject the observed and simulated vectors are decorrelated with the empirical mean and Cholesky factor of the simulated covariance (the Brendel/Comets procedure) before the empirical CDF and inverse- normal transform. This removes the within-subject correlation that NPD ignores.

Empirical-CDF probabilities at 0 or 1 are clamped to [1/(2K), 1 − 1/(2K)] before Φ⁻¹, matching the npde-package convention, so the scores stay finite.

Under a correctly specified model, NPDE follows a standard normal distribution. On the warfarin example (examples/warfarin.ferx, data/warfarin.csv), a converged FOCEI fit with npde_nsim = 1000 gives the expected mean-zero, unit-variance distribution:

Diagnostic mean variance reference
NPDE 0.006 1.02 ≈ N(0, 1)
NPD −0.002 0.88 mean ≈ 0; variance < 1 because NPD retains within-subject correlation

This matches the npde R package’s autonpde() (same Cholesky decorrelation and edge-clamping) to Monte-Carlo noise. See tests/npde_validation.rs for the engine-side check and an R reproduction recipe.

NONMEM cross-validation

NONMEM has no native NPDE output; the reference pipeline is a NONMEM $SIMULATION post-processed by the npde R package. The same warfarin data was fit with NONMEM 7.5.1 (ADVAN2 TRANS2, METHOD=1 INTER, proportional error) and the two fits agree, so any difference in the diagnostics is attributable to the NPDE computation rather than the fit:

Quantity ferx (FOCEI) NONMEM 7.5.1
TVCL 0.1328 0.1327
TVV 7.738 7.738
TVKA 0.817 0.811
OFV −286.0015 −286.0042

Feeding a 1000-replicate NONMEM $SIMULATION (under the NONMEM estimates) and the observed data to npde::autonpde() (Cholesky decorrelation, edge-clamping) reproduces ferx’s NPDE/NPD on the same 110 observations:

Diagnostic ferx mean ferx var npde-pkg mean npde-pkg var
NPDE 0.005 1.09 0.036 1.04
NPD −0.009 0.91 −0.004 0.91

The NPD variance (the quantity that should sit below 1 because NPD retains the within-subject correlation) matches to three figures (0.91 vs 0.91); the small NPDE-variance gap is Monte-Carlo noise between two independent simulation draws. Both engines give NPDE ≈ N(0, 1) and a mean-zero NPD, as expected under this well-specified model.

Limitations. M3/BLQ censored observations need the predictive-CDF variant and are out of scope: censored rows (CENS != 0) are emitted as empty/NaN (matching the CWRES/IWRES convention), and a subject’s NPDE is NaN as a whole when it has any censored row, since the within-subject decorrelation would otherwise mix the LLOQ into the uncensored rows. NPDE also requires more replicates than a subject has observations (npde_nsim > n_obs) for a full-rank simulated covariance; otherwise that subject’s NPDE is NaN (NPD is still computed). For IOV (kappa) models the NPDE/NPD simulation currently holds kappas at zero, so NPDE/NPD omit the inter-occasion component of the predictive variance. (This is specific to the NPDE path; simulate() itself samples per-occasion kappas — see Simulation. Aligning NPDE is a follow-up.)

Estimation Methods

FOCEI (default)

method = focei

FOCE with Interaction. Includes the dependence of the residual error on random effects. More accurate than FOCE when the error model depends on individual predictions, but slightly slower.

FOCE

method = foce

First-Order Conditional Estimation. Linearizes the model around the empirical Bayes estimates. Fast and reliable for most models.

Laplace

method = laplace

The Laplace approximation with the exact Hessian — NONMEM $EST METHOD=1 LAPLACIAN (alias laplacian).

This is not the same estimator as focei. Both approximate the marginal likelihood with a single Gaussian at the empirical-Bayes mode, but FOCEI builds that Gaussian from the Gauss-Newton Hessian CᵀC + Ω⁻¹ (which drops ∂²f/∂η²), while laplace differentiates the true conditional likelihood. It therefore carries the curvature of the η-dependent residual variance that the Gauss-Newton form discards — which is why it reproduces NONMEM’s LAPLACIAN to six significant figures where FOCEI lands on a different value.

At the default n_agq = 1 it is a single Gaussian at the mode; the n_agq argument generalizes it to adaptive Gauss-Hermite quadrature (below). Same objective, same analytic gradient, same covariance step. It is the cheapest configuration (one node, no grid) and on warfarin it is faster than FOCEI.

Use n_agq = 1 when you want a NONMEM-LAPLACIAN-comparable fit, or a Laplace baseline to raise n_agq from. See the AGQ page.

n_agq — adaptive Gauss-Hermite quadrature

method = laplace
n_agq  = 3

n_agq is an argument, not a separate method (there is no method = agq). It sets the Gauss-Hermite nodes per random effect. n_agq = 1 (default) is the classical single-point method; n_agq > 1 evaluates the exact conditional likelihood on a grid around the empirical-Bayes mode, refining the marginal toward the exact integral. It makes no Gaussian-residual assumption, so it is the deterministic option for non-Gaussian endpoints (TTE, categorical). Cost is n_agq ^ n_eta per subject per iteration, so it suits models with few random effects (n_eta ≤ 4).

The method name selects the Hessian anchor: method = laplace places the grid with the exact Hessian, method = focei with the Gauss-Newton H̃. So method = focei with n_agq > 1 is the Gauss-Newton-anchored quadrature — FOCEI refined toward the exact marginal. See the AGQ page.

SAEM

method = saem

Stochastic Approximation EM. Uses Metropolis-Hastings sampling instead of MAP optimization for random effects. More robust to local minima; recommended for complex models with many random effects.

AGQ-Specific Options

Key Default Description
n_agq 1 Gauss-Hermite nodes per random effect. 1 reproduces the Laplace approximation exactly; odd values are conventional (they keep a node at the mode). Max 21. The tensor grid costs n_agq ^ n_eta likelihood evaluations per subject per outer iteration, and a fit whose grid exceeds 100 000 nodes is rejected at check time — lower n_agq, use fewer random effects, or switch to saem / imp.

SAEM-Specific Options

Key Default Description
n_exploration 150 Phase 1 iterations (step size = 1)
n_convergence 250 Phase 2 iterations (step size = 1/k)
n_mh_steps 20 Block Metropolis-Hastings steps per subject per iteration. Also sizes the componentwise decorrelating kernel that prevents block-Ω collapse (max(2, n_mh_steps / n_eta) sweeps; multi-η models only — skipped when n_eta < 2). When n_leapfrog > 0, this applies to subjects that fall back to MH (see below); HMC subjects use one proposal per iteration regardless.
n_leapfrog 0 Leapfrog steps per HMC proposal (0 = use MH; see below). When > 0, subjects for which HMC is unavailable (ODE model, missing analytical PK path, non-finite Ω, unsupported TV-cov path) fall back to MH using n_mh_steps proposals.
adapt_interval 50 Iterations between step-size adaptation
omega_burnin 20 Initial exploration iterations during which Ω (and ΩIOV) are held at their starting values while the MH chain warms up. Clamped to n_exploration; set 0 to disable. Prevents the Ω collapse described in the SAEM page.
mstep_damping 0.03 Exploration-phase cap on the stochastic-approximation step for the numerical θ/σ M-step, the θ-side counterpart of the per-iteration SA step cap SAEM already applies to Ω. The M-step result is blended in as θ ← θ + γ_θ·(θ* − θ) rather than assigned, because assigning it outright is argmax of a single MCMC η draw and biases a theta that carries no ETA. Sigma comes out of the same joint NLopt solve and takes the same step. Must be in (0, 1]; smaller damps harder. 1.0 is a sentinel rather than a cap value: it disables the damping entirely — in the convergence phase as well as exploration, restoring the pre-damping outright assignment throughout — so the option is discontinuous there (0.999 still gives the full convergence schedule). Ignored (with a warning naming the reason) when every estimated theta is mu-referenced or FIX, when mu_referencing = false, and for mixture models. See the SAEM page.
seed 12345 RNG seed for reproducibility
conddist false Run a post-fit conditional-distribution pass estimating each subject’s p(η_i \| y_i) by MCMC. Reports per-subject conditional mean/SD and distribution-based η-shrinkage (the shrinkage-unbiased basis for η diagnostics). See the SAEM page.
conddist_nsamp 200 Retained MCMC draws per subject in the conditional-distribution pass. Larger tightens the mean/SD estimates at linear cost. Production diagnostics want values in the thousands. Only used when conddist = true.
conddist_burnin 20 Burn-in draws discarded before accumulation, to forget the EBE-mode warm start. Only used when conddist = true.
conddist_keep_samples false Retain the raw per-subject draws (written to {model}-conddist-samples.csv), not just the mean/SD. Only used when conddist = true.

SIR (Sampling Importance Resampling)

SIR provides non-parametric parameter uncertainty estimates as an optional post-estimation step. Requires covariance = true.

Key Default Description
sir false Enable SIR uncertainty estimation
sir_samples 1000 Number of proposal samples (M)
sir_resamples 250 Number of resampled vectors (m)
sir_seed 12345 RNG seed for reproducibility
sir_keep_samples false Retain resampled parameter vectors for simulate_with_uncertainty()
sir_df 5.0 Degrees of freedom for the Student-t proposal; higher values approach a normal proposal

When the covariance step fails

The standard SIR path draws from the inverted FD Hessian, so it needs a successful covariance step. When that step fails because the Hessian is not positive definite, sir = true automatically falls back to the same |eigenvalue|-rectified proposal that covariance_fallback = sir uses, and the fit reports SIR-based 95% CIs (covariance_status: sir_fallback) instead of Hessian SEs. Setting covariance_fallback = sir explicitly is therefore only needed when you want the fallback without asking for SIR on a successful covariance step.

SIR is skipped entirely — with a warning saying so — when the covariance step produces neither a covariance matrix nor a fallback proposal. That happens when the step was never run (covariance = false, or a Bayesian fit, which reports posterior credible intervals instead), and when it ran but failed for a reason other than a non-PD Hessian: a divergent eigendecomposition (NaN/Inf entries), a flat or non-finite FD stencil, a non-finite objective at the estimates, or a singular score cross-product under covariance_method = s/rsr. The warning names the covariance-step message that carries the actual cause. A cancelled fit skips SIR without an additional warning.

See SIR documentation for details.

Importance Sampling (IMP)

By default imp is a Monte-Carlo EM estimator (NONMEM METHOD=IMP): it maximises the importance-sampled marginal likelihood, updating θ/Ω/σ each iteration. Set imp_eval_only = true (NONMEM EONLY=1) to instead evaluate the marginal −2 log L at the fixed input parameters — a lower-bias −2 log L than the FOCE/Laplace OFV when subject posteriors of η are non-Gaussian (e.g. sparsely-sampled PK).

Behaviour change: imp previously only evaluated −2 log L. It now estimates by default. Add imp_eval_only = true to recover the old behaviour.

[fit_options]
  method        = imp            # estimate (NONMEM METHOD=IMP)
  imp_iterations = 200
  imp_samples    = 1000
  imp_averaging  = 50
  imp_proposal_df = 5             # or `normal` for a multivariate-normal proposal
  imp_seed       = 12345
[fit_options]
  method        = [focei, imp]   # evaluate FOCEI's fit (NONMEM EONLY=1)
  imp_eval_only  = true

On rich data prefer impmap or warm-start with [focei, imp]; plain imp’s one-iteration-lagged proposal can collapse the ESS on a sharp posterior. See the IMP documentation.

Key Default Description
imp_eval_only false true ⇒ evaluate −2 log L at fixed parameters (NONMEM EONLY=1); must be the terminal chain stage. false ⇒ estimate (NONMEM METHOD=IMP).
imp_iterations 200 MCEM iterations (estimator only).
imp_averaging 50 Terminal iterations averaged into the reported estimate (estimator only).
imp_samples 1000 Importance samples K per subject. 2000–5000 recommended for publication-quality MC SE.
imp_proposal_df 5.0 Student-t proposal degrees of freedom (≥ 1), or normal/mvn for a multivariate-normal proposal. Lower = heavier tails.
imp_auto true Adaptive sample count (NONMEM AUTO). When true, imp_samples is the starting count and is ramped up (×2/iteration, cap 10000) while the objective’s Monte-Carlo SE exceeds 1.0. Recommended for high-dimensional / FREM models, where a fixed count biases the M-step.
imp_seed 12345 RNG seed. Same seed → identical result.
imp_low_ess_threshold 0.1 Subjects with normalized ESS below this fraction get flagged in the result. Set 0 to silence.
imp_defensive_alpha 0.0 Defensive-mixture weight (issue #528), opt-in. Each subject draws this fraction of its samples from the prior N(0, Ω) rather than the mode-centred proposal, and every sample is scored under the mixture density. This bounds the importance weights so a weakly-identified subject — e.g. an analytical [initial_conditions] baseline whose V cancels in the amplitude — can’t hijack the weighted M-step and walk θ to the bounds. Must be in [0, 1); the default 0 is the legacy single-proposal sampler (bit-comparable with NONMEM), set a small positive value (e.g. 0.1) to enable the rescue. Applies to imp and impmap (including the FREM Rao-Blackwell path); for an impmap stage you may also write impmap_defensive_alpha. Enabling it disables Sobol QMC and raises the per-subject ESS floor (so imp_low_ess_threshold flags fewer subjects). See Importance Sampling → Defensive mixture.

See Importance Sampling documentation for the algorithm, the NONMEM mapping, IOV caveats, and tuning guidance.

IMPMAP (importance_sampling_map)

Like the estimating imp, impmap is a Monte-Carlo EM estimator — equivalent to NONMEM METHOD=IMPMAP. The difference is the proposal: impmap re-centers at the freshly-computed conditional mode (MAP) every iteration (robust on rich data), whereas imp re-centers from the previous iteration’s sample moments (cheaper, but fragile on rich data). Both update θ/Ω/σ from the importance-weighted posterior moments. It runs standalone or as a chain stage:

[fit_options]
  method             = importance_sampling_map   # alias: impmap
  impmap_iterations  = 200
  impmap_samples     = 300
  impmap_proposal_df = 4          # Student-t (default); `normal` for MVN
  impmap_seed        = 12345
Key Default Description
impmap_iterations 200 Number of MCEM iterations (parameter updates).
impmap_samples 300 Importance samples K per subject per iteration. Larger K reduces Monte-Carlo noise at linear cost.
impmap_proposal_df 4 Proposal degrees of freedom. A finite value ≥ 1 gives a heavier-tailed Student-t (default 4); normal (or mvn) gives a multivariate-normal proposal (NONMEM’s default). The Gaussian’s lighter tails under-cover the posterior of weakly-identified parameters and bias the M-step moments, so ferx defaults to a Student-t.
impmap_auto true Adaptive sample count (NONMEM AUTO). When true, impmap_samples is the starting count and is ramped up (×2/iteration, cap 10000) while the objective’s Monte-Carlo SE exceeds 1.0 (NONMEM STDOBJ). Strongly recommended for FREM / high-dimensional models — a fixed count leaves a sample-count-dependent bias in the typical-value and Ω estimates.
impmap_averaging 50 Final iterations whose parameters are averaged into the reported estimate (Monte-Carlo variance reduction).
impmap_seed 12345 RNG seed. Same seed → identical estimates.
impmap_low_ess_threshold 0.1 Subjects with normalized ESS below this fraction are flagged as poorly sampled.
impmap_trace false When true, collect per-iteration parameter values into FitResult.impmap_trace — analogous to NONMEM .ext output for traceplots.
impmap_mceta 0 Number of additional random starting points for per-subject MAP optimization (analogous to NONMEM MCETA). Each start draws η from N(0, Ω). The start with the lowest individual NLL wins. 0 = single warm-start (default). 3 is a good choice for high-dimensional models (e.g. FREM with ≥ 5 ETAs).
impmap_sobol false Use Sobol quasi-random sequences (with Cranley-Patterson randomization) for IS draws instead of pseudo-random. Gives more uniform posterior coverage with fewer samples. Only applies to MVN proposals (impmap_proposal_df = normal); Student-t falls back to pseudo-random.
iscale_min 0.1 Minimum proposal scaling factor for adaptive IS (NONMEM ISCALE_MIN). The IS proposal covariance is multiplied by s² where s is chosen from [iscale_min, iscale_max] to maximise per-subject ESS. Set both to 1.0 to disable.
iscale_max 10.0 Maximum proposal scaling factor (NONMEM ISCALE_MAX).
frem_rao_blackwell true FREM only: Rao-Blackwellise the covariate ETAs (integrate them analytically, importance-sample only the PK ETAs). Strongly recommended — brute-force sampling of the near-singular covariate dimensions has very poor ESS. Set false only to diagnose the RB path against the full-dimensional sampler.

impmap reuses inner_maxiter / inner_tol for the per-iteration MAP step. Inter-occasion variability ([iov] / kappa) and SDE ([diffusion]) models are not yet supported and are refused up front. See the IMPMAP documentation for the algorithm and the NONMEM comparison.

Optimizer Choices

Optimizer Description Recommended For
auto Picks nlopt_lbfgs (analytic gradient available) or bobyqa (FD only) per model Default; the recommended choice for most fits
bobyqa NLopt BOBYQA — derivative-free quadratic interpolation Robust on noisy / non-smooth FOCE surfaces (ODE/PD, sparse data, Hill-ridge models); what auto selects there
slsqp Sequential Least Squares Programming (NLopt) Smooth, well-conditioned analytical PK models where you want gradient-based convergence; pair with reconverge_gradient_interval if it stalls
nlopt_lbfgs NLopt L-BFGS Analytic-gradient FOCE/FOCEI models (fastest-to-optimum in #490 benchmarking); large parameter spaces; what auto selects there
lbfgs Deprecated alias → nlopt_lbfgs — (was the built-in L-BFGS; see #483)
bfgs Deprecated alias → nlopt_lbfgs — (was the built-in BFGS)
mma Method of Moving Asymptotes (NLopt) Constrained problems
trust_region Newton trust-region with Steihaug CG subproblem (argmin) Well-conditioned problems where second-order curvature helps convergence

Notes: - auto (the default) chooses between nlopt_lbfgs and bobyqa from the model: when the exact analytic FOCE/FOCEI gradient is available (analytical PK model, in the sensitivity provider’s scope, and gradient not forced to fd) it uses the gradient-based nlopt_lbfgs; otherwise — ODE/PD models, LTBS/SDE, or gradient = fd, where the outer loop runs on finite differences — it falls back to the derivative-free bobyqa. Benchmarking across ~10 real FOCEI datasets (#490) found nlopt_lbfgs fastest-to-optimum on every analytic-gradient problem and bobyqa both fastest and most reliable on the FD-only problems. The fit output records the resolved optimizer as auto (<resolved>). Set optimizer explicitly to override the choice. - bobyqa does not use gradients, so it is robust to small discontinuities in the FOCE surface caused by EBE re-estimation, but it converges more slowly than gradient-based methods on smooth problems. - trust_region uses the analytic outer gradient (same subject_nll_pop_grad as the outer FOCE optimizer) and a BHHH approximate Hessian (H ≈ 4 Σ gᵢgᵢᵀ). The BHHH matrix is always positive semi-definite, so the Steihaug subproblem is well-conditioned even near constraints. The Steihaug CG budget defaults to ceil(sqrt(n_params)) — typically 5 for standard NLME models, which is far cheaper than the previous FD-Hessian path (O(n²) OFV evaluations per Hessian).

Parameter Scaling and EBE Convergence

Key Default Description
parameter_scaling auto, none, abs, rescale2 auto
scale_params false Legacy boolean alias for parameter_scaling = abs. Divide each packed (log/Cholesky) coordinate by its initial magnitude before passing it to the optimizer. Off by default since issue #99: the scaling-enabled path only ever runs on log/Cholesky coordinates, where dividing by \|log value\| is counterproductive — e.g. ln(V)=ln(20)≈3 gets scale 3, turning the optimizer’s unit step into a ≈20× multiplicative jump in V, which overshoots and (via the uniform SLSQP gradient cap) starves the step in other dimensions such as OMEGA, halting short of the minimum. Prefer parameter_scaling = rescale2 for gradient-based optimizers. The OFV value is unchanged at any fixed point, but the optimizer trajectory and stop point are not.
max_unconverged_frac 0.1 Fraction of subjects (with at least min_obs_for_convergence_check observations) allowed to have unconverged EBEs before the outer optimizer rejects the step (returns OFV = ∞). Set to 1.0 to disable the guard.
min_obs_for_convergence_check 2 Subjects with fewer than this many observations are excluded from the max_unconverged_frac check (they still run normally).
stagnation_guard true Short-circuit the NLopt-based outer optimizers once recent evals show no OFV improvement above 1e-3 over a window of 3*(n+1).max(50) evals. This lets SLSQP / L-BFGS terminate quickly via their own xtol/ftol on numerically-flat plateaus (e.g. γ-bearing FOCEI problems) instead of grinding through the remaining outer_maxiter budget at full inner-loop cost. Set to false to let the optimizer run to its natural termination criterion — useful when debugging or for problems with very slow but real OFV improvements below the threshold.
ebe_warm_start false When an inner per-subject EBE solve fails its BFGS step and falls back to Nelder–Mead, warm-start the simplex from the BFGS partial η̂ rather than cold-starting from the prior mode η=0. A weakly-identified η (e.g. an unidentifiable peripheral volume) drives BFGS far out onto the steep prior slope, and NM slides down it to the mode in far fewer iterations than refining from 0 in the flat basin — substantially fewer prediction walks on fallback-heavy fits (≈1.7× on a 2-cpt unidentifiable-V2 benchmark). Opt-in: warm-starting moves the fallback subjects’ EBEs, which perturbs the outer optimiser’s trajectory — harmless for the derivative-free BOBYQA default, but it can derail a gradient-based outer optimiser (e.g. mma) into a worse basin on some models. Enable only after validating OFV/estimates on your model + optimizer; leave false for the historical cold-restart behaviour.

Options That Don’t Apply to the Selected Method

If you set an option that the chosen estimation method doesn’t consume (e.g. n_convergence with method = focei, or optimizer with method = saem), fit() emits a warning listing the option, the selected method, and the keys that are available for that method. The option is ignored — the fit still proceeds.

For a chained fit (method = [saem, focei]), an option is kept if it applies to any stage in the chain, so SAEM and FOCEI keys can be mixed without warnings.

Multi-Start Optimization

Key Default Description
n_starts 1 Number of independent optimization runs. 1 disables multi-start (no behaviour change). When > 1, all starts run in parallel via rayon; the converged run with the lowest OFV is returned. Start 0 always uses the exact initial values from the model file.
start_sigma 0.3 Log-space perturbation applied to initial theta values for starts 1..n. Log-packed thetas are multiplied by exp(N(0, start_sigma)); thetas with negative lower bounds are shifted additively.
multi_start_seed 42 RNG seed for the multi-start theta perturbations. Independent of seed (SAEM) so that changing the SAEM seed does not silently alter which perturbed starting points are used in FOCE multi-start runs.

Multi-start is most useful for models prone to local minima: nonlinear elimination (Michaelis-Menten), full-block omega, or many covariate parameters. On an 8-core machine n_starts = 8 costs the same wall-clock time as a single run but has ~8× lower probability of a local minimum.

Examples

Standard FOCEI with defaults:

[fit_options]
  method     = focei
  maxiter    = 300
  covariance = true

FOCEI with global search:

[fit_options]
  method        = focei
  maxiter       = 500
  covariance    = true
  global_search = true

SAEM with custom settings:

[fit_options]
  method        = saem
  n_exploration = 200
  n_convergence = 300
  n_mh_steps    = 5
  seed          = 42
  covariance    = true

FOCEI with SIR uncertainty:

[fit_options]
  method        = focei
  covariance    = true
  sir           = true
  sir_samples   = 1000
  sir_resamples = 250
  sir_seed      = 42

Derivative-free BOBYQA fit:

[fit_options]
  method        = foce
  optimizer     = bobyqa
  maxiter       = 300
  inner_maxiter = 100
  inner_tol     = 1e-6

Trust-region with tuned CG subproblem:

[fit_options]
  method             = foce
  optimizer          = trust_region
  maxiter            = 200
  steihaug_max_iters = 30

FOCE with Inter-Occasion Variability:

[fit_options]
  method     = foce
  iov_column = OCC
  covariance = true

Enable optimizer trace and EBE guard:

[fit_options]
  method                        = foce
  optimizer_trace               = true
  max_unconverged_frac          = 0.5
  min_obs_for_convergence_check = 3

Optimizer Trace

When optimizer_trace = true, a CSV is written to /tmp/ferx_trace_<pid>_<ts>.csv and the path is stored in FitResult::trace_path. Each row is one outer iteration.

Column Populated by Description
iter all Iteration number
method all foce, focei, laplace, gn, gn_hybrid, saem
phase gn_hybrid, saem focei (polish) or explore/converge
ofv all Objective function value
wall_ms all Wall time for this iteration (ms)
grad_norm BFGS, NLopt gradient-mode, GN Gradient norm
step_norm BFGS Step size
inner_iter_count (reserved) Reserved for future per-iteration inner-loop count; currently NA
optimizer FOCE/FOCEI Active NLopt algorithm
lm_lambda GN Levenberg-Marquardt damping factor
ofv_delta GN Change in OFV from previous iteration
step_accepted GN Whether the GN step was accepted
cond_nll SAEM Conditional observation NLL
gamma SAEM SAEM step-size (1 in exploration, 1/k in convergence)
mh_accept_rate SAEM Mean acceptance rate across all subjects (MH or HMC). In mixed HMC/MH runs (n_leapfrog > 0 with some MH-fallback subjects) this is an aggregate across both samplers.
n_ebe_unconverged FOCE/FOCEI Subjects whose inner optimizer did not converge
n_ebe_fallback FOCE/FOCEI Subjects that fell back to Nelder-Mead
val:<name> all Per-parameter estimate, one column per optimized coordinate, in natural / reporting scale (θ as values, Ω as variances/covariances, Σ as variances)
grad:<name> FOCE/FOCEI/GN Per-parameter gradient, in the optimizer’s transformed/scaled space (NA for SAEM); sqrt(sum(grad:<name>^2)) == grad_norm

After the fixed columns above, the trace appends two blocks of per-parameter columns — all val:<name> columns, then all grad:<name> columns — one pair per optimized coordinate (thetas, Ω diagonals and off-diagonal covariances, Σ, plus any IOV Ω). Headers use the declared parameter names where available (val:TVCL, val:ETA_CL, an Ω off-diagonal as val:ETA_V~ETA_CL, val:PROP_ERR), with THETA1 / OMEGA(2,1) / SIGMA(1) fallbacks (a name containing a comma is CSV-quoted). The column set is fixed across a method chain (e.g. saem -> focei), so one header serves the whole run. grad:* columns are NA on any iteration without an OFV gradient (SAEM, and derivative-free BOBYQA evals). The val:*/grad:* cells are written in scientific notation ({:.6e}) — unlike the fixed-decimal scalar columns — so small variances and near-converged gradient coordinates keep their precision instead of rounding to 0. This powers the per-parameter convergence view in ferx-r’s trace UI (#640).

Unused columns contain NA. The trace is buffered and flushed when the fit ends; if a run is killed mid-iteration the buffered tail may be lost.