Model workflow: inspect, edit, fit

Overview

This vignette covers two complementary workflows:

1. Path-based workflow — pass file paths directly to ferx_fit(). Straightforward for one-off fits.
2. Pipe-based workflow — bundle paths into a ferx_model object and chain steps with |>. Cleaner when you need to tweak options before fitting or chain post-fit steps.

Both workflows produce the same ferx_fit result.

Entry point: ferx_model(data, model) takes the data file first, which lets a data object flow naturally into a pipeline:

ex$data |> ferx_model(ex$model) |> ferx_fit() |> summary()

You can also call it positionally or by name:

ferx_model(ex$data, ex$model) |> ferx_fit() |> summary()
ferx_model(data = ex$data, model = ex$model) |> ferx_fit() |> summary()

Migration note: Earlier versions used the order ferx_model(model, data). Old-style positional calls (ferx_model("pk.ferx") or ferx_model("pk.ferx", "data.csv")) are detected by the .ferx extension and auto-corrected with a deprecation warning; this shim will be removed in a future release. Calls that pass data by name (ferx_model("pk.ferx", data = "data.csv")) keep working unchanged.

---

Path-based workflow

Step 1 — Get a model file

ferx_example() returns paths to the bundled example models and their data:

ex <- ferx_example("warfarin")
ex$model  # path to warfarin.ferx
#> [1] "/home/runner/work/_temp/Library/ferx/examples/models/warfarin.ferx"
ex$data   # path to warfarin.csv
#> [1] "/home/runner/work/_temp/Library/ferx/examples/data/warfarin.csv"

In a real workflow you would supply your own .ferx file here.

Step 2 — Inspect before fitting

ferx_model_inspect() parses the model file and prints a concise summary of its structure without running the optimizer. Use this to verify that ferx interprets your model the way you expect before committing to a long run.

ferx_model_inspect(ex$model)
#> Model structure (warfarin.ferx)
#>   Structural:  1-cpt oral  (TVCL, TVV, TVKA)
#>   IIV:         ETA_CL, ETA_V, ETA_KA
#>   IOV:         none
#>   Residual:    proportional

The printed summary has four lines:

Line	Meaning
Structural	Detected PK model type and population parameter (theta) names
IIV	Inter-individual variability (omega) parameter names
IOV	Inter-occasion variability (kappa) parameter names, or none
Residual	Error model type (proportional, additive, or combined)

ferx_model_inspect() returns the structure as a named list (invisibly), which is useful for programmatic checks:

s <- ferx_model_inspect(ex$model)
#> Model structure (warfarin.ferx)
#>   Structural:  1-cpt oral  (TVCL, TVV, TVKA)
#>   IIV:         ETA_CL, ETA_V, ETA_KA
#>   IOV:         none
#>   Residual:    proportional
s$theta_names   # population parameter names
#> [1] "TVCL" "TVV"  "TVKA"
s$model_type    # short label for the PK model
#> [1] "1-cpt oral"
s$iiv           # IIV parameter names
#> [1] "ETA_CL" "ETA_V"  "ETA_KA"
s$residual      # error model type
#> [1] "proportional"

If the printed summary does not match your expectations — for example the model type shows NULL or a parameter name is missing — fix the .ferx file before proceeding. This catches typos and structural mistakes cheaply, before the optimizer runs.

Step 3 — Fit

Once the structure looks correct, pass the same paths to ferx_fit():

fit <- ferx_fit(ex$model, ex$data)
fit

Step 4 — Re-inspect from the fitted result

After fitting you can call ferx_model_inspect() directly on the ferx_fit object returned by ferx_fit() — no need to supply the model path again:

ferx_model_inspect(fit)

The output is identical to the pre-fit inspection. This is convenient when you want to confirm the structure inside a script or report without keeping the model path in scope.

---

Pipe-based workflow

ferx_model() bundles model and data paths into a single object that flows through a |> chain. This is the preferred style when you need to adjust options before fitting or chain several post-fit steps.

Entry point

ex <- ferx_example("warfarin")
m <- ferx_model(ex$data, ex$model)
m   # prints path, data path, and structure summary
#> ferx_model
#>   Model: /home/runner/work/_temp/Library/ferx/examples/models/warfarin.ferx
#>   Data:  /home/runner/work/_temp/Library/ferx/examples/data/warfarin.csv
#>   ---
#>   Structural:  1-cpt oral  (TVCL, TVV, TVKA)
#>   IIV:         ETA_CL, ETA_V, ETA_KA
#>   IOV:         none
#>   Residual:    proportional

Minimal pipe — inspect, fit, summarise

ferx_fit() dispatches on the ferx_model object and picks up $data automatically.

ferx_model(ex$data, ex$model) |>
  ferx_fit(method = "focei", covariance = TRUE) |>
  summary()

Peek at a section mid-pipe

ferx_get_section() prints the named section to the console and passes the ferx_model object through unchanged, so the pipe continues:

ferx_model(ex$data, ex$model) |>
  ferx_get_section("parameters")
#> # [parameters]
#>   theta TVCL(0.134, 0.001, 10.0)
#>   theta TVV(8.1, 0.1, 500.0)
#>   theta TVKA(1.0, 0.01, 50.0)
#> 
#>   omega ETA_CL ~ 0.07
#>   omega ETA_V  ~ 0.02
#>   omega ETA_KA ~ 0.40
#> 
#>   sigma PROP_ERR ~ 0.01

# Print [parameters] then fit without interrupting the chain
fit <- ferx_model(ex$data, ex$model) |>
  ferx_get_section("parameters") |>
  ferx_fit(method = "focei")

Override fit options before fitting

ferx_set_section() rewrites a section on disk and returns the ferx_model unchanged. Important: it edits the file in place — copy first if you do not want to modify the original.

model_copy <- file.path(tempdir(), "warfarin.ferx")
file.copy(ex$model, model_copy)

fit <- ferx_model(ex$data, model_copy) |>
  ferx_set_section("fit_options", c(
    "  method     = focei",
    "  maxiter    = 500",
    "  covariance = true"
  )) |>
  ferx_fit()

summary(fit)
ferx_model_inspect(fit)    # no path needed after fitting
ferx_cor_matrix(fit)       # parameter correlation matrix
ferx_plot_trace(fit)       # OFV + gradient norm convergence plot

Validate before a long run

ferx_check_init() runs only 5 iterations and returns a trace and a summary table. Use it to confirm the OFV is dropping before committing to a long estimation:

chk <- ferx_check_init(ex$model, ex$data, method = "focei")
chk$summary      # ofv_start, ofv_end, ofv_drop — is OFV decreasing?
ferx_plot_trace(chk$fit)   # visual check: first few iterations

If ofv_drop is small or negative, the initial values are likely poor — adjust them in [parameters] or widen the parameter bounds before proceeding.

Multi-stage method chain

Pass a character vector to method to chain stages — each stage starts from the previous stage’s converged parameters:

fit <- ferx_model(ex$data, ex$model) |>
  ferx_fit(method = c("saem", "focei"), covariance = TRUE)

Simulate and predict from fitted parameters

# Simulate 100 replicates using individual parameters from the fit
sim <- ferx_simulate(ex$model, ex$data, n_sim = 100, seed = 42, fit = fit)

# Population predictions using fitted theta (eta = 0)
pred <- ferx_predict(ex$model, ex$data, fit = fit)

Post-fit diagnostics

All diagnostic outputs are fields on the ferx_fit object:

head(fit$sdtab)                # PRED, IPRED, CWRES, IWRES per observation
head(fit$ebe_etas)             # empirical Bayes ETAs per subject
head(fit$individual_estimates) # individual PK parameters (CL, V, KA, ...)
fit$eta_normality              # Shapiro-Wilk normality test per ETA

ferx_estimates() collects all estimated parameters (theta, omega, sigma) into a tidy data frame with SEs, %RSE, and 95 % CIs. Run it after any fit that used covariance = TRUE. The transform column reflects the parameterisation reported by the backend ("identity", "log", "logit", or "variance" for omega rows):

ferx_estimates(fit)
#   param           transform  estimate     se  rse_pct  lower_95  upper_95  ...
#   TVCL            identity     0.134  0.0108     8.1     0.113     0.155
#   TVV             identity     7.94   0.420      5.3     7.12      8.76
#   TVKA            identity     1.07   0.160     14.9     0.756     1.38
#   OMEGA(ETA_CL)   variance     0.042  0.012     28.6    0.019     0.066
#   OMEGA(ETA_V)    variance     0.020  0.009     44.9    0.002     0.037
#   OMEGA(ETA_KA)   variance     0.157  0.048     30.6    0.063     0.251
#   PROP_ERR        proportional 0.011  0.001      9.1    0.009     0.013

ferx_cor_matrix() converts the parameter covariance matrix to correlations. A value near ±1 between two parameters signals a structural identifiability problem worth investigating:

ferx_cor_matrix(fit)

Goodness-of-fit plots from fit$sdtab with ggplot2:

library(ggplot2)

sdtab <- fit$sdtab

# DV vs PRED
ggplot(sdtab, aes(PRED, DV)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, colour = "steelblue") +
  labs(title = "DV vs PRED", x = "Population prediction", y = "Observed")

# DV vs IPRED
ggplot(sdtab, aes(IPRED, DV)) +
  geom_point(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, colour = "steelblue") +
  labs(title = "DV vs IPRED", x = "Individual prediction", y = "Observed")

# CWRES vs TIME
ggplot(sdtab, aes(TIME, CWRES)) +
  geom_point(alpha = 0.4) +
  geom_hline(yintercept = c(-2, 0, 2), linetype = c("dashed", "solid", "dashed"),
             colour = "steelblue") +
  labs(title = "CWRES vs TIME")

# CWRES vs PRED
ggplot(sdtab, aes(PRED, CWRES)) +
  geom_point(alpha = 0.4) +
  geom_hline(yintercept = c(-2, 0, 2), linetype = c("dashed", "solid", "dashed"),
             colour = "steelblue") +
  labs(title = "CWRES vs PRED")

ETA-covariate correlations — identifies covariates worth testing in a formal covariate search. Pass the original data frame (not the path); any numeric column that is constant within a subject and not a standard NONMEM column (TIME, DV, AMT, EVID, MDV, CMT, RATE) is treated as a covariate:

# The warfarin dataset has no covariate columns. The following example uses
# two_cpt_oral_cov, which has body weight (WT) and creatinine clearance (CRCL).
ex_cov <- ferx_example("two_cpt_oral_cov")
fit_cov <- ferx_fit(ex_cov$model, ex_cov$data)
data_cov <- read.csv(ex_cov$data)
ferx_eta_cov(fit_cov, data_cov)
#   eta      covariate    r   p_val  flag
#   ETA_CL   WT        0.38  0.008   [!]
#   ETA_CL   CRCL      0.41  0.003   [!]
#   ETA_V1   WT        0.29  0.047
# [!] flags |r| >= 0.3; flag covariates worth formal testing

---

Parameter transforms

ferx supports log and logit transforms for both thetas and ETAs directly in the model DSL. When a transform is in use, print() and ferx_estimates() display results on the natural scale alongside the estimation-scale values.

# Log-normal CL (typical in PK — ETA_CL is additive on the log scale)
CL = TVCL * exp(ETA_CL)  # TVCL is on the natural scale; BSV is log-normal

# Logit-transformed bioavailability (constrained to (0, 1))
F = inv_logit(logit(THETA_F) + ETA_F)

See vignette("parameter-transforms") for worked examples with bundled models.

---

Saving and reloading fits

Long estimation runs can be saved to disk and reloaded in a later session with ferx_save_fit() / ferx_load_fit(). The saved file is a compressed RDS and preserves the full ferx_fit object including sdtab, ebe_etas, trace path, and all metadata.

# Save after a long run
fit <- ferx_fit(ex$model, ex$data, method = "focei", covariance = TRUE)
ferx_save_fit(fit, "warfarin_focei.rds")

# Later — restore and continue working
fit2 <- ferx_load_fit("warfarin_focei.rds")
identical(fit$theta, fit2$theta)   # TRUE
ferx_estimates(fit2)

The loaded object is indistinguishable from the original: all downstream functions (ferx_estimates(), ferx_cor_matrix(), ferx_plot_trace(), ferx_simulate(), etc.) work without modification.

Tip: save immediately after fitting so a crash or session restart does not lose results. Use a naming convention that records the run number, method, and date:

ferx_save_fit(fit, sprintf("run%02d_focei_%s.rds", run_number, Sys.Date()))

---

Summary

Path-based

ex <- ferx_example("warfarin")
ferx_model_inspect(ex$model)          # verify structure before fitting
fit <- ferx_fit(ex$model, ex$data)
ferx_model_inspect(fit)               # re-inspect post-fit (no path needed)

Pipe-based

ex <- ferx_example("warfarin")

# Inspect, fit, summarise in one chain
ferx_model(ex$data, ex$model) |>
  ferx_get_section("parameters") |>  # optional: peek before fitting
  ferx_fit(method = "focei", covariance = TRUE) |>
  summary()

# Edit options then fit (copy model file first)
model_copy <- file.path(tempdir(), "warfarin.ferx")
file.copy(ex$model, model_copy)

ferx_model(model_copy, data = ex$data) |>
  ferx_set_section("fit_options", c(
    "  method     = focei",
    "  maxiter    = 500",
    "  covariance = true"
  )) |>
  ferx_fit() |>
  summary()