Special features

Overview

This chapter is a tour of optional .ferx blocks and ferx_fit() modes that each cover a narrow but important case: unit conversion ([scaling]), post-fit derived columns ([derived] / [output]), data filtering ([data_selection]), steady-state and multiple additional doses (SS / II / ADDL columns), absorption lag time, persisting fits to disk for later analysis, log-transform-both-sides residual error, and stochastic differential equations ([diffusion] / EKF process noise).

library(ferx)

[scaling] — unit conversion and custom readouts

The [scaling] block performs a post-prediction transform before the likelihood is computed. It has three forms.

Form A — scalar obs_scale

The simplest case: divide every prediction by a constant, e.g. when AMT is recorded in micrograms but DV in mg/L. The bundled warfarin_scaled example uses this:

[scaling]
  obs_scale = 1000

ferx_model_show(ferx_example("warfarin_scaled")$model)

# model: warfarin_scaled.ferx 
# One-compartment oral PK model (warfarin) with [scaling] block (Form A)
#
# Demonstrates obs_scale for unit conversion.
#
# The dataset records AMT in micrograms (100 mg = 100000 ug) while DV is
# in mg/L.  Without [scaling] the model would predict in ug/L, a factor-1000
# mismatch.  obs_scale = 1000 divides every prediction by 1000 before the
# residual is computed, so the sigma is expressed on the mg/L (DV) scale.
#
# PK parameter values are identical to the standard warfarin example; the
# fit should converge to the same estimates.

[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)

[individual_parameters]
  CL = TVCL * exp(ETA_CL)
  V  = TVV  * exp(ETA_V)
  KA = TVKA * exp(ETA_KA)

[structural_model]
  pk one_cpt_oral(cl=CL, v=V, ka=KA)

[error_model]
  DV ~ proportional(PROP_ERR)

[scaling]
  obs_scale = 1000

[fit_options]
  method     = foce
  maxiter    = 300
  covariance = true
  gradient   = ad

The scalar form is gradient-safe, so the model can keep its default gradient setting.

Form B — per-CMT obs_scale[CMT=N]

Different scale per observation compartment:

[scaling]
  obs_scale[CMT=2] = 1       # plasma — already in obs units
  obs_scale[CMT=3] = 0.001   # effect — convert ng to ug

Form C — explicit y[CMT=N] = <expression>

A full custom readout. Used in emax_pkpd (see Overview) to drive a Hill-Emax response off the effect-compartment amount:

[scaling]
  y[CMT=2] = central / V
  y[CMT=3] = E0 + EMAX * effect^GAMMA / (EC50^GAMMA + effect^GAMMA)

Form C and the per-CMT variants of Form B move the model onto the finite-difference gradient path (gradient = fd).

[derived] and [output] — post-fit sdtab columns

After a fit, fit$sdtab contains prediction-level columns (IPRED, PRED, CWRES, …). The [derived] and [output] blocks add arbitrary extra columns to that table — computed from the fitted individual parameters, predicted concentrations, and time variables — without rerunning the likelihood.

Three computation kinds

Syntax	Kind	Evaluated
COL = <expr>	PerRow	At every observation row
COL = max(IPRED, TIME < 24)	Aggregate	One value per subject (max/min/tmax with optional filter)
COL = integral(IPRED, from=0, to=72, step=0.2)	Integral	Numerical trapezoid over a time window per subject

All three share the same expression namespace: thetas, individual parameters (CL, V, …), covariates (AGE, WT, …), TIME / TAFD / TAD, and IPRED.

Per-row: KE and half-life

[derived]
  KE     = CL / V
  T_HALF = 0.6931472 / KE

Columns defined earlier in the block can be reused in later expressions — KE is in scope when T_HALF is evaluated.

Aggregate: Cmax and Tmax

[derived]
  CMAX = max(IPRED)
  TMAX = tmax(IPRED)
  CMAX_48 = max(IPRED, TIME < 48)

max, min, and tmax aggregate over all observation rows for the subject. An optional filter expression (second argument) restricts which rows are included.

Integral: AUC over a defined window

[derived]
  AUC_72 = integral(IPRED, from=0, to=72, step=0.2)

The integrand is re-evaluated on a fine time grid (step in the same units as TIME) rather than just at the observation rows. Any per-row expression is valid as the integrand.

Referencing compartment states

Beyond IPRED, derived expressions and integral(...) integrands can reference the full compartment state vector via the compartments[i] subscript (or, for ODE models, the named ODE state variables). This lets a derived column or AUC target a compartment other than the observed one — for example the effect-site or a peripheral compartment.

[derived]
  AUC_Ce = integral(compartments[1], from=0, to=24, step=0.5)

When an integral(...) integrand references compartments[i] or a named ODE state, ferx re-runs the solver (or analytical formula) at the grid points rather than interpolating from the stored per-observation states — an exact result at the cost of a second solver pass per subject. Named-state access in integrals requires an [odes] block; analytical models must use the compartments[i] subscript form.

Time-above-threshold (TAT)

The constant integrand 1.0 with a filter measures duration above a threshold:

[derived]
  TAT_0_5 = integral(1.0, IPRED > 0.5, from=0, to=72, step=0.1)

This gives hours above 0.5 mg/L within the window. For per-interval TAT (using TAD) the anchor and window form should be used — see the Derived columns reference for the full syntax.

Trough detection with TAD < 1e-10

Because TAD (time after dose) resets at each dose event, the smallest positive TAD values are at dose times — but floating-point arithmetic can make them 1e-13 rather than 0. Use a tolerance:

[derived]
  CMIN_TAU = min(IPRED, TAD < 1e-10)

Do not use TAD == 0.

[output] — individual parameters and covariates in sdtab

The [output] block echoes individual parameters and covariates into sdtab. Without it, fit$sdtab contains only IPRED, PRED, residuals, and any [derived] columns — not CL, V, AGE, etc.

[output]
  CL V KA AGE WT

Only items that are individual parameters (defined in [individual_parameters]) or covariates (columns in the data) can appear in [output]. Thetas cannot.

Runnable example

ex_d <- ferx_example("warfarin_derived")
ferx_model_show(ex_d$model)

# model: warfarin_derived.ferx 
# One-compartment oral PK model (warfarin) -- [derived] and [output] demo

[parameters]
  theta TVCL(0.2, 0.001, 10.0)
  theta TVV(10.0, 0.1, 500.0)
  theta TVKA(1.5, 0.01, 50.0)

  omega ETA_CL ~ 0.09
  omega ETA_V  ~ 0.04
  omega ETA_KA ~ 0.30

  sigma PROP_ERR ~ 0.02 (sd)

[individual_parameters]
  CL = TVCL * exp(ETA_CL)
  V  = TVV  * exp(ETA_V)
  KA = TVKA * exp(ETA_KA)

[structural_model]
  pk one_cpt_oral(cl=CL, v=V, ka=KA)

[error_model]
  DV ~ proportional(PROP_ERR)

[derived]
  # Per-row: elimination rate and half-life
  KE     = CL / V
  T_HALF = 0.6931472 / KE

  # Which dosing day (tau = 24 h)
  DAY      = floor(TAFD / 24) + 1
  TAU_TIME = TAFD mod 24

  # Subject-level aggregates (one scalar per subject, repeated across rows)
  CMAX       = max(IPRED)
  TMAX       = tmax(IPRED)
  # 1e-10 absorbs floating-point residuals when TAD is mathematically 0
  CTROUGH    = min(IPRED, TAD < 1e-10)
  CMAX_D1    = max(IPRED, TAFD < 24)
  CMAX_D14   = max(IPRED, TAFD >= 312 && TAFD < 336)

  # AUC over first 72 h (fine internal grid, 500 steps)
  AUC_0_72   = integral(IPRED, from=0, to=72)

  # Periodic AUC: one value per 24-h dosing window
  AUC_TAU    = integral(IPRED, window=24, anchor=0, step=0.1)

  # DV-based AUC (observation times only -- no interpolation)
  AUC_DV_72  = integral(DV, from=0, to=72)

[output]
  # Echo individual PK parameters (ETAs already automatic)
  CL V KA

[fit_options]
  method     = foce
  maxiter    = 300
  covariance = true

fit_d <- ferx_fit(ex_d$model, ex_d$data, verbose = FALSE)
names(fit_d$sdtab)
head(unique(fit_d$sdtab[, c("ID", "KE", "T_HALF", "CMAX", "AUC_0_72")]))

Tip

For the full reference and more examples see Model DSL: Derived columns on the ferx site.

---

[data_selection] — excluding records without touching the CSV

The [data_selection] block filters out records at read time. It is the ferx equivalent of NONMEM’s $DATA IGNORE=. The CSV is never modified; the filter is part of the model specification.

Syntax

[data_selection]
  ignore = DV < 1.0
  ignore = EVID == 0 && TIME > 120   # multiple conditions: use && within a line
  accept = BW >= 30                  # keep only rows where BW >= 30
  ignore_subjects = [7, 8]           # exclude subjects 7 and 8 entirely

Multiple ignore lines are OR-combined (any fired line excludes the record). Use && to combine conditions within a single ignore or accept line; || within a single line is not supported — use separate lines instead.

ignore_subjects takes a list of subject IDs (integers matching the ID column). A single ID can be written without brackets: ignore_subjects = 7. It excludes all records for those subjects regardless of any other condition.

An accept condition keeps only rows where the expression is true.

R-side preview with ferx_selection()

Before fitting, preview which records would be dropped:

ex_sel <- ferx_example("warfarin_data_selection")
sel <- ferx_selection(ex_sel$data, ignore = "DV < 1.0")
cat("Excluded obs:", sel$exclusions$n_obs_excluded, "\n")

Excluded obs: 

ferx_selection_excluded(sel)[, c("ID", "TIME", "DV")]

   ID TIME     DV
84  7  120 0.9761
96  8  120 0.8700

ferx_selection() returns the filtered data frame (plus an $exclusions list); ferx_selection_excluded() returns only the dropped rows.

What fit$exclusions reports

After fitting, fit$exclusions carries:

fit_sel <- ferx_fit(ex_sel$model, ex_sel$data, verbose = FALSE)
fit_sel$exclusions$n_records_total   # total rows read
fit_sel$exclusions$n_obs_excluded    # excluded observations
fit_sel$exclusions$fired_ignore      # which conditions fired

The print(fit) header also shows a DATA SELECTION line when any condition fired.

NONMEM equivalents

NONMEM	ferx
$DATA IGNORE=@ (comment rows)	handled automatically
$DATA IGNORE=(column.EQ.value)	ignore = column == value
$DATA ACCEPT=(column.GE.value)	accept = column >= value

Tip

For the full reference see Model DSL: Data selection and the worked example.

---

Steady-state dosing (SS, II columns)

For studies recorded after the pre-study dosing has reached steady state, add two columns to the dose record:

• SS = 1 on the dose row: ferx computes analytical steady-state initial conditions from the dose, interval, and the current individual parameters.
• II (inter-dose interval) in the same time units as TIME.

The bundled warfarin_ss example is one-compartment oral at QD steady state:

ex_ss <- ferx_example("warfarin_ss")
head(read.csv(ex_ss$data))

  ID TIME      DV EVID   AMT CMT RATE MDV SS II
1  1  0.0       .    1 100.0   1    0   1  1 24
2  1  0.5 23.3448    0     .   1    0   0  0  0
3  1  1.0 25.9891    0     .   1    0   0  0  0
4  1  2.0 26.7357    0     .   1    0   0  0  0
5  1  4.0 30.6482    0     .   1    0   0  0  0
6  1  8.0 23.8695    0     .   1    0   0  0  0

fit_ss <- ferx_fit(ex_ss$model, ex_ss$data, verbose = FALSE)

The structural model is unchanged from the single-dose case — only the data columns and the dose row’s SS/II values differ. ferx resolves the steady-state initial condition from the dose, the interval, and the current individual parameters; the supported PK shortcuts use an analytical closed-form.

ADDL — multiple additional doses

Recording repeated identical doses one row at a time inflates the dataset. The ADDL column (NONMEM-compatible) collapses a train of doses into a single row:

ID,TIME,DV,EVID,AMT,CMT,RATE,MDV,II,ADDL
1,0,.,1,70,1,0,1,24,6     # day 1 dose + 6 more = 7 total, q24h
1,24,3.8,0,.,1,0,0,.,.
1,168,1.4,0,.,1,0,0,.,.

ADDL=6 with II=24 means: repeat the same dose 6 more times (so 7 doses total), each 24 h apart. TAD resets to zero at each expanded dose event.

No special DSL syntax is required. The structural model is identical to a single-dose model — multi-dose bookkeeping lives entirely in the data columns.

Combining ADDL with [derived]

Because TAD resets correctly after each expanded dose, trough detection and within-interval AUC work as expected:

[derived]
  CMIN_TAU = min(IPRED, TAD < 1e-10)   # trough
  AUC_TAU  = integral(IPRED, window=24, anchor=0, step=0.1)

Note

The warfarin_addl model file is bundled and demonstrates this design. However, the matching dataset is not yet bundled — supply your own ADDL-format CSV that matches the column schema above.

---

Absorption lag time

Lag time can be expressed in two equivalent ways depending on which structural model you use.

As an argument to an analytical PK shortcut

[structural_model]
  pk one_cpt_oral(cl=CL, v=V, ka=KA, lagtime=LAGTIME)

The named keyword lagtime= (preferred) is accepted by the oral analytical solvers; alag= is an accepted alias. Both route to the same internal PK slot (lagtime is ferx-core’s name for what NONMEM calls ALAG). There is no tlag keyword.

As an LAGTIME keyword in [individual_parameters] (ODE models)

For ODE models the convention is to assign a special-named individual parameter called LAGTIME. The solver detects this name and delays every dose event for that subject by the individual lag time; no change to the ODE equations themselves is needed.

ex_lag <- ferx_example("warfarin_ode_lagtime")
ferx_model_show(ex_lag$model)

# model: warfarin_ode_lagtime.ferx 
# Warfarin ODE with absorption lag time (LAGTIME keyword)
#
# Demonstrates the LAGTIME keyword in [individual_parameters] for ODE models.
# Assigning to the reserved name LAGTIME tells the event handler to delay
# every dose event for that subject by the individual lag time; no change to
# the ODE equations is needed.
#
# This is the ODE-model equivalent of passing lagtime= to an analytical solver:
#
#   analytical:  pk one_cpt_oral(cl=CL, v=V, ka=KA, lagtime=TLAG)
#   ODE:         LAGTIME = TVLAG * exp(ETA_LAG)   # in [individual_parameters]
#
# IIV on LAGTIME is log-normal. The individual LAGTIME is recovered per
# subject in fit$individual_estimates alongside CL, V, and KA.
#
# States (CMT numbers):
#   1  depot   — oral dose lands here (amount, mg)
#   2  central — drug in central compartment (amount, mg)
#
# Dataset: data/warfarin.csv (100 mg single oral dose, 32 subjects × 7 obs)
#
# True simulation parameters match warfarin.ferx:
#   TVCL=0.134, TVV=8.1, TVKA=1.0, TVLAG=0.5

[parameters]
  theta TVCL(0.134, 0.001, 10.0)    # typical clearance (L/h)
  theta TVV(8.1,    0.1,  500.0)    # typical volume (L)
  theta TVKA(1.0,   0.01,  50.0)    # typical absorption rate constant (/h)
  theta TVLAG(0.5,  0.01,   5.0)    # typical absorption lag time (h)

  omega ETA_CL  ~ 0.07              # BSV on CL  (log-normal, ~26% CV)
  omega ETA_V   ~ 0.02              # BSV on V   (log-normal, ~14% CV)
  omega ETA_KA  ~ 0.40              # BSV on KA  (log-normal, ~63% CV)
  omega ETA_LAG ~ 0.09              # BSV on LAGTIME (log-normal, ~30% CV)

  sigma PROP_ERR ~ 0.01 (sd)        # proportional residual error (~10% CV)

[individual_parameters]
  CL      = TVCL * exp(ETA_CL)
  V       = TVV  * exp(ETA_V)
  KA      = TVKA * exp(ETA_KA)
  LAGTIME = TVLAG * exp(ETA_LAG)    # dose event delayed by this many hours

[structural_model]
  ode(obs_cmt=central, states=[depot, central])

[odes]
  d/dt(depot)   = -KA * depot
  d/dt(central) =  KA * depot - CL/V * central

[scaling]
  obs_scale = V

[error_model]
  DV ~ proportional(PROP_ERR)

[fit_options]
  method     = focei
  maxiter    = 500
  covariance = true

The LAGTIME line in [individual_parameters] is the only addition relative to a lag-free ODE model. IIV on LAGTIME is log-normal here, but any expression is valid — including the additive form shown in warfarin_additive_eta.

In both forms the individual lag value flows through [individual_parameters] like any other PK parameter — meaning it can be ETA-driven and recovered per subject from fit$individual_estimates.

Time-to-event endpoints — [event_model]

Warning

Beta / experimental. The [event_model] block is a beta feature (ferx-core #191/#192). Syntax and behaviour may change in a future release.

The [event_model] block registers a CMT column as a time-to-event (TTE) endpoint. Observations on that compartment are routed to a parametric survival likelihood instead of the Gaussian residual-error model. Three hazard families are supported: exponential, weibull, and gompertz.

[event_model]
  cmt    = <integer>    # CMT value identifying TTE rows in the data
  family = exponential  # exponential | weibull | gompertz
  scale  = <expression> # rate λ — Exponential and Weibull (rate= is an Exponential alias)
  shape  = <expression> # Weibull only (required)
  alpha  = <expression> # Gompertz only: baseline hazard at t=0
  gamma  = <expression> # Gompertz only: hazard growth rate
  loghr  = <expression> # optional, all families: multiplies the hazard by exp(loghr)

Important

Expressions are evaluated in the theta / eta / covariate namespace only. Names defined in [individual_parameters] are not in scope — write the full theta/eta expression directly (e.g. TVSCALE * exp(ETA_SCALE), not SCALE).

A standalone Weibull TTE model needs no [structural_model], [error_model], or [individual_parameters] block when [event_model] is the sole endpoint:

[parameters]
  theta TVSCALE(20.0, 0.1, 500.0)
  theta TVSHAPE(1.5,  0.1, 10.0)

  omega ETA_SCALE ~ 0.04

[event_model]
  cmt    = 2
  family = weibull
  scale  = TVSCALE * exp(ETA_SCALE)
  shape  = TVSHAPE

[fit_options]
  method  = focei
  maxiter = 300

DV is coded 0 (right-censored), 1 (exact event), or 2 (interval-censored right bound, paired with a preceding DV=0 row on the same CMT). Named blocks ([event_model DROPOUT], [event_model DEATH]) allow multiple TTE endpoints in one model.

Saving and reloading fits

Long-running fits should be persisted so you do not pay the estimation cost twice. ferx_save_fit() writes a .fitrx file with the full fit object plus SHA-256 hashes of the model file and dataset; ferx_load_fit() restores it and warns if either has changed.

ferx_save_fit(fit, "warfarin.fitrx")                     # data not embedded
ferx_save_fit(fit, "warfarin.fitrx", include_data = TRUE)

fit_reloaded <- ferx_load_fit("warfarin.fitrx")

include_data = TRUE embeds the dataset as a base64 blob in the .fitrx file. The file is then self-contained — useful for sharing across machines, attaching to a regulatory submission, or moving a fit between branches of an analysis without the dataset’s path drifting.

Bioavailability

Estimating bioavailability F as a population parameter benefits from a logit transform so it stays in (0, 1). The bundled bioavailability example uses inv_logit(logit(THETA_F) + ETA_F); see Logit-normal for bounded parameters for the parameterisation and reporting conventions.

Log-transform-both-sides (LTBS) residual error

The familiar NONMEM idiom Y = LOG(F) + EPS(1) — fit log-transformed observations against log-predictions with additive error on the log scale — is written in .ferx as one of two equivalent forms:

[error_model]
  log(DV) ~ additive(SIGMA_LOG)   # DV on the natural scale; engine logs it

[error_model]
  DV ~ log_additive(SIGMA_LOG)    # DV already log-transformed in the data

The bundled warfarin_ltbs example uses the first form:

ferx_model_show(ferx_example("warfarin_ltbs")$model)

# model: warfarin_ltbs.ferx 
# One-compartment oral PK (warfarin) with log-transform-both-sides (LTBS)
# residual error. Writing `log(DV) ~ additive(SIGMA_LOG)` tells the engine to
# log-transform observations and predictions and fit additive (normal) error
# on the log scale. The sigma estimate is on the log scale; under LTBS
# `sigma` approximates the CV on the natural scale.
#
# Equivalent alternative form: `DV ~ log_additive(SIGMA_LOG)` (when DV is
# already supplied as log-transformed).

[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 SIGMA_LOG ~ 0.1 (sd)

[individual_parameters]
  CL = TVCL * exp(ETA_CL)
  V  = TVV  * exp(ETA_V)
  KA = TVKA * exp(ETA_KA)

[structural_model]
  pk one_cpt_oral(cl=CL, v=V, ka=KA)

[error_model]
  log(DV) ~ additive(SIGMA_LOG)

[fit_options]
  method     = focei
  maxiter    = 300
  covariance = true

Under LTBS the diagnostic outputs (IPRED, PRED, CWRES, IWRES in sdtab, and DV_SIM from ferx_simulate()) are on the log scale. ferx_model_inspect() reports the residual type as additive (log-transformed).

Stochastic differential equations (SDE / [diffusion])

When IWRES shows autocorrelation that cannot be absorbed by the residual error alone — residual structural misspecification, unmodelled absorption variability, within-subject biological noise — the standard ODE treats the system as deterministic and forces all unexplained variation into ETA and EPS. The [diffusion] block turns the ODE into a stochastic differential equation (SDE): the structural dynamics are themselves subject to continuous Gaussian process noise.

ferx uses the Extended Kalman Filter (EKF) to propagate both the state and its covariance forward in time. At each observation the integrated covariance inflates the residual variance:

\[V_\text{total}(t) = P_\text{central}(t) + \sigma^2(\hat{y})\]

The diffusion parameter is fitted on the variance scale (not SD) alongside the structural thetas.

DSL syntax

Add a [diffusion] block after [odes]:

[diffusion]
  central ~ 0.01

The number is the initial estimate for DIFF_CENTRAL, the process-noise variance on the central compartment state.

The bundled warfarin_sde example is a one-compartment oral ODE with SDE process noise:

ferx_model_show(ferx_example("warfarin_sde")$model)

# model: warfarin_sde.ferx 
# One-compartment oral PK model (warfarin) with SDE process noise
#
# Identical to the standard warfarin model but written as an ODE so that a
# [diffusion] block can be added.  The DIFF_CENTRAL parameter captures
# within-subject system noise on the central compartment — useful when
# IWRES shows autocorrelation that cannot be absorbed by the residual error
# alone (structural misspecification, unmodelled absorption variability, etc.).
#
# The EKF inflates the residual variance at each observation by the
# integrated covariance: V_total = P[central] + sigma^2(PRED).
# DIFF_CENTRAL is estimated on the variance scale (not SD).
#
# Note: SDE models use finite differences for the ODE+EKF evaluation, which is
# significantly slower than an equivalent analytical PK model.

[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)

[individual_parameters]
  CL = TVCL * exp(ETA_CL)
  V  = TVV  * exp(ETA_V)
  KA = TVKA * exp(ETA_KA)

[structural_model]
  ode(obs_cmt=central, states=[depot, central])

[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

[diffusion]
  central ~ 0.01

[error_model]
  DV ~ proportional(PROP_ERR)

[fit_options]
  method          = foce
  maxiter         = 300
  covariance      = true
  gradient_method = fd

Fitting an SDE model

ex_sde <- ferx_example("warfarin_sde")
fit_sde <- ferx_fit(ex_sde$model, ex_sde$data)

Key outputs that differ from a standard ODE fit:

fit_sde$uses_sde               # TRUE
fit_sde$theta["DIFF_CENTRAL"]  # fitted diffusion variance
ferx_estimates(fit_sde)        # DIFF_CENTRAL appears in the theta block

SDE-specific notes

Important

Gradient method: SDE models require gradient_method = fd. gradient_method = fd is set automatically in the model file — no manual override is needed. Finite-difference gradients make each function evaluation significantly slower than an equivalent analytical PK model; this slowdown is inherent to SDE fitting.

Property	SDE model	Standard ODE
fit$uses_sde	TRUE	FALSE
Gradient method	fd (required)	fd
Additional theta	DIFF_CENTRAL (variance)	—
CWRES / IWRES	Available	Available
EKF covariance propagation	Yes	No

Standard diagnostics (sdtab, CWRES, IWRES) remain available; they reflect the EKF-inflated total variance.

Tip

For a full worked example see Example: SDE on the ferx website.

Beyond this chapter

• [covariate_nn NAME] — Deep Compartment Models, a neural-network map from covariates to PK parameters. Feature-flagged (cargo --features nn); see the ferx-site DCM example.
• [derived] / [output] full reference — Model DSL: Derived columns covers all aggregate forms, anchor/window integrals, naming rules, and sequential scoping rules.
• [data_selection] full reference — Model DSL: Data selection documents evaluation order, ignore vs accept vs ignore_subjects, and supported column types.
• Worked examples — additional examples on the ferx site: Derived columns, Derived PK/PD, Data selection, ODE with TIME/TAFD/TAD.