ferx_fit.RdFits a NLME model using FOCE or FOCEI estimation with a Rust backend powered by automatic differentiation (Enzyme).
ferx_fit(model, data = NULL, ...)
# S3 method for class 'ferx_model'
ferx_fit(model, data = model$data, ...)Path to a .ferx model file, or a ferx_model object.
Path to a NONMEM-format CSV file (ID, TIME, DV, EVID, AMT, CMT, ...)
Estimation method(s). Either a single string or a character
vector of methods to run in sequence. Each stage is seeded with the
previous stage's converged parameters, and only the final stage produces
the reported covariance/diagnostics. Supported methods: "foce",
"focei", "saem", "gn" (Gauss-Newton / BHHH), or
"gn_hybrid" (Gauss-Newton followed by a FOCEI polish step).
Example chain: c("saem", "focei").
Logical; compute the covariance step for standard errors
Logical; print progress during estimation
Handling of observations below the lower limit of
quantification. NULL (default) keeps whatever the model file
specified; "m3" enables Beal's M3 likelihood (requires a
CENS column in the data, with DV carrying the LLOQ value
on CENS=1 rows); "drop" disables M3.
Number of worker threads for the per-subject parallel loops
in the Rust backend (inner EBE search, SAEM, SIR). NULL (default)
leaves the backend's thread pool at its default of one worker per logical
CPU. Pass an integer to cap the pool — useful on shared machines, in CI,
or to avoid SMT-induced contention (try
parallel::detectCores(logical = FALSE)). The setting is per-call,
so successive fits in the same R session can use different values.
Logical. If TRUE (default), automatically
detects mu-referencing from the model structure for faster and more
accurate convergence. Applies to all estimation methods. Set to
FALSE to disable for comparison purposes. Detection works
automatically for standard parameterizations such as
PARAM = THETA * exp(ETA); unusual parameterizations fall back
silently to zero-centred ETA initialisation with no error. No changes
to the .ferx model file are needed. Check fit$warnings
to see which ETAs were detected.
Inner-loop (per-subject EBE) gradient method.
One of "auto" (default), "ad", or "fd".
The inner optimizer is BFGS; what differs is how the gradient of the individual NLL w.r.t. ETA is computed:
"ad": reverse-mode automatic differentiation via
Enzyme. One forward + one reverse pass per gradient call,
regardless of the number of etas. Requires the package to have
been compiled with autodiff support (see
ferx_rust_autodiff_enabled()) and the model to have an
analytical PK path. ODE-based models currently have no AD
implementation and will silently fall back to FD even if
"ad" is requested.
"fd": central finite differences, 2 * n_eta
forward evaluations per gradient call. Always available.
"auto" (default): use AD whenever the package was
compiled with it and the model has an analytical PK path,
otherwise FD. This is the right choice for almost all uses.
When AD wins. AD's per-gradient cost is roughly independent
of the number of etas, while FD scales linearly. For typical
analytical PK fits we have measured AD ~5x faster per gradient call
on 1-cpt oral (n_eta = 3) and ~1.5-1.8x faster per call on 2-
and 3-cpt infusion. The advantage grows with more random effects,
more observations per subject, and (when implemented) ODE forward
models. On small problems the wall-clock gap is often small
because non-gradient work (NLopt steps, population likelihood
reductions, parallel scheduling) dominates the total fit time.
When to set "fd" explicitly. Mainly for validation
(cross-check that AD and FD converge to the same OFV on a new model),
for reproducibility against a known FD baseline, or to sidestep an
Enzyme codegen pathology on an unusual model structure. Both methods
converge to the same optimum within line-search tolerance on
well-conditioned problems.
Set FERX_TIME_GRADIENTS=1 in the environment to print
per-gradient-call timings at the end of a fit, which is the easiest
way to check which method is faster on your specific model/data.
Logical; run Sampling Importance Resampling after the fit to
produce non-parametric parameter uncertainty intervals. Requires
covariance = TRUE. Tuning knobs (sir_samples,
sir_resamples, sir_seed) still flow through
settings. Default FALSE.
Logical. If TRUE, write a per-iteration CSV
trace to a temporary file and store its path in fit$trace_path.
Pass the result to ferx_read_trace or
ferx_plot_trace to inspect optimizer progress. Default
FALSE.
Logical. If TRUE (default), apply a per-coordinate
scaling layer on top of the existing log/Cholesky parameterization so
that every transformed parameter the outer optimizer sees is O(1).
Why it helps. Population PK parameters often differ by orders
of magnitude on the natural scale (e.g. CL = 0.0015, V =
200, Ka = 0.8). Even after the log transform, the packed
optimizer vector mixes large negatives (log(0.0015) ≈ -6.5) with
large positives (log(200) ≈ 5.3), which leaves the Hessian
poorly conditioned and slows down NLopt / BFGS / Gauss-Newton. With
scale_params = TRUE each coordinate is divided by
|x0[i]| (when |x0[i]| > 0.1, otherwise by 1.0), so
the optimizer works in a near-unit-magnitude space. The scale is
computed once from the initial point and held fixed for the entire run.
When to set FALSE. Mathematically the scaling is
transparent — the OFV, estimates, standard errors, and diagnostics
should match the unscaled fit to numerical tolerance. Set
scale_params = FALSE to:
reproduce the bit-exact iteration trajectory of a pre-scaling
fit (e.g. when comparing against a baseline produced by an older
ferx build),
debug suspected scaling-related issues by toggling the layer off without changing anything else,
validate that scaled and unscaled fits land on the same optimum on a new model.
Applies to all outer optimizers: NLopt FOCE/FOCEI (BOBYQA, SLSQP, L-BFGS, MMA), the hand-rolled BFGS, Gauss-Newton / BHHH, and the SAEM M-step.
Optional named list of estimation-method-specific options
forwarded to the Rust FitOptions. Use this to tune knobs that do
not have a dedicated ferx_fit() argument, without needing a new
wrapper release for each option. Recognized keys include SAEM
(n_exploration, n_convergence, n_mh_steps,
adapt_interval, seed), SIR tuning (sir_samples,
sir_resamples, sir_seed), Gauss-Newton (gn_lambda),
optimizer selection (optimizer — one of "bobyqa" (default),
"slsqp", "lbfgs", "nlopt_lbfgs", "mma",
"bfgs", "trust_region"; also global_search,
global_maxeval), the outer optimizer iteration cap
(maxiter), the inner (per-subject EBE) loop
(inner_maxiter, inner_tol), and the Steihaug CG budget
for optimizer = "trust_region" (steihaug_max_iters).
Values that duplicate a dedicated argument
(method, covariance, verbose,
bloq_method, threads, sir) are rejected — pass
them via the dedicated argument. Unknown keys and malformed values
also raise an error.
Precedence: dedicated ferx_fit() arguments win over
settings, which in turn win over the model file's
[fit_options] block. A warning() is issued whenever a
call-time value overrides a different value from [fit_options].
Inspect fit$model_file_settings and fit$call_settings to
audit the full picture.
Optional path to a .fitrx file. When non-NULL,
ferx_save_fit is invoked on the result so the fit is
persisted to disk in a portable, cross-language bundle (zip of JSON +
CSV) before the function returns. Equivalent to calling
ferx_save_fit(fit, output) after the fit. See the format
reference in the ferx-core docs for the on-disk schema.
Logical. Only meaningful with output. When
TRUE, embeds the input data CSV verbatim inside the
.fitrx bundle so the file is self-contained. Default
FALSE.
A list with components:
Logical; did the optimizer converge
Estimation method used
Number of outer iterations run
Objective function value (-2 log-likelihood)
Akaike Information Criterion
Bayesian Information Criterion
Named numeric vector of fixed effect estimates
Between-subject variability covariance matrix. Row and column
names are the declared ETA names (e.g. "ETA_CL"); fallback is
"OMEGA(1,1)" when names are absent.
Residual error parameter estimates on the
standard-deviation scale. For a proportional component
CV% = sigma * 100; for an additive component the value is in
observation units. Variance is sigma^2 for either type.
Names of the sigma parameters as declared in the
[parameters] block of the .ferx file. Parallel to
sigma.
Per-component error type label, parallel to
sigma: "proportional" or "additive". Combined
error models report c("proportional", "additive") in that order.
Standard errors for theta (NULL if covariance step failed)
Standard errors for omega diagonal
Standard errors for sigma (on the SD scale, like
sigma itself)
Data frame with ID, TIME, DV, PRED, IPRED, CWRES, IWRES,
EBE_OFV, N_OBS; OCC if any subject carries an occasion column; CENS
if any rows are below LLOQ. Per-subject ETAs live in
ebe_etas; per-subject parameter values in
individual_estimates.
Data frame with one row per subject containing the BSV
empirical Bayes estimates: ID plus one column per eta named
after the model's eta declarations (e.g. ETA_CL, ETA_V).
NULL when the model declares no etas.
Data frame with one row per subject containing
each subject's individual parameter values: ID plus one column
per parameter declared in the [individual_parameters] block
(e.g. CL, V, KA). Computed by evaluating the
individual-parameter expressions at the subject's EBE eta with kappa
fixed at zero, so the result reflects each subject's typical value
under the model's covariate effects. For inter-occasion variation see
ebe_kappas.
Character vector of warnings
SIR effective sample size (NULL if SIR not run)
SIR 95\
with columns lower and upper (NULL if SIR not run)
Path to the optimizer trace CSV, or NULL when
optimizer_trace = FALSE. Pass to ferx_read_trace
or ferx_plot_trace.
Number of outer iterations in which too many EBEs were unconverged (step was rejected by the guard).
Worst-case number of unconverged subjects in a single outer iteration.
Total Nelder-Mead fallback invocations across all subjects and iterations.
String: "computed", "failed", or
"not_requested".
Numeric vector of ETA shrinkage per random effect
(1 - SD(eta_hat_k) / sqrt(omega_kk)). NA when
omega_kk = 0 or fewer than 2 subjects.
EPS shrinkage: 1 - SD(IWRES). NA when
fewer than 2 valid residuals.
Total wall-clock time for the fit in seconds.
Model name from the .ferx file. Falls back to
the model file's basename (without extension) when the file declares
no name.
Dataset name, derived as the basename of the data
path with its extension stripped.
The inner-loop gradient method as requested by the
caller (one of "auto", "ad", "fd"). The
resolved method may differ — see gradient_used and the
gradient argument.
The inner-loop gradient method the engine actually
used: "ad", "fd", or "N/A" (derivative-free /
sampling step). When gradient = "auto" this records which
branch resolved at fit time; for ODE models with gradient = "ad"
this shows the silent fallback to "fd". NA on older
engine binaries that did not populate this field.
Engine label for the inner-loop gradient
("Enzyme AD", "finite differences", or "N/A"). Prefer
gradient_used for display.
Engine label for the population-level optimizer's gradient ("Enzyme AD", "finite differences", or "N/A").
ferx-core library version string.
Full parameter covariance matrix as a named numeric
matrix (params × params). Row/column names use declared variable names
("TVCL", "ETA_CL", "EPS_PROP"); fallback is
"OMEGA(1,1)" / "SIGMA(1)". NULL when covariance
step was not run or failed. Use ferx_cor_matrix to
inspect correlations.
Numeric vector of eigenvalues of the correlation matrix
of estimated (non-fixed) parameters, sorted descending. Computed by the
Rust backend. NULL when the covariance step was not run, failed,
or fewer than two free parameters exist.
Ratio of the largest to smallest eigenvalue of the
correlation matrix of estimated parameters. Values above 1000 are flagged
as potentially ill-conditioned and also appear in warnings.
Inf when the smallest eigenvalue is non-positive. NULL
when eigenvalues is NULL.
Data frame with Shapiro-Wilk normality test for
each ETA: columns eta, W, p_val, flag.
A [!] flag appears when p < 0.05. NULL only
when ebe_etas is missing or has no eta columns. When an ETA
has fewer than 3 finite values (e.g. tiny populations or fixed
etas), its row is still returned but W and p_val are
NA.
IOV variance matrix for kappa parameters (NULL if no IOV).
Names of kappa (IOV) parameters (NULL if no IOV).
Standard errors for kappa parameters: length d
(diagonal kappa) or d*(d+1)/2 (block_kappa, lower-triangle order).
NULL if covariance step not run or no IOV.
Named list returned by the Rust engine, derived
from the parsed CompiledModel so it reflects exactly what
ferx-core ran. Fields: theta_names (character vector of
population parameter names), model_type (short label such as
"1-cpt oral", "ODE", or NULL when the structural
form is not one of the known analytical PK families), iiv
(omega names), iov (kappa names), residual (error type
string). Use ferx_model_inspect to view this before or
after fitting.
Shrinkage values for kappa EBEs (NULL if no IOV).
Data frame with columns ID, OCC, and one
column per kappa parameter containing per-subject per-occasion kappa
EBEs. ID carries the original subject identifier (matching
sdtab); OCC carries the labeled occasion in first-seen
order. NULL when no IOV.
Parameter-level correlation matrix for block omega.
Uses the bivariate lognormal formula for lognormal pairs and falls back
to the eta-level formula for additive or unknown parameterisations.
NULL when omega is diagonal (no off-diagonal elements to report).
Same as omega_param_corr but for the
IOV kappa block. NULL when kappa is diagonal or absent.
Character vector of ETA parameter names as declared in the
model file (e.g. "ETA_CL", "ETA_V"). Used to label OMEGA
rows in ferx_estimates() and print.ferx_fit().
Logical vector of length n_eta; TRUE
when the eta is lognormally parameterised (THETA * exp(ETA)),
FALSE for additive or unknown parameterisations. NULL when
the engine does not supply the information.
Named list of the effective settings passed to Rust,
with values typed back to their natural R types (logical, numeric, or
character). Includes merged defaults such as optimizer_trace and
scale_params. Empty list when no settings were supplied.
Named list of raw key = value pairs from
the model file's [fit_options] block, as character strings
(before Rust type coercion). Empty list when no [fit_options]
block is present. Use this alongside call_settings to audit which
values the model file requested and which were overridden at the R call
site. A warning() is issued automatically for each key that
appears in both sources with a different value.
There are two equivalent ways to supply the model and data:
1. Path strings (classic style):
2. ferx_model object (pipe style):
"data.csv" |> ferx_model("pk.ferx") |> ferx_fit()Both dispatch to the same Rust backend. The ferx_model form is
convenient when combined with ferx_set_section to modify model
options in the same chain (see Examples). The data path stored in the
ferx_model object can always be overridden by supplying data
explicitly to ferx_fit().
Estimation method — pass a single method or a vector to chain methods in sequence (each stage seeds the next with its converged parameters):
ferx_fit(m, d, method = "focei")
ferx_fit(m, d, method = c("saem", "focei")) # SAEM warm-start, FOCEI polishStandard errors — the covariance step is on by default:
ferx_fit(m, d, covariance = TRUE) # default — produces SE /
ferx_fit(m, d, covariance = FALSE) # skip for speed during developmentParallelism — cap the Rust thread pool:
ferx_fit(m, d, threads = parallel::detectCores(logical = FALSE))BLOQ handling:
ferx_fit(m, d, bloq_method = "m3") # Beal M3 likelihood (needs CENS column)
ferx_fit(m, d, bloq_method = "drop") # discard BLOQ rowsGradient method for the inner EBE loop:
ferx_fit(m, d, gradient = "auto") # default: AD when available, else FD
ferx_fit(m, d, gradient = "ad") # force automatic differentiation
ferx_fit(m, d, gradient = "fd") # force finite differencesOptimizer trace — write a per-iteration CSV for convergence diagnostics:
fit <- ferx_fit(m, d, optimizer_trace = TRUE)
ferx_plot_trace(fit)settingsThe settings argument forwards a named list of low-level options to
the Rust backend without requiring a new wrapper release. Arguments with
dedicated parameters (e.g. method, covariance) cannot be
duplicated in settings — pass them via the named argument.
Outer optimizer selection:
ferx_fit(m, d, settings = list(optimizer = "bobyqa")) # default
ferx_fit(m, d, settings = list(optimizer = "slsqp")) # gradient-based
ferx_fit(m, d, settings = list(optimizer = "lbfgs"))
ferx_fit(m, d, settings = list(optimizer = "bfgs"))
ferx_fit(m, d, settings = list(optimizer = "trust_region"))
ferx_fit(m, d, settings = list(optimizer = "mma"))Iteration cap and inner loop:
SAEM tuning:
ferx_fit(m, d, method = "saem", settings = list(
n_exploration = 200,
n_convergence = 400,
n_mh_steps = 3,
seed = 42L
))SIR uncertainty (requires sir = TRUE, covariance = TRUE):
ferx_fit(m, d, sir = TRUE, settings = list(
sir_samples = 2000L,
sir_resamples = 500L,
sir_seed = 1L
))Global search before local refinement:
Trust-region CG budget:
ferx_fit() returns a ferx_fit object. All of the following
work both as standalone calls and as the tail of a |> pipe:
fit |> print() # full parameter table
fit |> summary() # compact diagnostic summary
fit |> ferx_estimates() # tidy data frame: theta / omega / sigma + SE /
fit |> ferx_cor_matrix() # parameter correlation matrix (needs covariance = TRUE)
fit |> ferx_model_inspect() # model structure auto-derived by the engine
fit |> ferx_plot_trace() # convergence trace (needs optimizer_trace = TRUE)Diagnostics data frame (PRED, IPRED, CWRES, ETAs, …) lives in
fit$sdtab and can be used directly:
fit$sdtab
fit$ebe_etas
fit$individual_estimatesex <- ferx_example("warfarin")
fit <- ferx_fit(ex$model, ex$data, method = "gn", covariance = FALSE)
#> Warning: Model file [fit_options] sets `method = foce` but ferx_fit() argument overrides it with `gn`. The call-time value will be used.
#> Warning: Model file [fit_options] sets `covariance = true` but ferx_fit() argument overrides it with `false`. The call-time value will be used.
#> Mu-referencing detected for: ETA_CL, ETA_KA, ETA_V
#> Negative IWRES autocorrelation detected (Durbin-Watson = 2.61, lag-1 r = -0.37).
#> Possible over-parameterisation or misspecified residual error model.
fit$theta
#> TVCL TVV TVKA
#> 0.1327364 7.6905899 0.7572037
fit$ofv
#> [1] -280.1829
if (FALSE) { # \dontrun{
ex <- ferx_example("warfarin")
# ── Classic style (path strings) ─────────────────────────────────────────
# Minimal call — FOCEI with default optimizer, covariance step on
fit <- ferx_fit(ex$model, ex$data)
# Named arguments — fully equivalent
fit <- ferx_fit(model = ex$model, data = ex$data, method = "focei")
# All common options at once
fit <- ferx_fit(ex$model, ex$data,
method = "focei",
covariance = TRUE,
threads = 4L,
gradient = "auto",
verbose = TRUE,
scale_params = TRUE
)
# SAEM warm-start, then FOCEI polish
fit <- ferx_fit(ex$model, ex$data, method = c("saem", "focei"))
# Fine-tune via settings
fit <- ferx_fit(ex$model, ex$data,
method = "focei",
settings = list(
optimizer = "slsqp",
maxiter = 500L,
inner_maxiter = 100L,
inner_tol = 1e-6
)
)
# BLOQ (M3 method)
bloq <- ferx_example("warfarin_bloq")
fit <- ferx_fit(bloq$model, bloq$data, method = "focei", bloq_method = "m3")
# SIR parameter uncertainty
fit <- ferx_fit(ex$model, ex$data,
sir = TRUE,
settings = list(sir_samples = 2000L, sir_resamples = 500L)
)
# ── Pipe style (ferx_model object) ───────────────────────────────────────
# Basic pipe — data flows into ferx_model() which bundles the model path
ex$data |>
ferx_model(ex$model) |>
ferx_fit(method = "focei", covariance = TRUE)
# Modify a model section, then fit
ex$data |>
ferx_model(ex$model) |>
ferx_set_section("fit_options", c(
" method = focei",
" maxiter = 500",
" covariance = true"
)) |>
ferx_fit()
# Full pipeline including post-fit outputs
ex$data |>
ferx_model(ex$model) |>
ferx_set_section("fit_options", c(" method = focei")) |>
ferx_fit(covariance = TRUE, threads = 4L,
settings = list(optimizer = "slsqp")) |>
summary()
# Inspect then continue (ferx_get_section returns the ferx_model invisibly)
ex$data |>
ferx_model(ex$model) |>
ferx_get_section("parameters") |>
ferx_fit() |>
ferx_estimates()
# Override data stored in ferx_model at fit time
# (substitute the path to your own dataset for "other_cohort.csv")
m <- ferx_model(ex$data, ex$model)
ferx_fit(m, data = "other_cohort.csv", method = "focei")
# ── Post-fit outputs ─────────────────────────────────────────────────────
fit <- ferx_fit(ex$model, ex$data, covariance = TRUE, optimizer_trace = TRUE)
summary(fit) # compact diagnostic table
ferx_estimates(fit) # tidy data frame with SE and %RSE
ferx_cor_matrix(fit) # parameter correlation matrix
ferx_model_inspect(fit) # model structure auto-derived by the engine
ferx_plot_trace(fit) # convergence plot (optimizer_trace = TRUE required)
fit$sdtab # per-observation diagnostics (PRED, IPRED, CWRES, …)
fit$ebe_etas # per-subject empirical Bayes ETAs
fit$individual_estimates # per-subject individual PK parameters
fit$eigenvalues # sorted eigenvalues of parameter correlation matrix
fit$condition_number # > 1000 flags potential ill-conditioning
# Covariance diagnostics in a pipe
ferx_fit(ex$model, ex$data, covariance = TRUE) |> ferx_cor_matrix()
} # }