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 hazard may not depend on an inter-occasion (IOV/kappa) random effect — whether by referencing an [individual_parameters] value that uses one, or by naming a kappa directly in the hazard expression — because it is evaluated once per subject, with no occasion context; either 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

Covariates on a hazard must be constant within each subject. Time-varying covariates on a hazard are not yet supported. A covariate referenced by a hazard — in scale / shape / alpha / gamma / loghr, or in a drug-driven hazard = expression — would otherwise be read at its baseline value (the analytic hazard families take no time argument, and the joint PK-TTE ODE hazard integrates with the PK parameters frozen at t = 0), silently fitting the wrong hazard. fit() therefore rejects such a model up front (and predict() / simulate() error), naming the covariate and the subject (#741). A time-varying covariate used only elsewhere — e.g. on the PK parameters of a joint model whose hazard does not reference it — is unaffected.

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

Repeated events (RTTE)

A repeated TTE model (type = rtte in [event_model]) lets one subject experience the same event several times. Add type = rtte and give each subject multiple DV = 1 rows in time order, ending in a DV = 0 administrative-censoring row:

[parameters]
  theta TVLAMBDA(0.15, 0.001, 10.0)
  omega ETA_LAMBDA ~ 0.09

[event_model]
  cmt    = 2
  type   = rtte
  family = exponential
  scale  = TVLAMBDA * exp(ETA_LAMBDA)   # frailty shared across a subject's events

[fit_options]
  method = saem

clock = forward (the default, Andersen–Gill total time) integrates the cumulative hazard once across the whole window (Σ_k log h(t_k) − H(T)); clock = reset (gap time) restarts the hazard clock at each event (Σ_k log h(Δ_k) − Σ_k H(Δ_k) over inter-event gaps). For a time-homogeneous hazard the two coincide. See Repeated events for the DSL details.

Both clocks accept left truncation (a TENTRY column; see above). Under clock = forward the cumulative hazard integrates from H(TENTRY). Under clock = reset the first sojourn is measured from the time origin and conditioned on survival to entry — H(t₁) − H(TENTRY) − log h(t₁) — after which the renewal clock takes over for later gaps. This is the same delayed-entry conditioning across every survival endpoint: TENTRY always means “condition on survival past entry”, never “restart the clock at entry” (the two would differ for a time-varying hazard, and only the former is a standard, anchorable left-truncation model). Reset left truncation assumes the time origin t = 0 is the renewal origin of the first sojourn (the subject was event-free and at risk from 0); it coincides with clock-forward for a memoryless (exponential) hazard.

simulate() accepts left truncation for repeated events too (#740): the stream is drawn on the same time origin conditioned on survival to TENTRY — clock-forward seeds its conditioning clock at TENTRY, clock-reset draws its first sojourn conditional on survival to entry and renews from 0 thereafter — so a simulated left-truncated stream refits under the same convention it was drawn from. TENTRY = 0 is byte-identical to the non-truncated draw. (The simulated event rows carry the drawn times; as with single-event TTE simulation, TENTRY itself is a property of the input design, not re-emitted onto the output rows.)

WarningPrefer SAEM/IMP for RTTE

A Laplace-based fit (focei/foce) can underestimate the frailty ω² for RTTE at low event rates (Karlsson et al. 2009: −91 % to −96 % below a 43 % event rate). ferx warns and recommends method = saem or method = imp. For an event-rich design the methods coincide (see the anchor below).

Validation — constant-hazard RTTE

Dataset tests/reference/rtte_exponential/rtte_exp.csv (100 subjects; data-generating TVLAMBDA = 0.15, ω² = 0.09; censored at t = 20 → 305 events, 3.05/subject). Two of the anchors are exact and license-free — for a constant hazard the RTTE likelihood factorises, so the per-subject event counts are Poisson(λᵢ·T) and a Poisson-lognormal GLMM on the counts is the exact MLE (count_equivalence.R).

Method TVLAMBDA ω²(log rate) OFV
Analytic MLE D/Σt (fixed) 0.15250 — 1757.16
ferx fixed-effects 0.15250 — 1757.16
Poisson-LN GLMM (exact mixed MLE) 0.14062 0.16455 —
NONMEM LAPLACE 0.14062 0.16450 1748.49
ferx FOCEI 0.14060 0.16448 1748.49
ferx SAEM (seed 31) 0.14033 0.16993 1748.50

ferx FOCEI reproduces NONMEM LAPLACE to 4–5 significant figures and the OFV to the decimal (1748.49), and both land on the exact Poisson-LN GLMM MLE. The frailty ω² is well-identified here because the design is event-rich, so FOCEI, SAEM, NONMEM and the GLMM coincide — the Karlsson Laplace bias is a low-event-rate effect. Full write-up in tests/reference/rtte_exponential/expected.md.

Simulating repeated events

simulate() draws a recurrent event stream per subject up to the administrative [simulation] horizon (analytic hazard families; RTTE with a drug-driven ODE hazard is a later slice). Under clock = forward each next event is drawn conditional on survival past the previous one — the closed-form inverse of H(t) − H(t_{k−1}) = −log u — so a time-varying hazard keeps accumulating across the window; under clock = reset each inter-event gap is a fresh draw from zero. Every subject emits its events (DV = 1) followed by one administrative censor (DV = 0) at the horizon, so simulated data round-trips straight back through the data reader. A finite horizon is required (the stream has no implicit end); left truncation, multiple RTTE causes, an RTTE cause mixed with a competing single-event cause, and EVID=3/4 resets are rejected with a clear error rather than silently mis-sampled.

Unlike fitting, simulate() regenerates each subject’s stream from t = 0 out of the model and drawn frailty — it never reads the input event rows. The fit-boundary input-shape checks (out-of-order, interval-censored, or non-finite event times; see Repeated events) therefore do not apply to simulation: such rows are simply not consulted, not rejected.

NotePer-subject simulation diagnostics

A degenerate draw is handled per subject, not by aborting the run: a non-positive / non-finite analytic hazard rate produces no event (the subject is censored), and an RTTE hazard so extreme it would fire more than a million times over the window has its stream skipped (also censored). Rather than let these look like ordinary censoring, the simulate_with_options_diag API returns a SimulationOutput { results, warnings } naming the affected subjects (W_TTE_DEGENERATE_HAZARD / W_RTTE_DEGENERATE); ferx … --simulate echoes them alongside the fit warnings. The plain simulate() / simulate_with_seed() / simulate_with_options() entry points apply the same per-subject handling but return only the rows.

[simulation]
  n_subjects = 200
  horizon    = 20
  seed       = 33

The sampler is validated four ways (all in tests/tte_convergence.rs except the last):

  • Analytic count law (Tier-1, exact). For a constant hazard the stream is a homogeneous Poisson process, so N(T) ~ Poisson(λT); for clock-forward Weibull it is non-homogeneous with N(T) ~ Poisson(H(T)). The simulated counts match these exact laws.

  • Goodness-of-fit (PIT/KS, estimation-free). By the time-rescaling theorem the rescaled inter-event increments H(t_k) − H(t_{k−1}) are iid Exp(1) under correct sampling; a Kolmogorov–Smirnov test of V_k = exp(−ΔH_k) against Uniform(0,1) passes (D ≈ 0.018), confirming the conditional forward draws reproduce H(t).

  • Frailty variance (method of moments, estimation-free). The simulated event counts imply Var(log λ_i) = ω²; at a well-identified design (≈31 events/subject) the read is ω² = 0.092 against the injected 0.09.

  • Cross-tool anchor (NONMEM, external). ferx simulated tests/reference/rtte_exponential_sim/rtte_sim.csv (300 subjects; truth TVLAMBDA = 0.15, ω² = 0.09); ferx and NONMEM then both fit it:

    Engine Method TVLAMBDA ω²(log rate) OFV
    NONMEM LAPLACE 0.15617 0.12801 5557.266
    ferx FOCEI 0.15619 0.12794 5557.266

    An independent engine recovering the same parameters from ferx-simulated data is the external corroboration that the simulator is correct. Both land at ω² ≈ 0.128 (above the data-generating 0.09) because this single 300-subject realization is over-dispersed — a Monte-Carlo property of the drawn dataset, not the sampler, whose true injected variance the method-of-moments check pins at 0.09. Full write-up in tests/reference/rtte_exponential_sim/expected.md.

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