Model development and selection

Where you are

The base model is fitted (Initial estimates and a first fit) and evaluated (Diagnosing the model, Simulation-based evaluation: VPC). The diagnostics hinted at covariate relationships (ETA_CL correlated with CRCL and WT) and flagged autocorrelated residuals. Model development means fitting alternatives and choosing between them. You can do that by hand: write a candidate, fit it and compare. Or let ferx’s search tools generate, fit, gate and rank candidates from a search space. This chapter shows both.

The data

base <- ferx_example("two_cpt_oral_base")
cov <- ferx_example("two_cpt_oral_cov")

Comparing fitted models

two_cpt_oral_cov is a prespecified covariate model for the same data. Body weight scales CL and V1 through a shared exponent THETA_WT, and creatinine clearance scales CL:

invisible(ferx_model_get_section(cov$model, "individual_parameters"))
#> # [individual_parameters]
#>   CL = TVCL * (WT / 70)^THETA_WT * (CRCL / 100)^THETA_CRCL * exp(ETA_CL)
#>   V1 = TVV1 * (WT / 70)^THETA_WT * exp(ETA_V1)
#>   Q  = TVQ  * exp(ETA_Q)
#>   V2 = TVV2 * exp(ETA_V2)
#>   KA = TVKA * exp(ETA_KA)

fit_base <- ferx_fit(base$model, base$data, verbose = FALSE)
fit_cov <- ferx_fit(cov$model, cov$data, verbose = FALSE)

Likelihood ratio test and information criteria

Both are FOCEI fits to the same data, and the base model is nested in the covariate model (fix both exponents at 0). The difference in OFV is therefore a likelihood ratio statistic, with degrees of freedom equal to the number of added parameters:

d_ofv <- fit_base$ofv - fit_cov$ofv
df <- fit_cov$n_parameters - fit_base$n_parameters
c(d_ofv = d_ofv, df = df, p_value = pchisq(d_ofv, df, lower.tail = FALSE))
#>         d_ofv            df       p_value 
#> 13.9235928014  2.0000000000  0.0009473931

For non-nested models, use information criteria. fit$aic and fit$bic are the classical AIC and BIC. ferx_bic() computes four BIC variants that differ in how parameters are penalised:

`type`	Penalty added to the OFV
`"mixed"` (default)	random-effect parameters × log(subjects) + fixed-effect parameters × log(observations)
`"fixed"`	all parameters × log(observations); equals `fit$bic`
`"iiv"`	variance parameters × log(subjects), for comparing variability structures
`"random"`	all parameters × log(subjects)

Table 9.1: Base and covariate model compared on OFV and information criteria (lower is better).

criteria <- function(f) {
  c(ofv = f$ofv, parameters = f$n_parameters, aic = f$aic, bic = f$bic,
    bic_mixed = ferx_bic(f, "mixed"), bic_iiv = ferx_bic(f, type = "iiv"),
    bic_random = ferx_bic(f, "random"))
}
round(rbind(base = criteria(fit_base), covariate = criteria(fit_cov)), 2)
#>                ofv parameters      aic      bic bic_mixed  bic_iiv bic_random
#> base      -1185.40         11 -1163.40 -1122.66  -1145.69 -1168.40   -1147.99
#> covariate -1199.33         13 -1173.33 -1125.18  -1152.81 -1182.32   -1155.11

Here the covariate model is better on every criterion.

Is the fit eligible at all?

A converged fit can still be unfit for comparison: it may never have left its initial estimates, have a parameter at a bound, or have an ill-conditioned covariance matrix. check_strictness() applies such gates and returns the reasons for failure:

check_strictness(fit_cov, require_covariance = TRUE)
#> $passed
#> [1] TRUE
#> 
#> $failures
#> character(0)
#> 
#> $skipped
#> character(0)

The fit from NCA-based starting values in Initial estimates and a first fit converged to a much worse OFV. The gates reject it:

fit_nca <- ferx_fit(base$model, base$data, inits_from_nca = TRUE, verbose = FALSE)
check_strictness(fit_nca)
#> $passed
#> [1] FALSE
#> 
#> $failures
#> [1] "condition number 5.7691e+09 exceeds 1.0000e+03"   
#> [2] "parameter correlation |r| = 1.0000 exceeds 0.9500"
#> 
#> $skipped
#> character(0)

The gates and their defaults: require_converged = TRUE, require_covariance = FALSE, max_condition_number = 1000, max_correlation = 0.95, reject_on_boundary = TRUE and reject_init_stall = TRUE. Set a threshold to NULL, or a flag to FALSE, to disable that gate. reject_init_stall assumes a cold start from the model file’s initial values, so turn it off when candidates start from a previous model’s estimates.

Search configurations

Maturity: alpha

The ferx-core pages for the search tools carry the maturity label alpha. See feature maturity.

The search tools read a .ferxsearch configuration file (TOML) that names the base model and data and states the search. Four bundled examples carry one as $search:

Example	Search file for
`two_cpt_oral_base`	covariate search (`[covsearch]`), plus an `[allometry]` section (Covariate modeling)
`two_cpt_oral_cov`	a configuration with a BIC ranking, used here to show loading
`warfarin`	structural search (`[modelsearch]`)
`one_cpt_transit`	residual error search (`[ruvsearch]`)

ferx_search_config() loads a file with the engine’s own strict loader, so an invalid space or a misspelt section fails before anything is fitted:

cfg <- ferx_search_config(base$search)
cfg
#> <ferx search configuration>
#>   file:  two_cpt_oral_base.ferxsearch
#>   base:  /home/runner/work/_temp/Library/ferx/examples/search/../models/two_cpt_oral_base.ferx
#>   data:  /home/runner/work/_temp/Library/ferx/examples/search/../data/two_cpt_oral_cov.csv
#> 
#> Search space (1 feature):
#>   mfl = COVARIATE?(@IIV, @CONTINUOUS, [pow, lin])
#>   COVARIATE?(@IIV,@CONTINUOUS,[pow,lin])         (exploratory)
#> 
#> Rank: (tool default)   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 = 2, resume = FALSE
#> Tool sections: allometry, covsearch

The file’s sections, as shown in the printout: base and data (paths relative to the file); [space] with the search space as mfl; [rank] with the ranking criterion; [strictness] with the gates from above; [run] with threads, retries and resume; and a section per tool. The object returned holds the same information as a list:

cfg_cov <- ferx_search_config(ferx_example("two_cpt_oral_cov")$search)
cfg_cov$rank$type
#> [1] "bic"
cfg$tools
#> [1] "allometry" "covsearch"

The search space

The space is written in the Model Feature Language (MFL). COVARIATE?(@IIV, @CONTINUOUS, [pow, lin]) means: every parameter with a random effect, crossed with every continuous covariate, as a power or linear effect, each optional (?). The ? is accepted on four keywords only – COVARIATE, IIV, IOV and COVARIANCE – and marks that statement as something the search may decline; written without it, every candidate carries the feature. The other keywords, ABSORPTION and PERIPHERALS among them, enumerate the alternatives to try instead. The @ symbols resolve against the model’s [covariates] block and its random effects, so one space travels between models: @CATEGORICAL selects the categorical covariates as @CONTINUOUS selects the continuous ones. ferx_search_space() parses MFL text and, given a model and data, expands the symbols into the concrete candidates:

ferx_search_space("COVARIATE?(@IIV, @CONTINUOUS, [pow, lin])")
#> <ferx search space> (as written)
#>   mfl = COVARIATE?(@IIV, @CONTINUOUS, [pow, lin])
#> 
#>                                  feature   keyword optional
#> 1 COVARIATE?(@IIV,@CONTINUOUS,[pow,lin]) COVARIATE     TRUE
space <- ferx_search_space(cfg$mfl, model = base$model, data = base$data)
nrow(attr(space, "covariate_effects"))
#> [1] 20
head(attr(space, "covariate_effects"))
#>   parameter covariate effect op optional
#> 1        CL        WT    pow  *     TRUE
#> 2        CL        WT    lin  *     TRUE
#> 3        CL      CRCL    pow  *     TRUE
#> 4        CL      CRCL    lin  *     TRUE
#> 5        V1        WT    pow  *     TRUE
#> 6        V1        WT    lin  *     TRUE

ferx_search_coverage() reports, without fitting, whether the engine can build candidates for each feature of a space:

ferx_search_coverage(space)
#>                                  feature covered reason
#> 1 COVARIATE?(@IIV,@CONTINUOUS,[pow,lin])    TRUE   <NA>
ferx_search_coverage("ABSORPTION(SEQ-ZO-FO)")
#>                 feature covered
#> 1 ABSORPTION(SEQ-ZO-FO)   FALSE
#>                                                                                                                                                                                                                                                                                                                       reason
#> 1 sequential zero-order-then-first-order absorption is not one input term on a standard disposition but a depot of its own, filled at a constant rate and emptied by `ka`; there is no template for it and the search does not generate one. Write it by hand — `examples/sequential_absorption.ferx` — and search around it

Sequential zero- then first-order absorption is not generated by the search. Absorption and bioavailability shows the bundled model for it, which you can use as a base model instead.

An uncovered feature is not quietly dropped. Asking for one in a configuration file is an error at load, before anything is fitted, because a search that silently narrowed its space would answer a different question than the one you posed:

uncovered <- file.path(book_tempdir("search-space"), "uncovered.ferxsearch")
writeLines(sub('mfl = .*', 'mfl = "ABSORPTION([FO,SEQ-ZO-FO])"', readLines(base$search)), uncovered)
tryCatch(ferx_search_config(uncovered), error = function(e) sub(".*mfl: ", "", conditionMessage(e)))
#> [1] "the search space asks for 1 feature ferx cannot build a candidate for:\n  - ABSORPTION(SEQ-ZO-FO): sequential zero-order-then-first-order absorption is not one input term on a standard disposition but a depot of its own, filled at a constant rate and emptied by `ka`; there is no template for it and the search does not generate one. Write it by hand — `examples/sequential_absorption.ferx` — and search around it\nSee the coverage table at https://ferx-nlme.github.io/ferx-core/tools/search.html#coverage"

That is what ferx_search_coverage() is for: check a space before you file it.

Covariate search

ferx_covsearch() runs a stepwise covariate search. In each forward step it fits every remaining candidate effect and adds the one with the largest significant OFV drop at p_forward. The backward phase then removes effects whose removal is not significant at the stricter p_backward. The bundled configuration uses p_forward = 0.01 and p_backward = 0.001:

cs_dir <- book_tempdir("covsearch")
cs <- ferx_covsearch(config = base$search, directory = cs_dir, progress = FALSE)
cs
#> ferx covariate search (scm-forward-then-backward)
#>   Data:  /home/runner/work/_temp/Library/ferx/examples/search/../data/two_cpt_oral_cov.csv
#>   Wrote: /tmp/RtmpJuwo6H/covsearch
#>   OFV:   -1185 (base) -> -1185 (step 3)
#> 
#> Steps (the candidate selected at each step):
#>  step    phase parameter covariate   form  dofv  p_value converged passed
#>     1  forward        CL      CRCL linear 6.687 0.009712      TRUE   TRUE
#>     3 backward        CL      CRCL linear 6.675 0.009776      TRUE   TRUE
#> 
#> 4 candidates excluded by the strictness gate; summary() lists them.
#> 
#> Final covariate relations:
#>   (none - the base model was not improved on)

The forward phase added CL ~ CRCL (linear) at p = 0.0097, below p_forward. That p-value is above p_backward, so the backward phase removed the effect again. The search returns the base model. Candidates excluded by the strictness gate are kept in the step table with their reason:

cs$steps[!cs$steps$passed, c("step", "parameter", "covariate", "form", "failures")]
#>    step parameter covariate   form
#> 8     1        V1      CRCL linear
#> 20    1        KA      CRCL linear
#> 22    2        CL        WT linear
#> 38    2        KA      CRCL linear
#>                                                                                                                         failures
#> 8           condition number 8.1386e8 exceeds 1.0000e3; parameter correlation |r| = 1.0000 exceeds 0.9500 (TVKA ~ THETA_V1_CRCL)
#> 20     condition number 1.3933e10 exceeds 1.0000e3; parameter correlation |r| = 1.0000 exceeds 0.9500 (THETA_KA_CRCL ~ PROP_ERR)
#> 22          condition number 7.7525e8 exceeds 1.0000e3; parameter correlation |r| = 1.0000 exceeds 0.9500 (TVV1 ~ THETA_CL_CRCL)
#> 38 condition number 2.9288e9 exceeds 1.0000e3; parameter correlation |r| = 1.0000 exceeds 0.9500 (THETA_CL_CRCL ~ THETA_KA_CRCL)

print() shows the selected step of each phase. summary() on the result adds every candidate of every step; the same table is cs$steps. The run directory holds the per-step journals, steps.csv and final.ferx. ferx_search_results() reads the step table back from disk:

list.files(cs_dir)
#> [1] "backward-3" "base"       "final.ferx" "forward-1"  "forward-2" 
#> [6] "steps.csv"
steps <- ferx_search_results(cs_dir, type = "steps")
attr(steps, "tool")
#> [1] "covsearch"

The search can also be stated inline. Here it is restricted to CL, with a less strict backward level:

cs_cl <- ferx_covsearch(model = base$model, data = base$data,
                        search_space = "COVARIATE?(CL, [WT, CRCL], [pow, lin])",
                        p_forward = 0.01, p_backward = 0.01,
                        directory = book_tempdir("covsearch-cl"), progress = FALSE)
cs_cl$included
#>   parameter covariate   form         origin
#> 1        CL      CRCL linear forward step 1
c(base_ofv = cs_cl$base_ofv, final_ofv = cs_cl$final_ofv)
#>  base_ofv final_ofv 
#> -1185.423 -1192.110

The selected covariate model depends on the search space and the significance levels. A search result is evidence, not a decision. The rest of The analysis workflow continues with the prespecified two_cpt_oral_cov model, which the comparison above supports.

Structural search

ferx_modelsearch() searches structural features: absorption, elimination, peripheral and transit compartments, and lag time. It ranks the candidates on a criterion, the mixed BIC by default. The warfarin configuration searches a peripheral compartment and a lag time:

ms_dir <- book_tempdir("modelsearch")
warf <- ferx_example("warfarin")
ms <- ferx_modelsearch(config = warf$search, directory = ms_dir, progress = FALSE)
ms
#> ferx structural model search (reduced_stepwise, ranked on bic_mixed)
#>   Data:  /home/runner/work/_temp/Library/ferx/examples/search/../data/warfarin.csv
#>   Wrote: /tmp/RtmpJuwo6H/modelsearch
#>   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 run2       FO, 1 peripheral -282.5    -254.6       7.254      TRUE   TRUE
#>     4 run4  FO, 1 peripheral, lag -283.6    -251.1      10.730      TRUE   TRUE
#>     5 run3  FO, 1 peripheral, lag -282.7    -250.2      11.660      TRUE   TRUE

The base structure ranks best on the mixed BIC. The lag-time model has a lower OFV, but not by enough to pay for its extra parameters. Every candidate’s model text is kept, so the runner-up can be inspected or refitted:

runner_up <- file.path(ms_dir, "runner_up.ferx")
writeLines(ms$model_text[["run1"]], runner_up)
invisible(ferx_model_get_section(runner_up, "structural_model"))
#> # [structural_model]
#>   pk one_cpt_oral(cl=CL, v=V, ka=KA, lagtime=ALAG)

With resume = TRUE and the same directory, a search reuses the fits it already journalled:

ms_again <- ferx_modelsearch(config = warf$search, directory = ms_dir, resume = TRUE, progress = FALSE)
ms_again$models[, c("id", "structure", "reused")]
#>     id              structure reused
#> 1 base      FO, 0 peripherals   TRUE
#> 2 run1 FO, 0 peripherals, lag   TRUE
#> 3 run2       FO, 1 peripheral   TRUE
#> 4 run3  FO, 1 peripheral, lag   TRUE
#> 5 run4  FO, 1 peripheral, lag   TRUE
nrow(ferx_search_results(ms_dir, type = "models"))
#> [1] 5

Residual error search

ferx_ruvsearch() tests four families of residual error extensions against the current error model, using the likelihood ratio test at p_value:

IIV_on_RUV: a random effect on the residual error magnitude. It needs an estimation method with interaction.
power: a proportional error raised to an estimated power.
combined: an additive component added to the proportional one.
time_varying: a residual error magnitude that differs below a time-after-dose cutoff, set at quantiles given by groups.

Up to max_iter accepted features are stacked. The bundled configuration searches the analytic transit model one_cpt_transit:

transit <- ferx_example("one_cpt_transit")
rs <- ferx_ruvsearch(config = transit$search, directory = book_tempdir("ruvsearch"), progress = FALSE)
rs
#> ferx residual-error search (likelihood-ratio test at p = 0.05)
#>   Data:  /home/runner/work/_temp/Library/ferx/examples/search/../data/transit_oral.csv
#>   Wrote: /tmp/RtmpJuwo6H/ruvsearch
#>   OFV:   -1063.880 (input) -> -1070.083 (time_varying1-1, dOFV 6.203)
#>   2 iterations, 10 models fitted
#> 
#> Candidates by iteration:
#>  iteration       feature       family   ofv      dofv p_value significant
#>          1    IIV_on_RUV   IIV_on_RUV -1064 -0.001255 1.00000       FALSE
#>          1         power        power -1066  1.896000 0.16850       FALSE
#>          1      combined     combined -1059 -5.060000 1.00000       FALSE
#>          1 time_varying1 time_varying -1070  6.203000 0.01275        TRUE
#>          1 time_varying2 time_varying -1065  0.761800 0.38280       FALSE
#>          1 time_varying3 time_varying -1064  0.036270 0.84900       FALSE
#>          2    IIV_on_RUV   IIV_on_RUV -1070 -0.001305 1.00000       FALSE
#>          2         power        power -1071  0.986500 0.32060       FALSE
#>          2      combined     combined -1061 -9.261000 1.00000       FALSE
#>  selected converged passed
#>     FALSE      TRUE   TRUE
#>     FALSE      TRUE   TRUE
#>     FALSE      TRUE   TRUE
#>      TRUE      TRUE   TRUE
#>     FALSE      TRUE   TRUE
#>     FALSE      TRUE   TRUE
#>     FALSE      TRUE   TRUE
#>     FALSE      TRUE   TRUE
#>     FALSE      TRUE   TRUE
#> 
#> Selected residual-error model:
#>   time_varying1  (time_varying)
#>     DV ~ proportional(PROP_ERR * (if (TAD < 1.375) RUV_TV else 1.0))
rs$final_features
#> [1] "time_varying1"

cwres_prescreen changes what those numbers are. With it on, candidates are screened against the parent’s CWRES rather than fitted to the data, which is cheaper but puts the criterion on a different scale. Two runs that differ only in that flag, everything else held equal:

ruv <- function(prescreen) {
  ferx_ruvsearch(model = transit$model, data = transit$data, p_value = 0.05, max_iter = 1,
                 retries = 0, cwres_prescreen = prescreen,
                 directory = book_tempdir(paste0("ruvsearch-", prescreen)), progress = FALSE)
}
rs_data <- ruv(FALSE)
rs_screen <- ruv(TRUE)
screened <- startsWith(rs_screen$candidates$id, "cwres")
rbind(fitted_to_data = range(rs_data$candidates$criterion[rs_data$candidates$id != "input"]),
      screened_on_cwres = range(rs_screen$candidates$criterion[screened]))
#>                         [,1]       [,2]
#> fitted_to_data    -1070.0829 -1058.8200
#> screened_on_cwres   236.6697   240.6815
c(selected_data = length(rs_data$final_features), selected_screened = length(rs_screen$final_features))
#>     selected_data selected_screened 
#>                 1                 1

The two sets of numbers are not comparable, and the screening rows are the ones whose id starts with cwres. Note that the cheaper pass also selected differently here, so treat prescreening as a different search rather than a faster one. The engine default is cwres_prescreen = FALSE, which is what the configuration run above used.

The search accepted a time-varying residual error in the first iteration, and nothing further in the second. The base model of this workflow, which carried a residual autocorrelation warning in Diagnosing the model, can be searched the same way with the inline arguments (model, data, p_value, max_iter, …).

Variability search

ferx_iivsearch() decides the random-effect structure: which parameters carry an eta, and which etas share a covariance block. The space is written with IIV and COVARIANCE statements, and ? marks what the search may decide. The bundled warfarin_block_omega configuration makes every PK parameter’s eta exploratory and offers a block over whichever survive:

bo <- ferx_example("warfarin_block_omega")
iiv <- ferx_iivsearch(config = bo$search, directory = book_tempdir("iivsearch"), progress = FALSE)
#> Warning in .ferx_compute_cor_matrix(result$cov_matrix): One or more diagonal
#> elements are non-positive; correlation matrix may not be meaningful.
iiv
#> ferx variability-structure search (top_down_exhaustive + block structure, ranked on bic_iiv)
#>   Data:  /home/runner/work/_temp/Library/ferx/examples/search/../data/warfarin_block_omega.csv
#>   Wrote: /tmp/RtmpJuwo6H/iivsearch
#>   Input: [ETA_CL,ETA_V]+[ETA_KA]
#>   bic_iiv:   -271.3 (input) -> -272.5 (run10)
#>   12 models over 3 steps
#> 
#> Step 1 (no_of_etas), parent input ([ETA_CL,ETA_V]+[ETA_KA]):
#>  rank    id               structure criterion d_criterion  best
#>     1 input [ETA_CL,ETA_V]+[ETA_KA]   -271.30         0.0  TRUE
#>     2  run1       [ETA_CL]+[ETA_KA]    -14.16      -257.1 FALSE
#>     3  run5                [ETA_KA]     96.90      -368.2 FALSE
#>     4  run4                [ETA_CL]    108.60      -379.9 FALSE
#>     5  run6                 [ETA_V]    178.10      -449.4 FALSE
#>     6  run7                            181.70      -453.0 FALSE
#>    NA  run2          [ETA_CL,ETA_V]     86.80      -358.1 FALSE
#>    NA  run3        [ETA_KA]+[ETA_V]     27.38      -298.7 FALSE
#> 
#> Step 2 (block_structure), parent input ([ETA_CL,ETA_V]+[ETA_KA]):
#>  rank    id                 structure criterion d_criterion  best
#>     1 run10   [ETA_KA,ETA_V]+[ETA_CL]    -272.5      1.1930  TRUE
#>     2  run9   [ETA_CL,ETA_KA]+[ETA_V]    -272.2      0.9037 FALSE
#>     3 input   [ETA_CL,ETA_V]+[ETA_KA]    -271.3      0.0000 FALSE
#>     4 run11     [ETA_CL,ETA_KA,ETA_V]    -269.5     -1.7740 FALSE
#>    NA  run8 [ETA_CL]+[ETA_KA]+[ETA_V]    -273.5      2.1810 FALSE
#> 
#> Step 3 (compare_to_input), parent input ([ETA_CL,ETA_V]+[ETA_KA]):
#>  rank    id               structure criterion d_criterion  best
#>     1 run10 [ETA_KA,ETA_V]+[ETA_CL]    -272.5       1.193  TRUE
#>     2 input [ETA_CL,ETA_V]+[ETA_KA]    -271.3       0.000 FALSE
#> 
#> 3 models excluded by the strictness gate; summary() lists them.
#> 
#> Selected variability structure (run10):
#>   IIV:         ETA_CL, ETA_KA, ETA_V
#>   Correlated:  ETA_V ~ ETA_KA

The printed output has three steps. The eta-count stage removes etas from the input structure – top_down_exhaustive tries every subset, not one deletion at a time, so the candidates include structures with two or three etas gone. The block stage then offers covariance blocks over the etas that survived. The third step, compare_to_input, puts the winner of the block stage back against the input structure, and it is there that the selection is made. Candidates are ranked on the BIC(iiv), which is what rank = "bic" means for this tool.

Read the printed output carefully. The input model is a candidate in its own right: the step tables rank it, and in the eta-count stage it ranks first. The strictness gate is separate from that ranking. It rejected 3 of the fitted models, and one of them reached a better criterion than the winner did:

iiv$models[, c("id", "structure", "criterion", "rank", "passed")]
#>       id                 structure  criterion rank passed
#> 1  input   [ETA_CL,ETA_V]+[ETA_KA] -271.27544   NA   TRUE
#> 2   run1         [ETA_CL]+[ETA_KA]  -14.16006    2   TRUE
#> 3   run2            [ETA_CL,ETA_V]   86.79536   NA  FALSE
#> 4   run3          [ETA_KA]+[ETA_V]   27.37871   NA  FALSE
#> 5   run4                  [ETA_CL]  108.64292    4   TRUE
#> 6   run5                  [ETA_KA]   96.89883    3   TRUE
#> 7   run6                   [ETA_V]  178.10690    5   TRUE
#> 8   run7                            181.73650    6   TRUE
#> 9   run8 [ETA_CL]+[ETA_KA]+[ETA_V] -273.45621   NA  FALSE
#> 10  run9   [ETA_CL,ETA_KA]+[ETA_V] -272.17910    2   TRUE
#> 11 run10   [ETA_KA,ETA_V]+[ETA_CL] -272.46800    1   TRUE
#> 12 run11     [ETA_CL,ETA_KA,ETA_V] -269.50116    4   TRUE

run8 scores -273.46 against the selected -272.47, and is out because it did not pass. Note where that verdict lives: converged and passed are columns of the tables that report on fitted candidates – $models here, $candidates and $steps elsewhere – and not of the per-stage ranking summary, which ranks models those rows already judged:

names(iiv$steps)
#> [1] "step"        "kind"        "parent"      "id"          "criterion"  
#> [6] "d_criterion" "rank"        "best"
intersect(c("converged", "passed"), names(iiv$steps))
#> character(0)

The input carries no rank in the combined table above because it is the parent each stage starts from. A search result is a ranking of what passed, not proof that the winner is the best model fitted.

iiv$final_structure
#> [1] "[ETA_KA,ETA_V]+[ETA_CL]"
iiv$final_etas
#> [1] "ETA_CL" "ETA_KA" "ETA_V"

With algorithm = "skip" the eta-count stage is left out and only the block stage runs, which is the cheap way to ask one question: should CL and V share a covariance block? Skipping the eta-count stage means naming the block stage explicitly – asking for neither is an error rather than a search with nothing left to do. The statements name parameters, not etas:

iiv_block <- ferx_iivsearch(model = bo$model, data = bo$data,
                            search_space = "COVARIANCE?(IIV, [CL, V])",
                            algorithm = "skip",
                            correlation_algorithm = "top_down_exhaustive",
                            block_retries = 1,
                            directory = book_tempdir("iivsearch-block"), progress = FALSE)
#> Warning in .ferx_compute_cor_matrix(result$cov_matrix): One or more diagonal
#> elements are non-positive; correlation matrix may not be meaningful.
c(input = iiv_block$input_structure, final = iiv_block$final_structure)
#>                     input                     final 
#> "[ETA_CL,ETA_V]+[ETA_KA]" "[ETA_CL,ETA_V]+[ETA_KA]"
iiv_block$models[, c("id", "structure", "criterion", "rank", "passed")]
#>      id                 structure criterion rank passed
#> 1 input   [ETA_CL,ETA_V]+[ETA_KA] -271.2754   NA   TRUE
#> 2  run1 [ETA_CL]+[ETA_KA]+[ETA_V] -273.4562   NA  FALSE

Two structures are compared here: the input’s block_omega (ETA_CL, ETA_V) and the diagonal alternative. The input is kept – but read why. The diagonal candidate reaches the better criterion of the two, and it is out because it fails the strictness gate. The failures column gives the reason:

with(iiv_block$models, failures[!is.na(failures)])
#> [1] "stalled at the initial estimates: no free parameter moved more than 1% of its initial value (#751)"

No free parameter moved more than 1% of its initial value, so reject_init_stall refused the fit. Whether that is the right call here depends on where the candidate started, which the printed output does not say: a candidate handed its parent’s optimum can sit still because it is already converged, which is why ferx_amd() turns this gate off inside its steps (below). Read a gated-out candidate as unjudged rather than as rejected on merit, and refit it yourself if it matters. A selection is only as meaningful as the gate that produced it.

Inter-occasion variability search

ferx_iovsearch() decides which parameters carry an occasion random effect. It reads the occasion column from the base model’s iov_column, so a model without one has nothing to search. The bundled warfarin_iov configuration offers a kappa on every parameter that has an eta:

iov_ex <- ferx_example("warfarin_iov")
iov <- ferx_iovsearch(config = iov_ex$search, directory = book_tempdir("iovsearch"), progress = FALSE)
iov
#> ferx inter-occasion variability search (same-as-iiv, ranked on bic_random)
#>   Data:  /home/runner/work/_temp/Library/ferx/examples/search/../data/warfarin_iov.csv
#>   Occasions: OCC
#>   Wrote: /tmp/RtmpJuwo6H/iovsearch
#>   Input: IIV([ETA_CL]+[ETA_KA]+[ETA_V]);IOV([KAPPA_CL])
#>   bic_random:   222 (input) -> -313.9 (run5)
#>   7 models over 2 steps
#> 
#> Step 1 (no_of_etas), parent input:
#>  rank    id                                                           structure
#>     1  run1 IIV([ETA_CL]+[ETA_KA]+[ETA_V]);IOV([KAPPA_CL]+[KAPPA_KA]+[KAPPA_V])
#>     2  run2            IIV([ETA_CL]+[ETA_KA]+[ETA_V]);IOV([KAPPA_CL]+[KAPPA_V])
#>     3  run3           IIV([ETA_CL]+[ETA_KA]+[ETA_V]);IOV([KAPPA_CL]+[KAPPA_KA])
#>     4 input                      IIV([ETA_CL]+[ETA_KA]+[ETA_V]);IOV([KAPPA_CL])
#>  criterion d_criterion  best
#>    -311.70       533.7  TRUE
#>      32.34       189.6 FALSE
#>      69.06       152.9 FALSE
#>     222.00         0.0 FALSE
#> 
#> Step 2 (no_of_etas), parent run1:
#>  rank   id                                                           structure
#>     1 run5         IIV([ETA_CL]+[ETA_KA]);IOV([KAPPA_CL]+[KAPPA_KA]+[KAPPA_V])
#>     2 run1 IIV([ETA_CL]+[ETA_KA]+[ETA_V]);IOV([KAPPA_CL]+[KAPPA_KA]+[KAPPA_V])
#>     3 run6                  IIV([ETA_CL]);IOV([KAPPA_CL]+[KAPPA_KA]+[KAPPA_V])
#>     4 run4          IIV([ETA_CL]+[ETA_V]);IOV([KAPPA_CL]+[KAPPA_KA]+[KAPPA_V])
#>  criterion d_criterion  best
#>     -313.9       2.183  TRUE
#>     -311.7       0.000 FALSE
#>     -303.4      -8.302 FALSE
#>     -301.3     -10.480 FALSE
#> 
#> Selected variability structure (run5):
#>   IIV:  ETA_CL, ETA_KA
#>   IOV:  KAPPA_CL, KAPPA_KA, KAPPA_V
#> 
#> Notes:
#>   - CL already carries a κ in the base model and is left as it is

distribution = "same-as-iiv", the default, declares each added kappa the way its eta is declared. The first step added kappas on KA and V; the second then found that dropping ETA_V improved the BIC(random) further: an IOV search can change the IIV structure, because the two are fitted together. CL already carries a kappa in the base model and is left as it is.

iov$final_kappas
#> [1] "KAPPA_CL" "KAPPA_KA" "KAPPA_V"
iov$final_structure
#> [1] "IIV([ETA_CL]+[ETA_KA]);IOV([KAPPA_CL]+[KAPPA_KA]+[KAPPA_V])"

All axes at once

Every tool so far walks one axis with the others held fixed: ferx_modelsearch() decides the structure with the covariate model fixed, ferx_covsearch() the covariates with the structure fixed. That is what keeps them cheap, and it is blind to one thing – an effect that only earns its place after another decision has gone the other way. A clearance covariate that pays off only once 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. Each structural category in the space becomes an axis with its values as alleles; each optional COVARIATE? pair becomes an axis whose alleles are none and each of its forms. algorithm = "exhaustive" fits every point; algorithm = "ga", the default, runs a genetic algorithm 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 carrying WT and CRCL that were simulated from two compartments, so all four axes are live.

gs_ex <- ferx_example("two_cpt_oral_global")
gs <- ferx_globalsearch(config = gs_ex$search,
                        directory = book_tempdir("globalsearch"), progress = FALSE)
gs
#> ferx global model search (exhaustive, ranked on penalized)
#>   Data:  /home/runner/work/_temp/Library/ferx/examples/search/../data/two_cpt_oral_cov.csv
#>   Wrote: /tmp/RtmpJuwo6H/globalsearch
#>   Grid:  16 points over 4 axes
#>     LAGTIME: OFF | ON
#>     PERIPHERALS: 0 | 1
#>     CL-WT: none | power
#>     CL-CRCL: none | power
#>   Input: OFV -171.3, penalized -101.3
#>   Fitness: -101.3 (input) -> -799.5 (run5)
#>   17 models fitted of 17 evaluated
#> 
#> Ranked models (best first):
#>  rank    id                                              genome    ofv
#>     1  run5   LAGTIME=OFF;PERIPHERALS=1;CL-WT=none;CL-CRCL=none -889.5
#>     2  run6  LAGTIME=OFF;PERIPHERALS=1;CL-WT=none;CL-CRCL=power -896.1
#>     3  run7  LAGTIME=OFF;PERIPHERALS=1;CL-WT=power;CL-CRCL=none -895.6
#>     4  run8 LAGTIME=OFF;PERIPHERALS=1;CL-WT=power;CL-CRCL=power -898.8
#>     5 run13    LAGTIME=ON;PERIPHERALS=1;CL-WT=none;CL-CRCL=none -887.4
#>     6 run14   LAGTIME=ON;PERIPHERALS=1;CL-WT=none;CL-CRCL=power -895.7
#>     7 run15   LAGTIME=ON;PERIPHERALS=1;CL-WT=power;CL-CRCL=none -894.6
#>     8 run16  LAGTIME=ON;PERIPHERALS=1;CL-WT=power;CL-CRCL=power -900.9
#>     9  run1   LAGTIME=OFF;PERIPHERALS=0;CL-WT=none;CL-CRCL=none -172.2
#>    10 input                                                <NA> -171.3
#>    11  run2  LAGTIME=OFF;PERIPHERALS=0;CL-WT=none;CL-CRCL=power -179.3
#>    12  run3  LAGTIME=OFF;PERIPHERALS=0;CL-WT=power;CL-CRCL=none -179.2
#>    13  run4 LAGTIME=OFF;PERIPHERALS=0;CL-WT=power;CL-CRCL=power -185.4
#>  criterion fitness converged passed
#>    -799.50 -799.50      TRUE   TRUE
#>    -796.10 -796.10      TRUE   TRUE
#>    -795.60 -795.60      TRUE   TRUE
#>    -788.80 -788.80      TRUE   TRUE
#>    -777.40 -777.40      TRUE   TRUE
#>    -775.70 -775.70      TRUE   TRUE
#>    -774.60 -774.60      TRUE   TRUE
#>    -770.90 -770.90      TRUE   TRUE
#>    -102.20 -102.20      TRUE   TRUE
#>    -101.30 -101.30      TRUE   TRUE
#>     -99.32  -99.32      TRUE   TRUE
#>     -99.21  -99.21      TRUE   TRUE
#>     -95.43  -95.43      TRUE   TRUE
#> 
#> Models the search charged beyond the criterion:
#>     id                                             genome criterion
#>   run9   LAGTIME=ON;PERIPHERALS=0;CL-WT=none;CL-CRCL=none     16.69
#>  run10  LAGTIME=ON;PERIPHERALS=0;CL-WT=none;CL-CRCL=power     19.31
#>  run11  LAGTIME=ON;PERIPHERALS=0;CL-WT=power;CL-CRCL=none     19.53
#>  run12 LAGTIME=ON;PERIPHERALS=0;CL-WT=power;CL-CRCL=power     23.38
#>  charge_non_influential charge_gate charge_crash fitness
#>                       0         100            0   116.7
#>                       0         100            0   119.3
#>                       0         100            0   119.5
#>                       0         100            0   123.4

Read the top two rows before anything else. The runner-up has the better OFV and still ranks second, because it carries one more covariate relation and the criterion charges it for the parameter. A table showing only the OFV would name the wrong winner, which is why the criterion and the OFV are separate columns:

head(gs$models[order(gs$models$fitness), c("id", "genome", "ofv", "criterion", "fitness")], 3)
#>     id                                             genome       ofv criterion
#> 6 run5  LAGTIME=OFF;PERIPHERALS=1;CL-WT=none;CL-CRCL=none -889.4714 -799.4714
#> 7 run6 LAGTIME=OFF;PERIPHERALS=1;CL-WT=none;CL-CRCL=power -896.1283 -796.1283
#> 8 run7 LAGTIME=OFF;PERIPHERALS=1;CL-WT=power;CL-CRCL=none -895.5592 -795.5592
#>     fitness
#> 6 -799.4714
#> 7 -796.1283
#> 8 -795.5592

The criterion is not the whole of the ranking

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

Charge	Paid by
`charge_non_influential`	a gene that changed nothing in the rendered model, such as 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. Large, and finite, so the algorithm can still rank it.

They are columns beside the criterion and the fitness they add up to, which is the “never a silent drop” rule of the other tools taken one step further: a gated genome does not merely appear in the table, it says what the gate cost it.

charged <- gs$models[gs$models$charge_gate > 0, ]
charged[, c("id", "genome", "criterion", "charge_gate", "fitness", "passed")]
#>       id                                             genome criterion
#> 10  run9   LAGTIME=ON;PERIPHERALS=0;CL-WT=none;CL-CRCL=none  16.69030
#> 11 run10  LAGTIME=ON;PERIPHERALS=0;CL-WT=none;CL-CRCL=power  19.30548
#> 12 run11  LAGTIME=ON;PERIPHERALS=0;CL-WT=power;CL-CRCL=none  19.52656
#> 13 run12 LAGTIME=ON;PERIPHERALS=0;CL-WT=power;CL-CRCL=power  23.38225
#>    charge_gate  fitness passed
#> 10         100 116.6903  FALSE
#> 11         100 119.3055  FALSE
#> 12         100 119.5266  FALSE
#> 13         100 123.3823  FALSE
gs$penalties[c("non_influential", "gate", "crash")]
#> non_influential            gate           crash 
#>         0.00001       100.00000  99999999.00000

Without that column the first of those rows reads as a model that simply fitted worse than the winner, rather than one the gate refused and that sits 100 units away because of it.

The genetic algorithm

Sixteen points is small enough to enumerate, which is what makes the bundled file a reference: on a grid this size the algorithm has to find the same winner. On a grid of thousands, enumerating is not an option and finding it at all is the point.

ga_config <- file.path(book_tempdir("globalsearch-ga-config"), "global.ferxsearch")
ga_lines <- readLines(gs_ex$search)
ga_lines <- sub("^base = .*", sprintf('base = "%s"', normalizePath(gs_ex$model, winslash = "/")), ga_lines)
ga_lines <- sub("^data = .*", sprintf('data = "%s"', normalizePath(gs_ex$data, winslash = "/")), ga_lines)
ga_lines <- sub('algorithm = "exhaustive"', 'algorithm = "ga"', ga_lines, fixed = TRUE)
writeLines(ga_lines, ga_config)
gs_ga <- ferx_globalsearch(config = ga_config,
                           directory = book_tempdir("globalsearch-ga"), progress = FALSE)
gs_ga$generations
#>   index best best_fitness mean_fitness polished
#> 1     0 run5    -799.4714    -384.6412        0
#> 2     1 run5    -799.4714    -499.5705        0
#> 3     2 run5    -799.4714    -504.5793        0
#> 4     3 run5    -799.4714    -617.1745        1
data.frame(
  algorithm = c("exhaustive", "ga"),
  genome = c(gs$models$genome[gs$models$selected],
             gs_ga$models$genome[gs_ga$models$selected]),
  fitness = c(gs$final_fitness, gs_ga$final_fitness)
)
#>    algorithm                                            genome   fitness
#> 1 exhaustive LAGTIME=OFF;PERIPHERALS=1;CL-WT=none;CL-CRCL=none -799.4714
#> 2         ga LAGTIME=OFF;PERIPHERALS=1;CL-WT=none;CL-CRCL=none -799.4714

ga takes a named list checked 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 decision to enumerate a grid nobody has sized – and an exhaustive grid above [globalsearch] max_models is an error naming the size, never a truncated search.

ferx_globalsearch() is not 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 is shared by the engine, the .ferxsearch file ([globalsearch]) and the ferx globalsearch command, so a search stays portable between R and the command line.

The whole pipeline at once

ferx_amd() runs the tools above in order as one automatic model-development pipeline: structural, IIV, residual, IOV, allometry, covariates. One search space describes the whole pipeline, and the engine hands each step only the statements its tool can read. A step the space says nothing about is skipped, and says so.

ferx_amd_plan() shows that plan without fitting anything, which is the cheap way to check a configuration before spending the run:

amd_ex <- ferx_example("warfarin_amd")
plan <- ferx_amd_plan(config = amd_ex$search)
plan
#>   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

3 of the 6 steps are already ruled out by the model and the space, and the skipped column says why: the model declares no iov_column, and the space carries no ALLOMETRY or COVARIATE statement because warfarin.csv has no covariate column. Running it fits the rest:

amd <- ferx_amd(config = amd_ex$search, directory = book_tempdir("amd"), progress = FALSE)
amd
#> ferx AMD pipeline (default strategy, retries all_final)
#>   Model: /home/runner/work/_temp/Library/ferx/examples/search/../models/warfarin_amd.ferx
#>   Data:  /home/runner/work/_temp/Library/ferx/examples/search/../data/warfarin.csv
#>   Wrote: /tmp/RtmpJuwo6H/amd
#>   OFV:   -286.004 (start) -> -286.004 (final, dOFV -0.000)
#>   3 of 6 steps ran, 28 models fitted, 8.6 s across the steps
#> 
#> Pipeline:
#>  index       step  status criterion value_before value_after          d_value
#>      1 structural     ran bic_mixed       -267.5      -267.5  0.0000000008589
#>      2  iivsearch     ran   bic_iiv       -279.1      -279.1 -0.0000000011280
#>      3   residual     ran       ofv       -286.0      -286.0  0.0000000001716
#>      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
#>   0.0000000008589          5      TRUE   TRUE   1.765
#>  -0.0000000011280         12      TRUE   TRUE   1.648
#>   0.0000000001716          7        NA     NA   5.142
#>                NA          0        NA     NA   0.000
#>                NA          0        NA     NA   0.000
#>                NA          0        NA     NA   0.000
#> 
#> 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
#> 
#> Notes:
#>   - [strictness] reject_init_stall is off inside every step: the model a step is handed starts at its own optimum, where "did not leave the initial estimates" (#751) is the expected outcome rather than a failed fit. It still applies to the pipeline's own start fit, and every other gate applies throughout
#>   - the retries pass on `03-ruvsearch-retries` did not improve the fit (OFV -286.004 against -286.004); the selected model was kept

The pipeline improved nothing: warfarin_amd.ferx is already the model these three steps would choose, so each step re-selected its own input and the final OFV equals the starting one. That is a useful result to see printed rather than hidden – the step table carries each step’s criterion before and after, its candidate count, the strictness verdict and its wall clock, so a step that “won” by having its siblings rejected looks different from one that found something. summary() adds every candidate of every step, and the same two tables are written to the run directory as steps.csv and candidates.csv.

Getting the result out of an AMD run takes a little care, because two fields both look final and are not the same thing. $final_model is the model text the pipeline ended on, written to $final_model_path; the fitted object is $fit. $final_ofv and $d_ofv give the end point and the change from the start, negative being an improvement:

c(final_model = class(amd$final_model), fit = class(amd$fit))
#> final_model         fit 
#> "character"  "ferx_fit"
c(input_ofv = amd$input_ofv, final_ofv = amd$final_ofv, d_ofv = amd$d_ofv)
#>              input_ofv              final_ofv                  d_ofv 
#> -286.00421950648029679 -286.00421950657749903   -0.00000000009720225
basename(amd$final_model_path)
#> [1] "final.ferx"

Each step that ran also has an entry under $tools, named by its position and pipeline step. Watch the naming: the entry is keyed on the pipeline step, while its $directory is the tool that implemented it, and the two need not agree – 01-structural was run by modelsearch:

names(amd$tools)
#> [1] "01-structural" "02-iivsearch"  "03-residual"
vapply(amd$tools, function(step) step$directory, "")
#>    01-structural     02-iivsearch      03-residual 
#> "01-modelsearch"   "02-iivsearch"   "03-ruvsearch"
names(amd$tools[["01-structural"]])
#>  [1] "index"            "step"             "tool"             "directory"       
#>  [5] "status"           "reason"           "criterion"        "selected"        
#>  [9] "seconds"          "candidates"       "path"             "models"          
#> [13] "input_model_path" "final_model_path" "model_paths"

An entry is a plain list, not the object the tool would return on its own: it carries what the step decided (status, reason, criterion, selected, seconds), its slice of the candidates, and the engine’s own table read back from the run directory. Run a tool directly when you want its fit and its own options.

The pipeline can also be stated inline, with a shorter order and steps left out. strategy chooses which steps run and in what order, and ferx_amd_plan() will show you each without fitting anything:

strategies <- c("default", "reevaluation", "SIR", "SRI", "RSI")
data.frame(strategy = strategies, order = vapply(strategies, function(s) {
  paste(ferx_amd_plan(model = amd_ex$model, data = amd_ex$data,
                      search_space = "PERIPHERALS(0..1); IIV?(@PK, exp)",
                      strategy = s)$step, collapse = " -> ")
}, ""), row.names = NULL)
#>       strategy
#> 1      default
#> 2 reevaluation
#> 3          SIR
#> 4          SRI
#> 5          RSI
#>                                                                                                  order
#> 1                          structural -> iivsearch -> residual -> iovsearch -> allometry -> covariates
#> 2 structural -> iivsearch -> residual -> iovsearch -> allometry -> covariates -> iivsearch -> residual
#> 3                                                                  structural -> iivsearch -> residual
#> 4                                                                  structural -> residual -> iivsearch
#> 5                                                                  residual -> structural -> iivsearch

reevaluation is the default order with a second iivsearch and residual pass at the end. The three-letter names are shorter pipelines, not reorderings of the same six: they run the three steps their initials spell and leave out iovsearch, allometry and covariates entirely. skip drops a step from whichever order was chosen; retries_on says when the extra perturbed-restart pass runs:

amd_sir <- ferx_amd(model = amd_ex$model, data = amd_ex$data,
                    search_space = "PERIPHERALS(0..1); IIV?(@PK, exp)",
                    strategy = "SIR", skip = "residual", retries_on = "final",
                    retries = 0, directory = book_tempdir("amd-sir"), progress = FALSE)
amd_sir$steps[, c("step", "status", "criterion", "candidates", "seconds")]
#>         step  status criterion candidates   seconds
#> 1 structural     ran bic_mixed          2 0.2819569
#> 2  iivsearch     ran   bic_iiv          8 0.1606043
#> 3   residual skipped      <NA>          0 0.0000000

Configuration sections with no R tool

A .ferxsearch file may carry a section addressed to a tool this package does not bind. The loader validates section names against a fixed vocabulary, which runs ahead of what ferx-r binds and, in the case of [structsearch], ahead of the engine itself; a name outside that list is an error. The file loads, and ferx warns rather than running a search that quietly ignores part of its own configuration:

sections_dir <- book_tempdir("search-sections")
unconsumed <- function(section, body) {
  path <- file.path(sections_dir, paste0(gsub("\\[|\\]", "", section), ".ferxsearch"))
  writeLines(c(readLines(base$search), "", section, body), path)
  tryCatch(ferx_search_config(path), warning = conditionMessage)
}
unconsumed("[structsearch]", "  x = 1")
#> [1] "ferx_search_config: [structsearch] has no R tool in this package, so it is ignored - a search run from this file does what its other sections say."

[structsearch] is accepted vocabulary with no engine module yet, so it is reported without a remediation that would name a command nobody can run.

What the warning reports is the complement of the sections ferx-r has a tool for, not a list of known-bad names. 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 – which is what happened to [globalsearch], warned about until ferx_globalsearch() shipped and silent since. See the ferx-core tools pages for what the engine knows.

Options that matter

Function	Arguments
`ferx_bic()`	`fit`, `type` (`"mixed"`, `"fixed"`, `"iiv"`, `"random"`)
`check_strictness()`	`fit`, `require_converged`, `require_covariance`, `max_condition_number`, `max_correlation`, `reject_on_boundary`, `reject_init_stall`
`ferx_search_config()`	`path`
`ferx_search_space()`	`mfl` (text, config or space), `model`, `data`
`ferx_search_coverage()`	`space` (space, config or MFL text)
`ferx_search_results()`	`directory`, `partial` (`NULL` prefers the complete table), `type` (`"candidates"`, `"models"`, `"steps"`)

The search functions share the entry arguments model, data and config; the run arguments threads, retries, directory, resume and progress; and, wherever candidates are ranked rather than tested, search_space, rank and cutoff. A search is stated either by config or by the inline arguments, not both:

Function	Search arguments
`ferx_covsearch()`	`algorithm` (`"scm-forward-then-backward"` or `"scm-forward"`), `p_forward`, `p_backward`, `max_steps`, `adaptive_scope_reduction`
`ferx_modelsearch()`	`algorithm` (`"reduced_stepwise"`, `"exhaustive_stepwise"`, `"exhaustive"`), `iiv_strategy` (`"absorption_delay"`, `"add_diagonal"`, `"no_add"`)
`ferx_ruvsearch()`	`groups`, `p_value`, `skip` (families not to test), `max_iter` (1 to 3), `cwres_prescreen` (screen candidates on the parent’s CWRES before refitting)
`ferx_iivsearch()`	`algorithm` (`"top_down_exhaustive"`, `"bottom_up_stepwise"`, `"simultaneous_stepwise"`, `"skip"`), `correlation_algorithm` (the block stage: `"top_down_exhaustive"` or `"skip"`), `as_fullblock` (a bottom-up candidate blocks every eta it carries), `block_retries` (extra starts per eta in the largest block)
`ferx_iovsearch()`	`column` (the occasion column; defaults to the model’s `iov_column`, and disagreeing with it is an error), `distribution` (`"same-as-iiv"`, `"disjoint"`, `"joint"`, `"explicit"`), `groups` (the blocks `"explicit"` names), `block_retries`
`ferx_globalsearch()`	`algorithm` (`"ga"` or `"exhaustive"`, matched exactly), `iiv_strategy`, `max_models` (the cap on an exhaustive grid), `ga` (a named list of `[globalsearch.ga]` knobs), `penalties` (a named list of `[rank.penalties]` charges)
`ferx_amd()`, `ferx_amd_plan()`	`strategy` (`"default"`, `"reevaluation"`, `"SIR"`, `"SRI"`, `"RSI"`), `skip` (steps to leave out), `retries_on` (`"all_final"`, `"final"`, `"skip"`)

Only a .ferxsearch file can carry a per-tool section ([modelsearch], [iivsearch], [covsearch], …); the inline form runs each step at its own defaults.

print() and summary() exist for the results of all of them: ferx_covsearch, ferx_modelsearch, ferx_ruvsearch, ferx_iivsearch, ferx_iovsearch, ferx_globalsearch and ferx_amd objects. ferx_search_config and ferx_search_space objects have print() methods.

Pitfalls

Keep search directories. The book writes them to a temporary directory. In a project, write them next to the analysis: they hold every candidate model, the tables and final.ferx, and they make resume possible.
ferx_search_results() needs the right type. The covariate and residual error searches write steps.csv, the structural search models.csv. Neither wrote a candidates.csv here, so the default type = "candidates" finds no table for these runs.
Compare like with like. OFV differences and likelihood ratio tests need the same data, the same estimation method and nested models.
Thresholds change the answer. Record the space, the significance levels and the criterion together with the selected model. The .ferxsearch file is that record.

Summary

Compare fitted models by likelihood ratio test (nested) or with fit$aic, fit$bic and ferx_bic(). Gate candidates with check_strictness().
Describe searches in .ferxsearch files; check them with ferx_search_config(), ferx_search_space() and ferx_search_coverage().
Run covariate, structural and residual error searches with ferx_covsearch(), ferx_modelsearch() and ferx_ruvsearch(), and read their tables with ferx_search_results().
Search the random-effect structure with ferx_iivsearch() and the occasion structure with ferx_iovsearch(); both can change the IIV structure, so read what they selected rather than only the criterion.
Decide structure and covariates together with ferx_globalsearch(), which searches them as one grid rather than one axis at a time. It ranks on penalized fitness by default, and the charges it adds to the criterion are columns beside it.
Run the whole order as one pipeline with ferx_amd(), after checking it with ferx_amd_plan(). A step the space says nothing about is skipped, and the step table says why.

Next: Parameter uncertainty quantifies the uncertainty of the selected covariate model.

Reference

R help: ?ferx_bic, ?check_strictness, ?ferx_search_config, ?ferx_search_space, ?ferx_search_coverage, ?ferx_search_results, ?ferx_covsearch, ?ferx_modelsearch, ?ferx_ruvsearch, ?ferx_iivsearch, ?ferx_iovsearch, ?ferx_globalsearch, ?ferx_amd, ?ferx_amd_plan
ferx-core: search configuration, covsearch, modelsearch, ruvsearch, iivsearch, iovsearch, global search, amd, tools