随附材料 · 实验代码
experiment.py
#!/usr/bin/env python3
"""单细胞患者泄漏：细胞随机划分 vs 患者分组划分，性能被高估多少。
判定规则写在 ../claims.json，跑之前定死。零模型是标签打乱。
置换次数按实测单次拟合耗时自适应，卡住 40 分钟墙钟预算 —— 免费层只烧 CPU。
"""
import gzip, json, time
import numpy as np
from joblib import Parallel, delayed
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold, GroupKFold
from sklearn.metrics import f1_score, accuracy_score
N_GENES, SEED = 2000, 20260827
NULL_BUDGET_S, NULL_MIN, NULL_MAX, N_JOBS = 2400, 40, 200, 3
TYPES = {1: "T", 2: "B", 3: "Macro", 4: "Endo", 5: "CAF", 6: "NK"}
def load():
with gzip.open("gse72056.txt.gz", "rt") as fh:
fh.readline()
tumor = np.array(fh.readline().rstrip("\n").split("\t")[1:], dtype=int)
malig = np.array(fh.readline().rstrip("\n").split("\t")[1:], dtype=int)
ctype = np.array(fh.readline().rstrip("\n").split("\t")[1:], dtype=int)
rows = [line.rstrip("\n").split("\t")[1:] for line in fh]
X = np.array(rows, dtype=np.float32).T # 细胞 × 基因
keep = (malig == 1) & (ctype >= 1) & (ctype <= 6) # 只用非恶性、类型明确的细胞
return X[keep], ctype[keep], tumor[keep]
def fit_eval(X, y, splits, seed):
accs, f1s = [], []
for tr, te in splits:
clf = LogisticRegression(max_iter=300, C=0.1, random_state=seed)
clf.fit(X[tr], y[tr])
p = clf.predict(X[te])
accs.append(accuracy_score(y[te], p))
f1s.append(f1_score(y[te], p, average="macro", zero_division=0))
return float(np.mean(accs)), float(np.mean(f1s))
def one_perm(X, y, rnd2, grp2, k):
"""一次置换：同一套打乱标签，分别在随机划分和患者划分下评估，返回两者之差。"""
ys = y.copy()
np.random.default_rng(SEED + 1000 + k).shuffle(ys)
return fit_eval(X, ys, rnd2, SEED)[1] - fit_eval(X, ys, grp2, SEED)[1]
def main():
t0 = time.time()
X, y, g = load()
counts = {TYPES[t]: int((y == t).sum()) for t in sorted(set(y))}
print(f"[data] 非恶性细胞 {len(y)}｜患者 {len(set(g))} 位｜类型 {counts}", flush=True)
idx = np.argsort(X.var(axis=0))[-N_GENES:]
X = X[:, idx]
X = (X - X.mean(0)) / (X.std(0) + 1e-8)
print(f"[data] 取方差最大的 {N_GENES} 个基因｜载入耗时 {time.time()-t0:.0f}s", flush=True)
rnd = list(KFold(5, shuffle=True, random_state=SEED).split(X))
grp = list(GroupKFold(5).split(X, y, groups=g))
t1 = time.time()
a_acc, a_f1 = fit_eval(X, y, rnd, SEED)
b_acc, b_f1 = fit_eval(X, y, grp, SEED)
per_fit = (time.time() - t1) / 10
gap_f1, gap_acc = a_f1 - b_f1, a_acc - b_acc
print(f"\n[A] 细胞随机划分 准确率 {a_acc:.4f} macro-F1 {a_f1:.4f}", flush=True)
print(f"[B] 患者分组划分 准确率 {b_acc:.4f} macro-F1 {b_f1:.4f}", flush=True)
print(f"[差] 随机划分高估 准确率 +{gap_acc:.4f} macro-F1 +{gap_f1:.4f}"
f"（相对高估 {100*gap_f1/max(b_f1,1e-9):.1f}%）｜单次拟合 {per_fit:.1f}s", flush=True)
n_perm = int(np.clip(NULL_BUDGET_S * N_JOBS / max(per_fit * 4, 1e-6), NULL_MIN, NULL_MAX))
print(f"[零模型] 预算 {NULL_BUDGET_S}s × {N_JOBS} 并发 → 置换 {n_perm} 次", flush=True)
rnd2, grp2 = rnd[:2], grp[:2]
null = np.array(Parallel(n_jobs=N_JOBS, verbose=1)(
delayed(one_perm)(X, y, rnd2, grp2, k) for k in range(n_perm)))
p95 = float(np.percentile(null, 95))
emp_p = float((null >= gap_f1).mean())
h1, h2 = bool(gap_f1 > 0), bool(gap_f1 > p95)
print(f"\n[零模型] 标签打乱 {n_perm} 次的 F1 差：均值 {null.mean():+.4f}｜"
f"95分位 {p95:+.4f}｜最大 {null.max():+.4f}", flush=True)
print(f"[H1] 随机划分高于患者划分 → {'成立' if h1 else '被推翻'}", flush=True)
print(f"[H2] 差距超出零模型 95 分位 → {'成立' if h2 else '被推翻'}（经验 p={emp_p:.4f}）", flush=True)
json.dump({"seed": SEED, "n_cells": int(len(y)), "n_patients": int(len(set(g))),
"n_genes": N_GENES, "n_perm": int(n_perm), "sec_per_fit": per_fit,
"random_split": {"acc": a_acc, "macro_f1": a_f1},
"patient_split": {"acc": b_acc, "macro_f1": b_f1},
"gap_f1": gap_f1, "gap_acc": gap_acc,
"relative_overestimate": gap_f1 / max(b_f1, 1e-9),
"null_mean": float(null.mean()), "null_p95": p95, "empirical_p": emp_p,
"H1_passes": h1, "H2_passes": h2, "class_counts": counts,
"null_distribution": null.tolist(),
"wall_clock_s": round(time.time() - t0)},
open("results.json", "w"), ensure_ascii=False, indent=1)
print(f"[done] 总耗时 {round(time.time()-t0)}s", flush=True)
if __name__ == "__main__":
main()
check1.py
#!/usr/bin/env python3
"""核心复核 1 —— 这把尺子到底测的是不是"患者身份"。
主实验只给出 +0.0237 的差，太小，必须先证明测量本身有效：
a) 假患者（打乱患者标签）：分组划分退化成随机划分，差应该塌到 ~0
b) 注入患者特异偏移：差应该明显变大
c) 折大小匹配的随机划分：排除"KFold 与 GroupKFold 折几何不同"这个混杂
a) 或 b) 不成立，说明主实验那个差不能解释成泄漏。
"""
import gzip, json
import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold, GroupKFold
from sklearn.metrics import f1_score
SEED, N_GENES = 20260827, 2000
def load():
with gzip.open("gse72056.txt.gz", "rt") as fh:
fh.readline()
tumor = np.array(fh.readline().rstrip("\n").split("\t")[1:], dtype=int)
malig = np.array(fh.readline().rstrip("\n").split("\t")[1:], dtype=int)
ctype = np.array(fh.readline().rstrip("\n").split("\t")[1:], dtype=int)
rows = [line.rstrip("\n").split("\t")[1:] for line in fh]
X = np.array(rows, dtype=np.float32).T
keep = (malig == 1) & (ctype >= 1) & (ctype <= 6)
X, y, g = X[keep], ctype[keep], tumor[keep]
idx = np.argsort(X.var(axis=0))[-N_GENES:]
X = X[:, idx]
return (X - X.mean(0)) / (X.std(0) + 1e-8), y, g
def f1_of(X, y, splits):
out = []
for tr, te in splits:
clf = LogisticRegression(max_iter=300, C=0.1, random_state=SEED).fit(X[tr], y[tr])
out.append(f1_score(y[te], clf.predict(X[te]), average="macro", zero_division=0))
return float(np.mean(out))
def gap(X, y, g, seed=SEED):
rnd = list(KFold(5, shuffle=True, random_state=seed).split(X))
grp = list(GroupKFold(5).split(X, y, groups=g))
return f1_of(X, y, rnd) - f1_of(X, y, grp)
def matched_random(X, y, g, seed=SEED):
"""把患者折的大小照搬给随机划分 —— 折几何一样，只差"同患者是否跨侧"。"""
grp = list(GroupKFold(5).split(X, y, groups=g))
sizes = [len(te) for _, te in grp]
order = np.random.default_rng(seed).permutation(len(y))
splits, start = [], 0
for size in sizes:
te = order[start:start + size]
splits.append((np.setdiff1d(order, te), te))
start += size
return f1_of(X, y, splits) - f1_of(X, y, grp)
def main():
X, y, g = load()
out = {"observed_gap": gap(X, y, g)}
print(f"[主实验复现] 差 = {out['observed_gap']:+.4f}", flush=True)
rng = np.random.default_rng(SEED)
fake = [gap(X, y, rng.permutation(g)) for _ in range(3)]
out["fake_patient_gaps"] = fake
out["fake_patient_mean"] = float(np.mean(fake))
print(f"[a 假患者] 差 = {np.mean(fake):+.4f}（3 次 {['%+.4f' % v for v in fake]}）"
f" → 应≈0，实际{'塌了 ✅' if abs(np.mean(fake)) < abs(out['observed_gap']) / 2 else '没塌 ❌'}",
flush=True)
spiked = X.copy()
for pid in set(g):
m = g == pid
spiked[m] += rng.normal(0, 1.0, size=(1, X.shape[1])).astype(np.float32)
out["spiked_gap"] = gap(spiked, y, g)
print(f"[b 注入患者偏移] 差 = {out['spiked_gap']:+.4f} → 应明显变大，"
f"实际{'变大 ✅' if out['spiked_gap'] > out['observed_gap'] else '没变大 ❌'}", flush=True)
out["matched_geometry_gap"] = matched_random(X, y, g)
print(f"[c 折几何匹配] 差 = {out['matched_geometry_gap']:+.4f} → "
f"与主实验同号同量级即说明不是折几何造成的", flush=True)
out["verdict_instrument_valid"] = bool(
abs(out["fake_patient_mean"]) < abs(out["observed_gap"]) / 2
and out["spiked_gap"] > out["observed_gap"])
json.dump(out, open("check1_instrument.json", "w"), ensure_ascii=False, indent=1)
print(f"\n[结论] 尺子有效 = {out['verdict_instrument_valid']}", flush=True)
if __name__ == "__main__":
main()
check2.py
#!/usr/bin/env python3
"""核心复核 2 —— 2.6% 这个数字是不是被"任务太容易"压住了。
主实验准确率 98.8%，离 1.0 只剩 1.2 个点，泄漏再大也没地方涨（天花板效应）。
如果差随任务变难而变大，那 2.6% 就不是普适结论，而是这个易任务的下界 ——
这是最需要知道的一条，也是把这套做法搬到别的数据上时的判断依据。
把任务逐级变难：减基因数、减每类训练细胞数。
"""
import gzip, json
import numpy as np
from joblib import Parallel, delayed
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold, GroupKFold
from sklearn.metrics import f1_score
SEED = 20260827
def load(n_genes):
with gzip.open("gse72056.txt.gz", "rt") as fh:
fh.readline()
tumor = np.array(fh.readline().rstrip("\n").split("\t")[1:], dtype=int)
malig = np.array(fh.readline().rstrip("\n").split("\t")[1:], dtype=int)
ctype = np.array(fh.readline().rstrip("\n").split("\t")[1:], dtype=int)
rows = [line.rstrip("\n").split("\t")[1:] for line in fh]
X = np.array(rows, dtype=np.float32).T
keep = (malig == 1) & (ctype >= 1) & (ctype <= 6)
X, y, g = X[keep], ctype[keep], tumor[keep]
idx = np.argsort(X.var(axis=0))[-n_genes:]
X = X[:, idx]
return (X - X.mean(0)) / (X.std(0) + 1e-8), y, g
def score(X, y, splits, cap, seed):
"""cap = 每类最多用多少训练细胞，用来把任务调难。"""
rng = np.random.default_rng(seed)
accs = []
for tr, te in splits:
if cap:
sub = np.concatenate([rng.permutation(tr[y[tr] == c])[:cap] for c in np.unique(y[tr])])
tr = sub if len(np.unique(y[sub])) > 1 else tr
clf = LogisticRegression(max_iter=300, C=0.1, random_state=seed).fit(X[tr], y[tr])
accs.append(f1_score(y[te], clf.predict(X[te]), average="macro", zero_division=0))
return float(np.mean(accs))
def one(n_genes, cap):
X, y, g = load(n_genes)
rnd = list(KFold(5, shuffle=True, random_state=SEED).split(X))
grp = list(GroupKFold(5).split(X, y, groups=g))
a, b = score(X, y, rnd, cap, SEED), score(X, y, grp, cap, SEED)
return {"n_genes": n_genes, "cap_per_class": cap, "random_f1": a, "patient_f1": b,
"gap": a - b, "relative": (a - b) / max(b, 1e-9)}
def main():
grid = [(2000, None), (2000, 40), (2000, 15), (200, None), (200, 15), (50, 15)]
rows = Parallel(n_jobs=3, verbose=1)(delayed(one)(n, c) for n, c in grid)
rows.sort(key=lambda r: r["patient_f1"], reverse=True)
print(f"\n{'基因数':>6} {'每类训练细胞':>12} {'患者划分F1':>11} {'随机划分F1':>11} "
f"{'绝对高估':>9} {'相对高估':>9}")
for r in rows:
print(f"{r['n_genes']:>6} {str(r['cap_per_class'] or '全部'):>12} "
f"{r['patient_f1']:>11.4f} {r['random_f1']:>11.4f} "
f"{r['gap']:>+9.4f} {100*r['relative']:>8.1f}%")
easiest, hardest = rows[0], rows[-1]
grows = hardest["gap"] > easiest["gap"]
print(f"\n[结论] 任务从 F1 {easiest['patient_f1']:.3f} 变难到 {hardest['patient_f1']:.3f} 时，"
f"高估幅度 {easiest['gap']:+.4f} → {hardest['gap']:+.4f}"
f"（{'随难度放大 ✅ 2.6% 只是易任务下界' if grows else '未随难度放大'}）")
json.dump({"grid": rows, "gap_grows_with_difficulty": bool(grows)},
open("check2_difficulty.json", "w"), ensure_ascii=False, indent=1)
if __name__ == "__main__":
main()