Skip to content

Commit f993ed2

Browse files
fix(nca): filter mixed-endpoint data + data-driven fallback before defaults
Closes the second class of "Suite-C Phase-1 hits per-fit timeout" failures (the first class — FOCEI eta-drift on NMID-column data — was fixed in c721c55). Cause: NCAEstimator was not applying the same DVID PK-row filter the runner-side adapter uses, so warfarin's ``DVID="cp"`` PK rows + ``DVID="pca"`` PD rows produced terminal- slope regressions on a mixed-scale population. 25/26 subjects failed QC, the >50% exclusion gate fired, the estimator fell through to ``_default_estimates()`` (CL=5/V=70/ka=1), and SAEM+FOCEI were asked to traverse 4 OoMs (literature warfarin CL=0.106) inside the 600s per-fit budget — which they could not. Two changes: - ``apmode.data.adapters`` promotes ``_PK_DVID_ALLOWLIST`` to a public ``PK_DVID_ALLOWLIST`` constant (back-compat alias kept) so the adapter and the NCA estimator share a single source of truth for "what counts as a PK observation row". Drift between the two sets is now structurally impossible. - ``NCAEstimator.__init__`` filters its observation slice through ``PK_DVID_ALLOWLIST`` when DVID is present (fail-open when DVID is absent or the filter would empty the frame). On the warfarin fold02 train subset this recovers per-subject NCA from "25/26 excluded → defaults" to "12/26 excluded → CL=0.17, V=8.5, ka=0.11" — within 2x of literature on every parameter. - ``NCAEstimator._apply_fallback`` gains a data-driven layer in front of ``_default_estimates``. When per-subject NCA still fails QC and no dataset_card prior is supplied, it derives population estimates from observed ``Cmax_geomean`` / ``AUC_observed_geomean`` / ``Tmax_geomean`` with log-normal floors/caps ([1,1000] L, [0.01,500] L/h, [0.05,12] /h). The unit-scale heuristic is reused so a mg-dose / ng-mL CSV still gets the 1000x correction. ``fallback_source`` cascade is now nca → dataset_card → data_driven → defaults; defaults remains as the absolute last resort. End-to-end: - Theo + warfarin + mavoglurant nca regression: theo unchanged (CL=2.97, V=32.6, ka=0.95 vs literature 2.83/32.0/1.5), warfarin fold02 recovered (was: defaults 5/70/1; now nca 0.17/8.5/0.11; literature 0.106/10.9/0.888). - 2408 tests pass. Two new test classes pin the contract: TestNCADVIDFilter (filter recovers NCA on PK+PD synthetic frame; no-op on DVID-absent frame) and TestDataDrivenFallback (sparse data triggers data-driven over hard-coded defaults; dataset-card prior still wins over data-driven).
1 parent c721c55 commit f993ed2

4 files changed

Lines changed: 451 additions & 9 deletions

File tree

CHANGELOG.md

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,48 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
77

88
## [Unreleased]
99

10+
### Fixed — `NCAEstimator` filters mixed-endpoint datasets and uses a data-driven fallback when QC fails
11+
12+
A second, independent class of "Phase-1 hits per-fit timeout" failures —
13+
caused not by the FOCEI eta-drift loop (closed in the previous entry)
14+
but by initial-estimate quality — is closed here.
15+
16+
- **Mixed-endpoint datasets pass through the same `DVID` PK-row filter
17+
the runner-side adapter already applies.** Warfarin's canonical
18+
NONMEM CSV interleaves `DVID="cp"` PK rows with `DVID="pca"`
19+
prothrombin-complex-activity PD rows. `NCAEstimator.__init__` now
20+
filters its observation slice via
21+
`apmode.data.adapters.PK_DVID_ALLOWLIST` (the constant was promoted
22+
from `_PK_DVID_ALLOWLIST` to a public name so the NCA module and the
23+
runner-side adapter share one source of truth). Before the filter,
24+
the per-subject terminal-slope regression saw PK and PD values mixed
25+
together, lambda_z fits collapsed for 25/26 subjects, and the
26+
estimator fell through to `_default_estimates()` (CL=5/V=70/ka=1) —
27+
values 47x off for warfarin's literature CL=0.106. After the filter,
28+
warfarin fold02 NCA returns CL=0.17 / V=8.5 / ka=0.11 (within 2x of
29+
literature on every parameter; `fallback_source="nca"`). Fail-open:
30+
if `DVID` is absent (theo, mavoglurant) or the allowlist would empty
31+
the frame (custom DVID schemes), all observation rows are kept and
32+
behaviour is unchanged.
33+
- **New data-driven fallback layer in `NCAEstimator._apply_fallback`
34+
precedes the hard-coded defaults.** When per-subject NCA QC still
35+
fails *and* no `fallback_estimates` (dataset-card prior) is
36+
available, the estimator now derives initial estimates from the
37+
observed Cmax/AUC directly: `V = Dose_geo / Cmax_geo`,
38+
`CL = Dose_geo / AUC_obs_geo`, `ka = 2.5 / Tmax_geo`, with
39+
log-normal floors/caps `[1, 1000] L`, `[0.01, 500] L/h`,
40+
`[0.05, 12] 1/h`. The unit-scale heuristic (mg-dose / ng-mL detector)
41+
is reused from the NCA happy-path so a 1000x mass-units mismatch
42+
still gets corrected. The textual `fallback_source` cascade is now
43+
`nca → dataset_card → data_driven → defaults`. The hard-coded
44+
`_default_estimates()` remains as the final last-resort.
45+
- **Tests.** Two new test classes in `tests/unit/test_initial_estimates
46+
.py`: `TestNCADVIDFilter` (PK+PD synthetic frame recovers via filter;
47+
no-op when DVID is absent) and `TestDataDrivenFallback` (sparse
48+
2-obs-per-subject frame triggers data-driven over hard-coded defaults;
49+
dataset-card prior still wins over data-driven). Full non-live sweep
50+
stays green at 2408 tests.
51+
1052
### Fixed — `Nlmixr2Runner` pre-adapts the on-disk CSV; `r/harness.R` rename is the safety net
1153

1254
Two interlocking layers close a class of indefinite-hang bugs that the

src/apmode/data/adapters.py

Lines changed: 13 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,19 @@
2424
# PK concentration. Mirrors ``apmode.data.profiler._PK_DVIDS`` but is kept
2525
# in-sync manually because ``adapters`` must not import from ``profiler``
2626
# (profiler depends on adapters transitively via bundle artifacts).
27-
_PK_DVID_ALLOWLIST: frozenset[str] = frozenset({"1", "conc", "concentration", "cp"})
27+
#
28+
# Exported (no leading underscore) so ``apmode.data.initial_estimates``
29+
# can apply the same PK-row filter inside its NCA estimator. Without the
30+
# shared filter, NCA on mixed-endpoint datasets (e.g. warfarin's
31+
# ``DVID="cp"`` PK + ``DVID="pca"`` PD rows) sees PD concentrations as
32+
# PK noise, the per-subject lambda_z fits collapse for >50% of subjects,
33+
# and the estimator falls all the way to ``_default_estimates()``.
34+
PK_DVID_ALLOWLIST: frozenset[str] = frozenset({"1", "conc", "concentration", "cp"})
35+
36+
# Backwards-compatible alias for any in-tree consumer that grew an
37+
# import on the underscored name before it was promoted. New call sites
38+
# should use the public name.
39+
_PK_DVID_ALLOWLIST = PK_DVID_ALLOWLIST
2840

2941
# Canonical PK columns — never remap these even if their dtype looks
3042
# categorical. Everything else is treated as a candidate covariate.

src/apmode/data/initial_estimates.py

Lines changed: 184 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@
3737
InitialEstimates,
3838
NCASubjectDiagnostic,
3939
)
40+
from apmode.data.adapters import PK_DVID_ALLOWLIST
4041

4142
_logger = logging.getLogger(__name__)
4243

@@ -155,11 +156,53 @@ def __init__(
155156
"""
156157
self._df = df
157158
self._manifest = manifest
158-
self._obs: pd.DataFrame = cast("pd.DataFrame", df[df["EVID"] == 0].copy())
159+
obs: pd.DataFrame = cast("pd.DataFrame", df[df["EVID"] == 0].copy())
160+
# Filter mixed-endpoint datasets to PK rows only. Warfarin's
161+
# canonical NONMEM-style CSV interleaves DVID="cp" PK rows with
162+
# DVID="pca" prothrombin-complex-activity PD rows; before this
163+
# filter, the per-subject terminal-slope regression saw PK and
164+
# PD values mixed together, lambda_z fits collapsed for 25/26
165+
# subjects, and the estimator fell through to
166+
# ``_default_estimates()`` (CL=5/V=70/ka=1) — values 47x off
167+
# for warfarin, leaving SAEM+FOCEI to traverse 4 OoMs in the
168+
# outer loop. The allowlist is shared with
169+
# ``apmode.data.adapters.to_nlmixr2_format`` so the runner-side
170+
# nlmixr2-ready CSV and the NCA estimator agree on what
171+
# constitutes a PK observation. Fail-open: if DVID is absent
172+
# (theo, mavoglurant) or the allowlist filter would empty the
173+
# frame (custom DVID schemes from a research dataset), keep all
174+
# observation rows so behaviour is unchanged on those datasets.
175+
if "DVID" in obs.columns:
176+
dvid_str = obs["DVID"].astype(str).str.strip().str.lower()
177+
keep = dvid_str.isin(PK_DVID_ALLOWLIST)
178+
n_filtered = int((~keep).sum())
179+
if keep.any():
180+
obs = cast("pd.DataFrame", obs.loc[keep].reset_index(drop=True))
181+
if n_filtered > 0:
182+
_logger.info(
183+
"nca_dvid_filter_applied",
184+
extra={
185+
"kept_rows": int(keep.sum()),
186+
"dropped_rows": n_filtered,
187+
"allowlist": sorted(PK_DVID_ALLOWLIST),
188+
},
189+
)
190+
else:
191+
_logger.warning(
192+
"nca_dvid_filter_would_empty_frame_keeping_all_rows",
193+
extra={
194+
"n_obs": len(obs),
195+
"allowlist": sorted(PK_DVID_ALLOWLIST),
196+
"observed_dvid_values": sorted(
197+
{str(v).strip().lower() for v in obs["DVID"].dropna().unique()}
198+
),
199+
},
200+
)
201+
self._obs: pd.DataFrame = obs
159202
self._doses: pd.DataFrame = cast("pd.DataFrame", df[df["EVID"] == 1].copy())
160203
self._fallback_estimates = fallback_estimates
161204
self.diagnostics: list[NCASubjectDiagnostic] = []
162-
self.fallback_source: str = "nca" # nca | dataset_card | defaults
205+
self.fallback_source: str = "nca" # nca | dataset_card | defaults | data_driven
163206

164207
def estimate_per_subject(self) -> dict[str, float]:
165208
"""Derive population-median NCA estimates from per-subject NCA.
@@ -370,11 +413,22 @@ def _nca_for_subject(self, subj: object) -> NCAResult:
370413
return result
371414

372415
def _apply_fallback(self, excluded_fraction: float) -> dict[str, float]:
373-
"""Populate initial estimates from a literature prior or conservative defaults.
374-
375-
The textual source (``"dataset_card"`` vs ``"defaults"``) is tracked on
376-
``self.fallback_source``. Underscore-prefixed dict keys carry numeric
377-
metadata only.
416+
"""Populate initial estimates from a literature prior, data-driven
417+
heuristic, or conservative defaults — in that order of preference.
418+
419+
The textual source (``"dataset_card"`` vs ``"data_driven"`` vs
420+
``"defaults"``) is tracked on ``self.fallback_source``.
421+
Underscore-prefixed dict keys carry numeric metadata only.
422+
423+
Cascade rationale: a dataset card with explicit literature
424+
priors is always the highest-quality fallback. Absent that, the
425+
observed data itself carries usable signal — population
426+
``Dose / Cmax_geomean`` recovers V to within ~10% of literature
427+
on warfarin even when the per-subject lambda_z fits collapse,
428+
which is much better than the previous CL=5/V=70/ka=1 hard-
429+
coded defaults that were 47x off for warfarin (and triggered
430+
the SAEM+FOCEI 4-OoM traversal that exceeded the 600s per-fit
431+
budget on Suite-C Phase-1).
378432
379433
#27: when the dataset card does not carry a ``ka``, the rc8
380434
path silently defaulted to 1.0 /h - a 10-to-100x warm-start error
@@ -398,11 +452,24 @@ def _apply_fallback(self, excluded_fraction: float) -> dict[str, float]:
398452
est["_ka_defaulted"] = 1.0
399453
est["_excluded_fraction"] = round(excluded_fraction, 4)
400454
return est
455+
data_driven = self._data_driven_fallback()
456+
if data_driven is not None:
457+
self.fallback_source = "data_driven"
458+
data_driven["_excluded_fraction"] = round(excluded_fraction, 4)
459+
_logger.info(
460+
"initial_estimates_data_driven_fallback",
461+
extra={
462+
"reason": "per-subject NCA QC failed; derived from population Cmax/AUC",
463+
"excluded_fraction": excluded_fraction,
464+
**{k: v for k, v in data_driven.items() if not k.startswith("_")},
465+
},
466+
)
467+
return data_driven
401468
self.fallback_source = "defaults"
402469
_logger.warning(
403470
"initial_estimates_using_defaults",
404471
extra={
405-
"reason": "no NCA and no dataset_card prior available",
472+
"reason": "no NCA, no dataset_card prior, and data-driven heuristic failed",
406473
"excluded_fraction": excluded_fraction,
407474
},
408475
)
@@ -413,6 +480,115 @@ def _apply_fallback(self, excluded_fraction: float) -> dict[str, float]:
413480
est["_excluded_fraction"] = round(excluded_fraction, 4)
414481
return est
415482

483+
def _data_driven_fallback(self) -> dict[str, float] | None:
484+
"""Best-effort initial estimates from observed Cmax/AUC, no QC gates.
485+
486+
This bypasses the per-subject lambda_z + extrapolation gates
487+
that ``estimate_per_subject`` enforces. It is intentionally
488+
loose: the goal is order-of-magnitude correctness so SAEM/FOCEI
489+
starts inside its convergence basin, not BE-grade NCA.
490+
491+
Per-subject:
492+
- ``Cmax_obs`` = max DV
493+
- ``AUC_obs`` = trapezoid AUC over the observed window
494+
(no extrapolation; observed-window-only)
495+
- ``Tmax_obs`` = TIME of Cmax_obs (skipped if 0 — that
496+
row is the dose itself or first sample
497+
before absorption resolved)
498+
499+
Population aggregate uses the geometric mean (PK params are
500+
log-normally distributed; the geo-mean is the standard
501+
central-tendency measure for pop-PK initial estimates).
502+
503+
V = Dose_geomean / Cmax_geomean [floor 1, cap 1000 L]
504+
CL = Dose_geomean / AUC_obs_geomean [floor 0.01, cap 500 L/h]
505+
ka = 2.5 / Tmax_geomean [floor 0.05, cap 12 1/h]
506+
507+
Returns ``None`` if any input is empty or non-positive
508+
(e.g. data has no observations, only zero-DV samples, or only
509+
a single time point per subject so trapezoid AUC is undefined).
510+
``ka`` falls back to 1.0 with a flag if Tmax cannot be derived
511+
(all subject Tmax==0, IV-bolus-like data).
512+
"""
513+
if self._obs.empty or self._doses.empty:
514+
return None
515+
516+
per_subj_dose = cast("pd.Series", self._doses.groupby("NMID")["AMT"].sum().astype(float))
517+
positive_doses = per_subj_dose[per_subj_dose > 0]
518+
if positive_doses.empty:
519+
return None
520+
# Geo-mean dose: typical study has a single typical AMT, so
521+
# this collapses to that AMT. Robust to mixed-dose studies too.
522+
dose_geo = float(np.exp(np.log(positive_doses.to_numpy()).mean()))
523+
524+
cmax_per_subj: list[float] = []
525+
auc_per_subj: list[float] = []
526+
tmax_per_subj: list[float] = []
527+
for subj, sdf in self._obs.groupby("NMID"):
528+
sdf_sorted = sdf.sort_values("TIME")
529+
t = sdf_sorted["TIME"].to_numpy(dtype=float)
530+
c = sdf_sorted["DV"].to_numpy(dtype=float)
531+
pos = c > 0
532+
if pos.sum() < 2:
533+
continue
534+
t_pos = t[pos]
535+
c_pos = c[pos]
536+
cmax = float(c_pos.max())
537+
tmax_idx = int(np.argmax(c_pos))
538+
tmax = float(t_pos[tmax_idx])
539+
# Linear trapezoid over observed window only — no
540+
# extrapolation. Conservative: under-estimates AUC for fast
541+
# eliminators with truncated terminal phase, which biases
542+
# CL slightly high — that direction is benign for
543+
# SAEM/FOCEI numerical stability (avoids near-zero CL
544+
# divisions in the inner ODE solve).
545+
auc = float(np.trapezoid(c_pos, t_pos))
546+
if cmax <= 0 or auc <= 0:
547+
continue
548+
_ = subj
549+
cmax_per_subj.append(cmax)
550+
auc_per_subj.append(auc)
551+
if tmax > 0:
552+
tmax_per_subj.append(tmax)
553+
554+
if not cmax_per_subj or not auc_per_subj:
555+
return None
556+
557+
cmax_geo = float(np.exp(np.log(np.asarray(cmax_per_subj)).mean()))
558+
auc_geo = float(np.exp(np.log(np.asarray(auc_per_subj)).mean()))
559+
if cmax_geo <= 0 or auc_geo <= 0:
560+
return None
561+
562+
v_raw = dose_geo / cmax_geo
563+
cl_raw = dose_geo / auc_geo
564+
v = float(min(max(v_raw, 1.0), 1000.0))
565+
cl = float(min(max(cl_raw, 0.01), 500.0))
566+
567+
ka_defaulted = False
568+
if tmax_per_subj:
569+
tmax_geo = float(np.exp(np.log(np.asarray(tmax_per_subj)).mean()))
570+
ka_raw = 2.5 / tmax_geo
571+
ka = float(min(max(ka_raw, 0.05), 12.0))
572+
else:
573+
# All Tmax == 0 (e.g. IV bolus, or first-sample-Cmax data).
574+
# 1.0 /h is the historical default; flag for audit.
575+
ka = 1.0
576+
ka_defaulted = True
577+
578+
out: dict[str, float] = {"CL": cl, "V": v, "ka": ka}
579+
if ka_defaulted:
580+
out["_ka_defaulted"] = 1.0
581+
# Apply the same unit-scale heuristic the NCA happy-path uses
582+
# so a mg-dose / ng-mL CSV does not silently produce CL/V that
583+
# are 1000x off (the heuristic detects the 1e-9-vs-1e-3 mismatch
584+
# from the mass-units check in adapters.detect_unit_scale).
585+
scale, _reason = _detect_unit_scale_factor(self._doses, self._obs, cl_raw)
586+
if scale != 1.0:
587+
out["CL"] = float(min(max(cl_raw * scale, 0.01), 500.0))
588+
out["V"] = float(min(max(v_raw * scale, 1.0), 1000.0))
589+
out["_unit_scale_applied"] = scale
590+
return out
591+
416592

417593
# ---------------------------------------------------------------------------
418594
# Warm-start + bundle construction

0 commit comments

Comments
 (0)