Variability: random effects, residual error and IOV

Where you are

A population model splits the spread in the data into levels:

Between subjects: random effects (etas) with variances omega.
Between occasions within a subject: inter-occasion variability, random effects (kappas) with variances in omega_iov.
Within an occasion: residual error, with the sigma parameters of the error model.

This chapter shows how each level is declared in a ferx model, and how to read the estimates. It also covers two estimation options that act on the random effects: mu-referencing and parameter scaling.

Maturity: stable, beta and experimental parts

omega/block_omega, inter-occasion variability and the additive, proportional and combined error models are stable in ferx-core. Individual-parameter expressions are beta. Error magnitudes that depend on time or covariates, residual weighting and mixture models are experimental. See feature maturity.

The data

Most examples use the warfarin data (10 subjects, one oral dose each). warfarin_block_omega, warfarin_additive_eta and warfarin_ltbs ship their own copies of it, with the same concentrations:

w <- ferx_example("warfarin")
variability_dir <- book_tempdir("variability")
warfarin_dv <- read.csv(w$data, na.strings = ".")$DV
sapply(c("warfarin_block_omega", "warfarin_additive_eta", "warfarin_ltbs"), function(name) {
  isTRUE(all.equal(read.csv(ferx_example(name)$data, na.strings = ".")$DV, warfarin_dv))
})
#>  warfarin_block_omega warfarin_additive_eta         warfarin_ltbs 
#>                  TRUE                  TRUE                  TRUE

The IOV examples use warfarin_iov.csv: the same 10 subjects dosed twice, 120 hours apart, with an OCC column that numbers the occasions:

iov <- ferx_example("warfarin_iov")
iov_data <- read.csv(iov$data, na.strings = ".")
filter(iov_data, ID == 1, EVID == 1)
#>   ID TIME DV EVID AMT CMT RATE MDV OCC
#> 1  1    0 NA    1 100   1    0   1   1
#> 2  1  120 NA    1 100   1    0   1   2
count(filter(iov_data, EVID == 0), OCC)
#>   OCC   n
#> 1   1 110
#> 2   2 110

Minimal runnable call: between-subject variability

The warfarin model gives each individual parameter a log-normal random effect: CL = TVCL * exp(ETA_CL). Each eta has a variance declared with omega, and the residual error is proportional:

invisible(ferx_model_get_section(w$model, "parameters"))
#> # [parameters]
#>   theta TVCL(0.134, 0.001, 10.0)
#>   theta TVV(8.1, 0.1, 500.0)
#>   theta TVKA(1.0, 0.01, 50.0)
#> 
#>   omega ETA_CL ~ 0.07
#>   omega ETA_V  ~ 0.02
#>   omega ETA_KA ~ 0.40
#> 
#>   sigma PROP_ERR ~ 0.01 (sd)
invisible(ferx_model_get_section(w$model, "individual_parameters"))
#> # [individual_parameters]
#>   CL = TVCL * exp(ETA_CL)
#>   V  = TVV  * exp(ETA_V)
#>   KA = TVKA * exp(ETA_KA)
invisible(ferx_model_get_section(w$model, "error_model"))
#> # [error_model]
#>   DV ~ proportional(PROP_ERR)
fit_w <- ferx_fit(w$model, w$data, verbose = FALSE)

Reading the result

The variability part of the printout shows each omega as a variance with its coefficient of variation, the sigma on the standard-deviation scale, and shrinkage:

printout <- capture.output(print(fit_w))
printout[grep("^OMEGA", printout):(grep("^DIAGNOSTICS", printout) - 1)]
#>  [1] "OMEGA  (between-subject variability)"                                                                               
#>  [2] "------------------------------------------------------------"                                                       
#>  [3] "  ETA_CL                   [log-normal]  = 0.028595  CV% = 17.0  SE = 0.012798"                                     
#>  [4] "  ETA_V                    [log-normal]  = 0.009577  CV% = 9.8  SE = 0.004297"                                      
#>  [5] "  ETA_KA                   [log-normal]  = 0.348964  CV% = 64.6  SE = 0.160819"                                     
#>  [6] ""                                                                                                                   
#>  [7] "SIGMA  (residual error)"                                                                                            
#>  [8] "------------------------------------------------------------"                                                       
#>  [9] "  PROP_ERR         [proportional] = 0.010749  (var = 0.000116, CV% = 1.1)  SE = 0.000942  [initial specified as SD]"
#> [10] ""                                                                                                                   
#> [11] "SHRINKAGE"                                                                                                          
#> [12] "------------------------------------------------------------"                                                       
#> [13] " ETA_CL: 0.0%   ETA_V: 0.1%   ETA_KA: 0.2%   EPS: 16.1%"                                                            
#> [14] ""

As data, the same numbers are fit$omega (the covariance matrix), fit$sigma, fit$shrinkage_eta and the individual random effects fit$ebe_etas:

fit_w$omega
#>            ETA_CL       ETA_V    ETA_KA
#> ETA_CL 0.02859495 0.000000000 0.0000000
#> ETA_V  0.00000000 0.009576924 0.0000000
#> ETA_KA 0.00000000 0.000000000 0.3489642
fit_w$sigma
#> [1] 0.01074853
fit_w$shrinkage_eta
#> [1] 0.0003815202 0.0010733922 0.0016059600
head(fit_w$ebe_etas, 3)
#>   ID      ETA_CL       ETA_V     ETA_KA
#> 1  1  0.02794924  0.06598698  0.4815703
#> 2  2  0.01537051 -0.06325551 -0.3010812
#> 3  3 -0.02470858  0.09101450  0.2638643

Tables and figures turns these into a report table. Diagnosing the model covers shrinkage and eta diagnostics.

Options that matter

Level	Declaration	Meaning	Example
Between subjects	`omega ETA_CL ~ 0.07`	One random effect with its variance	`warfarin`
	`omega ETA_CL ~ 0.2645751 (sd)`	The same variance declared as a standard deviation	(below)
	`omega ETA_CL ~ 0.0 FIX`	A variance held fixed (at 0: no variability on that parameter)	(below)
	`block_omega (ETA_CL, ETA_V) = [0.07, 0.02, 0.02]`	Correlated random effects (lower triangle)	`warfarin_block_omega`
Between occasions	`kappa KAPPA_CL ~ 0.04`	One IOV random effect per occasion	`warfarin_iov`
	`block_kappa (KAPPA_CL, KAPPA_V) = [...]`	Correlated IOV random effects	(variant below)
	`[fit_options]` `iov_column`, `iov_occasion`	Where occasions come from	`warfarin_iov`
Residual	`DV ~ additive(S)`, `proportional(S)`, `combined(S1, S2)`	Standard error models	(variants below)
	`DV ~ power(S, THETA)`	Proportional error raised to an estimated power	(variant below)
	`log(DV) ~ additive(S)`, `DV ~ log_additive(S)`	Additive error on the log scale	`warfarin_ltbs`
	`block_sigma (S1, S2) = [...]`	Correlated residual components	(variant below)
	`iiv_on_ruv = ETA_RUV`	A random effect that scales each subject’s residual error	(variant below)
	`CMT=2: DV ~ ...`	One error model per observed compartment	PK/PD and multiple endpoints

A random effect can enter an individual parameter in any expression, and the expression is what decides the parameter’s domain. Which form to reach for:

The parameter is	Write it as	`transform` of the theta
Always positive, varying multiplicatively	`TVCL * exp(ETA_CL)`	`identity`
Near zero, or able to go negative	`TVTLAG + ETA_TLAG`	`identity`
Bounded to (0, 1), typical value on that scale	`inv_logit(logit(THETA_F) + ETA_F)`	`logit_probability`
Bounded to (0, 1), typical value on the logit scale	`inv_logit(THETA_F + ETA_F)`	`logit`
Typical value declared on the log scale	`exp(TVCL + ETA_CL)`	`log`

The last three change how the theta is parameterised, which Initial estimates and a first fit covers – log and logit are also reported on the transformed scale, while a logit_probability theta is reported on (0, 1). The first two do not, and are the forms this chapter uses. The three transformed rows are recognised by their exact shape: wrapping one in anything else, as in 2.0 * inv_logit(logit(TVKA) + ETA_KA), takes the theta back to identity and the eta to custom. The log-normal row is more forgiving, and tolerates the extra factors a covariate model puts around it.

Estimation options that act on these parameters:

Where	Name	Meaning
`ferx_fit()`	`mu_referencing =`	Re-centre the individual random-effect search on the population value (default from the model file, engine default `TRUE`)
`ferx_fit()`	`scale_params =`	Legacy switch for `parameter_scaling = "abs"`

book_settings_table(c("iov_column", "iov_occasion", "mu_referencing", "parameter_scaling", "scale_params"))

Table 16.1

Setting	Values	Default	Description
`iov_column`	string	—	Name of the occasion column in the dataset (e.g. `OCC`). Supplies occasions when the model uses `kappa` or `block_kappa` declarations.
`iov_occasion`	`column`, `dose`, `time(t₁, …)`	`column`	Derive the IOV occasion partition from the model instead of a data column.
`mu_referencing`	`true`, `false`	`true`	Re-centre inner-loop ETA estimates on the current population mean (auto-detected from `[individual_parameters]`).
`parameter_scaling`	`auto`, `none`, `abs`, `rescale2`	`auto`	Parameter-scaling strategy for the outer optimizer (supersedes `scale_params` when non-`none`).
`scale_params`		`false`	Legacy boolean alias for `parameter_scaling = abs`. Divide each packed (log/Cholesky) coordinate by its initial magnitude before passing it to the optimizer.

Variance or standard deviation?

An initial value for omega, kappa or sigma is a variance unless you say otherwise. Appending (sd) declares a standard deviation instead, which the parser squares; (variance) states the default explicitly. So these declare the same model, and the warfarin example’s omega ETA_CL ~ 0.07 could equally be written as an SD of 0.2645751:

as_sd <- file.path(variability_dir, "warfarin_omega_sd.ferx")
writeLines(sub("  omega ETA_CL ~ 0.07", "  omega ETA_CL ~ 0.2645751 (sd)",
               readLines(w$model), fixed = TRUE), as_sd)
sigma_as_variance <- file.path(variability_dir, "warfarin_sigma_variance.ferx")
writeLines(sub("  sigma PROP_ERR ~ 0.01 (sd)", "  sigma PROP_ERR ~ 0.0001",
               readLines(w$model), fixed = TRUE), sigma_as_variance)
explicit_variance <- file.path(variability_dir, "warfarin_omega_variance_tag.ferx")
writeLines(sub("  omega ETA_CL ~ 0.07", "  omega ETA_CL ~ 0.07 (variance)",
               readLines(w$model), fixed = TRUE), explicit_variance)
fit_as_sd <- ferx_fit(as_sd, w$data, verbose = FALSE)
fit_sigma_var <- ferx_fit(sigma_as_variance, w$data, verbose = FALSE)
fit_explicit_var <- ferx_fit(explicit_variance, w$data, verbose = FALSE)
c(bundled = fit_w$ofv, omega_as_sd = fit_as_sd$ofv,
  sigma_as_variance = fit_sigma_var$ofv, omega_variance_tag = fit_explicit_var$ofv)
#>            bundled        omega_as_sd  sigma_as_variance omega_variance_tag 
#>           -280.364           -280.364           -280.364           -280.364

The two scales differ only in what you type. What comes back is reported on a fixed scale, whichever form you declared: an omega row is a variance, and a sigma row is a standard deviation. print() on a fit spells the residual error out all three ways, and records which form the file used:

grep("PROP_ERR", capture.output(print(fit_w)), value = TRUE)
#> [1] "  PROP_ERR         [proportional] = 0.010749  (var = 0.000116, CV% = 1.1)  SE = 0.000942  [initial specified as SD]"
fit_w$estimates[, c("transform", "estimate", "init_as_sd")]
#>             transform    estimate init_as_sd
#> TVCL         identity 0.132969030      FALSE
#> TVV          identity 7.730699888      FALSE
#> TVKA         identity 0.725207009      FALSE
#> ETA_CL       variance 0.028594955      FALSE
#> ETA_V        variance 0.009576924      FALSE
#> ETA_KA       variance 0.348964240      FALSE
#> PROP_ERR proportional 0.010748528       TRUE

init_as_sd is the declaration, not the reported scale (Initial estimates and a first fit): it is TRUE for the bundled sigma PROP_ERR ~ 0.01 (sd) and FALSE for every variance-scale row.

A scale tag on a block is rejected

block_omega, block_kappa and block_sigma have no scale to choose: every entry of the lower triangle is a variance or a covariance, always. A single (sd) cannot say which of those entries it applies to, so a scale tag on a block is refused rather than interpreted. ferx_model_validate() reports the file invalid, with a code and the repair:

block_sd <- file.path(variability_dir, "warfarin_block_sd.ferx")
block_ex <- ferx_example("warfarin_block_omega")
writeLines(sub("block_omega (ETA_CL, ETA_V) = [0.07, 0.02, 0.02]",
               "block_omega (ETA_CL, ETA_V) = [0.07, 0.02, 0.02] (sd)",
               readLines(block_ex$model), fixed = TRUE), block_sd)
block_sd_check <- ferx_model_validate(block_sd)
#> Validating: warfarin_block_sd.ferx 
#> 
#> Sections present:
#>   parameters                     [ok]
#>   individual_parameters          [ok]
#>   structural_model               [ok]
#>   error_model                    [ok]
#>   fit_options                    [ok] (optional)
#> 
#> Result: INVALID
#>   * ERROR E_BLOCK_VARIANCE_ONLY: [parameters]: `block_omega` is variance-only — the scale tag `(sd)` is not accepted on a block declaration, since the lower triangle mixes variances and covariances and one tag cannot say which entry is on which scale. Square each SD into a variance and write the off-diagonals as covariances. Offending line: `block_omega (ETA_CL, ETA_V) = [0.07, 0.02, 0.02] (sd)`.
#>       hint: square each SD into a variance and write the off-diagonals as covariances
block_sd_check$ok
#> [1] FALSE
block_sd_check$diagnostics[, c("code", "suggestion")]
#>                    code
#> 1 E_BLOCK_VARIANCE_ONLY
#>                                                                  suggestion
#> 1 square each SD into a variance and write the off-diagonals as covariances

ferx_fit() refuses the same file, with the message but not the code (Writing and managing model files):

try(ferx_fit(block_sd, block_ex$data, verbose = FALSE))
#> Error in ferx_rust_fit(model_path = normalizePath(model), data_path = normalizePath(data),  : 
#>   Error parsing model: [parameters]: `block_omega` is variance-only — the scale tag `(sd)` is not accepted on a block declaration, since the lower triangle mixes variances and covariances and one tag cannot say which entry is on which scale. Square each SD into a variance and write the off-diagonals as covariances. Offending line: `block_omega (ETA_CL, ETA_V) = [0.07, 0.02, 0.02] (sd)`. [E_BLOCK_VARIANCE_ONLY]

The suggestion depends on which tag you wrote, because the two repairs are not the same edit. After (sd) the numbers are wrong and have to change: square each SD into a variance, and write the off-diagonals as covariances rather than as correlations or SD products. After (variance) or (var) only the tag is wrong – the lower triangle was already on that scale, so deleting the tag leaves the same model:

do.call(rbind, lapply(c("variance", "var"), function(tag) {
  tagged <- file.path(variability_dir, paste0("warfarin_block_", tag, ".ferx"))
  writeLines(sub("block_omega (ETA_CL, ETA_V) = [0.07, 0.02, 0.02]",
                 paste0("block_omega (ETA_CL, ETA_V) = [0.07, 0.02, 0.02] (", tag, ")"),
                 readLines(block_ex$model), fixed = TRUE), tagged)
  # The full check report is printed above for (sd); here only the repair differs,
  # so the report itself is captured away and just the diagnostic is shown.
  invisible(capture.output(check <- ferx_model_validate(tagged)))
  data.frame(tag = paste0("(", tag, ")"), check$diagnostics[, c("code", "suggestion")])
}))
#>          tag                  code
#> 1 (variance) E_BLOCK_VARIANCE_ONLY
#> 2      (var) E_BLOCK_VARIANCE_ONLY
#>                                                                                              suggestion
#> 1 delete the tag: the lower triangle is already variances and covariances, so the numbers do not change
#> 2 delete the tag: the lower triangle is already variances and covariances, so the numbers do not change

E_BLOCK_VARIANCE_ONLY is raised only when deleting the tag leaves a line that parses, so the suggestion is always sufficient on its own; a tag sitting next to other stray text is an ordinary parse error instead. The rule is the same for block_kappa and block_sigma, which are variance-only for the same reason – see the ferx-core parameters page. The diagonal omega, kappa and sigma forms above are the ones that take a scale tag.

Fixing a parameter

FIX after a declaration holds a parameter at its initial value. It works on a theta, on an omega or kappa variance, and on a sigma. Two reasons to reach for it: the data cannot identify the parameter, so you hold it at a value from elsewhere (Absorption and bioavailability fixes a bioavailability that way), or you want no variability at all on a parameter.

A variance of zero is the second case, and it has to be written FIX:

no_iiv_on_ka <- file.path(variability_dir, "warfarin_no_iiv_ka.ferx")
writeLines(sub("  omega ETA_KA ~ 0.40", "  omega ETA_KA ~ 0.0 FIX",
               readLines(w$model), fixed = TRUE), no_iiv_on_ka)
fit_no_iiv_ka <- ferx_fit(no_iiv_on_ka, w$data, verbose = 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(all_three = fit_w$n_parameters, ka_fixed = fit_no_iiv_ka$n_parameters)
#> all_three  ka_fixed 
#>         7         6
c(all_three = fit_w$ofv, ka_fixed = fit_no_iiv_ka$ofv)
#>  all_three   ka_fixed 
#> -280.36396   80.40003

Fixing the KA variance saves one parameter and costs 360.8 OFV on these data, so that random effect is clearly worth keeping. The variability it carried does not vanish, it moves: the CL variance goes from 0.029 to 1.27, and the residual error grows 21-fold, from 0.0107 to 0.231. Absorption variability the model can no longer express has to come out somewhere.

A fixed parameter still has a row in fit$estimates. It carries the value the engine stores for a declared zero, with a standard error of 0, and it no longer counts towards n_parameters:

fit_no_iiv_ka$estimates[c("ETA_CL", "ETA_KA"), c("estimate", "se", "rse_pct")]
#>          estimate        se rse_pct
#> ETA_CL 1.27417294 0.5930408 46.5432
#> ETA_KA 0.00000001 0.0000000  0.0000

Leaving the same zero free is refused rather than fitted: the optimizer would start on its lower rail and could not move off it (Dosing regimens and exposure metrics shows the message).

Variants

Correlated random effects

warfarin_block_omega estimates a covariance between ETA_CL and ETA_V, and keeps ETA_KA as a separate diagonal omega:

block <- ferx_example("warfarin_block_omega")
invisible(ferx_model_get_section(block$model, "parameters"))
#> # [parameters]
#>   theta TVCL(0.134, 0.001, 10.0)
#>   theta TVV(8.1, 0.1, 500.0)
#>   theta TVKA(1.0, 0.01, 50.0)
#> 
#>   block_omega (ETA_CL, ETA_V) = [0.07, 0.02, 0.02]
#>   omega ETA_KA ~ 0.40
#> 
#>   sigma PROP_ERR ~ 0.01 (sd)
fit_block <- ferx_fit(block$model, block$data, verbose = FALSE)
#> Warning in .ferx_compute_cor_matrix(result$cov_matrix): One or more diagonal
#> elements are non-positive; correlation matrix may not be meaningful.
fit_block$omega
#>             ETA_CL       ETA_V    ETA_KA
#> ETA_CL 0.028598598 0.001825408 0.0000000
#> ETA_V  0.001825408 0.009577901 0.0000000
#> ETA_KA 0.000000000 0.000000000 0.3488793

The printout reports each covariance with the correlation of the parameters themselves. For two log-normal random effects that is (exp(ω_ij) − 1) / √((exp(ω_ii) − 1)(exp(ω_jj) − 1)). fit$omega_param_corr holds the matrix:

printout <- capture.output(print(fit_block))
grep("~", printout, value = TRUE)
#> [1] "  ETA_V ~ ETA_CL : cov = 0.001825  (param corr = 0.1093)  SE = 0.005278"
round(fit_block$omega_param_corr, 3)
#>       [,1]  [,2] [,3]
#> [1,] 1.000 0.109    0
#> [2,] 0.109 1.000    0
#> [3,] 0.000 0.000    1

A partial block estimates only what it declares. ETA_KA keeps a variance of its own, but the model gives it no covariance with the other two, so those two off-diagonal entries stay exactly 0, and so do their standard errors. fit$se_omega is the lower triangle read column by column — (1,1), (2,1), (3,1), (2,2), (3,2), (3,3) — so the third and fifth entries are the held covariances, while the sixth is the standard error of the estimated ETA_KA variance:

fit_block$se_omega
#> [1] 0.012799626 0.005278263 0.000000000 0.004297678 0.000000000 0.160753224

Those two structural zeros are also the reason the fit prints an R warning that one or more diagonal elements of the correlation matrix are non-positive: ferx-r derives that matrix from fit$cov_matrix, which carries a zero variance for every held entry.

To estimate those covariances, declare them. A copy of the model with one full 3 × 3 block, and no separate diagonal omega:

full_block <- file.path(variability_dir, "warfarin_full_block.ferx")
full_lines <- sub("block_omega (ETA_CL, ETA_V) = [0.07, 0.02, 0.02]",
                  "block_omega (ETA_CL, ETA_V, ETA_KA) = [0.07, 0.02, 0.02, 0, 0, 0.40]",
                  readLines(block$model), fixed = TRUE)
writeLines(full_lines[!grepl("^  omega ETA_KA", full_lines)], full_block)
fit_full_block <- ferx_fit(full_block, block$data, verbose = FALSE)
fit_full_block$omega
#>              ETA_CL       ETA_V      ETA_KA
#> ETA_CL  0.028695094 0.001753564 -0.03165909
#> ETA_V   0.001753564 0.009629482  0.02068634
#> ETA_KA -0.031659091 0.020686341  0.35018945
data.frame(model = c("partial block", "full block"),
           n_parameters = c(fit_block$n_parameters, fit_full_block$n_parameters),
           ofv = c(fit_block$ofv, fit_full_block$ofv),
           aic = c(fit_block$aic, fit_full_block$aic),
           bic = c(fit_block$bic, fit_full_block$bic))
#>           model n_parameters       ofv       aic       bic
#> 1 partial block            8 -280.4858 -264.4858 -242.8819
#> 2    full block           10 -283.3167 -263.3167 -236.3119

The two extra covariances lower the OFV by 2.83 for two degrees of freedom, which a likelihood ratio test does not call significant (p = 0.243; Model development and selection). Both information criteria prefer the partial block. On these data ETA_KA is best left on its own, which is what the bundled model declares.

Additive random effects

warfarin_additive_eta gives the lag time an additive random effect, TLAG = TVTLAG + ETA_TLAG, so its variance is in squared hours:

additive <- ferx_example("warfarin_additive_eta")
invisible(ferx_model_get_section(additive$model, "individual_parameters"))
#> # [individual_parameters]
#>   CL   = TVCL  * exp(ETA_CL)
#>   V    = TVV   * exp(ETA_V)
#>   KA   = TVKA
#>   TLAG = TVTLAG + ETA_TLAG
fit_additive <- ferx_fit(additive$model, additive$data, verbose = FALSE)
fit_additive$estimates[c("TVTLAG", "ETA_TLAG"), c("estimate", "rse_pct")]
#>                 estimate     rse_pct
#> TVTLAG   0.0000001979303 83184.64655
#> ETA_TLAG 0.4764398317029    53.73498
summary(fit_additive$individual_estimates$TLAG)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -1.84258 -0.47863 -0.20734 -0.33790  0.07737  0.15992

These data carry no lag time: the typical value goes to its lower bound near 0. The additive random effect then produces negative individual lag times. A log-normal random effect keeps a positive parameter positive.

Residual error models

The residual error of the warfarin model can be written in several ways. The loop below replaces the sigma declarations and the [error_model] line of a copy of the model file, and fits each version:

Table 16.2: Residual error models for the warfarin model (FOCE).

warfarin_lines <- readLines(w$model)

# Replace the sigma lines and the [error_model] statement of the warfarin model.
set_error_model <- function(sigmas, error, extra = NULL) {
  lines <- warfarin_lines[!grepl("^  sigma ", warfarin_lines)]
  lines <- append(lines, c(extra, sigmas), after = grep("^  omega ETA_KA", lines))
  at <- grep("^\\[error_model\\]", lines)
  append(lines[-(at + 1)], error, after = at)
}
error_models <- list(
  proportional = warfarin_lines,
  additive = set_error_model("  sigma ADD_ERR ~ 1.0 (sd)", "  DV ~ additive(ADD_ERR)"),
  combined = set_error_model(c("  sigma PROP_ERR ~ 0.1 (sd)", "  sigma ADD_ERR ~ 0.5 (sd)"),
                             "  DV ~ combined(PROP_ERR, ADD_ERR)"),
  power = set_error_model("  sigma PROP_ERR ~ 0.1 (sd)", "  DV ~ power(PROP_ERR, RUV_POW)",
                          extra = "  theta RUV_POW(1.0, 0.01, 10.0)"),
  correlated_combined = set_error_model("  block_sigma (PROP_ERR, ADD_ERR) = [0.01, 0.0, 0.25]",
                                        "  DV ~ combined(PROP_ERR, ADD_ERR)")
)
error_fits <- lapply(names(error_models), function(name) {
  path <- file.path(variability_dir, paste0("warfarin_", name, ".ferx"))
  writeLines(error_models[[name]], path)
  ferx_fit(path, w$data, verbose = FALSE)
})
names(error_fits) <- names(error_models)
do.call(rbind, lapply(names(error_fits), function(name) {
  f <- error_fits[[name]]
  lines <- error_models[[name]]
  data.frame(model = name, error_model = trimws(lines[grep("^\\[error_model\\]", lines) + 1]),
             parameters = f$n_parameters, ofv = f$ofv, aic = f$aic,
             warnings = paste(ferx_get_warnings(f, as_df = TRUE)$category, collapse = ", "))
}))
#>                 model                      error_model parameters       ofv
#> 1        proportional      DV ~ proportional(PROP_ERR)          7 -280.3640
#> 2            additive           DV ~ additive(ADD_ERR)          7 -234.2689
#> 3            combined DV ~ combined(PROP_ERR, ADD_ERR)          8 -280.3634
#> 4               power    DV ~ power(PROP_ERR, RUV_POW)          8 -280.6655
#> 5 correlated_combined DV ~ combined(PROP_ERR, ADD_ERR)          9 -280.7510
#>         aic                                       warnings
#> 1 -266.3640                             dw_autocorrelation
#> 2 -220.2689                             dw_autocorrelation
#> 3 -264.3634 parameter_at_runaway_guard, dw_autocorrelation
#> 4 -264.6655                             dw_autocorrelation
#> 5 -262.7510                             dw_autocorrelation

The additive model is clearly worse. The other forms reach nearly the same OFV as the proportional model, so their extra parameters are not supported. In the combined model the additive component collapses:

error_fits$combined$estimates[c("PROP_ERR", "ADD_ERR"), c("estimate", "rse_pct")]
#>              estimate     rse_pct
#> PROP_ERR 0.0107483303    8.762397
#> ADD_ERR  0.0003354626 2925.132115
error_fits$power$estimates["RUV_POW", c("estimate", "rse_pct")]
#>         estimate  rse_pct
#> RUV_POW 1.084875 14.76553

All sigmas are on the standard-deviation scale: a proportional sigma of 0.1 is a CV of 10%. The additive start of a combined model should be on the scale of the data. ferx warns before fitting when it is below 1% of the median observation, as with 0.001 here (the warning has the category general):

small_start <- file.path(variability_dir, "warfarin_combined_small_start.ferx")
writeLines(set_error_model(c("  sigma PROP_ERR ~ 0.1 (sd)", "  sigma ADD_ERR ~ 0.001 (sd)"),
                           "  DV ~ combined(PROP_ERR, ADD_ERR)"), small_start)
fit_small_start <- ferx_fit(small_start, w$data, verbose = FALSE)
substr(ferx_get_warnings(fit_small_start, as_df = TRUE)$message[1], 1, 160)
#> [1] "Combined error model on all observations: additive SD initial estimate 'ADD_ERR' = 0.0010 is negligible (< 1%) relative to the data scale (median |DV| = 7.0715)"

In a single DV ~ line, the sigmas are taken in the order they are declared in [parameters], and the names in the call must match that order:

wrong_order <- file.path(variability_dir, "warfarin_sigma_order.ferx")
writeLines(set_error_model(c("  sigma ADD_ERR ~ 0.5 (sd)", "  sigma PROP_ERR ~ 0.1 (sd)"),
                           "  DV ~ combined(PROP_ERR, ADD_ERR)"), wrong_order)
try(ferx_fit(wrong_order, w$data, verbose = FALSE))
#> Error in ferx_rust_fit(model_path = normalizePath(model), data_path = normalizePath(data),  : 
#>   Error parsing model: [error_model] argument 1 names sigma 'PROP_ERR', but a single-endpoint [error_model] has its sigmas consumed positionally from the [parameters] declaration order, which supplies 'ADD_ERR' in that position. Reorder the sigma declarations to match the [error_model] argument order — if a block_sigma supplies them, permute its lower triangle to match the new name order. Reordering the arguments instead also parses, but for `combined` it swaps which sigma is the proportional component, so it changes the model rather than fixing the spelling. [E_SIGMA_ORDER_MISMATCH]

Log-transform-both-sides

warfarin_ltbs fits additive error to the logarithms of the observations and predictions:

ltbs <- ferx_example("warfarin_ltbs")
invisible(ferx_model_get_section(ltbs$model, "error_model"))
#> # [error_model]
#>   log(DV) ~ additive(SIGMA_LOG)
fit_ltbs <- ferx_fit(ltbs$model, ltbs$data, verbose = FALSE)
rbind(proportional = fit_w$theta, ltbs = fit_ltbs$theta)
#>                  TVCL      TVV      TVKA
#> proportional 0.132969 7.730700 0.7252070
#> ltbs         0.132697 7.738255 0.8109482
c(proportional_sigma = fit_w$sigma, ltbs_sigma = fit_ltbs$sigma)
#> proportional_sigma         ltbs_sigma 
#>         0.01074853         0.01056404

The structural estimates are close to the proportional model’s, and the log-scale sigma is close to the proportional CV. The OFV is on the log scale, so it is not comparable with the OFVs above. Predictions, residuals and simulated values are also on the log scale. Back-transform them with exp():

head(select(fit_ltbs$sdtab, ID, TIME, DV, PRED, IPRED, IWRES), 3)
#>   ID TIME      DV     PRED    IPRED       IWRES
#> 1  1  0.5  5.3653 1.455827 1.677672  0.21587411
#> 2  1  1.0  8.2578 1.961513 2.114356 -0.30265802
#> 3  1  2.0 10.7063 2.317392 2.370938 -0.01003281
head(mutate(select(fit_ltbs$sdtab, ID, TIME, DV, IPRED), IPRED_natural = exp(IPRED)), 3)
#>   ID TIME      DV    IPRED IPRED_natural
#> 1  1  0.5  5.3653 1.677672      5.353078
#> 2  1  1.0  8.2578 2.114356      8.284245
#> 3  1  2.0 10.7063 2.370938     10.707435

Under LTBS only additive error is allowed:

ltbs_proportional <- file.path(variability_dir, "warfarin_ltbs_proportional.ferx")
writeLines(sub("log(DV) ~ additive(SIGMA_LOG)", "log(DV) ~ proportional(SIGMA_LOG)",
               readLines(ltbs$model), fixed = TRUE), ltbs_proportional)
try(ferx_fit(ltbs_proportional, ltbs$data, verbose = FALSE))
#> Error in ferx_rust_fit(model_path = normalizePath(model), data_path = normalizePath(data),  : 
#>   Error parsing model: [error_model] log-transform-both-sides supports only additive error on the log scale; use `log(DV) ~ additive(...)` or `DV ~ log_additive(...)`: log(DV) ~ proportional(SIGMA_LOG) [E_PARSE]

Use DV ~ log_additive(...) instead when the DV column already holds logarithms.

Random effects on the residual error

iiv_on_ruv names a dedicated omega that scales the residual standard deviation of each subject by exp(ETA_RUV). It needs an estimation method with interaction, such as FOCEI:

ruv_model <- file.path(variability_dir, "warfarin_iiv_on_ruv.ferx")
writeLines(set_error_model(c("  omega ETA_RUV ~ 0.05", "  sigma PROP_ERR ~ 0.1 (sd)"),
                           c("  DV ~ proportional(PROP_ERR)", "  iiv_on_ruv = ETA_RUV")), ruv_model)
invisible(ferx_model_get_section(ruv_model, "error_model"))
#> # [error_model]
#>   DV ~ proportional(PROP_ERR)
#>   iiv_on_ruv = ETA_RUV
# the model file uses FOCE, without interaction
sub("\\)\\..*", ")", tryCatch(ferx_fit(ruv_model, w$data, verbose = FALSE), error = conditionMessage))  # first sentence
#> [1] "Fit error: IIV on residual error (iiv_on_ruv) requires an interaction or Monte-Carlo method: use method = focei, imp, impmap, or saem (got Foce with interaction = false)"
fit_ruv <- ferx_fit(ruv_model, w$data, method = "focei", verbose = FALSE)
#> Warning: Model file [fit_options] sets `method = foce` but ferx_fit() argument
#> overrides it with `focei`. The call-time value will be used.
fit_ruv$estimates[c("ETA_RUV", "PROP_ERR"), c("estimate", "rse_pct")]
#>                estimate     rse_pct
#> ETA_RUV  0.000006144212 4442.608422
#> PROP_ERR 0.010564986655    7.903861
ferx_get_warnings(fit_ruv, as_df = TRUE)[, c("severity", "category")]
#>   severity                   category
#> 1  warning              eta_shrinkage
#> 2  warning parameter_at_runaway_guard
#> 3  warning         dw_autocorrelation

With 10 subjects the variance collapses to the lower guard: these data do not show subjects with different residual error. The ferx-core error model page suggests Monte-Carlo methods (imp, impmap, saem) when the estimate collapses.

The error model can do more than shown here. The ferx-core error model page documents:

error models selected by a covariate with if/else;
magnitudes that depend on TIME or covariates (experimental);
residual weighting with weight = for meta-analysis data (experimental).

No bundled example uses them.

Inter-occasion variability

warfarin_iov adds a kappa to clearance: each subject gets one KAPPA_CL per occasion, on top of its ETA_CL. iov_column = OCC tells ferx where the occasions are:

invisible(ferx_model_get_section(iov$model, "parameters"))
#> # [parameters]
#>   theta TVCL(0.134, 0.001, 10.0)
#>   theta TVV(8.1, 0.1, 500.0)
#>   theta TVKA(1.0, 0.01, 50.0)
#>   omega ETA_CL ~ 0.07
#>   omega ETA_V  ~ 0.02
#>   omega ETA_KA ~ 0.40
#>   kappa KAPPA_CL ~ 0.04
#>   sigma PROP_ERR ~ 0.01 (sd)
invisible(ferx_model_get_section(iov$model, "individual_parameters"))
#> # [individual_parameters]
#>   CL = TVCL * exp(ETA_CL + KAPPA_CL)
#>   V  = TVV  * exp(ETA_V)
#>   KA = TVKA * exp(ETA_KA)
invisible(ferx_model_get_section(iov$model, "fit_options"))
#> # [fit_options]
#>   method     = foce
#>   iov_column = OCC
#>   covariance = false

The bundled model uses FOCE, and warfarin_iov_saem is the same model estimated with SAEM. The table compares both with FOCEI and with a SAEM → FOCEI chain:

Table 16.3: The warfarin IOV model with different estimation methods.

iov_saem <- ferx_example("warfarin_iov_saem")
iov_fits <- list(
  foce = ferx_fit(iov$model, iov$data, verbose = FALSE),
  focei = ferx_fit(iov$model, iov$data, method = "focei", verbose = FALSE),
  saem = ferx_fit(iov_saem$model, iov_saem$data, verbose = FALSE),
  saem_then_focei = ferx_fit(iov_saem$model, iov_saem$data, method = c("saem", "focei"), verbose = FALSE)
)
#> Warning: Model file [fit_options] sets `method = foce` but ferx_fit() argument
#> overrides it with `focei`. The call-time value will be used.
#> Warning: Model file [fit_options] sets `method = saem` but ferx_fit() argument
#> overrides it with `saem, focei`. The call-time value will be used.
do.call(rbind, lapply(names(iov_fits), function(name) {
  f <- iov_fits[[name]]
  data.frame(run = name, method = f$method, ofv = f$ofv, TVCL = f$theta[["TVCL"]], TVKA = f$theta[["TVKA"]],
             omega_CL = f$omega["ETA_CL", "ETA_CL"], kappa_CL = f$omega_iov[1, 1], sigma = f$sigma[[1]])
}))
#>               run method      ofv      TVCL     TVKA   omega_CL   kappa_CL
#> 1            foce   FOCE 203.5575 0.3226566 2.476479 0.47204747 0.03903501
#> 2           focei  FOCEI 308.8305 0.1727767 1.178559 0.03993581 0.03570809
#> 3            saem   SAEM 310.6889 0.1709340 1.200470 0.01185184 0.05152043
#> 4 saem_then_focei  FOCEI 308.8305 0.1727767 1.178559 0.03993581 0.03570809
#>       sigma
#> 1 0.2006693
#> 2 0.1881154
#> 3 0.1906719
#> 4 0.1881154

FOCEI, SAEM and the chain agree. FOCE does not: with a proportional error model, the method without interaction gives different estimates, including a clearance twice as high and a much larger between-subject variance. Its OFV is a different objective, so the numbers cannot be compared. Use an interaction method for models with proportional error (Estimation methods and controlling the fit).

The IOV results of the FOCEI fit:

fit_iov <- iov_fits$focei
fit_iov$kappa_names
#> [1] "KAPPA_CL"
fit_iov$omega_iov
#>            KAPPA_CL
#> KAPPA_CL 0.03570809
fit_iov$shrinkage_kappa
#>  KAPPA_CL 
#> 0.2275971
fit_iov$shrinkage_kappa_by_occ
#>   occ  KAPPA_CL
#> 1   1 0.1357611
#> 2   2 0.3319405
head(fit_iov$ebe_kappas, 4)
#>   ID OCC    KAPPA_CL
#> 1  1   1 -0.12885568
#> 2  1   2 -0.08893739
#> 3  2   1 -0.26486637
#> 4  2   2  0.08047400
head(select(fit_iov$sdtab, ID, TIME, OCC, DV, IPRED), 3)
#>   ID TIME OCC     DV    IPRED
#> 1  1  0.5   1 4.3571 4.466412
#> 2  1  1.0   1 6.9057 7.007888
#> 3  1  2.0   1 9.2120 9.239651

fit$ebe_kappas has one row per subject and occasion. sdtab gains the OCC column. Compare with the same model without the kappa, fitted to the same data:

fit_no_iov <- ferx_fit(w$model, iov$data, method = "focei", verbose = FALSE)
#> Warning: Model file [fit_options] sets `method = foce` but ferx_fit() argument
#> overrides it with `focei`. The call-time value will be used.
c(without_iov = fit_no_iov$ofv, with_iov = fit_iov$ofv)
#> without_iov    with_iov 
#>    366.5848    308.8305

The kappa lowers the OFV by 57.8 units for one extra parameter.

Simulation draws a new kappa for each occasion, so a VPC of an IOV model includes the occasion-to-occasion spread:

sim_iov <- ferx_simulate(iov$model, iov$data, n_sim = 2, seed = 1, fit = fit_iov)
head(sim_iov, 3)
#>   DRAW SIM ID TIME CMT    IPRED    DV_SIM OBSERVED
#> 1    1   1  1  0.5   1 4.678298  6.366642       NA
#> 2    1   1  1  1.0   1 7.250801  7.041806       NA
#> 3    1   1  1  2.0   1 9.370390 10.136965       NA

Occasions without an `OCC` column

iov_occasion derives the occasions from the dosing records or from time windows. dose starts an occasion at every dose time. time(120) makes one occasion before 120 hours and one from 120 hours on (interior boundaries only, no leading 0). The occasion column can also be given at call time with settings. All reproduce the fit with iov_column = OCC:

iov_lines <- readLines(iov$model)
occasion_variants <- list(
  "iov_occasion = dose" = sub("iov_column = OCC", "iov_occasion = dose", iov_lines, fixed = TRUE),
  "iov_occasion = time(120)" = sub("iov_column = OCC", "iov_occasion = time(120)", iov_lines, fixed = TRUE),
  "settings = list(iov_column = \"OCC\")" = iov_lines[!grepl("iov_column", iov_lines)]
)
sapply(names(occasion_variants), function(name) {
  path <- file.path(variability_dir, paste0(make.names(name), ".ferx"))
  writeLines(occasion_variants[[name]], path)
  settings <- if (grepl("^settings", name)) list(iov_column = "OCC") else NULL
  # suppress the R warning that reports overriding the model file's method
  suppressWarnings(ferx_fit(path, iov$data, method = "focei", settings = settings, verbose = FALSE))$ofv
})
#>                 iov_occasion = dose            iov_occasion = time(120) 
#>                            308.8305                            308.8305 
#> settings = list(iov_column = "OCC") 
#>                            308.8305

Without any occasion source, a model with kappas stops:

no_occasion_source <- file.path(variability_dir, paste0(make.names(names(occasion_variants)[3]), ".ferx"))
try(ferx_fit(no_occasion_source, iov$data, verbose = FALSE))
#> Error in ferx_rust_fit(model_path = normalizePath(model), data_path = normalizePath(data),  : 
#>   Fit error: Model declares kappa (IOV) parameters but no occasion labels were found. Either set `iov_column = "OCC"` (or the relevant column name) to read occasions from the dataset, or set `iov_occasion = dose` / `iov_occasion = time(...)` to derive them in the model — both in [fit_options] — so that per-occasion kappas can be estimated. [E_IOV_MISSING_OCC]

IOV on more parameters

A kappa on the volume as well, first as two independent kappas and then as a block_kappa:

Table 16.4: IOV on clearance and volume (FOCEI).

iov_v_lines <- sub("V  = TVV  * exp(ETA_V)", "V  = TVV  * exp(ETA_V + KAPPA_V)", iov_lines, fixed = TRUE)
iov_structures <- list(
  kappa_CL = iov_lines,
  kappa_CL_and_V = sub("kappa KAPPA_CL ~ 0.04", "kappa KAPPA_CL ~ 0.04\n  kappa KAPPA_V ~ 0.02", iov_v_lines, fixed = TRUE),
  block_kappa = sub("kappa KAPPA_CL ~ 0.04", "block_kappa (KAPPA_CL, KAPPA_V) = [0.04, 0.0, 0.02]", iov_v_lines, fixed = TRUE)
)
iov_structure_fits <- lapply(names(iov_structures), function(name) {
  path <- file.path(variability_dir, paste0("warfarin_iov_", name, ".ferx"))
  writeLines(iov_structures[[name]], path)
  # suppress the R warning that reports overriding the model file's method
  suppressWarnings(ferx_fit(path, iov$data, method = "focei", verbose = FALSE))
})
names(iov_structure_fits) <- names(iov_structures)
do.call(rbind, lapply(names(iov_structure_fits), function(name) {
  f <- iov_structure_fits[[name]]
  data.frame(model = name, parameters = f$n_parameters, ofv = f$ofv, omega_V = f$omega["ETA_V", "ETA_V"],
             sigma = f$sigma[[1]], warnings = paste(ferx_get_warnings(f, as_df = TRUE)$category, collapse = ", "))
}))
#>            model parameters      ofv        omega_V     sigma
#> 1       kappa_CL          8 308.8305 0.010777732526 0.1881154
#> 2 kappa_CL_and_V          9 206.2059 0.000006144212 0.1305205
#> 3    block_kappa         10 202.1861 0.000006144212 0.1304248
#>                                                        warnings
#> 1                             eta_shrinkage, dw_autocorrelation
#> 2 eta_shrinkage, parameter_at_runaway_guard, dw_autocorrelation
#> 3 eta_shrinkage, parameter_at_runaway_guard, dw_autocorrelation
iov_structure_fits$block_kappa$omega_iov
#>            KAPPA_CL    KAPPA_V
#> KAPPA_CL 0.04887553 0.02867844
#> KAPPA_V  0.02867844 0.04805688

Adding the kappa on the volume lowers the OFV a lot and the residual error drops. At the same time the between-subject variance of the volume collapses to the lower guard (parameter_at_runaway_guard). In these data the volume varies between occasions rather than between subjects, so the next model to try would drop ETA_V. The block adds a covariance for a small further decrease.

IOV with analytical absorption models

IOV works with every structural model. one_cpt_transit_iov puts a kappa on clearance in the analytical transit model of Absorption and bioavailability. ferx evaluates subjects with IOV through the equivalent ODE model. According to its model file, the data were simulated with TVCL 9, TVV 30, TVMTT 1, TVN 3 and a kappa variance of 0.04:

transit_iov <- ferx_example("one_cpt_transit_iov")
invisible(ferx_model_get_section(transit_iov$model, "individual_parameters"))
#> # [individual_parameters]
#>   CL  = TVCL * exp(KAPPA_CL)
#>   V   = TVV  * exp(ETA_V)
#>   MTT = TVMTT
#>   NTR = TVN
fit_transit_iov <- ferx_fit(transit_iov$model, transit_iov$data, verbose = FALSE)
fit_transit_iov$theta
#>       TVCL        TVV      TVMTT        TVN 
#>  9.9636094 25.8915920  0.9938112  2.9755746
fit_transit_iov$omega_iov
#>           KAPPA_CL
#> KAPPA_CL 0.0327846

Mu-referencing

With mu-referencing, ferx recognises parameters written as THETA * exp(ETA) (or THETA + ETA, or exp(log(THETA) + ETA)). It uses the mapping to centre the search for the individual random effects on the current population value. No model changes are needed. It is on by default, and mu_referencing = FALSE switches it off. For FOCEI on the two-compartment covariate model it makes no difference; for SAEM on warfarin it does:

Table 16.5: Fits with and without mu-referencing.

cov <- ferx_example("two_cpt_oral_cov")
saem <- ferx_example("warfarin_saem")
mu_runs <- expand.grid(model = c("two_cpt_oral_cov", "warfarin_saem"), mu_referencing = c(TRUE, FALSE),
                       stringsAsFactors = FALSE)
mu_runs$ofv <- mapply(function(model, mu) {
  ex <- if (model == "warfarin_saem") saem else cov
  # suppress the R warning that reports overriding the model file's covariance setting
  suppressWarnings(ferx_fit(ex$model, ex$data, mu_referencing = mu, covariance = FALSE, verbose = FALSE))$ofv
}, mu_runs$model, mu_runs$mu_referencing)
arrange(mu_runs, model)
#>              model mu_referencing        ofv
#> 1 two_cpt_oral_cov           TRUE -1199.3264
#> 2 two_cpt_oral_cov          FALSE -1199.3264
#> 3    warfarin_saem           TRUE  -285.9865
#> 4    warfarin_saem          FALSE  -277.8816

Mu-referencing is skipped for a parameter whose assignment sits inside an if block (individual parameters).

Parameter scaling

parameter_scaling rescales the internal parameters for the outer optimizer: auto (the default), none, abs or rescale2. scale_params = TRUE is the older switch for abs. On the two-compartment covariate model with the gradient-based nlopt_lbfgs optimizer, all choices reach the same optimum:

Table 16.6: Parameter scaling with the nlopt_lbfgs optimizer.

scaling_runs <- list(auto = list(), none = list(parameter_scaling = "none"), abs = list(parameter_scaling = "abs"),
                     rescale2 = list(parameter_scaling = "rescale2"))
scaling_table <- data.frame(
  parameter_scaling = names(scaling_runs),
  ofv = sapply(scaling_runs, function(s) {
    # suppress the R warning that reports overriding the model file's covariance setting
    suppressWarnings(ferx_fit(cov$model, cov$data, covariance = FALSE, verbose = FALSE,
                              settings = c(list(optimizer = "nlopt_lbfgs"), s)))$ofv
  }), row.names = NULL)
fit_scale_params <- suppressWarnings(ferx_fit(cov$model, cov$data, covariance = FALSE, scale_params = TRUE,
                                              verbose = FALSE, settings = list(optimizer = "nlopt_lbfgs")))
rbind(scaling_table, data.frame(parameter_scaling = "scale_params = TRUE", ofv = fit_scale_params$ofv))
#>     parameter_scaling       ofv
#> 1                auto -1199.326
#> 2                none -1199.326
#> 3                 abs -1199.326
#> 4            rescale2 -1199.326
#> 5 scale_params = TRUE -1199.326

The OFV at a given point does not depend on scaling, but the optimizer’s path can. Try it when a gradient-based fit stops early from a poor start.

Mixture models

A [mixture] block assigns each subject probabilistically to one of several subpopulations, each with its own typical values (selected with MIXNUM in [individual_parameters]). ferx-core marks it experimental, and no bundled ferx-r example uses it. See the ferx-core mixture models page.

Pitfalls

Proportional error needs interaction. FOCE without interaction gave clearly different IOV estimates above. Use FOCEI (or SAEM).
Extra variance parameters often collapse. An additive component, a residual random effect or a between-subject variance that goes to its lower guard is not supported by the data. Remove it rather than report it.
An additive random effect can make a positive parameter negative. Use exp() for parameters such as clearance, volume and lag time.
Keep sigma declarations in the order of the error model call. A single DV ~ line takes them by position.
A block estimates only the covariances it declares. Held entries stay 0 with a standard error of 0, and the correlation matrix of the estimates is then reported as not meaningful. Widen the block if a covariance is worth estimating.
LTBS reports on the log scale, and its OFV cannot be compared with natural-scale error models.

Warnings you may see here

omega_structure (warning): a block_omega combines a log-normal and an additive random effect. The parameter-level correlation then falls back to the correlation of the random effects themselves. Here ETA_CL (log-normal) and ETA_TLAG (additive) share a block:

mixed_block <- file.path(variability_dir, "warfarin_mixed_block.ferx")
mixed_lines <- sub("^  omega ETA_CL   ~ 0.07.*", "  block_omega (ETA_CL, ETA_TLAG) = [0.07, 0.01, 0.04]",
                   readLines(additive$model))
writeLines(mixed_lines[!grepl("^  omega ETA_TLAG", mixed_lines)], mixed_block)
fit_mixed_block <- ferx_fit(mixed_block, additive$data, verbose = FALSE)
#> Warning in .ferx_compute_cor_matrix(result$cov_matrix): One or more diagonal
#> elements are non-positive; correlation matrix may not be meaningful.
warnings_mixed <- ferx_get_warnings(fit_mixed_block, as_df = TRUE)
warnings_mixed[warnings_mixed$category == "omega_structure", c("severity", "message")]
#>   severity
#> 2  warning
#>                                                                                                                      message
#> 2 omega_param_corr: ETA_TLAG × ETA_CL have mixed lognormal/additive parameterizations; falling back to eta-level correlation

The warning concerns the reported parameter correlation (fit$omega_param_corr), not the fit itself. See the ferx-core section on parameter-level correlation.

Other warnings in this chapter have their home elsewhere: parameter_at_runaway_guard and inflated_rse (Parameter uncertainty), eta_shrinkage and dw_autocorrelation (Diagnosing the model).

Summary

Between-subject variability: omega and block_omega, entered as log-normal or additive random effects.
Residual error: additive, proportional, combined, power, log-scale, correlated components, and a random effect on the residual error. Sigmas are standard deviations.
Inter-occasion variability: kappa and block_kappa, with occasions from a column or from iov_occasion. Results are in omega_iov, ebe_kappas and the kappa shrinkage fields.
Mu-referencing is on by default and can matter for SAEM. parameter_scaling changes the optimizer’s path, not the objective.

Next: Covariate modeling explains part of the variability with covariates.

Reference

R help: ?ferx_fit (mu_referencing, scale_params, settings), ?ferx_simulate
ferx-core: parameters, individual parameters, error model, IOV, mixture models, packed parameter space