Fit Options

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

Syntax

[fit_options]
  key = value

Every line is a key = value pair, and a line that is not one is a parse error naming the repair — it is not ignored:

[fit_options]
  method saem       # rejected
[fit_options]: `method saem` is not a `key = value` pair — did you mean
`method = saem`? Every `[fit_options]` line assigns one key with `=`; a line
without it is rejected, not ignored.

Before 0.4.0 such a line was dropped without a diagnostic. That is worth spelling out because of what the dropped line usually is: method written without its = differs from the accepted spelling by one character, the surrounding DSL is not uniformly key = value ([parameters] and [error_model] lines are not), and the consequence is that the fit silently runs the default estimator. A run that falls back from SAEM to FOCEI is a different algorithm, not a perturbed start, and its objective is not comparable to the one you asked for. ferx check reported the file VALID, and only the {model}-fit.yaml recorded the method that actually ran. See #1390.

; does not start a comment in a .ferx file — use # or //. A line carrying one (method = saem ; run 3, or a whole line commented out with a leading ;) says so rather than reporting an unknown key.

General Options

Keys that apply to every estimation method, grouped by the machinery each one controls. Method-specific keys live in their own sections further down (AGQ, SAEM, VI, SIR, IMP, IMPMAP).

Run control

What runs, how long, and how it is resumed or traced.

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). A pinned count is an upper bound on the whole fit, n_starts > 1 included: a multi-start fit runs its start fan-out and the per-subject loops underneath it on the same pinned pool. The CLI’s --threads flag mirrors this key (see ferx --help) and overrides it when both are given: an explicit ferx model.ferx --threads 1 wins over a file that says threads = 8, and a warning names both counts. --threads 0 / --threads auto counts as explicit too — it asks for the engine default by name, so it also overrides a pinned count (#1416).
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 — except that the start is still clamped into the optimizer’s bounds, so a parameter that begins outside them is evaluated at the bound (e.g. a free omega ~ 0.0 at exp(-12), not 1e-8); FIX pins a coordinate at its own start and avoids this.
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. What is saved depends on the stage’s method. A deterministic stage (foce, focei, laplace, gn, gn_hybrid) saves the best point it has reached, not the last one evaluated (#1317): the objective is evaluated at every point the optimizer probes, including the trial steps a line search rejects, so saving the current evaluation could record a throwaway point far from the fit’s actual position — the file’s iter and ofv are then the evaluation index and the (unpenalized) objective at that best point. A saem stage instead saves its latest state, which is what a correct continuation of the stochastic-approximation chain resumes from; its ofv is that iteration’s conditional NLL, not a best-so-far and not comparable to a deterministic stage’s. method_chain and stage_idx in the file say which stage wrote it. 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.
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.

Worker reuse for repeated fits

Default fits share a persistent worker pool. Unpinned fits with the same fit-level ODE overrides also share a persistent pool keyed by those solver settings, so changing a tolerance does not multiply worker counts for concurrent calls. A positive threads value explicitly budgets the fit and exclusively leases a pool for the call; an idle pool with the same width and solver settings can be reused, while concurrent pinned fits keep independent budgets. This preserves the per-fit budgets assigned by PoolPlan. This isolates worker budgets and ODE settings; concurrent fits should still use consistent process-wide inner-optimizer and EBE warm-start settings.

Either way, a requested worker count sizes one pool: the one the fit runs on. The CLI’s --threads N used to size Rayon’s process-global pool as well, so a run held 2N worker threads to do N threads of work — the extra N measured 0.0% busy for the whole fit. Since #1460 the flag only declares the width, and ferx --threads N runs with N workers plus the main thread.

The same count applies to parallel work that is not a fit — GAM covariate screening, standalone run_sir / run_covariance, the NCA initial-estimate pass and npde. Those used to spread over every logical CPU regardless of the flag (ferx gam --threads 2 screened 15 wide on a 15-core host), because they run outside any ferx pool and a bare parallel loop there belongs to Rayon’s global pool. A tool built on ferx-core gets the same behaviour by wrapping its parallel loop in install_on_engine_pool.

The caches ordinarily retain at most twice the automatic thread budget across all settings (at most 16 workers each). Older entries are evicted. When a single pool is wider than that budget, it can remain as the sole most-recent entry so a repeated explicit wide configuration still avoids rebuilding its workers. Active fits keep their requested width; the idle limits are not process-wide limits on concurrent fitting. Worker stacks remain 32 MiB each, which is reserved capacity rather than necessarily resident memory.

During FOCE/FOCEI optimization, the same pool serves the inner EBE solves and outer gradient evaluations. The main NLopt and built-in BFGS/L-BFGS paths evaluate each subject’s marginal contribution immediately after its EBE solve. The mixed analytic/FD gradient path also finishes each subject’s fallback in the same task. This reduces synchronization between stages while preserving subject-order reductions. Workers still use Rayon’s normal idle policy; pool reuse does not mean keeping idle CPU cores busy continuously.

IOV inner solves also reuse per-event PK parameter buffer capacity between likelihood probes. Each probe still evaluates the current parameters and every event’s covariates, time, and occasion; stored numerical values are not reused.

Outer optimizer

The population-parameter search: which algorithm, when it stops, and whether a global pre-search runs first. See Optimizer Choices for how to pick one.

Algorithm and stopping rules

Key Values Default Description
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 (out-of-scope ODE/PD models, SDE, or gradient = fd) and for every [mixture] model — 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 — out-of-scope ODE/PD, SDE, gradient = fd, and TTE — and for every [mixture] model), 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.
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.

Reporting

Key Values Default Description
report_final_gradient true, false true Compute a finite-difference gradient at the reported estimates when the optimizer supplied none of its own, and report it as final_gradient with final_gradient_source = "finite_difference" (#997). Only bobyqa — the derivative-free auto pick for out-of-scope ODE/PD, SDE, gradient = fd and [mixture] models — has no gradient of its own; every other outer optimizer reports the one it used, labelled "optimizer". Without this a derivative-free run reports converged = true with nothing in the result to check the claim against, which is how a BOBYQA fit that stopped 1.54 OFV short of the L-BFGS optimum on the same model looked identical to one that had not. It is a reporting quantity: computed once, after the fit, never fed to the optimizer, so the estimates are bit-identical either way. Costs 2 × n_free objective evaluations — one gradient’s worth, the same as a single step of the optimizers that get it for free. Set false when that is not worth paying (a large ODE model whose every evaluation is seconds).

Inner loop (EBE)

The per-subject empirical-Bayes search that the outer objective is built from.

Key Values Default Description
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. Residual noise in the EBE solution propagates into the marginal objective the outer optimizer sees, so this is not a free speed knob in either direction — see Choosing inner_tol.
inner_restarts integer 1 Guarded multi-start count for the inner (per-subject EBE) optimizer, to escape a multimodal individual objective. A subject that a single warm-started BFGS would leave in the wrong basin re-solves from N Ω-scaled cold-start seeds and keeps the lowest-objective mode; unimodal subjects are bit-identical to 0 and pay only a flatness probe. Set 0 to disable — see Escaping a multimodal individual objective.
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, at roughly 5–6× the per-gradient cost; N reconverges on evaluations 0, N, 2N, … and uses the cheap gradient in between. Raise it when a gradient optimizer stalls above the derivative-free optimum. IOV models (kappa/block_kappa) always reconverge and ignore this setting. See Gradient accuracy vs cost.
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).

Choosing inner_tol

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.

Escaping a multimodal individual objective

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.

The initial BFGS metric also decides which basin a cold start lands in. Under focei, foce, gn, gn_hybrid and laplace (any n_agq), warm-started inner solves seed BFGS from the subject’s conditional η-Hessian, which speeds up re-convergence near the previous mode. A cold start (the first outer evaluation and the final cold EBE re-solve) runs twice, once with that seed and once with the historical metric, and keeps the lower objective. Neither metric picks the right basin on every subject (#1504). If both reach the same mode, the historical solve is kept, so unimodal subjects are unchanged. Models with a [covariate_nn] block never use the seed.

Covariance and standard errors

The post-fit covariance step. See Covariance Step for what each estimator means.

Which estimator to assemble

Key Values Default Description
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 three need a positive-definite R, analytic or FD (a non-PD R is not yet routed to the s cross-product). Requires covariance = true. The estimator that actually ran — not always the one requested — is reported back: see Which covariance estimator produced these numbers.
covariance_fallback none, sir none What to do when the covariance Hessian is non-PD, whether R is analytic or finite-difference. 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.

How the Hessian is computed

Key Values Default Description
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. A model-level exclusion (method = laplace, a mixture, gradient = fd, …) keeps the whole fit on the finite-difference stencil, which is correct for all of them. A subject-level one does not: since #1514 the covariance step assembles the exact term for every in-scope subject and finite-differences only the declining subjects’ own terms, adding them into the same information sum \(\sum_i R_i\). Each term is then an estimate of itself — the observed information is a per-subject sum, so this is not a blend of two approximations of one matrix — and the fit reports a W_COV_ANALYTIC_SALVAGE note naming the salvaged subjects. When at least half the population declines, the whole-population stencil is taken instead: there is no analytic majority left to carry the cost, and rebuilding it subject by subject measured +4% wall for a matrix identical to 7 significant figures. 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.
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. Must be positive and finite.

ODE solver tolerances

Accuracy and step budget for [odes] models, and for the auto-generated ODE twin that closed-form transit / inverse-Gaussian absorption models carry. See ODE Models.

Where the value comes from

These six keys (and the two stepper keys below) can also arrive from a caller, and precedence is per key either way: a key the caller sets wins over the model file for that fit, a key the caller does not set keeps the file’s value. What counts as “set” differs between the two front ends, so check which one you are on.

From R, ferx_fit(settings = list(ode_reltol = 1e-9)) merges your settings onto the model file’s own [fit_options] and stamps the result on the model. Every key you name wins — including one you name with its default value, so settings = list(ode_reltol = 1e-4) does pull a file that pinned 1e-10 back to 1e-4. Keys you leave out keep whatever the file said.

From Rust, fit(&model, &pop, &params, &options) — and equally run_covariance(), run_sir() and run_sir_core(), which integrate the same way — compares your FitOptions field by field against FitOptions::default() and carries only the moved fields for the duration of that call, so a hand-built settings struct cannot silently loosen a model that pinned ode_reltol = 1e-10. The consequence is the mirror image of the R rule: asking explicitly for a default value (ode_reltol = 1e-4) reads the same as not asking, and the model file still wins. That override also lasts exactly one fit() — a later predict() on the same &CompiledModel integrates at the model file’s values, so a fit’s IPREDs and a subsequent predict() can come from different tolerances. Call model.sync_ode_solver_opts(&options) on an owned model when both paths must match (that is what the R front end does).

From the CLI, a flag follows the R rule, not the Rust one: ferx is a caller, so a flag that appears on the command line wins over the model file even when its value is the default. --data has always worked that way (it beats the [data] block, with a warning when they differ), and since #1416 --threads does too — including --threads 0 / --threads auto, which name the engine’s own worker count rather than going unmentioned. A flag you do not type leaves the file in charge.

Put the key in [fit_options] when a model needs a specific accuracy regardless of who runs it.

Set one of these on a model that never integrates — an analytical PK model with no [odes] block and no closed-form absorption twin — and the fit warns that the key has no effect (#518), instead of dropping it silently. A model that does integrate, including a closed-form transit / inverse-Gaussian model reaching its twin, never draws that warning.

The keys

Key Values Default Description
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. Standard errors have their own sensitivity: when the covariance step falls back to the finite-difference Hessian (see covariance_regularized) it divides by h² and so amplifies integration noise, and SEs measured on the 3-cpt IV ODE fixture were 2-7x too large at the default, on the plateau from ode_reltol = 1e-6 / ode_abstol = 1e-8 (#520). The exact analytic covariance R-matrix, which most Gaussian [odes] fits now take, avoids that mechanism: it never second-differences the objective, so there is no 1/h² for integration error to be amplified through, and the #520 sweep measures it flat across eight orders of magnitude of ode_reltol. It is not step-free, though — third_order_fd_step finite-differences the Dual2 sensitivity jet to assemble the third-order blocks, so the analytic route keeps a step sensitivity of its own, across a smooth sensitivity rather than across the objective. Where that step would cross an event kink — a lagged-dose arrival, an infusion end, a zero_order window edge — the sweep bounds it by the observation-to-event gap, so no sample changes sides of the event between the two points of a difference pair; a subject whose mode sits on a moving event (a corner minimum of the inner objective, which a lagged arrival at a dense early sample does produce) has no third derivative there and is declined to the per-subject salvage instead. Before that bound, SE(TVLAG) on a lagged-dose fixture moved 27 % with the step (#1505).
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_stiff_abort_after integer ≥ 0, or off off Abandon an integration segment once this many of its steps have clamped at the minimum step size, instead of grinding on to ode_max_steps. A clamp means the step failed its local-error test and was accepted anyway because dt could not shrink further — the segment is stability-limited, not accuracy-limited, and the remaining thousands of steps buy wall time rather than accuracy. Off by default, and deliberately so: aborting freeze-pads the segment’s remaining output times with the last state, so it trades a slow-but-eventually-integrated segment for a padded one — a change to the objective, not only to its cost. On the time-to-event path (which returns an event time rather than a trajectory) an abort is reported as a failed segment instead of being padded. It applies to every likelihood evaluation, not just the final one, so a budget that actually bites makes the objective discontinuous in θ — a segment that clamps N-1 times at θ and N times at θ+h returns a differently-truncated trajectory — which stalls the outer optimizer and corrupts the finite-difference Hessian the covariance step builds. Treat estimates and standard errors from a fit whose ode_solver warning reports a non-zero abort count as diagnostic output, not as results. Reach for it when a fit is grinding and you want it to say so quickly; fix the grinding with ode_method = auto (the default) or a named stiff method, then turn the budget off. 0 / off / none / false disable it.

ODE stepper selection

Which integrator runs, and whether auto may change it mid-segment. See Stiff systems and Letting ferx pick the stepper.

Key Values Default Description
ode_method auto, rk45, vern7, rosenbrock23, rodas4, rodas5p auto Which stepper integrates the [odes] block. Two axes: vern7 (Verner 7(6)) buys order, for fits that are accuracy-limited at tight tolerances — not a blanket upgrade, since it is slower than rk45 at default tolerances. The three linearly implicit Rosenbrock methods buy stability, 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). rk45 is the explicit Dormand-Prince method and the fastest choice on a non-stiff model. Pick by tolerance: rosenbrock23 (aliases ros23/ode23s) for crude tolerances or rough right-hand sides, rodas4 at typical 1e-6–1e-9, rodas5p at ode_reltol ≤ 1e-9. A stiff step is not free — n + 1 extra right-hand-side evaluations for the finite-difference Jacobian, plus one n × n factorization — so switch only when stiffness, not accuracy, caps the step (see below for what n is). The default auto decides per integration segment; see Stiff systems for the full method table and Letting ferx pick the stepper for how it decides.
ode_auto_switch true, false true Let ode_method = auto change stepper inside an integration segment, not only at its start. The start-of-segment probe is blind to a model that turns stiff only once drug and target have accumulated — a depot-absorption binding model reads |Re λ|max = 10 when the dose lands and 7.6e3 an hour later. With this on, RK45 evaluates its dimensionless, stage-based stiffness estimate after every accepted step at no extra RHS cost; the verdict counters survive event boundaries. The Jacobian probe remains as a periodic backstop every 25 accepted steps. A switch keeps the accepted state and continues with the other stepper. Being a threshold test, it makes the objective discontinuous in θ: if a covariance step returns implausible standard errors on a model whose ode_solver warning reports mid-segment switches, re-run with ode_auto_switch = false (or a named stiff method, which is pinned and takes no such decisions) to see whether the switching is what moved them. Ignored under a named ode_method. See When a segment turns stiff halfway through.

Stiff-step cost follows the integrated system, not the state count

The n × n factorization a Rosenbrock step pays for is sized by 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 means 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.

Data handling and diagnostics

How the dataset is interpreted, and what diagnostics are produced alongside the fit.

Key Values Default Description
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.
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.
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

Two reference pipelines are used. The warfarin check below post-processes a NONMEM $SIMULATION with the npde R package; the IOV anchor further down uses NONMEM’s own $TABLE ... NPDE NPD ESAMPLE= output directly. 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: rows the fit scores as censored — CENS != 0 under bloq_method = m3 — are emitted as empty/NaN (matching the CWRES/IWRES convention), and a subject’s NPDE is NaN as a whole when it has any such row, since the within-subject decorrelation would otherwise mix the LLOQ into the uncensored rows. Under the default bloq_method = drop a CENS != 0 row is fitted as an ordinary observation at its DV, so it is also scored as one here and does not void its subject’s NPDE (#1499). 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).

Inter-occasion variability (IOV)

For an IOV (kappa) model the reference distribution draws one independent \(\kappa_{i,o} \sim N(0, \Omega_{\text{IOV}})\) per occasion, exactly as simulate() does, so the predictive variance NPDE/NPD are scored against carries its inter-occasion component. Earlier versions held every kappa at zero here, which left the reference too narrow and the scores correspondingly over-dispersed — a well-specified IOV model looked mis-specified. Non-IOV models are unaffected and their scores are unchanged.

The cross-tool check is tests/npde_iov_nonmem_anchor.rs: 60 subjects, three EVID=4 occasions each, all parameters fixed at the simulating values on both sides, ESAMPLE = nsim = 2000. Row by row, ferx and NONMEM agree to a worst |ΔNPD| of 0.226 and |ΔNPDE| of 0.280 — Monte-Carlo noise between two independent 2000-draw empirical CDFs. Forcing \(\Omega_{\text{IOV}}\) to zero (the old behaviour) breaks that agreement to 0.895 and 4.389.

NPDE decorrelation: ferx follows the npde package, not NONMEM’s default

NONMEM does emit NPDE and NPD natively — $TABLE ID TIME DV NPDE NPD ESAMPLE=2000 SEED=… — but its default decorrelation is not the one the npde R package (and ferx) use. Decorrelating a within-subject vector is not unique:

factor of the simulated covariance used by
symmetric square root eigen-based, the CWRES convention NONMEM, by default
Cholesky factor Brendel/Comets, npde::autonpde() ferx, and NONMEM with WRESCHOL

On the anchor above, the same NONMEM run writes a bit-identical NPD column either way (NPD carries no decorrelation) while its NPDE correlates 0.868 with ferx’s by default and 0.997 with WRESCHOL on. So compare NPDE against a WRESCHOL table; comparing against the default one is comparing two different statistics.

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.

agq_eval_only — AGQ as a likelihood evaluator

method        = vi, laplace
n_agq         = 21
agq_eval_only = true

With agq_eval_only = true a laplace stage stops being an estimator and becomes a pure evaluator: it reconverges the EBEs at the parameters it was handed, evaluates the adaptive-quadrature marginal there, and reports it as ofv. The preceding stage’s θ, Ω and σ are left exactly as they were. It must be the final stage — an evaluator mid-chain would report a likelihood at parameters the next stage overwrites.

This is the deterministic counterpart to methods = ..., imp with imp_eval_only = true. Both turn a fitted point into a −2 log L; the AGQ route carries no Monte-Carlo error, so two runs agree bit for bit, at the cost of n_agq ^ n_eta likelihood evaluations per subject.

Two uses:

  • Finishing a VI fit. The ELBO is a lower bound, not a likelihood, so ofv is NaN by default. This gives a real marginal likelihood at the VI estimate — and being deterministic, it is directly comparable across VI fits. vi_final_ofv = laplace is the same idea pinned at one node.
  • Checking a bound or an approximation. Because −2·ELBO ≥ −2 log p(y) must hold, evaluating the marginal at a VI point tells you how tight the variational approximation actually is. The same comparison against a FOCE/FOCEI OFV measures what the first-order approximation costs on your model.

The EBEs are recomputed rather than inherited, because adaptive quadrature centres its grid on the conditional mode and scales it by the posterior Hessian — and a VI stage reports variational means, which are not modes.

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.

VI

method = vi

Variational inference. Instead of profiling the random effects out at their mode (FOCE) or sampling them (SAEM), VI fits a tractable posterior q(η) per subject and optimizes its parameters jointly with θ/Ω/Σ by maximizing the evidence lower bound. There is no inner loop.

Suits deep compartment models ([covariate_nn]), models where the FOCE inner loop is unstable, and cases where you want a per-subject posterior covariance without paying for a Hessian — VI produces one as a by-product.

Its objective is a lower bound, not a likelihood. ofv is NaN after a VI fit by default rather than being filled with a number that is not comparable to a FOCE or SAEM OFV. Set vi_final_ofv = laplace, or chain method = vi, laplace with agq_eval_only = true (deterministic) or method = vi, imp with imp_eval_only = true (importance sampling), to evaluate a genuine marginal likelihood at the VI estimate. Inter-occasion variability ([iov]) is supported — q spans the stacked [η, κ₁ … κ_K] vector and the per-occasion posterior means are reported under vi$kappa_means. Non-Gaussian endpoints (TTE / categorical) are not: VI’s data term would omit their likelihood contribution, so the model is rejected rather than silently mis-fitted. See the VI page.

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

Schedule and E-step

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 auto Block Metropolis-Hastings steps per subject per iteration, or auto — the default — to size it from the dataset as 2.5 × observations per subject per η, clamped to [6, 20]. 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. method = bayes reads this option but not the auto rule: its η block keeps the historical fixed count of 20 unless you set one (the rule is calibrated on SAEM, whose componentwise kernel that sampler does not run). See Choosing n_mh_steps.
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, under scale_adaptation = interval. Also governs the κ (IOV) scales under either rule.
scale_adaptation robbins_monro How the MH step scales are adapted. interval is the historical rule (×1.1 / ×0.9 every adapt_interval); robbins_monro steps log δ every iteration by c·k^-0.6·(accept − target). The interval rule can only move δ by 1.1^n / 0.9^n over the n times it fires, which on a 400-iteration run is 8 — not enough reach for a model whose optimal step is an order of magnitude from the start — examples/warfarin_saem.ferx spends an entire run at 2–4 % acceptance against a 40 % target under it. robbins_monro became the default in #1449, together with scale_deadband; set scale_adaptation = interval for the historical rule. See the SAEM page.
scale_deadband 0.15,0.60 Acceptance window lo,hi in which the robbins_monro step is skipped, leaving the scale untouched; outside it the ordinary step applies. none steps on every iteration, which is the unbanded rule of #1444. Ignored under scale_adaptation = interval. A band wider than the spread of the per-iteration rate freezes the scales rather than disabling the band, so 0,1 is rejected — see the SAEM page.
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.
seed 12345 RNG seed for reproducibility

The numerical θ/σ M-step

The channel that estimates a theta the closed-form mu-reference update cannot reach — a covariate effect, or a structural parameter left without IIV. See Averaging a maximiser is not maximising the average.

Key Default Description
mstep_damping 1.0 (off); 0.03 with iiv_on_ruv Optional exploration-phase cap on the stochastic-approximation step for the numerical θ/σ M-step (the channel that moves a theta with no ETA). Below 1.0 the M-step result is blended in as θ ← θ + γ_θ·(θ* − θ) during exploration and averaged at γ = 1/(k−k1) in convergence. Sigma comes out of the same joint NLopt solve but takes its own step, min(γ_k, 0.2, γ_θ) on the variance scale, which is on whether or not you set this option (#1445). σ is therefore never given a larger step than θ, and setting this option at or below 0.2 makes the two equal throughout exploration, since σ never steps faster than θ. Off by default since #1415: a cap divides the number of EM steps the exploration phase amounts to and was measured to hold every no-ETA theta near its initial estimate (a covariate exponent at 0.32 from a 0.3 start against NONMEM’s 0.92 at the old default of 0.03). The one exception is a model with iiv_on_ruv, which keeps the #1011 default of 0.03: on that shape the undamped channel drifts (the same model without iiv_on_ruv does not), and the cap holds the no-ETA thetas near their start. A value you set wins either way. Must be in (0, 1]. 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.
mstep_solver bobyqa Which solver moves the numerical θ/σ M-step — the channel that estimates a theta with no ETA (#1458). bobyqa is the historical rule: a short derivative-free re-maximisation of the frozen-η objective at this iteration’s single η draw, whose maximiser is then adopted. Because a maximiser is a nonlinear function of the draw, that recursion converges to E[θ*(η)] rather than to the θ that maximises E[Q(θ, η)] — a Jensen-type bias that grows when the E-step mixes better. score_sa solves the score equation E[∇Q] = 0 by stochastic approximation instead — one Robbins-Monro step per M-step, preconditioned by the SA-averaged expected information — which has no Jensen term and costs one gradient sweep in place of mstep_maxiter·(n+1) population passes, so it runs every iteration rather than every third during exploration. Restricted to the plain Gaussian residual scope the expected information has a closed form on (no IOV, mixture, M3, TTE endpoint, block_sigma, [covariate_nn] θ, residual magnitude or FREM); a model outside it keeps bobyqa and says so by name in the fit’s warnings. σ takes the same Newton direction as every other coordinate but on the #1445 σ schedule and variance-scale blend, not at the shared γ (#1480 — before that fix it re-opened #1445’s additive-σ collapse on #1445’s own fixture). It stays opt-in: #1449 measured it tie-or-better than bobyqa on the importance-sampled −2 log L of all eight SAEM benchmarks at 25–54 % less CPU, but it converges more slowly than bobyqa on a no-ETA θ that starts far from its optimum — on the #619 covmuref_power anchor it reaches 0.746 against NONMEM’s 0.921 at the default 150/250 schedule and 0.9845 at 300/700 — see the SAEM page. See the SAEM page.
mstep_draws 1 Number of E-step draws the bobyqa M-step objective is averaged over (#1458). K > 1 evaluates the frozen-η objective as the mean over this iteration’s η draw and the previous K − 1 iterations’ draws, so what NLopt returns is the maximiser of an average rather than one draw’s maximiser; the Jensen bias then falls roughly as 1/K — it shrinks, it does not go away. Costs K times the M-step’s objective evaluations, and pays only to the extent that consecutive draws decorrelate, which under the default sticky random-walk E-step they largely do not. Ignored under mstep_solver = score_sa, whose Robbins-Monro average already spans every past draw, and restricted to the same scope.

Conditional distribution

Key Default Description
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.

VI-Specific Options

Only used when method = vi. See the VI page.

Key Default Description
vi_iters 25000 Adam ceiling, not a fixed budget — the run stops early once the objective and the estimates have both settled. vi$n_iterations reports what actually ran and vi$converged whether it settled; inspect vi$elbo_trace when it did not (or vi$objective_trace, the penalized objective convergence is actually judged on, when covariate-NN regularization is active). See Why there is no convergence tolerance.
vi_mc_samples 32 Monte-Carlo draws per subject per iteration. This sets the noise floor the settling test measures against, so it decides when the run stops as well as how noisy each step is; too few draws can settle early at the wrong answer. Cost is sublinear — a lower noise floor also settles sooner, so 4× the draws cost ~2.3× the wall time. See Why there is no convergence tolerance for the measured σ / OFV gap at 8 draws.
vi_lr 0.02 Adam learning rate. Lower it if the trace oscillates, or if σ settles above a FOCEI/AGQ fit of the same data.
vi_grad_clip 1e4 Global-L2 gradient clip applied to every Adam step; 0 disables it. See Gradient clipping, and the freeze it prevents.
vi_family full_rank full_rank fits a full posterior covariance per subject; mean_field fits a diagonal one — O(d) rather than O(d²). See Choosing a family.
vi_omega_update closed_form closed_form sets Ω to its exact ELBO maximizer each iteration, keeping it out of the stochastic optimization; adam learns it with everything else. See Ω structure.
vi_sigma_update closed_form As vi_omega_update, for the residual error σ. See How σ is updated.
vi_avg_last final 25% Polyak averaging window, in iterations — the reported estimate is the mean over the final window rather than the last iterate. See Why the estimate is averaged.
vi_eta_grad auto auto uses the analytic Dual2 η-gradient where available and central finite differences where not; analytic and fd pin one route.
vi_kl analytic How the KL(q ‖ N(0, Ω)) half is evaluated — analytic in closed form (only the data term is sampled) or mc. vi$kl reports the route actually taken. See The KL half.
vi_final_ofv none none leaves ofv as NaN, because the ELBO is a lower bound and is not comparable with a FOCEI OFV; laplace computes a comparable one. See The ELBO is not an OFV.
vi_seed 20240704 Seed for the common random numbers. See Reproducibility.

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 covariance 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 (the model is in the sensitivity provider’s scope and gradient is not forced to fd) it uses the gradient-based nlopt_lbfgs; otherwise — a model outside the analytic scope, or gradient = fd, where the outer loop runs on finite differences — it falls back to the derivative-free bobyqa. Models with more than 64 free parameters take nlopt_lbfgs regardless, and mixture models take bobyqa regardless (see the auto default). 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. - Because of that, setting gradient = fd and leaving optimizer at auto moves two factors, not one — on the model behind #1381 the optimizer swap was 4.6 of the 4.98 OFV difference and the gradient only 0.38. ferx emits W_AUTO_OPTIMIZER_FOLLOWS_GRADIENT when that happens; pin optimizer = bobyqa in both arms to measure the gradient alone. See gradient = fd moves the optimizer too. - 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 evals — max(3·(n+1), 50) in general, shortened to max(n+1, 10) for a FOCE/FOCEI gradient-based optimizer (L-BFGS, SLSQP, MMA; not Laplace/AGQ, which have reachable outer_xtol/outer_ftol stops) once the fit has improved on its first evaluation, unless its steps are growing — each of the last two at least 1.5× the one before (#1530). A gradient fit can otherwise spin on a flat objective whose outer gradient does not vanish; the growth veto keeps the short window off a fit still creeping off a saddle, where the OFV is flat but the steps are expanding. 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 derivative-free BOBYQA, 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.

Covariate-NN regularization (nn_l2, nn_smooth)

Two optional penalties regularize the small MLP behind a [covariate_nn] (DCM) block. Both are NN-only, take a non-negative strength, and default to 0.0 (off) — a strict no-op that leaves non-NN and unregularized-NN fits byte-identical. They are added to the population objective the optimizer minimizes, each with a matching analytic gradient term (and the penalty’s curvature where the optimizer builds a Hessian), and apply across the FOCE-family methods: foce / focei / laplace under every outer optimizer (bobyqa, slsqp, nlopt_lbfgs, mma, trust_region, and the built-in BFGS) and gn / gn_hybrid — and under vi, which folds the same penalty gradient into its Adam step (on the same scale, so a given nn_l2 means the same thing under vi as under focei). saem, imp, impmap and bayes do not apply them: setting either key with one of those as the final stage warns and fits unregularized. For a vi DCM the penalty does double duty — it breaks the network’s permutation/scale weight symmetry, which lets the fit settle and early-stop instead of running the full vi_iters ceiling.

Key Range Default Description
nn_l2 float ≥ 0 0.0 L2 (weight-decay) strength. Adds nn_l2 · Σ wᵢ² over the network’s weight matrices only (bias terms are left free). Shrinks all weights toward 0, damping spurious covariate-driven variation whose magnitude otherwise grows with network capacity.
nn_smooth float ≥ 0 0.0 Smoothness (curvature) strength. Penalizes the finite-difference 2nd derivative C = f(x+h) − 2·f(x) + f(x−h) of each output along every input’s marginal partial-dependence curve, i.e. nn_smooth · Σ C². The curve is swept in the network’s input space ((x − center) / scale, the values the MLP actually sees) in steps of a quarter of the input’s observed sd, other inputs held at their median; wide-range inputs get a coarser step over the full range rather than a truncated one. Punishes high-frequency wiggles in the learned covariate→modulator curves while leaving monotone slopes free — it targets curvature, not slope. The grid mirrors the partial-dependence curves used in downstream diagnostics, so it smooths exactly what those plots show.

Motivation: on a null covariate structure (no true effect) an unregularized DCM invents spurious covariate-driven variation that scales with net capacity, with sharp wiggles at higher capacity. nn_l2 shrinks all weights; nn_smooth targets the wiggles directly. The two are independently toggleable.

Reported OFV is unpenalized. The optimizer minimizes the penalized objective, but fit$ofv — and therefore AIC/BIC — report the unpenalized −2·log-likelihood, so a regularized DCM’s information criteria stay directly comparable to an analytic-covariate model’s. The optimizer_trace ofv column, the checkpoint and the verbose Eval/Iter lines carry that same unpenalized value. What is penalized: final_gradient (the gradient of the objective the fit converged on), the ranking of multi-start candidates and of the gn_hybrid phases (so a heavier λ cannot lose to a start that merely fits the data better), and the convergence gates. Standard errors are still the curvature of the unpenalized objective at the regularized optimum — see the covariate-NN page for how to read NN-weight SEs under a strong λ. (Set via the model file [fit_options] block or, in ferx-r, ferx_fit(settings = list(nn_l2 = ..., nn_smooth = ...)) — the same code path.)

[fit_options]
  method    = focei
  nn_l2     = 1e-2
  nn_smooth = 1e-1

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.