Time-to-Event (TTE) Models

Maturity: beta — see Feature Maturity for what this means.

Phase 1 (current): Parametric TTE with exponential, Weibull, and Gompertz hazard families; right-censored, exact, and interval-censored observations; left truncation via TENTRY; FOCEI and SAEM estimation; fixed-effects (n_eta = 0) and random-effects paths. [structural_model], [error_model], and [individual_parameters] are all optional for TTE-only models.

Overview

TTE models describe the probability distribution of the time until a clinical event (e.g. first adverse event, study dropout, disease progression). ferx uses the same FOCEI Laplace objective as for continuous PK/PD data, extended to the per-subject log-likelihood contribution

ℓᵢ = δᵢ · log h(Tᵢ) − H(Tᵢ)

where δᵢ is the event indicator, h(t) the hazard function, and H(t) = ∫₀ᵗ h(s) ds the cumulative hazard.

Model file syntax

Add one [event_model] block per TTE endpoint.

[event_model]
  cmt    = 2                        # data-file CMT value that carries TTE rows
  family = exponential              # exponential | weibull | gompertz
  scale  = TVLAMBDA * exp(ETA_LAMBDA)  # theta/eta/covariate expression

Note: [event_model] parameter expressions are evaluated in the theta / eta / covariate namespace, and may also reference names defined in [individual_parameters] (e.g. LAMBDA); those are resolved per subject at evaluation time. Writing the full expression in terms of theta and eta names directly also works. A referenced individual parameter may not depend on an inter-occasion (IOV/kappa) random effect — the hazard is evaluated once per subject, with no occasion context — and such a reference is rejected with an error.

For Weibull, a shape parameter is also required:

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

Gompertz uses alpha (baseline hazard) and gamma (growth rate):

[event_model]
  cmt    = 2
  family = gompertz
  alpha  = TVALPHA
  gamma  = TVGAMMA

Multiple TTE endpoints are supported by repeating the block with unique names:

[event_model DROPOUT]
  cmt    = 2
  family = exponential
  scale  = TVLAMBDA_DROPOUT * exp(ETA_LAMBDA)

[event_model DEATH]
  cmt    = 3
  family = weibull
  scale  = TVSCALE_DEATH
  shape  = TVSHAPE_DEATH

[structural_model] and [error_model] requirement

The parser currently requires both blocks even for TTE-only models. Include a dummy 1-cpt block with fixed parameters:

[parameters]
  theta TVLAMBDA(0.05, 0.001, 10.0)
  theta DUMMY_CL(1.0, FIX)
  theta DUMMY_V(1.0, FIX)
  omega ETA_LAMBDA ~ 0.09
  sigma SIGMA_DV ~ 0.01 FIX

[individual_parameters]
  LAMBDA = TVLAMBDA * exp(ETA_LAMBDA)
  CL     = DUMMY_CL
  V      = DUMMY_V

[structural_model]
  pk one_cpt_iv(cl=CL, v=V)     # inert: no CMT-1 observations in TTE-only data

[error_model]
  DV ~ additive(SIGMA_DV)       # inert: sigma carries zero gradient

[event_model]
  cmt    = 2
  family = exponential
  scale  = TVLAMBDA * exp(ETA_LAMBDA)  # theta/eta directly, or an [individual_parameters] name

This [structural_model] / [error_model] requirement will be lifted in a future release. (The separate [individual_parameters]-namespace restriction noted above has already been lifted — hazard expressions may reference individual parameters.)

Data format

TTE rows follow the NONMEM observation convention with EVID=0 and a CMT that matches the cmt = declaration in [event_model]:

DV Meaning
0 Right-censored (survived past T)
1 Exact event at T
2 Interval-censored: right bound (pair with preceding DV=0 row for same CMT)

For interval-censored observations, precede the DV=2 row with a DV=0 row carrying the left bound (entry into the at-risk window):

ID,TIME,DV,EVID,CMT
1, 5.0, 0, 0, 2   # left bound of interval (entry)
1,10.0, 2, 0, 2   # right bound  →  event occurred in (5, 10]

Left truncation (delayed entry)

If subjects enter the study at a time > 0 (e.g. they survived until enrolment), add a TENTRY column. ferx computes

H_eff(T) = H(T) − H(TENTRY)

conditioning the likelihood on survival past the entry time:

ID,TIME,DV,EVID,CMT,TENTRY
1,18.3, 1, 0, 2, 5.0   # subject entered at t=5, event at t=18.3
2,30.0, 0, 0, 2, 5.0   # entered at t=5, censored at t=30

Hazard families

Family Parameters h(t) H(t)
exponential scale = λ λ · exp(loghr) λ · exp(loghr) · t
weibull scale = α, shape = γ (γ/α)(t/α)^(γ−1) · exp(loghr) (t/α)^γ · exp(loghr)
gompertz alpha = α, gamma = γ α · exp(γ·t) · exp(loghr) (α/γ)(exp(γ·t) − 1) · exp(loghr)

All parameters must be positive. The optional loghr key adds a proportional-hazards covariate term: the entire hazard (and cumulative hazard) is multiplied by exp(loghr). When loghr is omitted it defaults to 0 (no effect).

# Exponential with PH covariate on sex (SEX=1 reference):
[event_model]
  cmt    = 2
  family = exponential
  scale  = TVLAMBDA * exp(ETA_LAMBDA)
  loghr  = BETA_SEX * SEX

# Weibull with PH covariate on weight:
[event_model]
  cmt    = 2
  family = weibull
  scale  = TVSCALE * exp(ETA_SCALE)
  shape  = TVSHAPE
  loghr  = BETA_WT * WT

Estimation

The FOCEI Laplace objective includes a TTE-specific FD Hessian contribution:

OFV = −2 · Σᵢ [ ℓᵢ(η̂ᵢ) − ½ log |Ω⁻¹ + ∇²ℓᵢ(η̂ᵢ)| ]

where ∇²ℓᵢ is computed via a 4-point central-difference finite-difference stencil (Shi 2021 step-size selection).

Set method = focei in [fit_options]. The gradient = fd setting is the default and the only supported path for TTE objectives.

Joint PK-TTE

When the hazard depends on drug exposure, write it as an expression over the PK state instead of a closed-form family. On an ODE model, the hazard key in [event_model] is accumulated as an extra ODE compartment and estimated jointly with the PK, with shared random effects:

[parameters]
  theta TVCL(1.0, 0.01, 100); theta TVV(10, 0.1, 500); theta TVKA(1, 0.01, 50)
  theta TVH0(0.01, 1e-5, 10); theta TVBETA(0.5, -10, 10)
  omega ETA_CL ~ 0.09
  sigma PROP_ERR ~ 0.05 (sd)

[individual_parameters]
  CL   = TVCL * exp(ETA_CL)
  V    = TVV
  KA   = TVKA
  H0   = TVH0
  BETA = TVBETA

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

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

[event_model]
  cmt    = 3                                  # TTE on its own CMT (PK observed on CMT 2)
  hazard = H0 * exp(BETA * (central / V))     # log-linear hazard in concentration

[error_model]
  DV ~ proportional(PROP_ERR)                  # Gaussian PK rows (CMT 2); TTE rows (CMT 3) route to the event likelihood

[fit_options]
  method = focei

The parser appends d/dt(__chz) = H0 * exp(BETA * (central / V)) (init 0) to the ODE system; the TTE likelihood reads H(T) from that state and h(T) from its derivative, combined with the Gaussian PK likelihood in one FOCEI objective (central / V is the concentration driving the hazard). Put the PK observation and the event on different compartments so the data reader can route Gaussian rows and TTE rows unambiguously.

predict_survival() returns S(t) / H(t) / h(t) for the endpoint by integrating the same augmented system at the typical values (its mean_survival is NaN — an ∫₀^∞ S summary for ODE hazards is a follow-up; median_survival is interpolated from the cumulative-hazard grid).

Simulating drug-driven event times

simulate() draws event times from a drug-driven hazard by integrating the augmented PK + cumulative-hazard system until H(T) = −log U — the ODE analogue of the analytic inverse-CDF — localizing the crossing inside the solver step. Because a drug-driven hazard can vanish (the drug clears, the hazard falls toward baseline, and an event may never fire), an ODE-TTE simulation requires a finite [simulation] horizon; a subject that has not crossed by the horizon is right-censored there. EVID-3/4 resets and left-truncated entry (entry_time > 0) are not yet supported for ODE-TTE simulation and raise a clear error rather than mis-sampling. A negative or non-finite hazard fails loudly during sampling — it is never silently treated as a censor.

NoteCurrent limitations

The hazard expression lives in the [odes] RHS namespace (ODE states + individual parameters — route thetas/covariates through [individual_parameters]), and inter-occasion variability (IOV) is rejected. Joint PK-TTE is the standard NONMEM $DES + cumulative-hazard construct; recovery against NONMEM / nlmixr2 is tracked in issue #564. See examples/pktte_joint.ferx.

Competing risks (cause-specific hazards)

When a subject is at risk of several mutually-exclusive event types, give each cause its own [event_model NAME] block on a distinct compartment. The blocks share the model’s random effects (a common frailty links the causes), so this is the standard cause-specific hazard approach:

[parameters]
  theta TVLAMBDA_A(0.10, 0.001, 10.0)
  theta TVLAMBDA_B(0.06, 0.001, 10.0)
  omega ETA_F ~ 0.25

[event_model cause_a]
  cmt    = 2
  family = exponential
  scale  = TVLAMBDA_A * exp(ETA_F)

[event_model cause_b]
  cmt    = 3
  family = exponential
  scale  = TVLAMBDA_B * exp(ETA_F)

Data format. Each subject has one row per cause CMT. The cause that occurred is DV=1 on its CMT; every other cause is right-censored (DV=0) at the same time, because experiencing one event removes the subject from risk of the others:

ID,TIME,DV,EVID,AMT,CMT,RATE,MDV
1,12.25,1,0,0,2,0,0     # subject 1: cause A (CMT 2) at t=12.25 …
1,12.25,0,0,0,3,0,0     #            … censored for cause B (CMT 3) at the same time
2,9.38,0,0,0,2,0,0      # subject 2: reached the horizon with no event —
2,9.38,0,0,0,3,0,0      #            right-censored for both causes

See examples/tte_competing_risks.ferx and data/tte_competing_risks.csv.

Prediction — cumulative incidence. predict_survival() reports, per cause, the cause-specific cumulative incidence F_k(t) = ∫₀ᵗ h_k(u)·S_all(u) du alongside the all-cause survival S_all(t) = exp(−Σ_j H_j(t)). In the presence of competing risks the naive 1 − exp(−H_k) over-states a cause’s incidence (it ignores that the other causes remove subjects from risk); the reported CIF is the correct quantity, and across all causes Σ_k F_k(t) + S_all(t) = 1 exactly. For a single endpoint it reduces to F(t) = 1 − S(t).

Estimation note. The cause-specific rates are well-identified and recover tightly. The shared frailty variance ω² is weakly identified for TTE and reads high under FOCEI-Laplace at higher event rates (the same approximation limit as single-endpoint frailty; see below) — prefer SAEM/IMP for the variance component. Subdistribution-hazard (Fine–Gray) modelling is not currently provided; cause-specific hazard covers the standard pharmacometric use case.

Validation of the cumulative incidence. The CIF is anchored at three levels:

  1. Closed form (constant hazards). For two exponential causes the actuarial allocation is exact: F_k(t) = (λ_k / Σλ)·(1 − e^{−Σλ·t}). Asserted to machine precision in cif_two_cause_exponential_closed_form (src/survival/mod.rs).

  2. Independent quadrature (time-varying hazards). For a Weibull + Gompertz pair the per-cause F_k(t) from cif_curves matches a high-resolution composite-Simpson integration of the defining integral F_k(t) = ∫₀ᵗ h_k(u)·S_all(u) du — a different numerical scheme — to within 1e-3 (cif_time_varying_matches_independent_quadrature).

  3. Nonparametric reference (R cmprsk::cuminc). Simulating the example’s two exponential causes (λ_A = 0.10, λ_B = 0.06, administrative censoring at t = 14, N = 200 000) and estimating the cumulative incidence with the license-free cmprsk::cuminc reproduces ferx’s analytic CIF (max absolute difference ≈ 2·10⁻³, i.e. Monte-Carlo noise):

    t cause A — cmprsk cause A — ferx cause B — cmprsk cause B — ferx
    2 0.1708 0.1712 0.1032 0.1027
    6 0.3848 0.3857 0.2324 0.2314
    10 0.4970 0.4988 0.3009 0.2993
    14 0.5565 0.5585 0.3370 0.3351 
    library(cmprsk)
set.seed(20260625)
lamA <- 0.10; lamB <- 0.06; tau <- 14; N <- 2e5
tA <- rexp(N, lamA); tB <- rexp(N, lamB)
tmin <- pmin(tA, tB); cause <- ifelse(tA < tB, 1L, 2L)
status <- ifelse(tmin <= tau, cause, 0L); ftime <- pmin(tmin, tau)
timepoints(cuminc(ftime, status, cencode = 0L), c(2, 6, 10, 14))$est
# ferx analytic: F_k(t) = (lam_k/(lamA+lamB)) * (1 - exp(-(lamA+lamB)*t))

Simulation and VPC. simulate() draws the competing causes with the same earliest-cause-wins rule: a latent time is drawn per cause, the earliest is the observed event, and every other cause is right-censored at that time. For a visual predictive check the simulation needs an administrative censoring horizon that is decoupled from the observed event times — otherwise a re-simulated subject that already carries an event re-draws unbounded, biasing the simulated event-time distribution. Supply it with the [simulation] horizon key (or SimulateOptions { horizon, .. } on the library API); when set, it overrides each record’s per-record window so every cause censors at the planned study end:

[simulation]
  n_subjects = 200
  horizon    = 14     # follow each subject to t = 14, censoring survivors there
  seed       = 12345

Running ferx examples/tte_competing_risks.ferx --simulate then generates one right-censored TTE row per cause compartment (CMT 2 and CMT 3) per synthetic subject and overwrites each with the drawn outcome. A TTE [simulation] block requires horizon (there are no continuous times to schedule).

Comparison with reference software

The exponential validation dataset is tests/reference/tte_exponential/tte_exp.csv (100 subjects; data-generating λ=0.1 h⁻¹, ω²=0.25 on log-rate; right-censored at t=24 → 82 events, 18% censored). All tools fit the same file.

Fixed-effects (no random effect) — exact, license-free anchor. A censored exponential MLE is the closed form λ = events/Σtime, so a ferx n_eta=0 fit must match base-R survival::survreg exactly:

Quantity ferx (fixed) survreg agreement
λ (rate) 0.074506 0.074506 exact (<1e‑4)
OFV / −2logLik 589.888 589.888 exact

This confirms ferx’s TTE hazard, censoring, and likelihood constants are correct.

Mixed-effects (frailty) — FOCEI. ferx recovers the data-generating truth under a clean large-N simulation–estimation check (N=2000): λ_pop = 0.099 (−1%), ω² = 0.232 (−7%). The small ω² shortfall is the expected mild FOCEI-Laplace bias for TTE (it grows at low event rates; SAEM/IMP reduce it). On the 100-subject file itself:

Parameter True ferx FOCEI nlmixr2 FOCEI NONMEM (LAPLACIAN)
TVLAMBDA (rate) 0.10 0.0768 0.0770 0.0771
ω²(ETA_LAMBDA) 0.25 0.290 0.293 0.286
OFV (−2LL) — 588.93 588.94 588.93

ferx and nlmixr2 FOCEI agree to ~3 digits on all parameters and the −2LL for all three families (Weibull and Gompertz likewise — see tests/reference/tte_*/expected.md), including the nonlinear-frailty Weibull omega^2: ferx 0.176, nlmixr2 0.173, NONMEM LAPLACIAN 0.175. Earlier ferx read 0.204 there — not a likelihood difference but a derivative-free outer-optimizer artifact: on the near-flat ω² ridge (the whole 0.175→0.204 span is <0.01 OFV) BOBYQA’s ftol_rel default stopped it short of ferx’s own profile minimum, which already sat at ~0.175. Tightening the TTE outer ftol_rel to 1e-8 (#469; auto-applied to pure-TTE models, whose hazard objective is exact) lands it on the consensus. This is separate from the FOCEI-Laplace method bias above (the large-N Weibull SSE still over-estimates the shape frailty, ~0.34 vs a true 0.20 — that is #440, unchanged by the optimizer fix).

NONMEM reports THETA(1)=log λ; compare exp(THETA(1)) to ferx TVLAMBDA. The NONMEM frailty column is from nonmem_frailty.ctl (CONDITIONAL LAPLACIAN INTERACTION). Two NONMEM gotchas for F_FLAG=1 likelihood models: a dummy $SIGMA 1 FIX is required and must be referenced (DUMMY = EPS(1)) or NM-TRAN infers single-subject data and rejects CONDITIONAL; and NUMERICAL SLOW are required with LAPLACIAN INTERACTION. All three tools agree on the linear-frailty exponential (λ within 0.4%, ω² within ~2%, identical −2LL). Gompertz (fixed-effects) also matches NONMEM to OFV 3011.12. nlmixr2’s OFV constants differ (§11.3) — compare parameters.

nlmixr2 equivalent

m <- function() {
  ini({
    TVLAMBDA <- 0.05
    eta.lambda ~ 0.09
  })
  model({
    lambda <- TVLAMBDA * exp(eta.lambda)
    ll(time.to.event) ~ event * log(lambda) - lambda * time.to.event
  })
}
fit <- nlmixr(m, data, est = "focei")

NONMEM equivalent

$ESTIMATION METHOD=LAPLACIAN INTERACTION MAXEVAL=500
$THETA (0.001, 0.05, 10)
$OMEGA 0.09

Weibull and Gompertz

The same validation kit exists for Weibull (tests/reference/tte_weibull/) and Gompertz (tests/reference/tte_gompertz/). Highlights:

  • Weibull fixed-effects fit matches survreg exactly (scale 22.18, shape 2.12); the 300-subject Gompertz RCT fixed-effects fit recovers log_alpha, log_gamma and the treatment log_hr essentially exactly, exercising the [event_model] covariate + loghr term.

Choosing an estimator for frailty TTE. When a random effect sits on a nonlinear hazard parameter (e.g. the Weibull shape or Gompertz growth rate, not the exponential rate), FOCEI-Laplace over-estimates the frailty variance omega^2 — a clean large-N simulation–estimation check recovers the structural parameters but inflates omega^2 (Weibull shape: 0.34 vs a true 0.20), while a SAEM fit of the same data reads lower (~0.13). This is an approximation limitation, not a likelihood error (fixed-effects fits match survreg exactly), and it is a property of the method, not of any one tool: a three-tool cross-check on the n=100 Weibull file has ferx FOCEI (ω² 0.176), nlmixr2 FOCEI (0.173) and NONMEM LAPLACIAN (0.175) all in agreement (structural params and −2LL match too, ≤0.5% / identical OFV). (ferx read 0.204 here before #469 — a derivative-free outer-optimizer convergence artifact on the flat ω² ridge, since fixed; it is unrelated to the method bias, which shows in the large-N SSE above.) Prefer SAEM/IMP for frailty TTE on nonlinear parameters, as in the standard pharmacometric guidance; FOCEI on a linear rate (exponential) is near-unbiased and all three tools agree there (ω² 0.286–0.293).

Joint PK-TTE — fit and simulation

The drug-driven joint model is anchored against NONMEM and nlmixr2 on a shared dataset (tests/reference/pktte_joint/, N=120; oral 1-cpt PK + h = H0·exp(BETA·Cc)).

Fit. All three tools recover the PK fixed effects and variance components closely:

Parameter Truth ferx FOCEI nlmixr2 FOCEI NONMEM LAPLACE
CL 1.0 1.039 1.026 1.031
V 10.0 9.879 9.946 9.897
KA 1.0 0.978 0.994 0.980
ω²(CL) 0.09 0.083 0.080 0.081
H0 0.015 0.0145 0.0210 0.0149
BETA 0.25 0.256 0.201 0.251

The baseline hazard H0 and the exposure slope BETA are collinear — H0·exp(BETA·Cc) trades baseline against slope over the narrow single-occasion concentration range. ferx and NONMEM (both truth-started) sit near truth; nlmixr2 lands at the other end of the same flat ridge — all valid optima of an under-determined likelihood, with corr(H0, BETA) = −0.93 (ferx) / −0.91 (NONMEM). So the PK block is the cross-tool agreement check; H0 and BETA are jointly, not individually, identified. (The objective value is not cross-tool comparable here — different F_FLAG likelihood normalization — so compare ΔOFV between nested models, not absolutes.)

Simulation. ferx’s simulate() event-time distribution is validated two ways (tests/reference/pktte_joint_sim/):

  • Goodness-of-fit (in-tree, license-free): simulated event times, probability-integral- transformed against the model’s own survival via an independent integrator, pass a Kolmogorov–Smirnov test against Uniform(0,1) — D = 0.040 vs the 5% critical 0.061 at N=500, ≈ the value expected for an exactly-correct sampler. This is immune to the H0/BETA ridge (no estimation), so it is the decisive sampler check.
  • Cross-tool (NONMEM $SIM): ferx’s Kaplan–Meier survival matches NONMEM’s analytic marginal survival mean_i exp(−CHZ_i(t)) to max |ΔS(t)| = 0.044 across the horizon (N=500).

See also

  • examples/tte_exponential.ferx — minimal worked example
  • data/tte_exponential.csv — small dataset (30 subjects) used by the smoke tests
  • tests/reference/tte_exponential/ — 100-subject validation dataset (tte_exp.csv), reference scripts (simulate.R, survreg.R, nlmixr2.R, nonmem.ctl), README.md (NONMEM/nlmixr2 hand-off) and expected.md (filled comparison)
  • examples/pktte_joint.ferx — joint PK-TTE worked example (drug-driven ODE hazard)
  • tests/reference/pktte_joint/ — joint PK-TTE fit anchor (ferx / NONMEM / nlmixr2)
  • tests/reference/pktte_joint_sim/ — joint PK-TTE simulation anchor (ferx vs NONMEM $SIM: make_template.py, sim.ctl, compare.R, expected.md)
  • tests/tte_smoke.rs — Tier-2 parse and short-run smoke tests
  • tests/tte_convergence.rs — Tier-3 convergence + SSE tests (--features survival,slow-tests)
  • plans/tte-survival-markov.md — full multi-phase roadmap