diff --git a/data/validation/clint_uncertainty_gate_2026-07-08.json b/data/validation/clint_uncertainty_gate_2026-07-08.json new file mode 100644 index 0000000..556bf1c --- /dev/null +++ b/data/validation/clint_uncertainty_gate_2026-07-08.json @@ -0,0 +1,188 @@ +{ + "seed": 20260422, + "K": 30, + "n": 107, + "boot_resamples": 10000, + "note": "gate: |rho|>=0.30 with CI excluding 0 = PASS", + "corr": { + "U_boot": { + "r_meta": { + "rho": -0.255, + "ci": [ + -0.423, + -0.069 + ] + }, + "e_eng": { + "rho": -0.006, + "ci": [ + -0.2, + 0.184 + ] + }, + "e_ml": { + "rho": -0.122, + "ci": [ + -0.303, + 0.071 + ] + }, + "diff": { + "rho": 0.119, + "ci": [ + -0.073, + 0.301 + ] + }, + "obs_mag": { + "rho": 0.002, + "ci": [ + -0.182, + 0.188 + ] + } + }, + "U_ad": { + "r_meta": { + "rho": 0.032, + "ci": [ + -0.159, + 0.227 + ] + }, + "e_eng": { + "rho": 0.145, + "ci": [ + -0.059, + 0.342 + ] + }, + "e_ml": { + "rho": 0.005, + "ci": [ + -0.177, + 0.186 + ] + }, + "diff": { + "rho": 0.121, + "ci": [ + -0.069, + 0.311 + ] + }, + "obs_mag": { + "rho": -0.018, + "ci": [ + -0.208, + 0.176 + ] + } + }, + "clint_mean": { + "r_meta": { + "rho": 0.026, + "ci": [ + -0.163, + 0.21 + ] + }, + "e_eng": { + "rho": -0.007, + "ci": [ + -0.181, + 0.171 + ] + }, + "e_ml": { + "rho": 0.012, + "ci": [ + -0.184, + 0.202 + ] + }, + "diff": { + "rho": -0.008, + "ci": [ + -0.199, + 0.178 + ] + }, + "obs_mag": { + "rho": 0.222, + "ci": [ + 0.034, + 0.397 + ] + } + } + }, + "trend": { + "U_boot": { + "r_meta": [ + 0.56, + 0.475, + 0.368, + 0.351 + ], + "e_eng": [ + 0.625, + 0.673, + 0.495, + 0.733 + ], + "diff": [ + 0.068, + 0.18, + 0.049, + 0.321 + ] + }, + "U_ad": { + "r_meta": [ + 0.449, + 0.347, + 0.515, + 0.439 + ], + "e_eng": [ + 0.584, + 0.518, + 0.644, + 0.775 + ], + "diff": [ + 0.108, + 0.083, + 0.108, + 0.316 + ] + }, + "clint_mean": { + "r_meta": [ + 0.489, + 0.357, + 0.463, + 0.441 + ], + "e_eng": [ + 0.649, + 0.588, + 0.68, + 0.607 + ], + "diff": [ + 0.149, + 0.208, + 0.056, + 0.206 + ] + } + }, + "u_boot_by_ad": { + "in_ad_mean": 0.156, + "out_ad_mean": 0.168, + "n_out_ad": 26 + }, + "verdict": "FAIL (null)" +} \ No newline at end of file diff --git a/docs/research/dead-ends.md b/docs/research/dead-ends.md index 4de173d..f15c2f2 100644 --- a/docs/research/dead-ends.md +++ b/docs/research/dead-ends.md @@ -494,6 +494,20 @@ Artifacts: `scripts/probe_liver_zonation.py`, `tests/integration/test_liver_zona --- +### DE-57 — Uncertainty-aware meta: no per-drug CLint epistemic-uncertainty signal exists to gate it (2026-07-08) + +**Date:** 2026-07-08 + +**The last §9 CLint-floor candidate able to sidestep error-cancellation, tested and foreclosed.** The idea does not improve the CLint point estimate (floored at R²≈0.24; every prior improvement tripped error-cancellation, §1/§2). Instead it exploits the floor's *heteroscedasticity*: if we knew **per-drug** which CLint predictions are unreliable, the meta could down-weight the engine track there — a *weighting* change that perturbs no track's residuals, so it cannot trip error-cancellation. **Distinct from DE-55**, which killed divergence/AD/magnitude and the *propagated MC-output* uncertainty (ρ=0.039); the shipped CLint `Distribution` CV is the fixed constant `_CLINT_CV=1.0` (`predict/adme.py`), which is *why* DE-55's propagated signal was flat. This test constructs a genuine per-drug epistemic uncertainty two independent ways: **U_boot** = SD of log10(CLint) across a K=30 XGBoost bootstrap ensemble on the 1,213-compound hepatocyte training set (range [0.069, 0.412], non-degenerate); **U_ad** = 1 − max Tanimoto(Morgan) to that set. + +**Result (gate = Spearman ρ + bootstrap 95% CI, seed 20260422, N=107).** Core null: **ρ(U_boot, e_eng) = −0.006** [−0.20, +0.18] — CLint epistemic uncertainty does not predict engine error, though the engine *is* the CLint-fed track. Actionable null: **ρ(U, e_eng−e_ml) ≈ +0.12** for **both** proxies (U_boot [−0.07, +0.30], U_ad [−0.07, +0.31]) — CI∋0, no reweighting signal. The one CI-significant term, ρ(U_boot, r_meta) = −0.255 [−0.42, −0.07], is sub-threshold **and wrong-signed** (higher CLint uncertainty ⇒ *lower* meta error — the meta already handles uncertain-CLint drugs well). Confound-clean: ρ(U_boot, obs-magnitude) = +0.002. + +**Why it failed — §10.** The binding wall is bioavailability-F blindness, shared across all four tracks and **orthogonal to CLint**; Cmax error is F-driven, so no CLint-derived signal can discriminate it. A perfect CLint-uncertainty router would have nothing to route on. Per §4 the gate precedes integration → the uncertainty-aware meta is **not built**. + +**Telltale if it returns:** "route the meta by how uncertain the CLint prediction is." No per-drug CLint uncertainty predicts Cmax error (ρ≈0 on the actionable target for two independent constructions); the discriminating error is F, not CLint. Artifacts: `scripts/clint_uncertainty_gate.py`, `data/validation/clint_uncertainty_gate_2026-07-08.json`. + +--- + ## 3. When to consult this list - Before writing a design spec for any accuracy improvement. diff --git a/docs/research/diagnosis.md b/docs/research/diagnosis.md index 0d13e8c..8ab486b 100644 --- a/docs/research/diagnosis.md +++ b/docs/research/diagnosis.md @@ -102,7 +102,7 @@ A verified literature review (IMI OrBiTo Part 4, Ahmad 2020/2021, 58 APIs, >200 **Correctness cross-link.** The post-FLUX-1 audit found the canonical well-stirred driving force is blood unbound `fu_b = fup/Rb` (Pang & Han 2019) — independently corroborating that the engine's RBP concentration-basis cluster (clearance sink applies `fup`, not `fup/RBP`) is a deviation from the textbook equation, not merely an internal inconsistency. Spec: `docs/_internal/specs/2026-06-04-rbp-concentration-basis-design.md`. Largely off the `realize_means` headline path (RBP clamped ~[0.5,1.5]); a correctness item, not a benchmark lever. -**Open follow-ups (research in flight 2026-06-04):** whether any *breakthrough* attacks the CLint floor without tripping error-cancellation (per-isoform rCYP direct prediction, kcat×Km×[E] structured decomposition, multi-assay target-denoising, calibrated-CLint-uncertainty→uncertainty-aware meta); and whether conformal/CQR + model-form-discrepancy recalibration fixes the §empirical PI under-coverage (29.9%@90%). Both must clear the end-to-end Cmax-holdout / error-decorrelation gate, not isolated metrics. +**Open follow-ups (research in flight 2026-06-04):** whether any *breakthrough* attacks the CLint floor without tripping error-cancellation (per-isoform rCYP direct prediction, kcat×Km×[E] structured decomposition, multi-assay target-denoising, ~~calibrated-CLint-uncertainty→uncertainty-aware meta~~ — **closed, DE-57**); and whether conformal/CQR + model-form-discrepancy recalibration fixes the §empirical PI under-coverage (29.9%@90%; the PI arm is **closed, DE-55/56**). Both must clear the end-to-end Cmax-holdout / error-decorrelation gate, not isolated metrics. **Update 2026-07-08 (DE-57):** the uncertainty-aware-meta candidate — the one that structurally sidesteps error-cancellation (a weighting signal, not a point-estimate change) — was gate-tested and failed: no per-drug CLint *epistemic* uncertainty (bootstrap-ensemble SD or AD distance) predicts Cmax error (ρ≈0 on `e_eng` and on `e_eng−e_ml`), because the discriminating error is bioavailability-F, orthogonal to CLint (§10). The remaining §9 breakthroughs (rCYP, kcat×Km, multi-assay denoising) all require external enzyme-kinetics data not in the repo. --- diff --git a/scripts/clint_uncertainty_gate.py b/scripts/clint_uncertainty_gate.py new file mode 100644 index 0000000..32711ae --- /dev/null +++ b/scripts/clint_uncertainty_gate.py @@ -0,0 +1,197 @@ +"""DE-57 decorrelation-gate: does per-drug CLint *epistemic* uncertainty predict +per-drug Cmax error — i.e. is there a signal an "uncertainty-aware meta" could +exploit to down-weight the engine track where CLint is unreliable? + +This is the last diagnosis.md §9 CLint-floor candidate that (a) is sourceable from +repo data, (b) does not touch the production headline, and (c) structurally sidesteps +error-cancellation (it is a *weighting* signal, not a point-estimate change). + +Distinct from DE-55, which killed divergence/AD/magnitude and the *propagated +MC-output* uncertainty (rho=0.039). The shipped CLint Distribution CV is a FIXED +constant (_CLINT_CV=1.0 in predict/adme.py) -> carries zero per-drug info, which is +exactly why the propagated signal was flat. Here we CONSTRUCT a genuine per-drug +epistemic uncertainty two independent ways: + U_boot : SD of log10(CLint) across a K-model bootstrap ensemble (model/data + disagreement -> high where the CLint training set is sparse). + U_ad : 1 - max Tanimoto similarity of the drug's Morgan FP to the CLint + training set (applicability-domain distance). + +Gate (DE-55 style): Spearman rho of each signal against per-drug residuals, with +bootstrap 95% CI (10000 resamples, seed 20260422 = project convention). + r_meta = |log10(meta/obs)| (overall meta error) + e_eng = |log10(eng /obs)| (engine error; engine is the CLint-fed track) + diff = e_eng - e_ml (>0 => engine worse than ml => an + uncertainty-aware meta SHOULD shift weight to + ml; this is THE actionable correlation) +Verdict: PASS iff any |rho| >= 0.30 with CI excluding 0. Else FAIL. + +Reproduce: PYTHONPATH=src python scripts/clint_uncertainty_gate.py +Deterministic; uses only repo data. Absolute residuals come from the canonical cache; +U is computed on the local stack, but the gate uses only *rank* correlation, which is +robust to CLint point-prediction stack drift. +""" +from __future__ import annotations + +import csv +import json + +import numpy as np +import xgboost as xgb +from scipy import stats + +from sisyphus.descriptors import compute_features + +SEED = 20260422 +K = 30 +rng = np.random.default_rng(SEED) + +CLINT_TAB = "data/training/clearance_hepatocyte_az.tab" +CACHE = "data/training/4track_holdout_predictions.json" +CPK = "data/reference/clinical_pk.json" +OUT = "data/validation/clint_uncertainty_gate_2026-07-08.json" + + +def feats(smiles: str): + try: + return compute_features(smiles) + except Exception: + return None + + +# ---- 1. CLint training set ---------------------------------------------------- +tr_smi, tr_y = [], [] +with open(CLINT_TAB) as f: + r = csv.reader(f, delimiter="\t") + next(r) + for row in r: + tr_smi.append(row[1].strip('"')) + tr_y.append(float(row[2])) +Xtr, ytt = [], [] +for s, y in zip(tr_smi, tr_y): + fv = feats(s) + if fv is None: + continue + Xtr.append(fv) + ytt.append(np.log10(max(y, 0.1))) # mirror shipped model target transform +Xtr = np.vstack(Xtr) +ytt = np.asarray(ytt) +print(f"CLint train: {Xtr.shape[0]} compounds, {Xtr.shape[1]} feats, " + f"y(log10) range [{ytt.min():.2f},{ytt.max():.2f}]") + +# ---- 2. Holdout 107: residuals (cache) + SMILES (clinical_pk) ----------------- +cache = json.load(open(CACHE)) +drugs = cache["drugs"] +cpk = json.load(open(CPK))["drugs"] +smap = {k.strip().lower(): v["smiles"] for k, v in cpk.items() if v.get("smiles")} + +names, obs, eng, ml, meta, in_ad = [], [], [], [], [], [] +Xho = [] +for d in drugs: + sm = smap.get(d["name"].strip().lower()) + fv = feats(sm) if sm else None + if fv is None: + print(" !! no features for", d["name"]) + continue + names.append(d["name"]) + obs.append(d["obs"]) + eng.append(d["eng"]) + ml.append(d["ml"]) + meta.append(d["meta"]) + in_ad.append(bool(d.get("in_ad", True))) + Xho.append(fv) +Xho = np.vstack(Xho) +obs = np.asarray(obs) +eng = np.asarray(eng) +ml = np.asarray(ml) +meta = np.asarray(meta) +n = len(names) +print(f"Holdout: {n} drugs with features") + +# ---- 3. Bootstrap ensemble -> U_boot ----------------------------------------- +preds = np.zeros((K, n)) +for k in range(K): + idx = rng.integers(0, len(ytt), len(ytt)) + m = xgb.XGBRegressor( + n_estimators=200, max_depth=6, learning_rate=0.05, + subsample=0.8, colsample_bytree=0.8, random_state=SEED + k, n_jobs=4, + ) + m.fit(Xtr[idx], ytt[idx]) + preds[k] = m.predict(Xho) +U_boot = preds.std(axis=0) # per-drug epistemic SD, log10 units +clint_mean = preds.mean(axis=0) # ensemble-mean log10(CLint) +print(f"U_boot range [{U_boot.min():.3f},{U_boot.max():.3f}] " + f"median {np.median(U_boot):.3f}") + +# ---- 4. AD distance -> U_ad --------------------------------------------------- +# Morgan bits = first 2048 features (0/1). Tanimoto on bit sets. +Btr = (Xtr[:, :2048] > 0.5) +Bho = (Xho[:, :2048] > 0.5) +tr_pop = Btr.sum(1).astype(float) +U_ad = np.zeros(n) +for i in range(n): + inter = (Bho[i] & Btr).sum(1).astype(float) + union = Bho[i].sum() + tr_pop - inter + sim = np.where(union > 0, inter / union, 0.0) + U_ad[i] = 1.0 - sim.max() # distance to nearest training neighbour +print(f"U_ad range [{U_ad.min():.3f},{U_ad.max():.3f}] median {np.median(U_ad):.3f}") + +# ---- 5. Residuals ------------------------------------------------------------- +eps = 1e-12 +r_meta = np.abs(np.log10(meta / obs + eps)) +e_eng = np.abs(np.log10(eng / obs + eps)) +e_ml = np.abs(np.log10(ml / obs + eps)) +diff = e_eng - e_ml +mag = np.abs(np.log10(obs + eps)) # confound: obs magnitude + + +def boot_spearman(x, y, nb=10000): + rho = stats.spearmanr(x, y).correlation + idxs = rng.integers(0, len(x), (nb, len(x))) + bs = np.array([stats.spearmanr(x[i], y[i]).correlation for i in idxs]) + lo, hi = np.percentile(bs, [2.5, 97.5]) + return float(rho), float(lo), float(hi) + + +def quartile_trend(u, resid): + q = np.quantile(u, [0.25, 0.5, 0.75]) + bins = np.digitize(u, q) + return [round(float(resid[bins == b].mean()), 3) if (bins == b).any() else None + for b in range(4)] + + +targets = {"r_meta": r_meta, "e_eng": e_eng, "e_ml": e_ml, "diff": diff, "obs_mag": mag} +signals = {"U_boot": U_boot, "U_ad": U_ad, "clint_mean": clint_mean} + +result = {"seed": SEED, "K": K, "n": n, "boot_resamples": 10000, + "note": "gate: |rho|>=0.30 with CI excluding 0 = PASS", "corr": {}, "trend": {}} +print("\n=== Spearman rho (95% CI) ===") +for sname, sig in signals.items(): + result["corr"][sname] = {} + for tname, t in targets.items(): + rho, lo, hi = boot_spearman(sig, t) + result["corr"][sname][tname] = {"rho": round(rho, 3), "ci": [round(lo, 3), round(hi, 3)]} + star = " <-- PASS" if abs(rho) >= 0.30 and (lo > 0 or hi < 0) else "" + print(f" {sname:11s} vs {tname:8s}: rho={rho:+.3f} [{lo:+.3f},{hi:+.3f}]{star}") + result["trend"][sname] = { + t: quartile_trend(sig, targets[t]) for t in ("r_meta", "e_eng", "diff") + } + +# in_ad split of U (relation to DE-55 AD signal) +in_ad = np.asarray(in_ad) +result["u_boot_by_ad"] = { + "in_ad_mean": round(float(U_boot[in_ad].mean()), 3), + "out_ad_mean": round(float(U_boot[~in_ad].mean()), 3) if (~in_ad).any() else None, + "n_out_ad": int((~in_ad).sum()), +} + +# verdict +passed = any( + abs(result["corr"][s][t]["rho"]) >= 0.30 + and (result["corr"][s][t]["ci"][0] > 0 or result["corr"][s][t]["ci"][1] < 0) + for s in ("U_boot", "U_ad") for t in ("r_meta", "e_eng", "diff") +) +result["verdict"] = "PASS" if passed else "FAIL (null)" +print(f"\nVERDICT: {result['verdict']}") + +json.dump(result, open(OUT, "w"), indent=2) +print("wrote", OUT)