#!/usr/bin/env python3 """按 claims.json 的判据做统计判定，并把每个计数式说法绑上连续量（p 值 + 效应量）。""" import json, os import numpy as np from scipy import stats HERE = os.path.dirname(os.path.abspath(__file__)) RES = os.path.join(HERE, "..", "results") def exact_wilcoxon_p(diffs): """精确 Wilcoxon 符号秩双侧 p，纯枚举（n<=20）。定义写死在这里， audit_delivery.py 用同一定义独立重写一遍，两边必须逐位一致。""" d = [x for x in diffs if x != 0.0] n = len(d) if n == 0 or n > 20: return None a = sorted((abs(x), i) for i, x in enumerate(d)) ranks = [0.0] * n i = 0 while i < n: j = i while j + 1 < n and a[j + 1][0] == a[i][0]: j += 1 r = (i + j) / 2.0 + 1.0 for k in range(i, j + 1): ranks[a[k][1]] = r i = j + 1 w_plus = sum(ranks[i] for i in range(n) if d[i] > 0) counts = {0.0: 1} for r in ranks: nxt = {} for v, c in counts.items(): nxt[v] = nxt.get(v, 0) + c nxt[v + r] = nxt.get(v + r, 0) + c counts = nxt total = 2 ** n le = sum(c for v, c in counts.items() if v <= w_plus) ge = sum(c for v, c in counts.items() if v >= w_plus) return min(1.0, 2.0 * min(le, ge) / total), w_plus, n def paired_block(rows, key_a, key_b, label): a = np.array([r[key_a] for r in rows], dtype=float) b = np.array([r[key_b] for r in rows], dtype=float) d = a - b w = stats.wilcoxon(a, b, alternative="two-sided", zero_method="wilcox") t = stats.ttest_rel(a, b) n = len(d) wins = int((d > 0).sum()) # 配对秩双列相关（rank-biserial）：Wilcoxon 的标准效应量 nz = d[d != 0] r_rb = float("nan") if len(nz): ranks = stats.rankdata(np.abs(nz)) rplus = ranks[nz > 0].sum() total = ranks.sum() r_rb = float(2 * rplus / total - 1) return { "comparison": label, "n_pairs": n, "count_claim": f"{wins}/{n} 折 {key_a} 胜过 {key_b}", "wins": wins, "median_difference": float(np.median(d)), "mean_difference": float(d.mean()), "difference_ci95_bootstrap": None, "wilcoxon_p_two_sided": float(w.pvalue), "wilcoxon_p_exact_enumerated": (exact_wilcoxon_p(d.tolist()) or [None])[0], "wilcoxon_w_plus": (exact_wilcoxon_p(d.tolist()) or [None, None])[1], "wilcoxon_statistic": float(w.statistic), "paired_t_p_two_sided": float(t.pvalue), "rank_biserial_effect_size": r_rb, "mean_a": float(a.mean()), "mean_b": float(b.mean()), } def boot_ci(rows, key_a, key_b, seed=20260912, n=10000): a = np.array([r[key_a] for r in rows], dtype=float) b = np.array([r[key_b] for r in rows], dtype=float) d = a - b rng = np.random.default_rng(seed) idx = rng.integers(0, len(d), size=(n, len(d))) means = d[idx].mean(axis=1) return [float(np.quantile(means, 0.025)), float(np.quantile(means, 0.975))] def main(): an = json.load(open(os.path.join(RES, "analysis.json"), encoding="utf-8")) z = np.load(os.path.join(RES, "scores.npz")) out = {"params": an["params"], "universe": an["universe"]} # ---- H4 负对照（尺子自检）先判：它不过，别的都不算数 ---- lp = an["label_permutation"] out["H4_label_permutation"] = { "mean_auroc": lp["mean"], "ci95": lp["ci95"], "n": lp["n"], "covers_0.5": bool(lp["ci95"][0] <= 0.5 <= lp["ci95"][1]), "within_0.5_pm_0.02": bool(abs(lp["mean"] - 0.5) <= 0.02), "verdict": "passed" if (lp["ci95"][0] <= 0.5 <= lp["ci95"][1] and abs(lp["mean"] - 0.5) <= 0.02) else "FAILED", } # ---- H0 正对照 ---- pc = an["positive_control_h0"] h0 = paired_block(pc, "auroc_gout_rwr", "auroc_matched_null", "正对照：痛风自预测 vs 同折度数配对随机种子") h0["difference_ci95_bootstrap"] = boot_ci(pc, "auroc_gout_rwr", "auroc_matched_null") h0["mean_auroc_gout_rwr"] = float(np.mean([r["auroc_gout_rwr"] for r in pc])) h0["mean_auroc_degree_only"] = float(np.mean([r["auroc_degree_only"] for r in pc])) h0["verdict"] = "passed" if (h0["wilcoxon_p_two_sided"] < 0.05 and h0["median_difference"] > 0) else "REFUTED" out["H0_positive_control"] = h0 # ---- H1 主判据 ---- m = an["main"] z_score = (m["auroc_gout_rwr"] - m["null_mean"]) / m["null_sd"] out["H1_main"] = { "auroc_gout_rwr": m["auroc_gout_rwr"], "auprc_gout_rwr": m["auprc_gout_rwr"], "auroc_degree_only": m["auroc_degree_only"], "auprc_degree_only": m["auprc_degree_only"], "prevalence_baseline_auprc": an["universe"]["prevalence"], "null_mean": m["null_mean"], "null_sd": m["null_sd"], "null_p95": m["null_p95"], "null_max": m["null_max"], "n_null": m["n_null"], "empirical_p": m["empirical_p"], "z_vs_null": float(z_score), "delta_vs_null_mean": float(m["auroc_gout_rwr"] - m["null_mean"]), "delta_vs_degree_only": float(m["auroc_gout_rwr"] - m["auroc_degree_only"]), "beats_degree_only": bool(m["auroc_gout_rwr"] > m["auroc_degree_only"]), "null_by_seed": m["null_by_seed"], "count_claim_bound": f"3 个种子的零模型分布全部低于实测值；对应连续量：经验 p={m['empirical_p']:.4g}，" f"z={z_score:.2f}，AUROC 差 {m['auroc_gout_rwr'] - m['null_mean']:+.4f}", "verdict": "passed" if (m["empirical_p"] < 0.05 and m["auroc_gout_rwr"] > m["auroc_degree_only"]) else "REFUTED", } # ---- H2 配对留出 ---- fr = an["paired_folds_h2"] h2 = paired_block(fr, "auroc_gout_rwr", "auroc_matched_null", "五折留出：4/5 痛风种子 vs 同折度数配对随机种子") h2["difference_ci95_bootstrap"] = boot_ci(fr, "auroc_gout_rwr", "auroc_matched_null") h2b = paired_block(fr, "auroc_gout_rwr", "auroc_degree_only", "五折留出：4/5 痛风种子 vs 只用度数") h2b["difference_ci95_bootstrap"] = boot_ci(fr, "auroc_gout_rwr", "auroc_degree_only") h2["verdict"] = "passed" if (h2["wilcoxon_p_two_sided"] < 0.05 and h2["median_difference"] > 0) else "REFUTED" out["H2_paired_holdout"] = h2 out["H2b_vs_degree_only"] = h2b out["H2_per_seed"] = {} for s in an["params"]["seeds"]: rows = [r for r in fr if r["seed"] == s] d = np.array([r["auroc_gout_rwr"] - r["auroc_matched_null"] for r in rows]) out["H2_per_seed"][str(s)] = {"n_folds": len(rows), "wins": int((d > 0).sum()), "median_difference": float(np.median(d)), "mean_auroc": float(np.mean([r["auroc_gout_rwr"] for r in rows]))} # ---- H3 时间留出 ---- t = an["temporal_h3"] zt = (t["auroc_gout_rwr"] - t["null_mean"]) / t["null_sd"] out["H3_temporal"] = dict(t) out["H3_temporal"].update({ "z_vs_null": float(zt), "delta_vs_null_mean": float(t["auroc_gout_rwr"] - t["null_mean"]), "delta_vs_degree_only": float(t["auroc_gout_rwr"] - t["auroc_degree_only"]), "beats_degree_only": bool(t["auroc_gout_rwr"] > t["auroc_degree_only"]), "verdict": "passed" if (t["empirical_p"] < 0.05 and t["auroc_gout_rwr"] > t["auroc_degree_only"]) else "REFUTED", }) # ---- 敏感性：去掉 textmining 通道 ---- sp_path = os.path.join(RES, "sensitivity_channels.json") if os.path.isfile(sp_path): sc = json.load(open(sp_path, encoding="utf-8")) sc["z_vs_null"] = float((sc["auroc_gout_rwr"] - sc["null_mean"]) / sc["null_sd"]) sc["delta_vs_null_mean"] = float(sc["auroc_gout_rwr"] - sc["null_mean"]) sc["beats_degree_only"] = bool(sc["auroc_gout_rwr"] > sc["auroc_degree_only"]) sc["verdict"] = "passed" if (sc["empirical_p"] < 0.05 and sc["beats_degree_only"]) else "REFUTED" out["sensitivity_no_textmining"] = sc json.dump(out, open(os.path.join(RES, "final_stats.json"), "w"), ensure_ascii=False, indent=1) print(json.dumps({k: (v.get("verdict") if isinstance(v, dict) else v) for k, v in out.items() if isinstance(v, dict) and "verdict" in v}, ensure_ascii=False, indent=1)) print(json.dumps({"H1_auroc": out["H1_main"]["auroc_gout_rwr"], "H1_emp_p": out["H1_main"]["empirical_p"], "H2_wilcoxon_p": out["H2_paired_holdout"]["wilcoxon_p_two_sided"], "H2_median_diff": out["H2_paired_holdout"]["median_difference"], "H2_count": out["H2_paired_holdout"]["count_claim"], "H3_auroc": out["H3_temporal"]["auroc_gout_rwr"], "H3_emp_p": out["H3_temporal"]["empirical_p"]}, ensure_ascii=False, indent=1)) if __name__ == "__main__": main()