Covariate modeling

Where you are

Model development and selection selected the covariates of the analysis thread with a formal search. This chapter covers the rest of covariate modeling:

declaring covariates;
screening them against the random effects;
writing covariate effects by hand, as conditions, or with the allometric-scaling tool;
the full random effects model (FREM);
covariates that change over time.

It starts again from the covariate-free base model.

Maturity: stable, experimental and alpha parts

The [covariates] block is stable in ferx-core. The [covariate_model] block is experimental. GAM screening and the search tools, including allometric scaling, are alpha. See feature maturity.

The data

two_cpt_oral_base and two_cpt_oral_cov share one dataset: 30 subjects with a single oral dose of 250, body weight WT and creatinine clearance CRCL. Both covariates are constant within each subject:

base <- ferx_example("two_cpt_oral_base")
cov <- ferx_example("two_cpt_oral_cov")
covariate_dir <- book_tempdir("covariates")
subjects <- read.csv(base$data, na.strings = ".") |> distinct(ID, WT, CRCL)
nrow(subjects)
#> [1] 30
summary(subjects[, c("WT", "CRCL")])
#>        WT             CRCL       
#>  Min.   :45.00   Min.   : 46.50  
#>  1st Qu.:60.73   1st Qu.: 73.47  
#>  Median :70.85   Median : 82.90  
#>  Mean   :70.43   Mean   : 86.11  
#>  3rd Qu.:80.55   3rd Qu.:102.17  
#>  Max.   :93.70   Max.   :150.00

Minimal runnable call: declaring covariates

The [covariates] block lists the dataset columns that are covariates, and whether each is continuous or categorical. Declaring a covariate does not put it in the model. two_cpt_oral_base declares both columns but uses neither:

invisible(ferx_model_get_section(base$model, "covariates"))
#> # [covariates]
#>   WT   continuous
#>   CRCL continuous
invisible(ferx_model_get_section(base$model, "individual_parameters"))
#> # [individual_parameters]
#>   CL = TVCL * exp(ETA_CL)
#>   V1 = TVV1 * exp(ETA_V1)
#>   Q  = TVQ  * exp(ETA_Q)
#>   V2 = TVV2 * exp(ETA_V2)
#>   KA = TVKA * exp(ETA_KA)
fit_base <- ferx_fit(base$model, base$data, verbose = FALSE)

Reading the result: the covariate table

A fit of a model with a [covariates] block carries the covariate table fit$covtab, one row per data record, and the declared types:

head(fit_base$covtab, 3)
#>   ID TIME EVID   WT CRCL
#> 1  1  0.0    1 70.6 73.7
#> 2  1  0.5    0 70.6 73.7
#> 3  1  1.0    0 70.6 73.7
fit_base$covariate_types
#>           WT         CRCL 
#> "continuous" "continuous"

The screening tools below read these. Categorical covariates must be coded as numbers in the data (SEX as 0 and 1). No bundled example has one; see the ferx-core covariates page.

Screening against the random effects

ferx_cov_screen() correlates each covariate, per subject, with each parameter that has a random effect. It uses both the individual parameter and its eta, and reports the pairs whose association reaches threshold:

screen <- ferx_cov_screen(fit_base)
#>   parameter covariate       type   ebe   eta
#> 1        CL      CRCL continuous 0.478 0.469
#> 2        CL        WT continuous 0.360 0.382
#> 3        V1        WT continuous 0.326 0.305
ferx_cov_screen(fit_base, threshold = 0.5)
#> No covariate correlations with |r| >= 0.50.

Clearance is associated with CRCL and WT, and the central volume with WT. ferx_gam_screen() goes one step further. For each eta and covariate it compares a regression of the eta on the covariate (linear, or a natural spline with spline_df degrees of freedom) with an intercept-only model, and ranks the pairs by the AIC improvement:

gam <- ferx_gam_screen(fit_base)
#>    eta_name covariate  delta_aic    best_form        aic  aic_null   r_squared
#> 1    ETA_CL      CRCL  5.4348439       Linear  -59.25279 -53.81795 0.219505924
#> 2    ETA_CL        WT  2.8451085 Spline(df=2)  -56.66306 -53.81795 0.204011295
#> 3    ETA_V1        WT  0.9340646       Linear  -64.80154 -63.86747 0.093171701
#> 4    ETA_V1      CRCL  0.2941567 Spline(df=2)  -64.16163 -63.86747 0.133366019
#> 5     ETA_Q        WT -0.9187237 Spline(df=2)  -87.23562 -88.15434 0.097610660
#> 6     ETA_Q      CRCL -1.9393902       Linear  -86.21495 -88.15434 0.002018287
#> 7    ETA_V2        WT  2.7840622 Spline(df=2) -100.00225 -97.21819 0.202389908
#> 8    ETA_V2      CRCL -1.7757908       Linear  -95.44240 -97.21819 0.007445782
#> 9    ETA_KA        WT -1.1772312 Spline(df=3)  -51.53800 -52.71523 0.148502710
#> 10   ETA_KA      CRCL -1.7289480 Spline(df=2)  -50.98628 -52.71523 0.072907313
#>       shrinkage
#> 1  -0.001744168
#> 2  -0.001744168
#> 3   0.036475626
#> 4   0.036475626
#> 5   0.029265476
#> 6   0.029265476
#> 7   0.024072075
#> 8   0.024072075
#> 9   0.041527717
#> 10  0.041527717

The strongest pair is again ETA_CL with CRCL, as a linear effect. The etas and covariates arguments restrict the screen, include_linear = FALSE keeps only the spline forms, and shrinkage_warn sets the shrinkage above which a warning is given:

ferx_gam_screen(fit_base, etas = "ETA_CL", covariates = "CRCL", spline_df = 2L,
                include_linear = FALSE, shrinkage_warn = 0.3)
#>   eta_name covariate delta_aic    best_form       aic  aic_null r_squared
#> 1   ETA_CL      CRCL  3.510319 Spline(df=2) -57.32827 -53.81795  0.221467
#>      shrinkage
#> 1 -0.001744168

Screening uses the individual random effects, so it is only informative when shrinkage is low. It is here (Diagnosing the model):

setNames(round(100 * fit_base$shrinkage_eta, 1), fit_base$eta_names)
#> ETA_CL ETA_V1  ETA_Q ETA_V2 ETA_KA 
#>   -0.2    3.6    2.9    2.4    4.2

A screen suggests candidates. The decision needs fits: a covariate search (Model development and selection) or the comparison of hand-written models below.

Options that matter

Function	Argument	Meaning
`ferx_cov_screen()`	`fit`	Fit of a model with a `[covariates]` block
	`threshold`	Minimum absolute association to report (default 0.2)
`ferx_gam_screen()`	`fit`	As above
	`etas`, `covariates`	Restrict the pairs (default: all)
	`spline_df`	Spline degrees of freedom to try (default `c(2L, 3L)`)
	`include_linear`	Also try the linear form (default `TRUE`)
	`shrinkage_warn`	Warn for etas with shrinkage above this fraction (default 0.3)
`ferx_allometry()`	`model`, `data`	The model to scale and its dataset
	`config`	A `.ferxsearch` file with an `ALLOMETRY(WT, 70)` statement, instead of the arguments
	`covariate`, `reference`	Size covariate and reference value (default `WT` and 70)
	`parameters`, `exponents`	Parameters to scale and their exponents (default: all clearances with 0.75, all volumes with 1)
	`estimate`, `lower`, `upper`	Estimate the exponents within bounds instead of fixing them
	`fit`	Fit the base and scaled models (default `TRUE`), or only build the scaled model
	`threads`, `retries`, `directory`	Worker threads, restarts per fit, and a directory for the run journal
`ferx_model_to_frem()`	`model`, `data`	The base model (with a `[covariates]` block) and its dataset
	`covariates`	A subset of the declared covariates (default: all)
	`output_dir`, `output_model`, `output_data`	Where the FREM model and data are written
	`fit`	A fit of the base model, whose estimates become the initial values

Variants

Covariate effects in `[individual_parameters]`

The usual way to add a covariate effect is to write it into the individual parameter. two_cpt_oral_cov scales clearance and central volume by weight with an estimated exponent, and clearance by creatinine clearance:

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_cov <- ferx_fit(cov$model, cov$data, verbose = FALSE)
fit_cov$estimates[c("THETA_WT", "THETA_CRCL"), c("estimate", "rse_pct")]
#>             estimate  rse_pct
#> THETA_WT   0.6534563 36.13571
#> THETA_CRCL 0.5646128 39.23061

Conditional effects

An if/else block in [individual_parameters] lets a covariate choose the expression. warfarin_if uses the same dataset. Subjects above 70 kg get allometric clearance and lighter subjects the typical clearance, with volume proportional to weight:

cond <- ferx_example("warfarin_if")
invisible(ferx_model_get_section(cond$model, "individual_parameters"))
#> # [individual_parameters]
#>   if (WT > 70) {
#>     CL = TVCL * (WT / 70)^0.75 * exp(ETA_CL)
#>   } else {
#>     CL = TVCL * exp(ETA_CL)
#>   }
#>   V1 = TVV1 * (WT / 70) * exp(ETA_V1)
#>   Q  = TVQ  * exp(ETA_Q)
#>   V2 = TVV2 * exp(ETA_V2)
#>   KA = TVKA * exp(ETA_KA)
fit_if <- ferx_fit(cond$model, cond$data, verbose = FALSE)
ferx_get_warnings(fit_if, as_df = TRUE)[, c("severity", "category", "message")]
#>   severity           category
#> 1     info     mu_referencing
#> 2  warning dw_autocorrelation
#>                                                                                                                                                                           message
#> 1 Mu-referencing disabled for conditional parameter(s): CL. Assign TV* unconditionally and apply the if-block to the individual parameter expression to re-enable mu-referencing.
#> 2                                                     Negative IWRES autocorrelation detected (Durbin-Watson = 2.91). Possible over-parameterization or misspecified error model.

A parameter assigned inside an if block loses mu-referencing (Variability: random effects, residual error and IOV), and the fit says so with an info note. On the same data, the three models compare as follows:

Table 17.1: Covariate models on the two_cpt_oral_cov dataset.

data.frame(model = c("two_cpt_oral_base", "warfarin_if", "two_cpt_oral_cov"),
           covariate_effects = c("none", "WT (conditional, fixed exponents)", "WT and CRCL (estimated exponents)"),
           parameters = c(fit_base$n_parameters, fit_if$n_parameters, fit_cov$n_parameters),
           ofv = c(fit_base$ofv, fit_if$ofv, fit_cov$ofv))
#>               model                 covariate_effects parameters       ofv
#> 1 two_cpt_oral_base                              none         11 -1185.403
#> 2       warfarin_if WT (conditional, fixed exponents)         11 -1184.600
#> 3  two_cpt_oral_cov WT and CRCL (estimated exponents)         13 -1199.326

Allometric scaling

ferx_allometry() adds body-size scaling to a model: (WT/70)^0.75 on every clearance and (WT/70)^1 on every volume of the pk line. With fit = FALSE it only builds the scaled model. The scaling is written as a [covariate_model] block:

scaled <- ferx_allometry(base$model, base$data, fit = FALSE)
scaled$scalings
#>   parameter exponent fixed theta
#> 1        CL     0.75  TRUE  <NA>
#> 2         Q     0.75  TRUE  <NA>
#> 3        V1     1.00  TRUE  <NA>
#> 4        V2     1.00  TRUE  <NA>
invisible(ferx_model_get_section(scaled$model_path, "covariate_model"))
#> # [covariate_model]
#>   CL ~ WT power(center = 70, fix = 0.75)
#>   Q ~ WT power(center = 70, fix = 0.75)
#>   V1 ~ WT power(center = 70, fix = 1.0)
#>   V2 ~ WT power(center = 70, fix = 1.0)

Each [covariate_model] line states a relation (PARAM ~ COV form(...)), and ferx rewrites it into the expression you would otherwise write by hand. The ferx-core covariate model page lists the forms (linear, power, exponential, hockey, categorical, …). No bundled ferx-r example writes the block by hand.

With fit = TRUE, the default, the base and the scaled model are both fitted. print() on the ferx_allometry result shows the relations and both fits (its digits argument sets the precision):

allometry <- ferx_allometry(base$model, base$data)
#> Warning in .ferx_compute_cor_matrix(result$cov_matrix): One or more diagonal
#> elements are non-positive; correlation matrix may not be meaningful.
allometry
#> ferx allometric scaling
#>   Data:  /home/runner/work/_temp/Library/ferx/examples/data/two_cpt_oral_cov.csv
#>   Model: /tmp/Rtmp5Kzcce/ferx-allometry-4cdf2ae5e4b8.ferx
#> 
#> Scaling on WT (reference 70):
#>   CL ~ WT power(center = 70, fix = 0.75)
#>   Q ~ WT power(center = 70, fix = 0.75)
#>   V1 ~ WT power(center = 70, fix = 1)
#>   V2 ~ WT power(center = 70, fix = 1)
#> 
#> Fits:
#>   model   ofv converged passed failures
#>    base -1185      TRUE   TRUE     <NA>
#>  scaled -1173      TRUE   TRUE     <NA>
#> 
#>   dOFV (base - scaled): -12.74
allometry$dofv
#> [1] -12.73983

Fixed allometric exponents make this model worse. With estimate = TRUE the exponents are estimated between lower and upper, and parameters limits the scaling to clearance and central volume. directory keeps a journal of the run:

allometry_dir <- file.path(covariate_dir, "allometry")
estimated <- ferx_allometry(base$model, base$data, parameters = c("CL", "V1"), exponents = c(0.75, 1),
                            estimate = TRUE, lower = 0, upper = 2, directory = allometry_dir)
estimated$comparison
#>    model       ofv converged passed failures
#> 1   base -1185.423      TRUE   TRUE     <NA>
#> 2 scaled -1193.593      TRUE   TRUE     <NA>
estimated$fit$estimates[c("THETA_CL_WT", "THETA_V1_WT"), c("estimate", "rse_pct")]
#>              estimate  rse_pct
#> THETA_CL_WT 0.7496625 53.00103
#> THETA_V1_WT 0.5567606 60.76338
list.files(allometry_dir)
#> [1] "allometric.ferx"      "candidates.csv"       "fits"                
#> [4] "search_journal.jsonl" "search_run.json"

With estimated exponents the scaled model improves on the base by 8.2 OFV units. Both exponents are poorly determined by these 30 subjects.

The same scaling can come from a .ferxsearch file with an ALLOMETRY statement. The bundled two_cpt_oral_base.ferxsearch holds a covariate search, so a copy with the statement is written here:

config_lines <- readLines(base$search)
config_lines <- sub("^base = .*", sprintf('base = "%s"', base$model), config_lines)
config_lines <- sub("^data = .*", sprintf('data = "%s"', base$data), config_lines)
config_lines <- sub("^mfl = .*", 'mfl = "ALLOMETRY(WT, 70)"', config_lines)
config_lines <- config_lines[!grepl("^\\[covsearch\\]|^algorithm|^p_forward|^p_backward", config_lines)]
allometry_config <- file.path(covariate_dir, "allometry.ferxsearch")
writeLines(config_lines, allometry_config)
ferx_allometry(config = allometry_config, fit = FALSE)$scalings
#>   parameter exponent fixed theta
#> 1        CL     0.75  TRUE  <NA>
#> 2         Q     0.75  TRUE  <NA>
#> 3        V1     1.00  TRUE  <NA>
#> 4        V2     1.00  TRUE  <NA>

threads and retries control the fits as in the other search tools (Model development and selection).

Full random effects modeling (FREM)

FREM treats each covariate as an additional observation. The covariate gets its own random effect in one large omega block, and the covariances between the parameter and covariate random effects describe the covariate relations, without a stepwise search. ferx_model_to_frem() builds the FREM model and dataset from a base model with a [covariates] block.

By default the generated files go next to the model file, which for a ferx_example() model is the installed package. output_dir writes them to a temporary directory instead:

frem <- ferx_model_to_frem(base$model, base$data, fit = fit_base,
                           output_dir = file.path(covariate_dir, "frem"))
frem
#> ferx_model
#>   Model: /tmp/Rtmp5Kzcce/covariates/frem/two_cpt_oral_base_frem.ferx
#>   Data:  /tmp/Rtmp5Kzcce/covariates/frem/two_cpt_oral_base_frem_data.csv
#>   ---
#>   Structural:  2-cpt oral  (TVCL, TVV1, TVQ, TVV2, TVKA, TV_WT, TV_CRCL)
#>   IIV:         ETA_CL, ETA_V1, ETA_Q, ETA_V2, ETA_KA, ETA_WT_FREM, ETA_CRCL_FREM
#>   IOV:         none
#>   Residual:    proportional

The dataset gains one row per subject and covariate. Their FREMTYPE codes the covariate (100 for WT, 200 for CRCL), and DV holds the covariate value:

frem_data <- read.csv(frem$data, na.strings = ".")
count(frem_data, FREMTYPE)
#>   FREMTYPE   n
#> 1        0 330
#> 2      100  30
#> 3      200  30
head(filter(frem_data, FREMTYPE > 0), 4)
#>   ID TIME   DV EVID AMT CMT RATE MDV II SS CENS FREMTYPE   WT CRCL
#> 1  1  0.5 70.6    0   0   1    0   0  0  0    0      100 70.6 73.7
#> 2  1  0.5 73.7    0   0   1    0   0  0  0    0      200 70.6 73.7
#> 3  2  0.5 80.9    0   0   1    0   0  0  0    0      100 80.9 80.7
#> 4  2  0.5 80.7    0   0   1    0   0  0  0    0      200 80.9 80.7

The model gains a fixed typical value and a random effect for each covariate, one omega block, and a small fixed covariate sigma. fit = fit_base put the base estimates in as initial values. The [fit_options] entries frem_predictions and frem_sigma tell the engine which rows belong to which covariate and which sigma they use:

invisible(ferx_model_get_section(frem, "parameters"))
#> # [parameters]
#>   theta TVCL(4.464757578550209, 0.1, 100.0)
#>   theta TVV1(48.24224473329911, 1.0, 500.0)
#>   theta TVQ(9.584245594278087, 0.1, 100.0)
#>   theta TVV2(93.78773933299064, 1.0, 500.0)
#>   theta TVKA(1.2474671134949786, 0.01, 10.0)
#>   theta TV_WT(70.43333333333334, FIX)
#>   theta TV_CRCL(86.11000000000003, FIX)
#> 
#>   block_omega (ETA_CL, ETA_V1, ETA_Q, ETA_V2, ETA_KA, ETA_WT_FREM, ETA_CRCL_FREM) = [
#>     1.551248e-1,
#>     0.000000e0, 1.199053e-1,
#>     0.000000e0, 0.000000e0, 5.257442e-2,
#>     0.000000e0, 0.000000e0, 0.000000e0, 3.845706e-2,
#>     0.000000e0, 0.000000e0, 0.000000e0, 0.000000e0, 1.758346e-1,
#>     4.956079e-2, 4.357292e-2, 2.885259e-2, 2.467661e-2, 5.276546e-2, 1.583416e2,
#>     9.377998e-2, 8.244961e-2, 5.459548e-2, 4.669361e-2, 9.984392e-2, 0.000000e0, 5.669423e2
#>   ]
#> 
#>   sigma PROP_ERR ~ 0.02 (sd)
#>   sigma EPSCOV ~ 1e-6 FIX
invisible(ferx_model_get_section(frem, "individual_parameters"))
#> # [individual_parameters]
#>   CL = TVCL * exp(ETA_CL)
#>   V1 = TVV1 * exp(ETA_V1)
#>   Q  = TVQ  * exp(ETA_Q)
#>   V2 = TVV2 * exp(ETA_V2)
#>   KA = TVKA * exp(ETA_KA)
#>   COV_WT = TV_WT + ETA_WT_FREM
#>   COV_CRCL = TV_CRCL + ETA_CRCL_FREM
invisible(ferx_model_get_section(frem, "fit_options"))
#> # [fit_options]
#>   method     = focei
#>   maxiter    = 500
#>   covariance = true
#>   gradient   = fd
#>   frem_predictions = TV_WT/ETA_WT_FREM:100, TV_CRCL/ETA_CRCL_FREM:200
#>   frem_sigma = EPSCOV
#> 
#> # Data: /tmp/Rtmp5Kzcce/covariates/frem/two_cpt_oral_base_frem_data.csv

A ferx_model object fits directly. The ferx-core FREM page recommends SAEM for FREM models, as the large omega block is hard for gradient-based methods. The covariance step is skipped here to save time:

# suppress the R warnings that report overriding the model file's method and covariance setting
fit_frem_focei <- suppressWarnings(ferx_fit(frem, covariance = FALSE, verbose = FALSE))
fit_frem_saem <- suppressWarnings(ferx_fit(frem, method = "saem", covariance = FALSE, verbose = FALSE))
frem_correlations <- function(f) round(cov2cor(f$omega)[c("ETA_CL", "ETA_V1"), c("ETA_WT_FREM", "ETA_CRCL_FREM")], 2)
list(focei = frem_correlations(fit_frem_focei), saem = frem_correlations(fit_frem_saem))
#> $focei
#>        ETA_WT_FREM ETA_CRCL_FREM
#> ETA_CL        0.01          0.01
#> ETA_V1        0.01          0.01
#> 
#> $saem
#>        ETA_WT_FREM ETA_CRCL_FREM
#> ETA_CL        0.39          0.47
#> ETA_V1        0.33          0.07

With SAEM the correlations show what the covariate model found: clearance goes with creatinine clearance and weight, and central volume with weight. FOCEI leaves them near zero on this model. The FREM fits carry several warnings that belong to the construction:

count(ferx_get_warnings(fit_frem_saem, as_df = TRUE), severity, category)
#>   severity           category  n
#> 1  warning dw_autocorrelation  1
#> 2  warning            general  3
#> 3  warning    omega_structure 10

Two of the general warnings report that the covariate sigma and the covariate expressions are not used by the error or structural model, which is intended here. The third reports that the model file’s maxiter does not apply to SAEM (Estimation methods and controlling the fit). The omega_structure warnings concern the reported parameter correlations (Variability: random effects, residual error and IOV).

covariates builds a FREM model for a subset of the declared covariates:

frem_crcl <- ferx_model_to_frem(base$model, base$data, covariates = "CRCL",
                                output_dir = file.path(covariate_dir, "frem_crcl"))
invisible(ferx_model_get_section(frem_crcl, "fit_options"))
#> # [fit_options]
#>   method     = focei
#>   maxiter    = 500
#>   covariance = true
#>   gradient   = fd
#>   frem_predictions = TV_CRCL/ETA_CRCL_FREM:100
#>   frem_sigma = EPSCOV
#> 
#> # Data: /tmp/Rtmp5Kzcce/covariates/frem_crcl/two_cpt_oral_base_frem_data.csv

FREM options:

book_settings_table(c("frem_rao_blackwell", "frem_predictions", "frem_sigma"))

Table 17.2

Setting	Default	Description
`frem_rao_blackwell`	`true`	FREM only: Rao-Blackwellise the covariate ETAs (integrate them analytically, importance-sample only the PK ETAs).
`frem_predictions`		See `?ferx_fit` and the ferx-core fit options page.
`frem_sigma`		See `?ferx_fit` and the ferx-core fit options page.

frem_rao_blackwell applies to importance sampling (method = "imp"). On this model an IMP fit takes several minutes, so it is not run here.

Time-varying covariates

A covariate may change within a subject. ferx carries the last value forward from each record, and the model uses the value in force at each time. In a copy of the data, subject 1’s weight rises by 50% from 12 hours. Only that subject’s predictions change, and only from 12 hours on:

tv_data <- read.csv(cov$data, na.strings = ".")
heavier <- tv_data$ID == 1 & tv_data$TIME >= 12
tv_data$WT[heavier] <- 1.5 * tv_data$WT[heavier]
tv_file <- file.path(covariate_dir, "two_cpt_oral_cov_tv.csv")
write.csv(tv_data, tv_file, row.names = FALSE, na = ".")
constant_pred <- ferx_predict(cov$model, cov$data, fit = fit_cov)
tv_pred <- ferx_predict(cov$model, tv_file, fit = fit_cov)
filter(constant_pred, ID == 1) |> select(TIME, PRED) |>
  mutate(PRED_time_varying = tv_pred$PRED[tv_pred$ID == 1])
#>    TIME      PRED PRED_time_varying
#> 1   0.5 2.2114606         2.2114606
#> 2   1.0 3.1116226         3.1116226
#> 3   2.0 3.2905693         3.2905693
#> 4   4.0 2.3602311         2.3602311
#> 5   6.0 1.6778386         1.6778386
#> 6   8.0 1.3044925         1.3044925
#> 7  12.0 0.9755278         0.8291789
#> 8  24.0 0.6747188         0.5698892
#> 9  36.0 0.5035083         0.4059945
#> 10 48.0 0.3762776         0.2895959
isTRUE(all.equal(constant_pred$PRED[constant_pred$ID != 1], tv_pred$PRED[tv_pred$ID != 1]))
#> [1] TRUE

A fit uses the same records. Its covariate table shows the change:

fit_tv <- ferx_fit(cov$model, tv_file, verbose = FALSE)
filter(fit_tv$covtab, ID == 1, TIME %in% c(8, 12))
#>   ID TIME EVID    WT CRCL
#> 1  1    8    0  70.6 73.7
#> 2  1   12    0 105.9 73.7

Pitfalls

Declaring is not using. A covariate in [covariates] has no effect until an individual parameter or a [covariate_model] relation uses it.
Screens are not tests. Correlations and GAM results on individual random effects suggest candidates, and they weaken as shrinkage grows. Confirm with fits.

Allometry does not see effects written by hand. It leaves a parameter alone only when a [covariate_model] relation already uses the size covariate. On two_cpt_oral_cov, which writes (WT / 70)^THETA_WT into [individual_parameters], it adds weight a second time:

double <- ferx_allometry(cov$model, cov$data, fit = FALSE)
double$notes
#> character(0)
grep("WT", readLines(double$model_path), value = TRUE)
#> [1] "# Two-compartment oral PK model with covariates (WT, CRCL)"              
#> [2] "  theta THETA_WT(0.75, 0.01, 5.0)"                                       
#> [3] "  CL = TVCL * (WT / 70)^THETA_WT * (CRCL / 100)^THETA_CRCL * exp(ETA_CL)"
#> [4] "  V1 = TVV1 * (WT / 70)^THETA_WT * exp(ETA_V1)"                          
#> [5] "  WT   continuous"                                                       
#> [6] "  CL ~ WT power(center = 70, fix = 0.75)"                                
#> [7] "  Q ~ WT power(center = 70, fix = 0.75)"                                 
#> [8] "  V1 ~ WT power(center = 70, fix = 1.0)"                                 
#> [9] "  V2 ~ WT power(center = 70, fix = 1.0)"

FREM needs a suitable method. Check the covariate–parameter correlations with SAEM: FOCEI found almost none on these data. Give output_dir a writable directory of your own, or the generated files land next to the model file.

Warnings you may see here

mu_referencing (info): a conditional assignment switched off mu-referencing for a parameter, as for warfarin_if above.
general (warning) on FREM models: a sigma declared but not used in [error_model], or individual parameters that are computed but not used. For the generated covariate sigma and covariate expressions this is expected.

Summary

Declare covariates in [covariates]; the fit returns covtab and covariate_types.
Screen with ferx_cov_screen() and ferx_gam_screen() while shrinkage is low, then confirm with fits.
Write effects in [individual_parameters], as if/else conditions, or as [covariate_model] relations. ferx_allometry() adds body-size scaling, fixed or estimated.
ferx_model_to_frem() builds a FREM model from the declared covariates; fit it with SAEM.
Covariates that change within a subject are carried forward from each record.

Next: Dosing regimens and exposure metrics covers dosing records and regimens.

Reference

R help: ?ferx_cov_screen, ?ferx_gam_screen, ?ferx_allometry, ?ferx_model_to_frem
ferx-core: covariates, covariate model, individual parameters, FREM, GAM screening