Writing and managing model files

Where you are

The dataset from Preparing and checking the analysis dataset is checked. The next step is a model file: a plain-text .ferx file that describes the parameters, the structural model, the residual error and the default fit settings. Keeping the model in a file separates it from R code, makes it easy to version, and lets the same file run from R or from the ferx command line.

This chapter shows how to read, create, check and edit model files from R. It uses the base model of the workflow thread, two_cpt_oral_base: a two-compartment oral model without covariates, for the same dataset.

The data

The base model ships with the same dataset as Preparing and checking the analysis dataset:

base <- ferx_example("two_cpt_oral_base")
basename(unlist(base))
#> [1] "two_cpt_oral_base.ferx"       "two_cpt_oral_cov.csv"        
#> [3] "two_cpt_oral_base.ferxsearch"

Reading a model file

ferx_model_show() prints the file:

ferx_model_show(base$model)
#> # model: two_cpt_oral_base.ferx 
#> # Two-compartment oral PK model -- the covariate-free base of two_cpt_oral_cov.
#> #
#> # Same structure, same dataset (two_cpt_oral_cov.csv, which carries WT and
#> # CRCL), but no covariate effect on any parameter. It exists to be the *start*
#> # of a search: a covariate search whose base model already multiplies CL by
#> # (WT/70)^THETA_WT is asking whether to add an effect the model has, and
#> # allometric scaling on top of an inline (WT/70) term would count body size
#> # twice.
#> #
#> # Paired tools:
#> #   ferx_covsearch()  -- two_cpt_oral_base.ferxsearch searches WT and CRCL on
#> #                        every parameter carrying an eta
#> #   ferx_allometry()  -- (WT/70)^0.75 on CL and Q, (WT/70)^1.0 on V1 and V2
#> #
#> # The `[covariates]` block below is what makes `@CONTINUOUS` resolvable: a
#> # search space written as `COVARIATE?(@IIV, @CONTINUOUS, [pow, lin])` expands
#> # against the declared columns, not against every column in the dataset.
#> 
#> [parameters]
#>   theta TVCL(5.0, 0.1, 100.0)
#>   theta TVV1(50.0, 1.0, 500.0)
#>   theta TVQ(10.0, 0.1, 100.0)
#>   theta TVV2(100.0, 1.0, 500.0)
#>   theta TVKA(1.2, 0.01, 10.0)
#> 
#>   omega ETA_CL ~ 0.10
#>   omega ETA_V1 ~ 0.10
#>   omega ETA_Q  ~ 0.05
#>   omega ETA_V2 ~ 0.05
#>   omega ETA_KA ~ 0.15
#> 
#>   sigma PROP_ERR ~ 0.02 (sd)
#> 
#> [individual_parameters]
#>   CL = TVCL * exp(ETA_CL)
#>   V1 = TVV1 * exp(ETA_V1)
#>   Q  = TVQ  * exp(ETA_Q)
#>   V2 = TVV2 * exp(ETA_V2)
#>   KA = TVKA * exp(ETA_KA)
#> 
#> [structural_model]
#>   pk two_cpt_oral(cl=CL, v1=V1, q=Q, v2=V2, ka=KA)
#> 
#> [covariates]
#>   WT   continuous
#>   CRCL continuous
#> 
#> [error_model]
#>   DV ~ proportional(PROP_ERR)
#> 
#> [fit_options]
#>   method     = focei
#>   maxiter    = 500
#>   covariance = true
#>   # Finite differences rather than the default `auto`. This was a workaround
#>   # for ferx-core #1290: the analytic-gradient path used to stall at the
#>   # initial estimates on this model and dataset (`NLopt stopped: Failure` on
#>   # the first evaluation, no parameter moves), which a search then had to
#>   # reject as an init-stalled base fit. That is fixed, and both paths now
#>   # converge; `fd` is kept explicit so the example's published numbers do not
#>   # move with the gradient default.
#>   gradient   = fd

A model file is a sequence of [blocks]. These five appear in almost every model:

Block	What it declares	ferx-core reference
`[parameters]`	Fixed effects `theta NAME(initial, lower, upper)`, between-subject variances `omega NAME ~ value`, residual error `sigma` (`(sd)` marks an initial value given as a standard deviation)	parameters
`[individual_parameters]`	Each subject’s parameters as expressions of thetas, etas and covariates	individual parameters
`[structural_model]`	A built-in PK model (`pk two_cpt_oral(...)`), an ODE system, or a compartment-free model	structural model
`[error_model]`	The residual error model for `DV`	error model
`[fit_options]`	Default estimation settings; arguments to `ferx_fit()` override them	fit options

What an expression may contain

The expressions in [individual_parameters] and [odes] take a fixed set of functions. Rather than reproduce the list, ask the engine: an unrecognised name is a parse error that names the whole supported set, so this stays true at whatever build you are on.

fn_probe <- file.path(book_tempdir("model-files"), "unknown_function.ferx")
original <- readLines(ferx_example("warfarin")$model)
probed <- sub("  CL = TVCL * exp(ETA_CL)", "  CL = TVCL * tanh(ETA_CL)",
              original, fixed = TRUE)
stopifnot(!identical(probed, original))   # the line must actually have been replaced
writeLines(probed, fn_probe)
cat(tryCatch(ferx_fit(fn_probe, ferx_example("warfarin")$data, verbose = FALSE),
             error = function(e) sub("^Error parsing model: ", "", conditionMessage(e))),
    fill = TRUE)
#> unknown function `tanh`. Supported: exp, log, ln, sqrt, abs, floor, ceil, round, logit, inv_logit, expit, min(a, b), max(a, b), clamp(x, lo, hi), present(x) (conditions only) [E_PARSE]

logit() and inv_logit() are how a parameter is kept inside (0, 1) – Variability: random effects, residual error and IOV gives the forms and Initial estimates and a first fit shows how they are reported. Names match case-insensitively, so EXP(...) works as well as exp(...).

This model also has a [covariates] block, which declares WT and CRCL. The other blocks are introduced in the chapters where they are used:

Block	Chapter
`[data]`, `[data_selection]`	Preparing and checking the analysis dataset
`[odes]`, `[initial_conditions]`, `[scaling]`	Structural models: analytical and ODE
`[covariates]`, `[covariate_model]`	Covariate modeling
`[mixture]`	Variability: random effects, residual error and IOV
`[derived]`, `[output]`	Tables and figures
`[simulation]`	Simulating scenarios
`[binary_model]`	Binary endpoints
`[event_model]`	Time-to-event models
`[adaptive_dosing]`	Adaptive dosing and TDM strategies
`[diffusion]`, `[covariate_nn]`	Experimental features
`[markov_model]`	not available in ferx-r (Function, option and example index)

Checking a model before fitting

ferx_model_inspect() parses the file without fitting and reports how ferx interprets it:

structure_info <- ferx_model_inspect(base$model)
#> Model structure (two_cpt_oral_base.ferx)
#>   Structural:  2-cpt oral  (TVCL, TVV1, TVQ, TVV2, TVKA)
#>   IIV:         ETA_CL, ETA_V1, ETA_Q, ETA_V2, ETA_KA
#>   IOV:         none
#>   Residual:    proportional
structure_info$theta_names
#> [1] "TVCL" "TVV1" "TVQ"  "TVV2" "TVKA"

ferx_model_validate() runs the engine’s parser and structural checks. With data, it also runs the checks that need the dataset, such as whether declared covariates are present:

check <- ferx_model_validate(base$model, data = base$data)
#> Validating: two_cpt_oral_base.ferx 
#>        data: two_cpt_oral_cov.csv 
#> 
#> Sections present:
#>   parameters                     [ok]
#>   individual_parameters          [ok]
#>   structural_model               [ok]
#>   error_model                    [ok]
#>   covariates                     [ok] (optional)
#>   fit_options                    [ok] (optional)
#> 
#> Result: VALID
check$ok
#> [1] TRUE

The report is printed, and the invisible return value holds ok, the model and data names, and a diagnostics data frame. To see a failing check, misspell a block name in a copy of the model:

model_dir <- book_tempdir("model-files")
typo <- file.path(model_dir, "typo.ferx")
writeLines(sub("^\\[fit_options\\]", "[fit_option]", readLines(base$model)), typo)

typo_check <- ferx_model_validate(typo)
#> Validating: typo.ferx 
#> 
#> Sections present:
#>   parameters                     [ok]
#>   individual_parameters          [ok]
#>   structural_model               [ok]
#>   error_model                    [ok]
#>   covariates                     [ok] (optional)
#>   fit_option                     [unknown section]
#> 
#> Result: INVALID
#>   * ERROR E_UNKNOWN_BLOCK [fit_option:51]: Unknown block `[fit_option]` — did you mean `[fit_options]`. Valid blocks: adaptive_dosing, binary_model, covariate_model, covariate_nn, covariates, data, data_selection, derived, diffusion, error_model, event_model, fit_options, individual_parameters, initial_conditions, mixture, odes, output, parameters, priors, scaling, simulation, structural_model.
#>       hint: did you mean `[fit_options]`?
typo_check$ok
#> [1] FALSE
typo_check$diagnostics[, c("severity", "code", "block", "line", "suggestion")]
#>   severity            code      block line                    suggestion
#> 1    error E_UNKNOWN_BLOCK fit_option   51 did you mean `[fit_options]`?

Diagnostic codes are stable identifiers: E_* for errors, W_* for warnings. Function, option and example index lists them all with a one-line meaning, and the ferx-core check report page gives the full text of each.

Starting a new model from a template

ferx_model() can scaffold a new model file from a built-in template. The templates are "1cpt_oral", "1cpt_iv", "2cpt_oral", "2cpt_iv" and "ode". With print = TRUE it only prints the skeleton:

ferx_model(template = "2cpt_oral", print = TRUE)
#> # Two-compartment oral PK model
#> 
#> [parameters]
#>   theta TVCL(5.0, 0.1, 100.0)
#>   theta TVV1(50.0, 1.0, 500.0)
#>   theta TVQ(10.0, 0.1, 100.0)
#>   theta TVV2(100.0, 1.0, 1000.0)
#>   theta TVKA(1.2, 0.01, 10.0)
#> 
#>   omega ETA_CL ~ 0.10
#>   omega ETA_V1 ~ 0.10
#>   omega ETA_Q  ~ 0.05
#>   omega ETA_V2 ~ 0.05
#>   omega ETA_KA ~ 0.15
#> 
#>   sigma PROP_ERR ~ 0.02
#> 
#> [individual_parameters]
#>   CL = TVCL * exp(ETA_CL)
#>   V1 = TVV1 * exp(ETA_V1)
#>   Q  = TVQ  * exp(ETA_Q)
#>   V2 = TVV2 * exp(ETA_V2)
#>   KA = TVKA * exp(ETA_KA)
#> 
#> [structural_model]
#>   pk two_cpt_oral(cl=CL, v1=V1, q=Q, v2=V2, ka=KA)
#> 
#> [error_model]
#>   DV ~ proportional(PROP_ERR)
#> 
#> [fit_options]
#>   method     = focei
#>   maxiter    = 500
#>   covariance = true

With a path it writes the file and returns a ferx_model object. Use edit = FALSE to skip opening an editor, and overwrite = TRUE to replace an existing file:

scaffold <- ferx_model(template = "2cpt_oral", path = file.path(model_dir, "run1.ferx"),
                       edit = FALSE, overwrite = TRUE)
#> Created /tmp/Rtmpe9t14F/model-files/run1.ferx
scaffold
#> ferx_model
#>   Model: /tmp/Rtmpe9t14F/model-files/run1.ferx
#>   Data:  <none>
#>   ---
#>   Structural:  2-cpt oral  (TVCL, TVV1, TVQ, TVV2, TVKA)
#>   IIV:         ETA_CL, ETA_V1, ETA_Q, ETA_V2, ETA_KA
#>   IOV:         none
#>   Residual:    proportional
ferx_model_validate(scaffold$model, data = base$data)$ok
#> Validating: run1.ferx 
#>        data: two_cpt_oral_cov.csv 
#> 
#> Sections present:
#>   parameters                     [ok]
#>   individual_parameters          [ok]
#>   structural_model               [ok]
#>   error_model                    [ok]
#>   fit_options                    [ok] (optional)
#> 
#> Result: VALID
#> [1] TRUE

A template is a starting point. Adapt its initial estimates, bounds and blocks to your drug and data before you rely on it.

The `ferx_model` object

ferx_model(data, model) bundles a model file with a dataset. print() on a ferx_model shows both paths and the parsed structure:

m <- ferx_model(base$data, base$model)
m
#> ferx_model
#>   Model: /home/runner/work/_temp/Library/ferx/examples/models/two_cpt_oral_base.ferx
#>   Data:  /home/runner/work/_temp/Library/ferx/examples/data/two_cpt_oral_cov.csv
#>   ---
#>   Structural:  2-cpt oral  (TVCL, TVV1, TVQ, TVV2, TVKA)
#>   IIV:         ETA_CL, ETA_V1, ETA_Q, ETA_V2, ETA_KA
#>   IOV:         none
#>   Residual:    proportional

data comes first so that a dataset can be piped in. ferx_fit() accepts the object directly and takes the data from it:

fit <- base$data |> ferx_model(base$model) |> ferx_fit(verbose = FALSE)
fit$converged
#> [1] TRUE

Fitting itself is the topic of Initial estimates and a first fit.

Piping through an analysis

The argument order across the package is chosen so that a whole step reads as one pipeline: ferx_model() takes the data first, and the accessors that summarise a fit take the fit first. So a dataset can go to a fitted model and on to its warnings without naming an intermediate object:

base$data |>
  ferx_model(base$model) |>
  ferx_fit(verbose = FALSE) |>
  ferx_get_warnings(as_df = TRUE) |>
  select(severity, category)
#>   severity           category
#> 1  warning dw_autocorrelation

The same holds for the other fit accessors, and for simulation and prediction, where the model comes first and a fit is passed by name:

fit |> check_diagnostics() |> getElement("shrinkage") |> head(3)
#>    param type    shrinkage shrinkage_pct
#> 1 ETA_CL  eta -0.001744168    -0.1744168
#> 2 ETA_V1  eta  0.036475626     3.6475626
#> 3  ETA_Q  eta  0.029265476     2.9265476
base$model |> ferx_simulate(base$data, n_sim = 1, seed = 1) |> head(3)
#>   DRAW SIM ID TIME CMT    IPRED   DV_SIM OBSERVED
#> 1    1   1  1  0.5   1 1.262502 1.310942       NA
#> 2    1   1  1  1.0   1 1.926472 1.920568       NA
#> 3    1   1  1  2.0   1 2.302516 2.322542       NA

Which object pipes in is decided by the name of each function’s first argument, as the argument table at the end of this chapter records. The four ferx_model_* helpers whose first argument is path – ferx_model_show(), ferx_model_inspect(), ferx_model_validate() and ferx_model_edit() – read a model file, so a path pipes into them but a ferx_model object does not. The two whose first argument is x, ferx_model_get_section() and ferx_model_set_section(), take either.

try(m |> ferx_model_validate())
#> Error in file.exists(path) : invalid 'file' argument
m |> ferx_model_get_section("error_model", strip = TRUE)
#> # [error_model]
#> DV ~ proportional(PROP_ERR)

Editing model files from R

Reading and replacing a block

ferx_model_get_section() returns the lines of one block. It prints them and returns them invisibly; strip = TRUE trims leading whitespace:

param_lines <- ferx_model_get_section(m, "parameters", strip = TRUE)
#> # [parameters]
#> theta TVCL(5.0, 0.1, 100.0)
#> theta TVV1(50.0, 1.0, 500.0)
#> theta TVQ(10.0, 0.1, 100.0)
#> theta TVV2(100.0, 1.0, 500.0)
#> theta TVKA(1.2, 0.01, 10.0)
#> 
#> omega ETA_CL ~ 0.10
#> omega ETA_V1 ~ 0.10
#> omega ETA_Q  ~ 0.05
#> omega ETA_V2 ~ 0.05
#> omega ETA_KA ~ 0.15
#> 
#> sigma PROP_ERR ~ 0.02 (sd)
head(param_lines, 5)
#> [1] "theta TVCL(5.0, 0.1, 100.0)"   "theta TVV1(50.0, 1.0, 500.0)" 
#> [3] "theta TVQ(10.0, 0.1, 100.0)"   "theta TVV2(100.0, 1.0, 500.0)"
#> [5] "theta TVKA(1.2, 0.01, 10.0)"

ferx_model_set_section() replaces the body of a block with new lines. Given a ferx_model that points to a bundled example, it first copies the file to a temporary directory and edits the copy. The installed example stays unchanged:

m_focei <- ferx_model_set_section(m, "fit_options", c(
  "  method     = focei",
  "  maxiter    = 300",
  "  covariance = true"
))
#> ferx_model_set_section: copying read-only package model to /tmp/Rtmpe9t14F/two_cpt_oral_base_47f9788241c3.ferx before editing.
m_focei$model == m$model
#> [1] FALSE
ferx_model_get_section(m_focei, "fit_options")
#> # [fit_options]
#>   method     = focei
#>   maxiter    = 300
#>   covariance = true

Opening a model in an editor

ferx_model_edit() opens a model file in your editor. A bundled example is copied to dest first (default: the working directory), and overwrite controls whether an existing copy may be replaced. save_as copies the file after editing, which is a simple way to keep numbered versions of a model (run1.ferx, run2.ferx, …). The .editor argument replaces the editor function. The book passes a function that does nothing, so the example can run while the book is built:

no_editor <- function(...) invisible(NULL)
my_copy <- ferx_model_edit(base$model, dest = model_dir, overwrite = TRUE, .editor = no_editor)
#> Copied to /tmp/Rtmpe9t14F/model-files/two_cpt_oral_base.ferx; editing your copy.
run2 <- ferx_model_edit(my_copy, save_as = file.path(model_dir, "run2.ferx"), .editor = no_editor)
#> Saved copy to /tmp/Rtmpe9t14F/model-files/run2.ferx
basename(c(my_copy, run2))
#> [1] "two_cpt_oral_base.ferx" "run2.ferx"

In an interactive session, call ferx_model_edit(path) without .editor.

Options that matter

Function	Argument	Meaning
`ferx_model_show()`	`path`	Model file to print
`ferx_model_inspect()`	`path`	Model file, or a fit (reads its stored structure)
`ferx_model_validate()`	`path`	Model file (must exist and end in `.ferx`)
	`data`	Optional dataset for data-dependent checks
`ferx_model()`	`data`	Dataset path; if missing, the model’s `[data]` block is used at fit time
	`model`	Existing model file to wrap
	`template`	Scaffold from `"1cpt_oral"`, `"1cpt_iv"`, `"2cpt_oral"`, `"2cpt_iv"` or `"ode"`
	`path`	Where to write the scaffold
	`overwrite`	Replace an existing file at `path`
	`edit`	Open the new file in an editor (default `TRUE`)
	`print`	Print the skeleton only
`ferx_model_get_section()`	`x`, `section`, `strip`	Model object or path, block name without brackets, trim whitespace
`ferx_model_set_section()`	`x`, `section`, `lines`	Model object or path, block name, replacement body lines
`ferx_model_edit()`	`path`, `dest`, `overwrite`, `save_as`, `.editor`	File, copy destination, replace existing copy, save a copy after editing, editor function

Pitfalls

ferx_model_set_section() on a plain path edits that file in place. Copy-on-write only applies to a ferx_model object pointing into the installed package. Never pass ferx_example(...)$model itself as a path. Wrap it with ferx_model() or copy it first.
ferx_model_set_section() replaces an existing block. It cannot add a new one. To add a block such as [data], append the lines to a copy of the file (Preparing and checking the analysis dataset does this).
Argument order of ferx_model(). data is the first argument. An old-style call ferx_model("model.ferx") still works but warns about the changed order. Name model = to be explicit.

Summary

Read a model with ferx_model_show(). Check how it is parsed with ferx_model_inspect(), and validate it (optionally against the data) with ferx_model_validate().
Start new models from ferx_model(template = ...). Keep versions with ferx_model_edit(save_as = ...).
Edit blocks programmatically with ferx_model_get_section() and ferx_model_set_section(), on a ferx_model object or a copy.

Next: Initial estimates and a first fit fits this base model.

Reference

R help: ?ferx_model, ?ferx_model_show, ?ferx_model_inspect, ?ferx_model_validate, ?ferx_model_get_section, ?ferx_model_set_section, ?ferx_model_edit
ferx-core: model file reference, check report codes

Where you are

The data

Reading a model file

What an expression may contain

Checking a model before fitting

Starting a new model from a template

The ferx_model object

Piping through an analysis

Editing model files from R

Reading and replacing a block

Opening a model in an editor

Options that matter

Pitfalls

Summary

The `ferx_model` object