Experimental features
Where you are
Two model extensions in ferx are explicitly experimental: stochastic differential equations, which add random noise to the model states between observations, and neural networks for the covariate model, known as deep compartment models. Both fit from R with the bundled examples. This chapter runs them, compares them with their conventional counterparts, and shows what the fits report. It closes with parameter priors, which carry the same label without being a model extension.
ferx-core labels [diffusion] (SDE) and [covariate_nn] (neural networks) experimental: validated on a small set of examples, with syntax that may still change. Every fit that uses them carries an experimental warning. The feature maturity table’s third experimental entry is the Vern7 ODE stepper (Structural models: analytical and ODE). Parameter priors, at the end of this chapter, are labelled experimental on their own page but are absent from that table, and a priored fit raises no experimental warning.
The data
warfarin_sde uses the warfarin concentrations (10 subjects). warfarin_dcm uses the 30-subject dataset of the covariate thread, with weight and creatinine clearance:
Minimal runnable call: an SDE model
An SDE model is an ODE model with a [diffusion] block. Each line gives a state a diffusion variance, estimated as the parameter DIFF_<STATE>. ferx evaluates the likelihood with an extended Kalman filter, which adds the propagated state uncertainty to the residual variance at each observation:
invisible(ferx_model_get_section(sde$model, "odes"))
#> # [odes]
#> # `central` is concentration (mg/L): the absorption term carries `/V` so the
#> # ODE state matches the concentration-scale DV column. See the existing
#> # `examples/bioavailability_ode.ferx` in ferx-core for the same convention.
#> d/dt(depot) = -KA * depot
#> d/dt(central) = KA * depot / V - (CL / V) * central
invisible(ferx_model_get_section(sde$model, "diffusion"))
#> # [diffusion]
#> central ~ 0.01
fit_sde <- ferx_fit(sde$model, sde$data, verbose = FALSE)
fit_sde$uses_sde
#> [1] TRUE
fit_sde$estimates[, c("estimate", "rse_pct")]
#> estimate rse_pct
#> TVCL 0.13406927696 40780.260
#> TVV 7.67239278287 1636.661
#> TVKA 0.74065619542 272022.128
#> DIFF_CENTRAL 0.00003101962 246142.085
#> ETA_CL 0.02758213931 89501.362
#> ETA_V 0.00993763038 275367.023
#> ETA_KA 0.29894993012 115017.715
#> PROP_ERR 0.00993620496 38136.534Reading the result
The fit carries the experimental notice with its other warnings:
ferx_get_warnings(fit_sde, as_df = TRUE)[, c("severity", "category")]
#> severity category
#> 1 warning experimental
#> 2 warning covariance_regularized
#> 3 warning inflated_rse
#> 4 info ode_solver
#> 5 warning dw_autocorrelation
#> 6 warning high_correlation
#> 7 critical condition_number
substr(ferx_get_warnings(fit_sde, as_df = TRUE)$message[1], 1, 150)
#> [1] "Stochastic differential equations ([diffusion] / Extended Kalman Filter) are an EXPERIMENTAL feature: validated only on a small set of toy examples, w"The diffusion variance is small, and the covariance step does not support it: the relative standard error of every parameter, not only of DIFF_CENTRAL, runs to thousands of percent (the smallest is 1,637%), and the fit carries covariance_regularized, inflated_rse, high_correlation and a critical condition_number warning. The standard errors of this fit are not usable. Removing the [diffusion] block gives the ordinary ODE model, which reaches a similar OFV:
sde_lines <- readLines(sde$model)
start <- grep("^\\[diffusion\\]", sde_lines)
ode_model <- file.path(experimental_dir, "warfarin_ode_without_diffusion.ferx")
writeLines(sde_lines[-(start:(start + 2))], ode_model)
fit_ode <- ferx_fit(ode_model, sde$data, verbose = FALSE)
data.frame(model = c("with [diffusion]", "without"), ofv = c(fit_sde$ofv, fit_ode$ofv),
durbin_watson = c(fit_sde$dw_statistic, fit_ode$dw_statistic))
#> model ofv durbin_watson
#> 1 with [diffusion] -276.3050 2.569116
#> 2 without -280.0208 2.609533These data show no system noise. A [diffusion] term is not a remedy for residuals that are correlated in time within a subject (Diagnosing the model). The filter propagates the state variance and adds it to the residual variance; the predicted state itself stays the ODE solution and is never pulled towards a subject’s observations, so the term cannot follow a drift. Here the Durbin-Watson statistic barely moves, 2.57 with the block and 2.61 without. What [diffusion] does estimate is how much unexplained variance builds up between observations.
Options that matter
| Block or setting | Meaning |
|---|---|
[diffusion] STATE ~ variance
|
Diffusion variance for an ODE state, estimated as DIFF_STATE (FIX to hold it) |
[covariate_nn NAME] |
A neural network with inputs (covariates), outputs (parameter names), hidden layers, activation and output activation |
NAME.OUTPUT in [individual_parameters]
|
A network output, used like any other quantity |
book_settings_table(c("nn_l2", "nn_smooth"))Constraints from the ferx-core diffusion page: [diffusion] needs an ODE model, cannot be combined with SAEM or [scaling] (Observation scaling with [scaling]), and uses finite-difference gradients.
Variants
Estimation method for the SDE model
The bundled model uses FOCE. With FOCEI the diffusion variance is even smaller. The OFVs of the two methods are different objectives:
# suppress the R warning that reports overriding the model file's method
fit_sde_focei <- suppressWarnings(ferx_fit(sde$model, sde$data, method = "focei", verbose = FALSE))
rbind(foce = c(ofv = fit_sde$ofv, fit_sde$theta), focei = c(ofv = fit_sde_focei$ofv, fit_sde_focei$theta))
#> ofv TVCL TVV TVKA DIFF_CENTRAL
#> foce -276.3050 0.1340693 7.672393 0.7406562 0.00003101962
#> focei -272.2985 0.1333708 7.757245 0.8115071 0.00005155504Simulation from an SDE fit works like any other simulation (Simulating scenarios):
head(ferx_simulate(sde$model, sde$data, n_sim = 1, seed = 1, fit = fit_sde), 3)
#> DRAW SIM ID TIME CMT IPRED DV_SIM OBSERVED
#> 1 1 1 1 0.5 1 3.516303 3.494879 NA
#> 2 1 1 1 1.0 1 5.968058 5.914088 NA
#> 3 1 1 1 2.0 1 8.832364 9.000727 NAA neural-network covariate model
warfarin_dcm replaces the covariate expressions of two_cpt_oral_cov with a network. The network reads WT and CRCL and outputs a multiplier for each PK parameter. The multipliers start at 1, so the thetas carry the typical values:
# print the two blocks without their comment lines
section_lines <- function(block) {
capture.output(lines <- ferx_model_get_section(dcm$model, block, strip = TRUE))
trimws(lines[nzchar(trimws(lines)) & !grepl("^\\s*#", lines)])
}
section_lines("covariate_nn TYPICAL_PK")
#> [1] "inputs = [WT, CRCL]" "outputs = [CL, V1, Q, V2, KA]"
#> [3] "layers = [8, 8]" "activation = tanh"
#> [5] "output = exp"
section_lines("individual_parameters")
#> [1] "CL = TVCL * TYPICAL_PK.CL * exp(ETA_CL)"
#> [2] "V1 = TVV1 * TYPICAL_PK.V1 * exp(ETA_V1)"
#> [3] "Q = TVQ * TYPICAL_PK.Q * exp(ETA_Q)"
#> [4] "V2 = TVV2 * TYPICAL_PK.V2 * exp(ETA_V2)"
#> [5] "KA = TVKA * TYPICAL_PK.KA * exp(ETA_KA)"The network has 2 inputs, two hidden layers of 8 and 5 outputs: 141 weights, estimated as thetas. The covariance step is skipped here:
# suppress the R warning that reports overriding the model file's covariance setting
fit_dcm <- suppressWarnings(ferx_fit(dcm$model, dcm$data, covariance = FALSE, verbose = FALSE))
network <- fit_dcm$neural_networks$TYPICAL_PK
network[c("shape", "hidden_activation", "output_activation", "n_weights", "input_names", "output_names")]
#> $shape
#> [1] 2 8 8 5
#>
#> $hidden_activation
#> [1] "tanh"
#>
#> $output_activation
#> [1] "exp"
#>
#> $n_weights
#> [1] 141
#>
#> $input_names
#> [1] "WT" "CRCL"
#>
#> $output_names
#> [1] "CL" "V1" "Q" "V2" "KA"
ferx_get_warnings(fit_dcm, as_df = TRUE)[, c("severity", "category")]
#> severity category
#> 1 warning experimental
#> 2 warning dw_autocorrelationThe same data with the base model and with the explicit covariate model:
base <- ferx_example("two_cpt_oral_base")
cov <- ferx_example("two_cpt_oral_cov")
fit_base <- suppressWarnings(ferx_fit(base$model, base$data, covariance = FALSE, verbose = FALSE))
fit_cov <- suppressWarnings(ferx_fit(cov$model, cov$data, covariance = FALSE, verbose = FALSE))
data.frame(model = c("two_cpt_oral_base (no covariates)", "two_cpt_oral_cov (WT, CRCL)", "warfarin_dcm (network)"),
parameters = c(fit_base$n_parameters, fit_cov$n_parameters, fit_dcm$n_parameters),
ofv = c(fit_base$ofv, fit_cov$ofv, fit_dcm$ofv),
aic = c(fit_base$aic, fit_cov$aic, fit_dcm$aic))
#> model parameters ofv aic
#> 1 two_cpt_oral_base (no covariates) 11 -1185.403 -1163.4028
#> 2 two_cpt_oral_cov (WT, CRCL) 13 -1199.326 -1173.3264
#> 3 warfarin_dcm (network) 152 -1185.462 -881.4621With 30 subjects the network does not improve on the model without covariates, and its 141 weights make the AIC far worse than either alternative. The explicit covariate model captures the covariate effects with two parameters.
How much the network departs from multipliers of 1 can be read from the individual parameters. For clearance, CL / (TVCL * exp(ETA_CL)) is the network’s multiplier for each subject:
subjects <- read.csv(dcm$data) |> distinct(ID, WT, CRCL)
multipliers <- mutate(fit_dcm$individual_estimates, ID = as.integer(ID)) |>
select(ID, CL) |>
left_join(mutate(fit_dcm$ebe_etas, ID = as.integer(ID)) |> select(ID, ETA_CL), by = "ID") |>
left_join(subjects, by = "ID") |>
mutate(multiplier = CL / (fit_dcm$theta[["TVCL"]] * exp(ETA_CL)))
summary(multipliers$multiplier)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1.032 1.032 1.032 1.032 1.032 1.032
range(subjects$WT)
#> [1] 45.0 93.7
range(subjects$CRCL)
#> [1] 46.5 150.0The multiplier is the same for every subject to seven significant digits, across the whole range of weight and creatinine clearance: the fitted network learned no covariate dependence at all, and only rescaled the typical clearance.
Regularizing the network
nn_l2 (weight decay) and nn_smooth (a penalty on curvature) keep a network from fitting spurious covariate effects. They are set in [fit_options] or with settings, apply to the FOCE-family methods, and the reported OFV stays the unpenalized one, so regularized and unregularized fits remain comparable (regularization). On this example a fit with settings = list(nn_l2 = 0.01) did not finish within 10 minutes, with either the model’s lbfgs or the bobyqa optimizer, against about half a minute without the penalty. It is therefore not run in this book.
Parameter priors
A prior pins a parameter near a value you already believe and lets the data move it only as far as the data can pay for – penalized maximum likelihood, the mechanism behind updating a published model on local data. It is written inline on the [parameters] row it belongs to, as a central value and a relative standard error, both on the scale the parameter was declared on:
warf <- ferx_example("warfarin")
priored <- file.path(experimental_dir, "warfarin_prior.ferx")
writeLines(sub("theta TVCL(0.134, 0.001, 10.0)",
"theta TVCL(0.134, 0.001, 10.0) prior(0.25, rse = 10%)",
readLines(warf$model), fixed = TRUE), priored)
fit_plain <- ferx_fit(warf$model, warf$data, verbose = FALSE)
fit_prior <- ferx_fit(priored, warf$data, verbose = FALSE)
c(no_prior = fit_plain$theta[["TVCL"]], with_prior = fit_prior$theta[["TVCL"]])
#> no_prior with_prior
#> 0.1329690 0.1627248These data want a TVCL near 0.133; a prior centred at 0.25 with a 10% relative standard error pulls the estimate to 0.1627. The penalty is added to the objective, so fit$ofv is the penalized total and does not compare with the unpenalized fit. The fit reports the two halves separately: fit$ofv_data is the data’s share, fit$ofv_prior the penalty, and they sum to fit$ofv. The information criteria are built from the data half alone, so fit$aic is fit$ofv_data + 2 * fit$n_parameters and not fit$ofv + 2 * fit$n_parameters once a prior is in the model:
data.frame(fit = c("no prior", "prior on TVCL"),
ofv = c(fit_plain$ofv, fit_prior$ofv),
ofv_data = c(fit_plain$ofv_data, fit_prior$ofv_data),
ofv_prior = c(fit_plain$ofv_prior, fit_prior$ofv_prior),
aic = c(fit_plain$aic, fit_prior$aic),
ofv_data_plus_2k = c(fit_plain$ofv_data + 2 * fit_plain$n_parameters,
fit_prior$ofv_data + 2 * fit_prior$n_parameters))
#> fit ofv ofv_data ofv_prior aic ofv_data_plus_2k
#> 1 no prior -280.3640 -280.3640 0.0000 -266.3640 -266.3640
#> 2 prior on TVCL -250.4404 -268.9709 18.5305 -254.9709 -254.9709A fit without a prior has ofv_data equal to ofv and an ofv_prior of 0. ofv_data is the number to set beside the unpenalized fit: -268.97 against -280.36, so pulling TVCL towards 0.25 costs the data 11.39 points of objective, on top of a penalty of 18.53. print() and summary() state the split under the OFV line:
grep("OFV|prior", capture.output(print(fit_prior)), value = TRUE)
#> [1] " Model: warfarin_prior Dataset: warfarin"
#> [2] " OFV: -250.4404 AIC: -254.9709 BIC: -236.0676"
#> [3] " (prior on 1 parameter: data OFV -268.9709 + prior penalty 18.5305; AIC/BIC use the data half)"fit$prior_summary has one row per priored parameter, and is NULL for a model without one:
fit_prior$prior_summary
#> name prior_value estimate shift_in_prior_sds penalty family
#> 1 TVCL 0.25 0.1627248 -4.304707 18.5305 lognormal
#> prior_lower_95 prior_upper_95
#> 1 0.205604 0.3039824
is.null(fit_plain$prior_summary)
#> [1] TRUEshift_in_prior_sds is how far the estimate landed from the prior’s centre, with its sign, in prior standard deviations on the scale the prior is evaluated on. For a lognormal family that is the log scale, so it is not (estimate - 0.25) / 0.025. penalty is its square, and that parameter’s share of ofv_prior. prior_lower_95 and prior_upper_95 are the interval that rse = 10% implies for the family the engine chose (?ferx_fit describes every column). An estimate 4.3 prior standard deviations from the centre, outside that interval, says the prior and these data disagree. The standard errors of a priored fit are those of the penalized objective. What the prior cost in total is fit_prior$ofv - fit_plain$ofv, 29.92 points: the data’s share plus the penalty.
A [priors] block takes the priors from a previous fit instead of writing them out by hand. It has exactly one key, from_fit, naming a saved fit. That path does not work from R at this build: the only fit file ferx-r writes is the .fitrx bundle of ferx_save_fit(), and the engine’s reader rejects it.
saved <- file.path(experimental_dir, "warfarin.fitrx")
ferx_save_fit(fit_plain, saved)
from_fit_model <- file.path(experimental_dir, "warfarin_from_fit.ferx")
# normalizePath(): `//` starts a comment in a .ferx file, and R's tempdir()
# contains one on macOS, which would truncate the path before the engine saw it.
writeLines(c(readLines(warf$model), "", "[priors]",
paste0(" from_fit = ", normalizePath(saved))),
from_fit_model)
try(ferx_fit(from_fit_model, warf$data, verbose = FALSE))
#> Error in ferx_rust_fit(model_path = normalizePath(model), data_path = normalizePath(data), :
#> Error parsing model: [priors] from_fit: failed to read `/tmp/Rtmp6tL7So/experimental/warfarin.fitrx`: JSON error: invalid type: string "foce", expected a sequence at line 3 column 24 [E_PARSE]The bundle stores every one-element array as a plain value, so the first such field the reader meets – here the method chain, which it expects as a list of methods – ends the import. That one field is not the whole of it: a chained fit writes a correct list and then stops at the next single-element field instead (ferx-r #379). Until the writer and the format agree, write priors inline.
Pitfalls
- Experimental means validate yourself. Compare an SDE or network model with its conventional counterpart, as above, before relying on it.
-
[priors] from_fitdoes not work from R at this build. The.fitrxbundleferx_save_fit()writes is rejected by the engine’s reader (ferx-r #379). Write the prior inline on the[parameters]row instead, as above. - Extra flexibility is not extra information. A diffusion variance or a network with many weights needs data that support it; check standard errors, the OFV change and information criteria.
-
Standard errors of network weights are not reliable, and the covariance step with many weights is slow. Skip it with
covariance = FALSEwhile exploring. - SDE fits are slow. The extended Kalman filter uses finite differences; the warfarin SDE fit takes many times longer than the ODE fit.
Warnings you may see here
-
experimental(warning): the model uses an experimental component ([diffusion]or[covariate_nn]). The message states what has and has not been validated. It appears on every such fit.
Summary
-
[diffusion]adds system noise to ODE states, estimated with an extended Kalman filter; compare with the model without it. -
[covariate_nn]lets a neural network produce covariate multipliers; compare its OFV and AIC with an explicit covariate model. -
nn_l2andnn_smoothregularize a network. - An inline
prior(value, rse = …)regularizes a parameter towards a value you already believe; the reported OFV is then penalized, and the[priors]block’sfrom_fitimport is not reachable from R at this build. -
[diffusion]and[covariate_nn]fits carry anexperimentalwarning.
Next: Function, option and example index lists every function, option and example with the chapter that covers it.
- R help:
?ferx_fit,?ferx_save_fit - ferx-core: SDE and diffusion, covariate neural networks, neural networks overview, parameter priors