随附材料 · 实验代码
fetch_corpus.py
#!/usr/bin/env python3
"""从 PubMed 公开 E-utilities 取结直肠癌腹膜转移文献的 MeSH 标注。"""
import json, time, urllib.parse, urllib.request, pathlib
KEY = pathlib.Path.home().joinpath(".config/ncbi/api_key").read_text().strip()
EU = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"
TERM = ('("colorectal neoplasms"[MeSH] OR colorectal cancer[Title/Abstract]) AND '
'("peritoneal neoplasms"[MeSH] OR peritoneal metastas*[Title/Abstract] OR '
'peritoneal carcinomatosis[Title/Abstract]) AND 2015:2026[dp]')
def get(path, **kw):
kw["api_key"] = KEY
with urllib.request.urlopen(f"{EU}/{path}?{urllib.parse.urlencode(kw)}", timeout=180) as fh:
return fh.read()
ids = json.loads(get("esearch.fcgi", db="pubmed", term=TERM, retmax=10000,
retmode="json"))["esearchresult"]["idlist"]
print(f" 检索命中 {len(ids)} 篇")
recs = []
for i in range(0, len(ids), 200):
xml = get("efetch.fcgi", db="pubmed", id=",".join(ids[i:i + 200]), retmode="xml").decode("utf-8", "ignore")
for art in xml.split("<PubmedArticle>")[1:]:
yr = None
for tag in ("<Year>", "<MedlineDate>"):
if tag in art:
seg = art.split(tag, 1)[1][:12]
digits = "".join(ch for ch in seg if ch.isdigit())[:4]
if len(digits) == 4:
yr = int(digits); break
mesh = []
for chunk in art.split("<DescriptorName")[1:]:
if ">" in chunk:
mesh.append(chunk.split(">", 1)[1].split("<")[0].strip())
if yr and mesh:
recs.append({"year": yr, "mesh": sorted(set(mesh))})
print(f" 已解析 {len(recs)} 篇", end="\r")
time.sleep(0.15)
json.dump({"term": TERM, "n": len(recs), "records": recs}, open("corpus.json", "w"), ensure_ascii=False)
print(f"\n 有年份且有 MeSH 的: {len(recs)} 篇")
experiment.py
#!/usr/bin/env python3
"""结直肠癌腹膜转移文献：2015–2019 vs 2020–2026 的主题词占比变化。
判定规则写在 ../claims.json，跑之前定死。零模型是年份标签随机打乱 1000 次。
"""
import json, collections
import numpy as np
from scipy.stats import fisher_exact
N_PERM, MIN_DOCS, SEED = 1000, 20, 20260827
def bh(p):
p = np.asarray(p, float); n = len(p); o = np.argsort(p); adj = np.empty(n); prev = 1.0
for r, i in enumerate(o[::-1]):
prev = min(prev, p[i] * n / (n - r)); adj[i] = prev
return adj
def screen(recs, labels, terms):
"""labels: True=后期。返回每个词的 (p, 前期占比, 后期占比)。"""
n_late = int(labels.sum()); n_early = len(labels) - n_late
out = []
for t in terms:
has = np.array([t in r["mesh"] for r in recs])
a = int((has & labels).sum()) # 后期有
b = int((has & ~labels).sum()) # 前期有
_, p = fisher_exact([[a, n_late - a], [b, n_early - b]])
out.append((t, float(p), b / n_early, a / n_late))
return out
def main():
d = json.load(open("corpus.json"))
recs = [r for r in d["records"] if 2015 <= r["year"] <= 2026]
labels = np.array([r["year"] >= 2020 for r in recs])
print(f"[data] 文献 {len(recs)} 篇｜2015–2019 {int((~labels).sum())} 篇｜2020–2026 {int(labels.sum())} 篇")
cnt = collections.Counter(m for r in recs for m in r["mesh"])
terms = [t for t, c in cnt.items() if c >= MIN_DOCS]
print(f"[data] 至少出现在 {MIN_DOCS} 篇里的 MeSH 词: {len(terms)} 个")
obs = screen(recs, labels, terms)
q = bh([x[1] for x in obs])
sig = [(t, p, e, l, qq) for (t, p, e, l), qq in zip(obs, q) if qq < 0.05]
sig.sort(key=lambda x: -(x[3] - x[2]))
print(f"[实测] FDR<0.05 的主题词: {len(sig)} 个")
rng = np.random.default_rng(SEED)
null = []
for _ in range(N_PERM):
lab = labels.copy(); rng.shuffle(lab)
pp = [x[1] for x in screen(recs, lab, terms)]
null.append(int((bh(pp) < 0.05).sum()))
null = np.array(null)
p95 = float(np.percentile(null, 95))
h1 = bool(len(sig) > p95 and len(sig) >= 5)
print(f"[零模型] {N_PERM} 次打乱：均值 {null.mean():.2f}｜95分位 {p95:.1f}｜最大 {null.max()}")
print(f"[H1] {len(sig)} > {p95:.1f} 且 >=5 → {'✅成立' if h1 else '❌被推翻'}")
rising = [s for s in sig if s[3] > s[2]][:12]
falling = sorted([s for s in sig if s[3] < s[2]], key=lambda x: x[3] - x[2])[:12]
print("\n 升幅最大：")
for t, p, e, l, qq in rising[:6]:
print(f" {t[:44]:<46} {e*100:5.1f}% → {l*100:5.1f}% (FDR {qq:.1e})")
print(" 降幅最大：")
for t, p, e, l, qq in falling[:6]:
print(f" {t[:44]:<46} {e*100:5.1f}% → {l*100:5.1f}% (FDR {qq:.1e})")
HOT = ("Immunotherapy", "Single-Cell", "Molecular Targeted", "Immune Checkpoint",
"Antineoplastic Combined", "Tumor Microenvironment", "Biomarkers")
hot_hits = [s[0] for s in rising if any(h.lower() in s[0].lower() for h in HOT)]
h2 = len(hot_hits) > 0
print(f"\n[H2] 升幅词里含免疫/单细胞/靶向类新兴主题 → {'✅成立' if h2 else '❌被推翻'}｜命中 {hot_hits}")
json.dump({"seed": SEED, "n_docs": len(recs), "n_early": int((~labels).sum()),
"n_late": int(labels.sum()), "n_terms": len(terms), "min_docs": MIN_DOCS,
"n_significant": len(sig), "null_mean": float(null.mean()), "null_p95": p95,
"null_max": int(null.max()),
"empirical_p": float((null >= len(sig)).mean()),
"H1_passes": h1, "H2_passes": h2, "hot_hits": hot_hits,
"rising": [{"term": t, "early": e, "late": l, "fdr": qq} for t, p, e, l, qq in rising],
"falling": [{"term": t, "early": e, "late": l, "fdr": qq} for t, p, e, l, qq in falling],
"null_distribution": null.tolist()},
open("results.json", "w"), ensure_ascii=False, indent=1)
if __name__ == "__main__":
main()
adversarial_check.py
#!/usr/bin/env python3
"""敌意复核：这些"变化"里有多少是 MeSH 词表版本与索引完整度造成的假象？"""
import json, collections
import numpy as np
d = json.load(open("corpus.json")); recs = d["records"]
print("== 质疑 1：每篇的 MeSH 标注数逐年变了吗（索引完整度）==")
by = collections.defaultdict(list)
for r in recs:
by[r["year"]].append(len(r["mesh"]))
for y in sorted(by):
print(f" {y} 文献 {len(by[y]):>4} 篇 平均 MeSH 词数 {np.mean(by[y]):5.1f}")
early = [n for y, v in by.items() if y < 2020 for n in v]
late = [n for y, v in by.items() if y >= 2020 for n in v]
print(f" → 前期平均 {np.mean(early):.1f} 词/篇，后期 {np.mean(late):.1f} 词/篇"
f"（{'⚠️ 后期标注更少，会系统性压低所有词的占比' if np.mean(late) < np.mean(early) else '差异不大'}）")
print("\n== 质疑 2：HIPEC 从 0.0% 跳到 29.1%，是新概念还是新词条 ==")
first = {}
for r in sorted(recs, key=lambda x: x["year"]):
for m in r["mesh"]:
first.setdefault(m, r["year"])
res = json.load(open("results.json"))
for row in res["rising"][:8]:
t = row["term"]
print(f" {t[:46]:<48} 语料里首次出现于 {first.get(t)} 年")
print(" → 首次出现年份 >=2020 的词，极可能是 MeSH 当年新增词条，而不是研究方向新出现。")
print("\n== 质疑 3：把人口学检索标签剔掉，还剩多少真变化 ==")
CHECK_TAGS = {"Male", "Female", "Humans", "Adult", "Aged", "Middle Aged", "Young Adult",
"Adolescent", "Aged, 80 and over", "Child", "Animals", "Mice",
"Retrospective Studies", "Prospective Studies", "Treatment Outcome",
"Survival Rate", "Follow-Up Studies"}
sig_terms = [r["term"] for r in res["rising"]] + [r["term"] for r in res["falling"]]
real = [t for t in sig_terms if t not in CHECK_TAGS]
print(f" 显著词（示例集 {len(sig_terms)} 个）里，剔掉检索标签后剩 {len(real)} 个内容词")
print(f" 被剔掉的: {sorted(set(sig_terms) & CHECK_TAGS)}")
print("\n== 质疑 4：只看 2020 年后就存在的词，变化还在吗 ==")
stable = [t for t in real if first.get(t, 9999) < 2018]
print(f" 语料里 2018 年前就出现过的内容词: {len(stable)} 个 → {stable[:8]}")
print("\n== 质疑 5：零模型均值 0.04 是不是太干净了 ==")
null = np.array(res["null_distribution"])
print(f" 1000 次打乱：{int((null==0).sum())} 次给 0，最大 {null.max()}")
print(" → 零模型这么干净，说明 Fisher+BH 在这个规模下不会随便造出显著词；")
print(" 所以 55 个显著词是真的统计信号 —— 但信号来源可能是索引而非研究内容，见质疑 1–3。")
json.dump({"early_mesh_per_doc": float(np.mean(early)), "late_mesh_per_doc": float(np.mean(late)),
"first_year": {t: first.get(t) for t in [r["term"] for r in res["rising"][:8]]},
"content_terms_after_filter": len(real), "stable_content_terms": len(stable)},
open("adversarial.json", "w"), ensure_ascii=False, indent=1)
corrected.py
#!/usr/bin/env python3
"""校正版：扣掉 MeSH 词表版本与索引完整度两个混杂，再看还剩多少真变化。
两处校正：
1) 只保留 2018 年前就在语料里出现过的词 —— 排除 MeSH 新增词条
2) 用「该词占本篇 MeSH 总数的份额」而不是「出现率」—— 抵消索引完整度漂移
另加一条对照：把人口学检索标签单独拎出来当"内部对照"。若校正有效，
这些与研究内容无关的标签在校正后应当不再显著。
"""
import json, collections
import numpy as np
from scipy.stats import mannwhitneyu
CHECK_TAGS = {"Male", "Female", "Humans", "Adult", "Aged", "Middle Aged", "Young Adult",
"Adolescent", "Aged, 80 and over", "Child", "Animals", "Retrospective Studies",
"Prospective Studies", "Treatment Outcome", "Survival Rate", "Follow-Up Studies"}
def bh(p):
p = np.asarray(p, float); n = len(p); o = np.argsort(p); adj = np.empty(n); prev = 1.0
for r, i in enumerate(o[::-1]):
prev = min(prev, p[i] * n / (n - r)); adj[i] = prev
return adj
recs = [r for r in json.load(open("corpus.json"))["records"] if 2015 <= r["year"] <= 2026]
first = {}
for r in sorted(recs, key=lambda x: x["year"]):
for m in r["mesh"]:
first.setdefault(m, r["year"])
cnt = collections.Counter(m for r in recs for m in r["mesh"])
terms = [t for t, c in cnt.items() if c >= 20 and first[t] < 2018]
print(f" 校正后参与检验的词: {len(terms)} 个（原 115 个，剔除 MeSH 新增词条）")
late = np.array([r["year"] >= 2020 for r in recs])
share = {t: np.array([(1.0 / len(r["mesh"])) if t in r["mesh"] else 0.0 for r in recs])
for t in terms}
ps = []
for t in terms:
a, b = share[t][late], share[t][~late]
ps.append(mannwhitneyu(a, b, alternative="two-sided")[1])
qs = bh(ps)
sig = [(t, float(q), float(share[t][~late].mean()), float(share[t][late].mean()))
for t, q in zip(terms, qs) if q < 0.05]
sig.sort(key=lambda x: -(x[3] - x[2]))
tag_sig = [s for s in sig if s[0] in CHECK_TAGS]
content_sig = [s for s in sig if s[0] not in CHECK_TAGS]
print(f" 校正后显著: {len(sig)} 个｜其中内容词 {len(content_sig)}｜检索标签 {len(tag_sig)}")
print(f" 内部对照：检索标签{'仍有 %d 个显著 ⚠️ 校正不彻底' % len(tag_sig) if tag_sig else '全部不再显著 ✅ 校正有效'}")
print("\n 校正后真正上升的内容词：")
for t, q, e, l in content_sig[:8]:
if l > e: print(f" {t[:44]:<46} 份额 {e:.4f} → {l:.4f} (FDR {q:.1e})")
print(" 校正后真正下降的内容词：")
for t, q, e, l in sorted(content_sig, key=lambda x: x[3] - x[2])[:6]:
if l < e: print(f" {t[:44]:<46} 份额 {e:.4f} → {l:.4f} (FDR {q:.1e})")
json.dump({"n_terms_after_filter": len(terms), "n_sig": len(sig),
"n_content_sig": len(content_sig), "n_tag_sig": len(tag_sig),
"correction_worked": len(tag_sig) == 0,
"rising": [{"term": t, "early_share": e, "late_share": l, "fdr": q}
for t, q, e, l in content_sig if l > e][:12],
"falling": [{"term": t, "early_share": e, "late_share": l, "fdr": q}
for t, q, e, l in sorted(content_sig, key=lambda x: x[3] - x[2]) if l < e][:12]},
open("corrected.json", "w"), ensure_ascii=False, indent=1)
plot_results.py
#!/usr/bin/env python3
import json
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np, collections
d=json.load(open("corpus.json"))["records"]
by=collections.defaultdict(list)
for r in d:
if 2015<=r["year"]<=2026: by[r["year"]].append(len(r["mesh"]))
ys=sorted(by); m=[np.mean(by[y]) for y in ys]
fig,ax=plt.subplots(1,2,figsize=(12,4.8))
ax[0].plot(ys,m,"o-",color="#c0392b"); ax[0].axvline(2019.5,ls="--",c="k",lw=1)
ax[0].text(2019.6,max(m)*0.97,"period split",fontsize=8)
ax[0].set_xlabel("year"); ax[0].set_ylabel("MeSH terms per article")
ax[0].set_title("Indexing completeness drifts:\nthis alone deflates every term's rate in the later period")
ax[0].grid(alpha=.3)
r=json.load(open("results.json")); null=np.array(r["null_distribution"])
ax[1].hist(null,bins=np.arange(0,max(null.max(),6)+2)-.5,color="#95a5a6",edgecolor="white")
ax[1].axvline(r["n_significant"],color="#c0392b",lw=2.5)
ax[1].text(r["n_significant"]*0.55,ax[1].get_ylim()[1]*.6,f"observed = {r['n_significant']}",color="#c0392b",fontsize=10)
ax[1].set_xlabel("significant MeSH terms"); ax[1].set_ylabel("permutations")
ax[1].set_title(f"Year-shuffle null (n={len(null)}): the statistics are sound —\nthe problem is the data, not the test")
plt.tight_layout(); plt.savefig("figs/indexing_drift.png",dpi=150); print("figs/indexing_drift.png 已生成")