Example: Model search and automatic model development

Model development is a sequence of decisions — one compartment or two, which parameters carry between-subject variability, whether the residual error needs a time-varying term, which covariates earn their place. ferx can make those decisions by search rather than by hand, using the same tool chain as Pharmpy’s amd, driven from R.

There are seven tools and one pipeline that runs six of them in order:

Function Decides Ranks on
ferx_modelsearch() structural model — absorption, peripheral compartments, transits, lag time BIC (mixed)
ferx_iivsearch() which parameters carry an eta, and which etas are correlated BIC (IIV)
ferx_iovsearch() which parameters carry inter-occasion variability BIC (random)
ferx_ruvsearch() the residual-error model likelihood-ratio test
ferx_covsearch() which covariate effects to include (SCM, forward / forward-backward) likelihood-ratio test
ferx_allometry() allometric scaling on weight — a transform, not a search —
ferx_globalsearch() structure and covariates together, as one grid — exhaustively or by genetic algorithm penalized fitness
ferx_amd() the six stepwise tools, in order, each step starting from what the last one selected per step

Everything below runs against bundled examples, so it is reproducible without your own data.

The search space is MFL

A space is written in Pharmpy’s Model Feature Language, quoted verbatim — ferx never translates R arguments into MFL, which is what makes a space portable between the two tools.

ABSORPTION(FO); PERIPHERALS(0..1); LAGTIME([OFF,ON])
IIV?(@PK, exp); COVARIANCE?(IIV, *)

A statement without ? is structural: every candidate has it. A statement with ? is exploratory: the search decides. @PK, @IIV, @CONTINUOUS and @CATEGORICAL are resolved against your model and your dataset, so one space travels across models.

Check the space before paying for it

ferx_search_config() loads a .ferxsearch file through the engine’s own loader. An unknown section, an unparseable space, an empty space, an unrecognised [rank] type and a feature the engine cannot build are all errors here, before the first candidate is fitted.

library(ferx)

ex  <- ferx_example("warfarin_amd")
cfg <- ferx_search_config(ex$search)
cfg
<ferx search configuration>
  file:  warfarin_amd.ferxsearch
  base:  .../examples/models/warfarin_amd.ferx
  data:  .../examples/data/warfarin.csv

Search space (5 features):
  mfl = ABSORPTION(FO); PERIPHERALS(0..1); LAGTIME([OFF,ON])
IIV?(@PK, exp); COVARIANCE?(IIV, *)

  ABSORPTION(FO)                                 (structural)
  PERIPHERALS(0..1)                              (structural)
  LAGTIME([OFF,ON])                              (structural)
  IIV?(@PK,EXP)                                  (exploratory)
  COVARIANCE?(IIV,*)                             (exploratory)

Rank: bic   cutoff: (tool default)
Strictness (* = set by the file):
  * require_converged      TRUE
    require_covariance     FALSE
    max_condition_number   1000
    max_correlation        0.95
    reject_on_boundary     TRUE
    reject_init_stall      TRUE
Run: threads = auto, retries = 1, resume = FALSE
Tool sections: amd, iivsearch, modelsearch, ruvsearch

Two more functions make an interactive space debuggable. ferx_search_space() expands the @-symbols against a model and its data, so you can see what a one-line space actually means on your model:

ferx_search_space("COVARIATE?(@IIV, @CONTINUOUS, [pow, lin])",
                  model = ferx_example("two_cpt_oral_cov")$model,
                  data  = ferx_example("two_cpt_oral_cov")$data)
<ferx search space> (resolved against the base model)
  mfl = COVARIATE?(@IIV, @CONTINUOUS, [pow, lin])

                         feature   keyword optional
1    COVARIATE?(CL,WT,[pow,lin]) COVARIATE     TRUE
2  COVARIATE?(CL,CRCL,[pow,lin]) COVARIATE     TRUE
3    COVARIATE?(V1,WT,[pow,lin]) COVARIATE     TRUE
...
10 COVARIATE?(KA,CRCL,[pow,lin]) COVARIATE     TRUE

20 covariate effects (see attr(, "covariate_effects"))

And ferx_search_coverage() reports a feature the engine cannot express as a row rather than as an aborted run:

ferx_search_coverage("ABSORPTION([FO, ZO, SEQ-ZO-FO]); ELIMINATION(FO)")
                feature covered  reason
1 ABSORPTION(SEQ-ZO-FO)   FALSE  sequential zero-order-then-first-order absorption is
                                 not one input term on a standard disposition but a
                                 depot of its own ... write it by hand and search
                                 around it
2       ELIMINATION(FO)    TRUE  <NA>

The tools themselves refuse an uncovered feature outright — a search that quietly narrowed your space would answer a question you did not ask.

One tool at a time

Every tool takes either a .ferxsearch file or inline arguments, never a mix. The file is the reproducible artifact; the inline form is for a one-off.

ex <- ferx_example("warfarin")

res <- ferx_modelsearch(
  model        = ex$model,
  data         = ex$data,
  search_space = "ABSORPTION(FO); PERIPHERALS(0..1); LAGTIME([OFF,ON])",
  directory    = "modelsearch-run",
  progress     = interactive()
)
res
ferx structural model search (reduced_stepwise, ranked on bic_mixed)
  Data:  .../examples/data/warfarin.csv
  Wrote: modelsearch-run
  Base:  base - FO, 0 peripherals
  bic_mixed:   -261.8 (base) -> -261.8 (base)
  5 models over 2 layers

Ranked models (best first):
 rank   id              structure    ofv criterion d_criterion converged passed
    1 base      FO, 0 peripherals -280.4    -261.8       0.000      TRUE   TRUE
    2 run1 FO, 0 peripherals, lag -282.7    -259.6       2.258      TRUE   TRUE

3 models excluded by the strictness gate; summary() lists them.

Two things in that output are the point of the whole family:

  • The base model won. A search that cannot answer “none of the above” is a search that always finds something.

  • Three models were excluded, and they are still in the table. A candidate that failed to converge, hit a bound or came back ill-conditioned is not a result to be ranked. summary(res) prints each exclusion with the gate’s own words, and res$models$passed is a column you can filter on — never a silent drop:

    res$models[!res$models$passed, c("id", "structure", "failures")]
        id             structure
3 run2      FO, 1 peripheral
4 run3 FO, 1 peripheral, lag
5 run4 FO, 1 peripheral, lag
                                                                      failures
3 condition number 2.5654e5 exceeds 1.0000e3; parameter correlation |r| = 1.0000
  exceeds 0.9500 (TVKA ~ TVV2)
4 estimate pinned to a declared bound: TVQ
5 condition number 2.3613e5 exceeds 1.0000e3; parameter correlation |r| = 1.0000
  exceeds 0.9500 (TVKA ~ TVV2)

    Those three are the lowest-OFV models in the run. A winner-only view would have reported one of them.

The other tools read the same way. ferx_iivsearch() decides the variability structure in two stages — which parameters keep an eta, then which of the survivors are correlated:

ferx_iivsearch(model = ex$model, data = ex$data,
               search_space = "IIV?(@PK, exp); COVARIANCE?(IIV, *)",
               directory = "iivsearch-run")
ferx variability-structure search (top_down_exhaustive + block structure, ranked on bic_iiv)
  Input: [ETA_CL]+[ETA_KA]+[ETA_V]
  bic_iiv:   -273.5 (input) -> -273.5 (input)
  12 models over 2 steps

Step 1 (no_of_etas), parent input ([ETA_CL]+[ETA_KA]+[ETA_V]):
 rank    id                 structure criterion d_criterion  best
    1 input [ETA_CL]+[ETA_KA]+[ETA_V]   -273.50         0.0  TRUE
    2  run1         [ETA_CL]+[ETA_KA]    -14.16      -259.3 FALSE
    3  run2          [ETA_CL]+[ETA_V]     85.01      -358.5 FALSE
    ...

Step 2 (block_structure), parent input ([ETA_CL]+[ETA_KA]+[ETA_V]):
 rank    id                 structure criterion d_criterion  best
    1 input [ETA_CL]+[ETA_KA]+[ETA_V]    -273.5      0.0000  TRUE
    2 run11     [ETA_CL,ETA_KA,ETA_V]    -269.5     -3.9550 FALSE
   NA  run8   [ETA_CL,ETA_KA]+[ETA_V]    -274.1      0.6501 FALSE
   ...

Selected variability structure (input):
  IIV:         ETA_CL, ETA_KA, ETA_V
  Correlated:  (none - every eta is diagonal)

The etas are named, not numbered: a search over variability structure that reported OMEGA(2,1) would make you look the mapping up on every row. A block prints as ETA_V ~ ETA_CL and an IOV term as KAPPA_CL, with OMEGA(i,j) only where the model gave no name.

ferx_ruvsearch() tests residual-error features by likelihood-ratio test, one iteration at a time, and reports the p-value it decided on:

ferx_ruvsearch(model = ex$model, data = ex$data, directory = "ruvsearch-run")
ferx residual-error search (likelihood-ratio test at p = 0.001)
  OFV:   -280.364 (input) -> -280.364 (input, dOFV 0.000)
  1 iteration, 6 models fitted

Candidates by iteration:
 iteration       feature       family    ofv       dofv p_value significant passed
         1         power        power -280.7  0.3016000  0.5829       FALSE  FALSE
         1      combined     combined -280.4 -0.0005840  1.0000       FALSE   TRUE
         1 time_varying1 time_varying -280.9  0.5605000  0.4540       FALSE   TRUE
         1 time_varying2 time_varying -281.0  0.5963000  0.4400       FALSE   TRUE
         1 time_varying3 time_varying -280.4  0.0002039  0.9886       FALSE  FALSE

Selected residual-error model:
  (none - the input model's error model was not improved on)
    DV ~ proportional(PROP_ERR)

Notes:
  - IIV_on_RUV not tested: [fit_options] method has no eta-eps interaction, and
    `iiv_on_ruv` needs one (focei, imp, impmap or saem)

The remaining two need a model that declares its covariates, so they switch examples: two_cpt_oral_base is the covariate-free two-compartment oral model, and its [covariates] block is what makes @CONTINUOUS resolvable and WT available to scale on.

cov_ex <- ferx_example("two_cpt_oral_base")

# Stepwise covariate modelling (forward, then backward elimination)
ferx_covsearch(model = cov_ex$model, data = cov_ex$data,
               search_space = "COVARIATE?(@IIV, @CONTINUOUS, [pow, lin])",
               p_forward = 0.01, p_backward = 0.001,
               directory = "covsearch-run")

# Allometric scaling: a transform, so it can run without fitting at all.
# `parameters` is left at its default, which takes every clearance and volume
# the `pk` template line binds.
ferx_allometry(cov_ex$model, cov_ex$data, covariate = "WT",
               reference = 70, fit = FALSE)
ferx allometric scaling
  Data:  .../examples/data/two_cpt_oral_cov.csv
  Model: .../ferx-allometry-544f261810ae.ferx

Scaling on WT (reference 70):
  CL ~ WT power(center = 70, fix = 0.75)
  Q ~ WT power(center = 70, fix = 0.75)
  V1 ~ WT power(center = 70, fix = 1)
  V2 ~ WT power(center = 70, fix = 1)

Not fitted (`fit = FALSE`); the scaled model is in `$model`.

All axes at once: ferx_globalsearch()

Every stepwise search tool above keeps the other decision classes fixed. ferx_modelsearch() decides the structure with the covariate model fixed; ferx_covsearch() decides the covariates with the structure fixed. That is cheap, and it is blind to exactly one thing: an effect that only pays off once another decision has gone the other way. A clearance covariate that earns its place only after the distribution phase has a compartment of its own is invisible to a forward step taken on the one-compartment model.

ferx_globalsearch() lays both kinds of decision out as one grid and searches it whole. Every structural category in the space becomes an axis with its values as alleles; every optional COVARIATE? pair becomes an axis whose alleles are none and each of its forms. algorithm = "exhaustive" fits every point of the grid; algorithm = "ga" — the default — runs pyDarwin’s genetic algorithm over it and fits a fraction of them.

The bundled two_cpt_oral_global example is built to have something to find: a one-compartment, covariate-free model on data that carry WT and CRCL and were simulated from two compartments, so all four axes are live.

ex  <- ferx_example("two_cpt_oral_global")
res <- ferx_globalsearch(config = ex$search, directory = "globalsearch-run")
res
ferx global model search (exhaustive, ranked on penalized)
  Data:  .../examples/data/two_cpt_oral_cov.csv
  Wrote: globalsearch-run
  Grid:  16 points over 4 axes
    LAGTIME: OFF | ON
    PERIPHERALS: 0 | 1
    CL-WT: none | power
    CL-CRCL: none | power
  Input: OFV -172.2, penalized -102.2
  Fitness: -102.2 (input) -> -796 (run6)
  17 models fitted of 17 evaluated

Ranked models (best first):
 rank    id                                              genome    ofv criterion fitness converged passed
    1  run6  LAGTIME=OFF;PERIPHERALS=1;CL-WT=none;CL-CRCL=power -896.0   -796.00 -796.00      TRUE   TRUE
    2  run8 LAGTIME=OFF;PERIPHERALS=1;CL-WT=power;CL-CRCL=power -901.3   -791.30 -791.30      TRUE   TRUE
    3 run13    LAGTIME=ON;PERIPHERALS=1;CL-WT=none;CL-CRCL=none -887.5   -777.50 -777.50      TRUE   TRUE
    4 run14   LAGTIME=ON;PERIPHERALS=1;CL-WT=none;CL-CRCL=power -895.4   -775.40 -775.40      TRUE   TRUE
    5 run15   LAGTIME=ON;PERIPHERALS=1;CL-WT=power;CL-CRCL=none -893.4   -773.40 -773.40      TRUE   TRUE
    6 run16  LAGTIME=ON;PERIPHERALS=1;CL-WT=power;CL-CRCL=power -896.4   -766.40 -766.40      TRUE   TRUE
    7  run1   LAGTIME=OFF;PERIPHERALS=0;CL-WT=none;CL-CRCL=none -172.2   -102.20 -102.20      TRUE   TRUE
    8 input                                                <NA> -172.2   -102.20 -102.20      TRUE   TRUE
    9  run2  LAGTIME=OFF;PERIPHERALS=0;CL-WT=none;CL-CRCL=power -179.4    -99.36  -99.36      TRUE   TRUE
   10  run3  LAGTIME=OFF;PERIPHERALS=0;CL-WT=power;CL-CRCL=none -179.3    -99.27  -99.27      TRUE   TRUE
   11  run4 LAGTIME=OFF;PERIPHERALS=0;CL-WT=power;CL-CRCL=power -185.5    -95.49  -95.49      TRUE   TRUE

Read the first two rows before anything else. run8 has the better OFV (-901.3 against -896.0) and ranks second. It carries one more covariate relation than run6, and the penalized criterion charges it for the extra parameter. A table that showed only the OFV would name the wrong winner — which is the reason the criterion and the OFV are separate columns.

The criterion is not the whole of the ranking

[rank] type defaults to "penalized" here, where every other tool defaults to a BIC: pyDarwin’s fitness, the objective function plus a charge per estimated parameter and per unhealthy fit. Whatever criterion you name, the global search then charges three further things the criterion cannot see, and ranks on the sum — the fitness column, always finite:

charge paid by
charge_non_influential a gene that changed nothing in the rendered model — a covariate on a parameter the structural choice removed. A tie-break towards the simpler genotype among models that are otherwise identical.
charge_gate a fit the [strictness] gate refused. It can never be selected, but it still steers the genetic algorithm away from the region that produced it.
charge_crash a candidate that produced no fit at all — it does not compile, or the fit errored. Large, and finite, so the algorithm can still rank it.

They are columns of $models, beside the criterion and the fitness they add up to. print() lists the models that paid one:

Models the search charged beyond the criterion:
    id                                             genome criterion charge_non_influential charge_gate charge_crash fitness
  run5  LAGTIME=OFF;PERIPHERALS=1;CL-WT=none;CL-CRCL=none   -411.00                      0         100            0  -311.0
  run7 LAGTIME=OFF;PERIPHERALS=1;CL-WT=power;CL-CRCL=none   -404.20                      0         100            0  -304.2
  run9   LAGTIME=ON;PERIPHERALS=0;CL-WT=none;CL-CRCL=none     16.52                      0         100            0   116.5
 run10  LAGTIME=ON;PERIPHERALS=0;CL-WT=none;CL-CRCL=power     19.34                      0         100            0   119.3
 run11  LAGTIME=ON;PERIPHERALS=0;CL-WT=power;CL-CRCL=none     19.46                      0         100            0   119.5
 run12  LAGTIME=ON;PERIPHERALS=0;CL-WT=power;CL-CRCL=power     23.34                      0         100            0   123.3

This is the same “never a silent drop” rule the stepwise tools follow, one level further: a gated genome is not merely present in the table, it says what the gate cost it. Without that column, run5 at fitness -311.0 would look like a model that simply fitted worse than run6 — rather than one the gate refused, which happens to sit 100 units away because of it.

$penalties is the effective schedule the run charged, so the fitness column can be read against the numbers that built it.

The genetic algorithm

Sixteen points is small enough to enumerate, which makes the exhaustive result a reference for the GA. With the settings and seed below, the GA finds the same winner. On a grid of thousands — six binary optional covariate pairs alone already give 64 covariate combinations — enumerating is not an option and finding a strong candidate at all is the point.

res_ga <- ferx_globalsearch(
  model        = ex$model,
  data         = ex$data,
  search_space = c("PERIPHERALS(0..1); LAGTIME([OFF,ON])",
                   "COVARIATE?(CL, WT, pow); COVARIATE?(CL, CRCL, pow)"),
  algorithm    = "ga",
  ga           = list(population_size = 8, generations = 3, seed = 20250914),
  directory    = "globalsearch-ga"
)
res_ga$generations
  index best best_fitness mean_fitness polished
1     0 run1    -779.4346    -323.4043        0
2     1 run9    -796.0057    -438.3091        0
3     2 run9    -796.0057    -441.3968        0
4     3 run9    -796.0057    -497.5284        1

ga takes a named list validated against the engine’s own key list, so a mistyped knob is refused by name before anything is fitted. The seed is what makes resume = TRUE reuse a journal: the same seed on the same grid proposes the same genomes.

algorithm is matched exactly rather than by prefix — algorithm = "e" is an error naming the value, not a silent choice to enumerate a grid nobody has sized. For the same reason an exhaustive grid above [globalsearch] max_models (500 by default) is an error naming the size, never a truncated search.

NoteNot the same as ferx_fit(settings = list(global_search = TRUE))

That setting is a global optimizer phase inside the estimation of one model — a wider search for the parameter values of a fixed model. ferx_globalsearch() is a global search over models. The name here is the engine’s, the .ferxsearch file’s ([globalsearch]), the ferx globalsearch command’s and Pharmpy/pyDarwin’s, so a search stays portable between them.

The whole pipeline: ferx_amd()

ferx_amd() runs the tools in order, each step starting from the model the previous step selected. One space describes the whole pipeline; the engine partitions it by statement kind and hands each step only the statements its tool can read.

Plan first

ferx_amd_plan() answers “what would run, and what would not” without fitting anything:

ex <- ferx_example("warfarin_amd")
ferx_amd_plan(config = ex$search)
  index       step        tool rerun      directory
1     1 structural modelsearch FALSE 01-modelsearch
2     2  iivsearch   iivsearch FALSE   02-iivsearch
3     3   residual   ruvsearch FALSE   03-ruvsearch
4     4  iovsearch   iovsearch FALSE   04-iovsearch
5     5  allometry   allometry FALSE   05-allometry
6     6 covariates   covsearch FALSE   06-covsearch
                                                                        skipped
1                                                                          <NA>
2                                                                          <NA>
3                                                                          <NA>
4 the starting model declares no `iov_column` in [fit_options], so the dataset's
                                                     occasions were not read
5                                   the search space has no ALLOMETRY statement
6                                   the search space has no COVARIATE statement

A step is skipped for a reason you can read, not because it silently found nothing to do.

Then run it

res <- ferx_amd(config = ex$search, directory = "amd-run",
                progress = interactive())
res
ferx AMD pipeline (default strategy, retries all_final)
  Model: .../examples/models/warfarin_amd.ferx
  Data:  .../examples/data/warfarin.csv
  Wrote: amd-run
  OFV:   -286.004 (start) -> -286.004 (final, dOFV 0.000)
  3 of 6 steps ran, 28 models fitted, 8.8 s across the steps

Pipeline:
 index       step  status criterion value_before value_after    d_value
     1 structural     ran bic_mixed       -267.5      -267.5 -4.865e-10
     2  iivsearch     ran   bic_iiv       -279.1      -279.1  1.021e-09
     3   residual     ran       ofv       -286.0      -286.0 -4.380e-10
     4  iovsearch skipped      <NA>           NA          NA         NA
     5  allometry skipped      <NA>           NA          NA         NA
     6 covariates skipped      <NA>           NA          NA         NA
      d_ofv candidates converged passed seconds
 -4.865e-10          5      TRUE   TRUE  0.6229
  1.021e-09         12      TRUE   TRUE  0.6219
 -4.380e-10          7        NA     NA  7.5060
         NA          0        NA     NA  0.0000
         NA          0        NA     NA  0.0000
         NA          0        NA     NA  0.0000

Steps that did not run:
  4 iovsearch    skipped: the starting model declares no `iov_column` in [fit_options], so the dataset's occasions were not read
  5 allometry    skipped: the search space has no ALLOMETRY statement
  6 covariates   skipped: the search space has no COVARIATE statement

Selected:
  structural   FO, 0 peripherals
  iivsearch    [CL]+[KA]+[V]
  residual     no residual-error feature added

3 candidates excluded by the strictness gate; summary() lists them.

The strategy decides the order: default is structural → IIV → residual → IOV → allometry → covariates; reevaluation appends a second IIV and residual pass once the model has its covariates; SIR, SRI and RSI run three steps in the order their initials spell.

What the result carries

# One row per planned step, skipped ones included
res$steps[, c("index", "step", "status", "criterion", "d_value",
              "converged", "passed", "selected")]

# Every candidate of every step, with the tool that produced it
head(res$candidates[, c("step", "tool", "id", "description", "criterion",
                        "value", "rank", "passed", "selected")])

# Candidates the gate excluded, in the gate's own words
res$candidates[!res$candidates$passed, c("step", "tool", "id", "failures")]

# Each step's own result object, keyed by its directory
names(res$tools)
#> [1] "01-structural" "02-iivsearch"  "03-residual"

# The model the pipeline ended on: a fitted model, seeded with its own estimates
res$fit
res$final_model_path
res$d_ofv                  # final - start; negative is an improvement

# The engine's own report, as `ferx amd` prints it
cat(res$summary_text)

One trap worth naming: criterion says what value is on, and it is not always the step’s own scale. A ruvsearch candidate is fitted to its parent’s CWRES, so its number must not be compared with a data OFV. The column is there so you do not have to remember that.

Reading a run back

Every tool writes its run to directory, and the tables are the same ones the R object carries — so an interrupted run, or one driven from the ferx command-line tool, reads back the same way:

ferx_search_results("modelsearch-run", type = "models")     # models.csv
ferx_search_results("covsearch-run",   type = "steps")      # steps.csv
ferx_search_results("modelsearch-run")                      # candidates.csv

Columns come from the engine’s own column list, typed on the way in: logicals are logical, an empty cell is NA rather than NaN. A partial run (*.partial.csv, written when a search is cancelled) reads back too.

resume = TRUE picks a run up where it stopped, using the journal in directory; threads and retries mirror ferx_bootstrap(), the closest existing tool.

Notes

  • Some .ferxsearch sections have no R tool. The engine accepts a section for every tool it knows, which is more than this package binds — [structsearch] today. Such a file loads, and ferx warns that no R tool consumes the section, because the tool it is handed would otherwise run a different search without saying so. [structsearch] is accepted vocabulary with no engine module behind it yet, so it is reported without a remediation rather than with one that names a command that does not exist. What the warning reports is the complement of the sections this package has a tool for, so a section a later ferx-core adds is named from the day a file can carry it — and stops being named on the day a binding lands, as [globalsearch] did when ferx_globalsearch() shipped.
  • Cancelling a search from R interrupts it at the next candidate boundary; the run still writes its partial tables.
  • The strictness gate is shared. [strictness] in the configuration applies to every tool, and every table with one row per fitted model carries converged and passed columns — $models for modelsearch, iivsearch and iovsearch, $steps for covsearch and ruvsearch, $comparison for allometry. (The per-stage ranking $steps of iivsearch and iovsearch, and allometry’s $scalings, are not model tables and carry no gate columns; the models they rank are in $models.) Inside an AMD step the init-stall gate is turned off by design: each step starts from a model already at its own optimum, where “did not move from the initial estimates” is the expected outcome.

See also