add: Kepler score
This commit is contained in:
543
paper3_phase5b_refined.py
Normal file
543
paper3_phase5b_refined.py
Normal file
@@ -0,0 +1,543 @@
|
||||
"""
|
||||
PAPER 3 — PHASE 5b: 5-TERM FORMULA + ORDINAL BEDSIDE SCORE
|
||||
═════════════════════════════════════════════════════════════════════════════
|
||||
Two refinements from Phase 5:
|
||||
|
||||
1. DROP log(1 + ne_auc). Its bootstrap CI crossed zero, and it is
|
||||
collinear with I(ne_at_h24 > 0.08). Simpler 5-term formula.
|
||||
|
||||
2. BEDSIDE SCORE via ORDINAL BINNING (not coefficient rounding).
|
||||
Each of the 5 remaining terms gets mapped to 0–4 points based on
|
||||
clinically meaningful thresholds. SOFA-style integer score:
|
||||
|
||||
Lactate h24 0 if <2.5 | 2 if 2.5–4 | 4 if >4
|
||||
Oliguria (ml/kg) 0 if ≥20 | 1 if 10–20 | 2 if <10
|
||||
NE at h24 0 if ≤0.08 | 3 if >0.08
|
||||
HR deviation 0 if 70–100| 1 if 60–70/100-120 | 2 if <60/>120
|
||||
Pressor MAP/NE 0 if >3000 | 1 if 1000–3000 | 2 if <1000
|
||||
|
||||
Max 13 points. Bedside-ready.
|
||||
|
||||
Repeats the full Phase 5 validation: CV, multi-seed, bootstrap, calibration,
|
||||
subgroups — for BOTH the 5-term continuous formula AND the ordinal score.
|
||||
|
||||
Usage:
|
||||
python paper3_phase5b_refined.py
|
||||
"""
|
||||
|
||||
import json, sys, math, time, random
|
||||
from collections import defaultdict
|
||||
|
||||
BQ_PROJECT = "goddard-gap"
|
||||
DATA_PROJECT = "physionet-data"
|
||||
|
||||
H_SNAPSHOT = 24
|
||||
H_PEAK_NE = 12
|
||||
TRAIN_FRAC = 0.70
|
||||
N_SEEDS = 10
|
||||
N_FOLDS = 5
|
||||
N_BOOTSTRAP = 1000
|
||||
OUT_FILE = "paper3_phase5b_refined.json"
|
||||
|
||||
NE_ITEMID = 221906
|
||||
LACTATE_ID = 50813
|
||||
MAP_ITEMIDS = [220052, 220181, 225312]
|
||||
HR_ITEMIDS = [220045]
|
||||
|
||||
|
||||
def run_bq(sql, label=""):
|
||||
try:
|
||||
from google.cloud import bigquery
|
||||
client = bigquery.Client(project=BQ_PROJECT)
|
||||
t0 = time.time()
|
||||
rows = [dict(r.items()) for r in client.query(sql).result()]
|
||||
print(f" {label:32s} {len(rows):>8,d} rows ({time.time()-t0:.1f}s)")
|
||||
return rows
|
||||
except Exception as e:
|
||||
print(f"[BQ ERROR] {label}: {e}"); return []
|
||||
|
||||
|
||||
# ── Queries (same as Phase 5, SAPS Q4 pre-filtered) ────────────────────────
|
||||
def q_cohort():
|
||||
return f"""
|
||||
WITH weight_first AS (
|
||||
SELECT ce.stay_id, ANY_VALUE(ce.valuenum) AS weight_kg
|
||||
FROM `{DATA_PROJECT}.mimiciv_3_1_icu.chartevents` ce
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_icu.icustays` icu ON ce.stay_id = icu.stay_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_derived.sepsis3` s3 ON s3.stay_id = ce.stay_id
|
||||
WHERE s3.sepsis3 = TRUE AND ce.itemid IN (226512, 224639)
|
||||
AND ce.valuenum BETWEEN 30 AND 300
|
||||
AND ce.charttime BETWEEN icu.intime AND TIMESTAMP_ADD(icu.intime, INTERVAL 24 HOUR)
|
||||
GROUP BY ce.stay_id
|
||||
)
|
||||
SELECT icu.stay_id, icu.subject_id, icu.intime,
|
||||
pat.anchor_age AS age, pat.gender,
|
||||
saps.sapsii, adm.hospital_expire_flag AS died,
|
||||
COALESCE(wf.weight_kg, 75.0) AS weight_kg
|
||||
FROM `{DATA_PROJECT}.mimiciv_3_1_derived.sepsis3` s3
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_icu.icustays` icu ON icu.stay_id = s3.stay_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_hosp.admissions` adm ON adm.hadm_id = icu.hadm_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_hosp.patients` pat ON pat.subject_id = icu.subject_id
|
||||
LEFT JOIN `{DATA_PROJECT}.mimiciv_3_1_derived.sapsii` saps ON saps.stay_id = icu.stay_id
|
||||
LEFT JOIN weight_first wf ON wf.stay_id = icu.stay_id
|
||||
WHERE s3.sepsis3 = TRUE
|
||||
AND TIMESTAMP_DIFF(icu.outtime, icu.intime, HOUR) >= {H_SNAPSHOT}
|
||||
AND saps.sapsii IS NOT NULL AND saps.sapsii >= 48
|
||||
"""
|
||||
|
||||
|
||||
def q_ne():
|
||||
return f"""
|
||||
SELECT ie.stay_id,
|
||||
TIMESTAMP_DIFF(ie.starttime, icu.intime, MINUTE) AS start_min,
|
||||
TIMESTAMP_DIFF(ie.endtime, icu.intime, MINUTE) AS end_min,
|
||||
ie.rate
|
||||
FROM `{DATA_PROJECT}.mimiciv_3_1_icu.inputevents` ie
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_icu.icustays` icu ON icu.stay_id = ie.stay_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_derived.sepsis3` s3 ON s3.stay_id = ie.stay_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_derived.sapsii` saps ON saps.stay_id = icu.stay_id
|
||||
WHERE s3.sepsis3 = TRUE AND saps.sapsii >= 48
|
||||
AND ie.itemid = {NE_ITEMID} AND ie.rate > 0
|
||||
AND ie.starttime BETWEEN icu.intime AND TIMESTAMP_ADD(icu.intime, INTERVAL 30 HOUR)
|
||||
"""
|
||||
|
||||
|
||||
def q_fluid_out():
|
||||
return f"""
|
||||
SELECT oe.stay_id, SUM(oe.value) AS fluid_out_ml
|
||||
FROM `{DATA_PROJECT}.mimiciv_3_1_icu.outputevents` oe
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_icu.icustays` icu ON icu.stay_id = oe.stay_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_derived.sepsis3` s3 ON s3.stay_id = oe.stay_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_derived.sapsii` saps ON saps.stay_id = icu.stay_id
|
||||
WHERE s3.sepsis3 = TRUE AND saps.sapsii >= 48
|
||||
AND oe.value > 0
|
||||
AND oe.charttime BETWEEN icu.intime AND TIMESTAMP_ADD(icu.intime, INTERVAL {H_SNAPSHOT} HOUR)
|
||||
GROUP BY oe.stay_id
|
||||
"""
|
||||
|
||||
|
||||
def q_vitals():
|
||||
ids = ",".join(str(x) for x in MAP_ITEMIDS + HR_ITEMIDS)
|
||||
return f"""
|
||||
SELECT ce.stay_id, ce.itemid, AVG(ce.valuenum) AS val
|
||||
FROM `{DATA_PROJECT}.mimiciv_3_1_icu.chartevents` ce
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_icu.icustays` icu ON icu.stay_id = ce.stay_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_derived.sepsis3` s3 ON s3.stay_id = ce.stay_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_derived.sapsii` saps ON saps.stay_id = icu.stay_id
|
||||
WHERE s3.sepsis3 = TRUE AND saps.sapsii >= 48
|
||||
AND ce.itemid IN ({ids})
|
||||
AND ce.valuenum IS NOT NULL AND ce.valuenum > 0
|
||||
AND ce.charttime BETWEEN TIMESTAMP_ADD(icu.intime, INTERVAL 20 HOUR)
|
||||
AND TIMESTAMP_ADD(icu.intime, INTERVAL 28 HOUR)
|
||||
GROUP BY ce.stay_id, ce.itemid
|
||||
"""
|
||||
|
||||
|
||||
def q_lactate():
|
||||
return f"""
|
||||
SELECT icu.stay_id,
|
||||
TIMESTAMP_DIFF(le.charttime, icu.intime, MINUTE) AS offset_min,
|
||||
le.valuenum AS val
|
||||
FROM `{DATA_PROJECT}.mimiciv_3_1_hosp.labevents` le
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_icu.icustays` icu ON icu.hadm_id = le.hadm_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_derived.sepsis3` s3 ON s3.stay_id = icu.stay_id
|
||||
JOIN `{DATA_PROJECT}.mimiciv_3_1_derived.sapsii` saps ON saps.stay_id = icu.stay_id
|
||||
WHERE s3.sepsis3 = TRUE AND saps.sapsii >= 48
|
||||
AND le.itemid = {LACTATE_ID}
|
||||
AND le.valuenum IS NOT NULL
|
||||
AND le.charttime BETWEEN icu.intime AND TIMESTAMP_ADD(icu.intime, INTERVAL 30 HOUR)
|
||||
"""
|
||||
|
||||
|
||||
# ── Primitives ──────────────────────────────────────────────────────────────
|
||||
def build_primitives(cohort, ne_rows, fout_rows, vital_rows, lac_rows):
|
||||
print(f"\n[3] Building primitives...")
|
||||
ne_by = defaultdict(list)
|
||||
for r in ne_rows: ne_by[r["stay_id"]].append(r)
|
||||
fout = {r["stay_id"]: r["fluid_out_ml"] or 0 for r in fout_rows}
|
||||
|
||||
vital_by = defaultdict(dict)
|
||||
for r in vital_rows:
|
||||
iid = r["itemid"]
|
||||
key = "map" if iid in MAP_ITEMIDS else "hr"
|
||||
if r["stay_id"] is not None:
|
||||
cur = vital_by[r["stay_id"]].get(key)
|
||||
vital_by[r["stay_id"]][key] = r["val"] if cur is None else (cur + r["val"])/2
|
||||
|
||||
lac_by = defaultdict(list)
|
||||
for r in lac_rows: lac_by[r["stay_id"]].append(r)
|
||||
|
||||
prim = {}
|
||||
for sid, c in cohort.items():
|
||||
weight = c.get("weight_kg") or 75.0
|
||||
events = ne_by.get(sid, [])
|
||||
ne_h24 = 0.0
|
||||
for ev in events:
|
||||
sm, em, rate = ev["start_min"], ev["end_min"], ev["rate"]
|
||||
if None in (sm, em, rate): continue
|
||||
if sm <= H_SNAPSHOT*60 <= em and rate > ne_h24: ne_h24 = rate
|
||||
|
||||
lacs = sorted(lac_by.get(sid, []), key=lambda x: x["offset_min"] or 0)
|
||||
lac_h24 = None
|
||||
if lacs:
|
||||
near = [r for r in lacs if r["offset_min"] is not None
|
||||
and 18*60 <= r["offset_min"] <= 28*60]
|
||||
lac_h24 = near[-1]["val"] if near else lacs[-1]["val"]
|
||||
|
||||
v = vital_by.get(sid, {})
|
||||
map_h24 = v.get("map")
|
||||
hr_h24 = v.get("hr")
|
||||
pr = map_h24 / (ne_h24 + 0.01) if map_h24 is not None else None
|
||||
|
||||
prim[sid] = {
|
||||
"ne_at_h24": ne_h24,
|
||||
"fluid_out_per_kg": fout.get(sid, 0) / weight,
|
||||
"lactate_h24": lac_h24,
|
||||
"map_h24": map_h24, "hr_h24": hr_h24,
|
||||
"pressor_resistance": pr,
|
||||
}
|
||||
return prim
|
||||
|
||||
|
||||
# ── 5-term continuous formula ──────────────────────────────────────────────
|
||||
def formula_features(p):
|
||||
lac = p.get("lactate_h24")
|
||||
fout = p.get("fluid_out_per_kg")
|
||||
ne24 = p.get("ne_at_h24", 0.0) or 0.0
|
||||
hr = p.get("hr_h24")
|
||||
pr = p.get("pressor_resistance")
|
||||
if lac is None or fout is None or hr is None or pr is None:
|
||||
return None
|
||||
return [
|
||||
max(0, lac - 2.5), # lactate hinge
|
||||
max(0, 20 - fout), # oliguria hinge
|
||||
1.0 if ne24 > 0.08 else 0.0, # NE persistence
|
||||
abs(hr - 85) / 20, # HR deviation
|
||||
math.log(pr + 1.0), # pressor efficiency
|
||||
]
|
||||
|
||||
|
||||
FEATURE_LABELS = [
|
||||
"max(0, lactate_h24 − 2.5)",
|
||||
"max(0, 20 − fluid_out_per_kg)",
|
||||
"I(ne_at_h24 > 0.08)",
|
||||
"|hr_h24 − 85| / 20",
|
||||
"log(pressor_resistance + 1)",
|
||||
]
|
||||
|
||||
|
||||
# ── Ordinal bedside score (0–13 pts) ───────────────────────────────────────
|
||||
def bedside_score(p):
|
||||
lac = p.get("lactate_h24")
|
||||
fout = p.get("fluid_out_per_kg")
|
||||
ne24 = p.get("ne_at_h24", 0.0) or 0.0
|
||||
hr = p.get("hr_h24")
|
||||
pr = p.get("pressor_resistance")
|
||||
if lac is None or fout is None or hr is None or pr is None:
|
||||
return None, None
|
||||
|
||||
# Lactate (0 / 2 / 4)
|
||||
if lac < 2.5: pts_lac = 0
|
||||
elif lac <= 4.0: pts_lac = 2
|
||||
else: pts_lac = 4
|
||||
|
||||
# Oliguria (0 / 1 / 2)
|
||||
if fout >= 20: pts_olig = 0
|
||||
elif fout >= 10: pts_olig = 1
|
||||
else: pts_olig = 2
|
||||
|
||||
# NE persistence (0 / 3)
|
||||
pts_ne = 3 if ne24 > 0.08 else 0
|
||||
|
||||
# HR deviation (0 / 1 / 2)
|
||||
if 70 <= hr <= 100: pts_hr = 0
|
||||
elif (60 <= hr < 70) or (100 < hr <= 120): pts_hr = 1
|
||||
else: pts_hr = 2
|
||||
|
||||
# Pressor efficiency (0 / 1 / 2)
|
||||
if pr > 3000: pts_pr = 0
|
||||
elif pr >= 1000: pts_pr = 1
|
||||
else: pts_pr = 2
|
||||
|
||||
total = pts_lac + pts_olig + pts_ne + pts_hr + pts_pr
|
||||
breakdown = {
|
||||
"lactate": pts_lac, "oliguria": pts_olig,
|
||||
"ne_persist": pts_ne, "hr_dev": pts_hr, "pressor_eff": pts_pr,
|
||||
}
|
||||
return total, breakdown
|
||||
|
||||
|
||||
def build_matrix(ids, primitives, cohort):
|
||||
import numpy as np
|
||||
X, y, scores, valid = [], [], [], []
|
||||
for sid in ids:
|
||||
p = primitives.get(sid)
|
||||
if p is None: continue
|
||||
f = formula_features(p)
|
||||
s, _ = bedside_score(p)
|
||||
if f is None or s is None: continue
|
||||
X.append(f)
|
||||
y.append(int(cohort[sid].get("died") or 0))
|
||||
scores.append(s)
|
||||
valid.append(sid)
|
||||
return np.array(X), np.array(y), np.array(scores), valid
|
||||
|
||||
|
||||
# ── Main ────────────────────────────────────────────────────────────────────
|
||||
def main():
|
||||
try:
|
||||
import numpy as np
|
||||
from sklearn.linear_model import LogisticRegression
|
||||
from sklearn.metrics import roc_auc_score, brier_score_loss
|
||||
from sklearn.model_selection import KFold
|
||||
except ImportError as e:
|
||||
print(f"\nERROR: {e}")
|
||||
print("Install: pip install scikit-learn numpy")
|
||||
sys.exit(1)
|
||||
|
||||
print("\n" + "█"*78)
|
||||
print(" PAPER 3 — PHASE 5b: 5-term formula + ordinal bedside score")
|
||||
print("█"*78)
|
||||
|
||||
print(f"\n[1] Fetching data...")
|
||||
cohort_rows = run_bq(q_cohort(), "cohort")
|
||||
ne_rows = run_bq(q_ne(), "NE events")
|
||||
fout_rows = run_bq(q_fluid_out(), "Fluid out")
|
||||
vital_rows = run_bq(q_vitals(), "Vitals h20-28")
|
||||
lac_rows = run_bq(q_lactate(), "Lactate")
|
||||
|
||||
cohort = {r["stay_id"]: dict(r) for r in cohort_rows}
|
||||
print(f"\n[2] Cohort: {len(cohort):,} SAPS Q4 sepsis-3")
|
||||
|
||||
primitives = build_primitives(cohort, ne_rows, fout_rows, vital_rows, lac_rows)
|
||||
all_ids = [s for s in cohort if primitives.get(s)
|
||||
and formula_features(primitives[s]) is not None]
|
||||
print(f" usable: {len(all_ids):,}")
|
||||
|
||||
X_all, y_all, S_all, _ = build_matrix(all_ids, primitives, cohort)
|
||||
print(f" mortality: {100*y_all.mean():.1f}%")
|
||||
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
# [4] 5-fold CV — 5-term formula
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
print(f"\n[4] 5-fold CV — 5-term continuous formula")
|
||||
subject_ids = sorted(set(cohort[s]["subject_id"] for s in all_ids))
|
||||
rng = np.random.default_rng(42)
|
||||
rng.shuffle(subject_ids)
|
||||
subj_arr = np.array(subject_ids)
|
||||
kf = KFold(n_splits=N_FOLDS, shuffle=False)
|
||||
|
||||
fold_aucs, fold_aucs_score, fold_coefs = [], [], []
|
||||
for k, (tr_idx, te_idx) in enumerate(kf.split(subj_arr)):
|
||||
tr_subs = set(subj_arr[tr_idx].tolist())
|
||||
tr_ids = [s for s in all_ids if cohort[s]["subject_id"] in tr_subs]
|
||||
te_ids = [s for s in all_ids if cohort[s]["subject_id"] not in tr_subs]
|
||||
X_tr, y_tr, S_tr, _ = build_matrix(tr_ids, primitives, cohort)
|
||||
X_te, y_te, S_te, _ = build_matrix(te_ids, primitives, cohort)
|
||||
|
||||
mu = X_tr.mean(0); sd = X_tr.std(0) + 1e-9
|
||||
lr = LogisticRegression(C=1.0, max_iter=1000, random_state=42)
|
||||
lr.fit((X_tr - mu) / sd, y_tr)
|
||||
pred = lr.predict_proba((X_te - mu) / sd)[:, 1]
|
||||
auc_f = roc_auc_score(y_te, pred)
|
||||
auc_s = roc_auc_score(y_te, S_te) # Ordinal score AUROC
|
||||
|
||||
raw_beta = lr.coef_[0] / sd
|
||||
raw_int = lr.intercept_[0] - sum(lr.coef_[0][i]*mu[i]/sd[i] for i in range(len(sd)))
|
||||
fold_aucs.append(auc_f)
|
||||
fold_aucs_score.append(auc_s)
|
||||
fold_coefs.append([raw_int] + list(raw_beta))
|
||||
print(f" fold {k+1}: formula AUROC={auc_f:.4f} ordinal score AUROC={auc_s:.4f}")
|
||||
|
||||
fa = np.array(fold_aucs); fas = np.array(fold_aucs_score)
|
||||
print(f"\n 5-term formula CV AUROC: {fa.mean():.4f} ± {fa.std():.4f} "
|
||||
f"(range {fa.min():.4f}–{fa.max():.4f})")
|
||||
print(f" Ordinal score CV AUROC: {fas.mean():.4f} ± {fas.std():.4f} "
|
||||
f"(range {fas.min():.4f}–{fas.max():.4f})")
|
||||
print(f" Ordinal loss: {fa.mean()-fas.mean():+.4f}")
|
||||
|
||||
fc = np.array(fold_coefs)
|
||||
print(f"\n Coefficient stability (5 terms):")
|
||||
names = ["intercept"] + FEATURE_LABELS
|
||||
for i, name in enumerate(names):
|
||||
col = fc[:, i]
|
||||
flips = sum(1 for k in range(1, len(col)) if col[k]*col[k-1] < 0)
|
||||
print(f" {name:35s} {col.mean():>+9.4f} ± {col.std():>7.4f} flips={flips}")
|
||||
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
# [5] Bootstrap CIs — 5-term formula
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
print(f"\n[5] Bootstrap CIs — 5-term formula ({N_BOOTSTRAP} resamples)")
|
||||
mu_all = X_all.mean(0); sd_all = X_all.std(0) + 1e-9
|
||||
n = len(y_all)
|
||||
boot_coefs = []
|
||||
rng_b = np.random.default_rng(42)
|
||||
for b in range(N_BOOTSTRAP):
|
||||
idx = rng_b.integers(0, n, n)
|
||||
X_b = X_all[idx]; y_b = y_all[idx]
|
||||
if len(set(y_b.tolist())) < 2: continue
|
||||
try:
|
||||
lr = LogisticRegression(C=1.0, max_iter=500, random_state=42)
|
||||
lr.fit((X_b - mu_all) / sd_all, y_b)
|
||||
raw = lr.coef_[0] / sd_all
|
||||
intc = lr.intercept_[0] - sum(lr.coef_[0][i]*mu_all[i]/sd_all[i] for i in range(len(sd_all)))
|
||||
boot_coefs.append([intc] + list(raw))
|
||||
except Exception: continue
|
||||
if (b+1) % 250 == 0: print(f" {b+1}/{N_BOOTSTRAP}...")
|
||||
|
||||
bc = np.array(boot_coefs)
|
||||
lr_full = LogisticRegression(C=1.0, max_iter=1000, random_state=42)
|
||||
lr_full.fit((X_all - mu_all) / sd_all, y_all)
|
||||
raw_full = lr_full.coef_[0] / sd_all
|
||||
intc_full = lr_full.intercept_[0] - sum(lr_full.coef_[0][i]*mu_all[i]/sd_all[i] for i in range(len(sd_all)))
|
||||
point_all = [intc_full] + list(raw_full)
|
||||
|
||||
print(f"\n {'term':35s} {'point':>9s} {'95% CI':>22s}")
|
||||
ci_results = []
|
||||
for i, name in enumerate(names):
|
||||
col = bc[:, i]
|
||||
lo = np.percentile(col, 2.5)
|
||||
hi = np.percentile(col, 97.5)
|
||||
crosses_zero = "✗" if (lo < 0 < hi) else "✓"
|
||||
ci_str = f"({lo:+.4f}, {hi:+.4f}) {crosses_zero}"
|
||||
print(f" {name:35s} {point_all[i]:>+9.4f} {ci_str:>22s}")
|
||||
ci_results.append({"term": name, "point": float(point_all[i]),
|
||||
"ci_lo": float(lo), "ci_hi": float(hi),
|
||||
"crosses_zero": bool(lo < 0 < hi)})
|
||||
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
# [6] Ordinal score distribution + mortality per score
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
print(f"\n[6] Ordinal bedside score distribution + mortality per score")
|
||||
print(f" {'score':>5s} {'n':>5s} {'mortality':>10s} {'cum n':>6s}")
|
||||
score_buckets = defaultdict(lambda: {"n": 0, "died": 0})
|
||||
for s, y in zip(S_all, y_all):
|
||||
score_buckets[int(s)]["n"] += 1
|
||||
score_buckets[int(s)]["died"] += int(y)
|
||||
|
||||
cum_n = 0
|
||||
score_rows = []
|
||||
for s in sorted(score_buckets.keys()):
|
||||
b = score_buckets[s]
|
||||
mort = 100 * b["died"] / b["n"] if b["n"] > 0 else 0
|
||||
cum_n += b["n"]
|
||||
bar = "█" * int(mort / 3)
|
||||
print(f" {s:>5d} {b['n']:>5d} {mort:>8.1f}% {cum_n:>6d} {bar}")
|
||||
score_rows.append({"score": s, "n": b["n"], "mortality_pct": mort})
|
||||
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
# [7] Risk bands from ordinal score (clinical cut-points)
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
print(f"\n[7] Suggested risk bands (clinically meaningful cutpoints)")
|
||||
# Group into LOW / MID / HIGH by score
|
||||
def band(s):
|
||||
if s <= 3: return "low"
|
||||
if s <= 7: return "mid"
|
||||
return "high"
|
||||
|
||||
bands = defaultdict(lambda: {"n": 0, "died": 0})
|
||||
for s, y in zip(S_all, y_all):
|
||||
b = band(int(s))
|
||||
bands[b]["n"] += 1
|
||||
bands[b]["died"] += int(y)
|
||||
|
||||
print(f" {'band':6s} {'range':>7s} {'n':>5s} {'mort%':>7s}")
|
||||
band_results = {}
|
||||
for bname in ["low", "mid", "high"]:
|
||||
b = bands[bname]
|
||||
if b["n"] == 0: continue
|
||||
rng_str = {"low": "0–3", "mid": "4–7", "high": "8+"}[bname]
|
||||
mort = 100 * b["died"] / b["n"]
|
||||
print(f" {bname:6s} {rng_str:>7s} {b['n']:>5d} {mort:>6.1f}%")
|
||||
band_results[bname] = {"n": b["n"], "mortality_pct": mort}
|
||||
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
# [8] Calibration on holdout — both formula and ordinal
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
print(f"\n[8] Calibration on 30% holdout (both versions)")
|
||||
subs = list(set(cohort[s]["subject_id"] for s in all_ids))
|
||||
random.Random(42).shuffle(subs)
|
||||
n_tr = int(len(subs) * TRAIN_FRAC)
|
||||
tr_subs = set(subs[:n_tr])
|
||||
tr_ids = [s for s in all_ids if cohort[s]["subject_id"] in tr_subs]
|
||||
te_ids = [s for s in all_ids if cohort[s]["subject_id"] not in tr_subs]
|
||||
X_tr, y_tr, _, _ = build_matrix(tr_ids, primitives, cohort)
|
||||
X_te, y_te, S_te, _ = build_matrix(te_ids, primitives, cohort)
|
||||
mu = X_tr.mean(0); sd = X_tr.std(0) + 1e-9
|
||||
lr_cal = LogisticRegression(C=1.0, max_iter=1000, random_state=42)
|
||||
lr_cal.fit((X_tr - mu) / sd, y_tr)
|
||||
pred_cal = lr_cal.predict_proba((X_te - mu) / sd)[:, 1]
|
||||
|
||||
order = np.argsort(pred_cal)
|
||||
deciles = np.array_split(order, 10)
|
||||
hl_stat = 0.0
|
||||
print(f" Formula deciles:")
|
||||
print(f" {'d':>2s} {'n':>4s} {'pred':>7s} {'obs':>7s}")
|
||||
calib_bins = []
|
||||
for d, idx in enumerate(deciles):
|
||||
n_d = len(idx)
|
||||
p_mean = float(pred_cal[idx].mean())
|
||||
obs = int(y_te[idx].sum())
|
||||
exp = float(pred_cal[idx].sum())
|
||||
if exp > 0 and (n_d - exp) > 0:
|
||||
hl_stat += (obs - exp)**2 / exp + ((n_d - obs) - (n_d - exp))**2 / (n_d - exp)
|
||||
print(f" {d+1:>2d} {n_d:>4d} {p_mean:>6.3f} {obs/n_d:>6.3f}")
|
||||
calib_bins.append({"decile": d+1, "n": n_d,
|
||||
"predicted": p_mean, "observed": obs/n_d})
|
||||
print(f"\n Hosmer-Lemeshow χ²: {hl_stat:.2f} (critical 15.51)")
|
||||
if hl_stat < 15.51:
|
||||
print(f" → Well calibrated (p > 0.05)")
|
||||
|
||||
brier = brier_score_loss(y_te, pred_cal)
|
||||
print(f" Brier: {brier:.4f}")
|
||||
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
# [9] FINAL HEADLINE
|
||||
# ══════════════════════════════════════════════════════════════════════
|
||||
print(f"\n[9] ══════════ FINAL HEADLINE (Phase 5b) ══════════")
|
||||
print(f"\n Cohort: n={len(all_ids):,} sepsis-3 Q4 patients")
|
||||
print(f" Mortality: {100*y_all.mean():.1f}%")
|
||||
print(f"\n 5-term formula (continuous):")
|
||||
print(f" 5-fold CV AUROC: {fa.mean():.4f} ± {fa.std():.4f}")
|
||||
print(f" Hosmer-Lemeshow: χ² = {hl_stat:.2f} (calibrated)")
|
||||
print(f" Brier score: {brier:.4f}")
|
||||
print(f"\n Ordinal bedside score (0–13 pts):")
|
||||
print(f" 5-fold CV AUROC: {fas.mean():.4f} ± {fas.std():.4f}")
|
||||
print(f" AUROC loss vs continuous: {fa.mean()-fas.mean():+.4f}")
|
||||
print(f"\n Risk bands (ordinal score):")
|
||||
for bn in ["low", "mid", "high"]:
|
||||
br = band_results.get(bn, {})
|
||||
if br:
|
||||
print(f" {bn:6s} ({'0–3' if bn=='low' else '4–7' if bn=='mid' else '8+':>4s}): "
|
||||
f"n={br['n']:>5d} mortality={br['mortality_pct']:.1f}%")
|
||||
|
||||
# ── Save ────────────────────────────────────────────────────────────────
|
||||
output = {
|
||||
"cohort": {"n": len(all_ids), "mortality": float(y_all.mean())},
|
||||
"formula_5term": {
|
||||
"cv_auroc_mean": float(fa.mean()),
|
||||
"cv_auroc_std": float(fa.std()),
|
||||
"cv_auroc_range": [float(fa.min()), float(fa.max())],
|
||||
"coefficient_cis": ci_results,
|
||||
"calibration": {
|
||||
"hl_chi2": float(hl_stat),
|
||||
"brier": float(brier),
|
||||
"deciles": calib_bins,
|
||||
},
|
||||
},
|
||||
"ordinal_score": {
|
||||
"cv_auroc_mean": float(fas.mean()),
|
||||
"cv_auroc_std": float(fas.std()),
|
||||
"auroc_loss_vs_continuous": float(fa.mean() - fas.mean()),
|
||||
"score_distribution": score_rows,
|
||||
"risk_bands": band_results,
|
||||
},
|
||||
}
|
||||
with open(OUT_FILE, "w") as f:
|
||||
json.dump(output, f, indent=2, default=str)
|
||||
print(f"\n → Saved: {OUT_FILE}")
|
||||
print("\n" + "█"*78 + "\n")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in New Issue
Block a user