Day 26 已在相同 1,267 筆紀錄與五次重複五折下執行三個簡單基準。B0 永遠預測訓練折的多數級數，平均準確率是 0.3844；B2 使用檢傷當下可得的結構化欄位建立序位 Ridge 迴歸（Ridge Regression）基準，平均準確率提高到 0.5376。這個線性模型把較大的輸出分數對應到較不急迫級數，再四捨五入並限制在第一至第五級。Day 27 則確認完整 P0 與六個消融版本仍沒有病患紀錄折外預測，不能加入同一場病患比較。

如果今天只把「0.5376 高於 0.3844」寫成結論，至少會漏掉三個問題：錯誤離參考級數多遠、錯誤朝較不急迫或較急迫方向移動，以及這個差值在重抽樣後是否穩定。更重要的是，準確率較高不等於已經證明臨床更安全。

下圖是概念示意。相同五級輸出要經過多種量測工具檢視，最後仍由臨床人員覆核；圖中的卡片、量尺與轉盤不代表真實病患、正式門檻或本篇實驗數字。

上圖要表達的不是「指標越多越專業」，而是每一個指標回答不同問題。準確率回答完全相符幾筆；平均絕對誤差回答平均差幾級；檢傷不足率回答錯誤方向；信賴區間則提醒我們點估計不是固定真理。

本篇會完成四件事：

1. 實作五級分類、序位距離、方向與高急迫級安全指標。
2. 使用 Day 26 的 6,335 列折外預測，但把同一 record_index 的五次結果綁在一起，不假裝成 6,335 個獨立病患。
3. 對 B2 與 B0 做 2,000 次配對的 record_index 分組重抽樣，報告描述性 95% 百分位區間。
4. 明確保留目前不能計算的候選覆蓋與病患風險－覆蓋曲線，不從最終級數倒推不存在的資料。

本篇的結果是已觀察 Day 26 折外預測的探索性再分析。Day 09 的多重比較策略與安全界值仍未在看結果前鎖定，因此本文不會把區間或 p 值包裝成確認性證據。

本篇會用到的名詞

先鎖定數字方向與資料角色

本篇使用韓國急診檢傷與急迫度分級量表（Korean Triage and Acuity Scale, KTAS）的五級數字。第一級最急迫，第五級最不急迫，因此：

• 預測數字比參考數字大，是檢傷不足方向。
• 預測數字比參考數字小，是檢傷過度方向。
• 預測數字相同，才是完全相符。

例如參考第二級、預測第四級，數字增加兩級，屬於檢傷不足；參考第四級、預測第二級，數字減少兩級，屬於檢傷過度。兩者的絕對距離都是 2，但方向與可能代價不同，所以 MAE 不能取代方向指標。

Day 28 只讀取下列逐列欄位：

公開資料沒有病患識別碼。record_index 只能證明同一資料列跨五次重複要綁在一起，不能證明不同資料列一定來自不同病患。因此原大綱中的「patient-level bootstrap」在本篇修正為「record_index 分組 bootstrap」，公開結果也固定 patient_level_independence_claim_allowed=false。

第一層：完全相符與類別平衡

Accuracy 回答「完全答對幾筆」

準確率（Accuracy）的公式是：

[
Accuracy = \frac{\sum_{i=1}^{n} I(\hat{y}_i = y_i)}{n}
]

其中：

• (n) 是已回答案例數。
• (y_i) 是第 (i) 筆參考級數。
• (\hat{y}_i) 是第 (i) 筆預測級數。
• (I(\cdot)) 是條件成立時為 1、否則為 0 的指示函數。

若十筆案例有六筆級數完全相同，準確率就是 (6/10=0.60)。這個數字不會告訴我們剩下四筆差一級還是差四級，也不會告訴我們錯誤方向。

Balanced Accuracy 讓五級各有一票

第 (k) 級召回率（Recall）是：

[
Recall_k = \frac{TP_k}{N_k}
]

(TP_k) 是參考與預測都為第 (k) 級的筆數，(N_k) 是參考為第 (k) 級的總筆數。平衡準確率（Balanced Accuracy）再把五級召回率平均：

[
Balanced\ Accuracy = \frac{1}{5}\sum_{k=1}^{5} Recall_k
]

假設五級召回率依序為 0.50、0.40、0.80、0.60、0.20，平衡準確率就是 ((0.50+0.40+0.80+0.60+0.20)/5=0.50)。即使第三級樣本最多，它仍只占平均中的五分之一。

巨觀 F1 分數（Macro-averaged F1 Score, Macro F1）則在每一級同時考慮精確率（Precision）與召回率，再把五級 F1 平均。精確率回答「所有被預測成這一級的案例中，有多少真的屬於這一級」；Macro F1 適合補充模型是否把太多其他級數塞進某一級，但一樣沒有直接表示級數距離。

第二層：把五級順序放回指標

MAE 直接回答平均差幾級

平均絕對誤差（Mean Absolute Error, MAE）定義為：

[
MAE = \frac{1}{n}\sum_{i=1}^{n}|\hat{y}_i-y_i|
]

四筆案例若參考為 [1, 2, 3, 4]、預測為 [1, 3, 5, 3]，絕對誤差是 [0, 1, 2, 1]，所以 MAE 為 ((0+1+2+1)/4=1.0) 級。MAE 越低，平均距離越小；但 +1 與 -1 都會變成絕對值 1，所以仍要另外報方向。

QWK 對較遠錯誤給更高權重

二次加權 Kappa（Quadratic Weighted Kappa, QWK）使用混淆矩陣、兩邊類別分布推導的期望矩陣，以及二次距離權重：

[
\kappa_w = 1-\frac{\sum_{i,j}w_{ij}O_{ij}}{\sum_{i,j}w_{ij}E_{ij}},
\qquad
w_{ij}=\frac{(i-j)^2}{(K-1)^2}
]

其中：

• (O_{ij}) 是參考第 (i) 級、預測第 (j) 級的觀察筆數。
• (E_{ij}) 是依參考與預測邊際分布計算的偶然期望筆數。
• (K=5) 是級數數量。
• (w_{ij}) 是錯開距離的二次權重。

五級量表中，差一級的權重是 (1^2/4^2=1/16)，差兩級是 (4/16)，第一級錯到第五級則是 (16/16=1)。因此跨越越多級，分子中的懲罰增加得更快。QWK 同時校正偶然一致，不能解讀成「答對比例」；本篇實作與 scikit-learn 官方文件所述的 quadratic weights 語意一致。scikit-learn：Cohen’s kappa

第三層：方向與高急迫級安全終點

本篇把三個方向指標的分母寫清楚：

[
Undertriage\ Rate = \frac{#(\hat{y}>y)}{n}
]

[
Overtriage\ Rate = \frac{#(\hat{y}<y)}{n}
]

[
Severe\ Undertriage\ Rate =
\frac{#(y\in{1,2}\ \land\ \hat{y}\in{3,4,5})}
{#(y\in{1,2})}
]

前兩式的分母是所有已回答案例；第三式的分母只包含參考第一、二級。假設十筆案例中有四筆參考第一、二級，其中一筆被預測為第三級，嚴重檢傷不足率就是 (1/4=0.25)，不是 (1/10=0.10)。

高急迫側敏感度（High-acuity-side Sensitivity）使用相同分母，但分子改成參考第一、二級且預測仍留在第一、二級。對每一筆高急迫案例而言，它和本篇嚴重檢傷不足事件互補，所以兩率相加為 1；這只是目前二分定義的數學關係，不代表第一級與第二級之間分錯沒有代價。

下圖使用 B2 的 repeat 1、共 1,267 筆 OOF 畫混淆矩陣。請先看對角線，再看右上方的檢傷不足方向與左下方的檢傷過度方向；每格同時保留筆數與該參考級的列百分比。

上圖顯示 repeat 1 的 B2 完全相符 685／1,267 筆，檢傷不足 297 筆，檢傷過度 285 筆，跨兩級以上 73 筆。參考第一、二級共 246 筆，其中 191 筆被分到第三至第五級。這張圖只使用一次 repeat，沒有把五次重複的 6,335 列假裝成獨立病患。

候選層與最終分類層不能互相代替

前 k 名候選覆蓋率（Top-k Candidate Coverage）檢查正確規則或級數是否仍在候選集合中；跨側錯誤（Straddling Error）則描述未覆蓋時，候選是否同時落在參考級數兩側。這兩個指標需要逐筆候選成員與排序。

Day 26 的 OOF 產物只有最終 B0、B2 級數，沒有檢索候選。因此 Day 28 的公開結果保留：

{
  "top_k_candidate_coverage": null,
  "straddling_error": null,
  "must_not_infer_candidates_from_final_prediction": true
}

一個最終預測不能還原曾經有哪些候選。把最終級數旁邊自行補出 {2, 3, 4}，不是分析，而是捏造中間產物。

為什麼要做配對的 record_index 分組 bootstrap

重抽樣（Bootstrap Resampling）從現有抽樣單位有放回抽取，重複計算統計量，觀察估計值如何波動。Efron 與 Tibshirani 的綜述介紹了 bootstrap 對標準誤、信賴區間與複雜統計量的應用；方法仍依賴抽樣單位與資料結構是否合理。Efron 與 Tibshirani：Bootstrap Methods

本篇比較相同紀錄上的 B2 與 B0，所以每次抽到某個 record_index 時，必須同時帶入：

1. 該紀錄的五次 reference_level。
2. 該紀錄的五次 B0_predicted_level。
3. 該紀錄的五次 B2_predicted_level。

每次重抽樣都先在每個 repeat 個別計算指標，再取五次平均。若某個越高越好的指標記為 (M)，原始差值是：

[
\Delta_M = M_{B2}-M_{B0}
]

若 MAE 或檢傷不足率是越低越好，圖中為了統一閱讀方向，會改畫：

[
Effect_{favorable}=M_{B0}-M_{B2}
]

例如 B2 MAE 為 0.519、B0 為 0.695，原始 B2-B0 差是 (-0.176)；轉成「MAE 降幅」後是 (+0.176)，正值表示方向有利於 B2。

2,000 次重抽樣的 95% 百分位區間，取差值分布的第 2.5 與第 97.5 百分位。這個區間描述本資料與本重抽樣設計下的變動，不會修復資料代表性、參考標籤誤差或缺少病患識別碼的限制。臨床預測模型報告指引 TRIPOD+AI 也要求效能估計搭配信賴區間；這是透明報告要求，不是安全認證。TRIPOD+AI statement

實際配對結果：多數指標改善，但方向故事不只一個

下圖把四個主要描述性效果轉成同一方向：正值代表 B2 相對 B0 較有利。請注意圖中刻意寫「探索性」，因為 Day 09 的決策界值尚未鎖定，而且 Day 26 分數已經看過。

上圖顯示四個主要效果的 2,000 次百分位區間都沒有跨 0；這只能稱為本次探索性重抽樣中的穩定方向。完整數字如下，差值一律使用原始 B2-B0，因此越低越好的指標若為負值，才代表 B2 較低。

這張表最值得停下來看的不是最大增益，而是兩個同時成立的方向：B2 的嚴重檢傷不足率比 B0 低 23.09 個百分點，但整體檢傷不足率反而高 4.25 個百分點。原因是 B0 永遠預測第三級：它會把所有第一、二級都分到較不急迫側，也會把大量第四、五級分到較急迫側；B2 修正了一部分高急迫案例，同時在其他級數產生新的較不急迫方向錯誤。

因此「嚴重不足下降」不能改寫成「所有不足都下降」，「Accuracy 上升」也不能改寫成「已證明更安全」。

McNemar、效應量與 Holm 校正放在哪裡

McNemar 檢定（McNemar Test）處理相同案例上的兩個二元結果，只使用兩方法不一致的配對。原始方法就是為相關比例差異提出；本篇使用精確雙尾二項版本，避免把相同案例當成兩組獨立樣本。McNemar 1947

以 repeat 1 的「是否完全相符」為例：

• B2 對、B0 錯：305 筆。
• B0 對、B2 錯：107 筆。
• 配對差值：B2-B0 = +0.1563。
• 不一致配對勝算比：(305/107=2.8505)。

檢定只使用 305 與 107 這兩格；兩方法都對或都錯不會影響 McNemar 的不一致方向。本文也對跨兩級以上錯誤與嚴重檢傷不足建立兩個二元檢查，再用 Holm 校正（Holm Correction）處理這三項探索性檢查。Holm 方法依排序後的 p 值逐步調整，用來控制同一家族的第一類錯誤。Holm 1979

這些 p 值仍不能把分析變成確認性研究。Day 09 沒有在結果前鎖定這三項家族、Holm 策略或臨床安全界值；本篇也只任意指定 repeat 1 作教學性二元檢查。機器學習演算法比較研究曾提醒，重複切分產生的相依性會讓某些常見檢定出現過高第一類錯誤，因此不能把多次交叉驗證分數當成一般獨立樣本直接做 paired t test。Dietterich：比較分類演算法的統計檢定

Coverage 與 selective risk 必須一起看

覆蓋率（Coverage）是回答數占所有符合評估條件案例數的比例：

[
Coverage=\frac{N_{answered}}{N_{eligible}}
]

選擇性風險（Selective Risk）只在已回答案例計算平均損失：

[
Selective\ Risk=
\frac{\sum_{i=1}^{n}\ell(\hat{y}i,y_i)g_i}
{\sum{i=1}^{n}g_i}
]

其中 (g_i=1) 代表系統回答，(g_i=0) 代表拒答；(\ell) 是損失函數。本篇教學曲線使用 0–1 error，所以答錯是 1、答對是 0。選擇性分類研究把這個 coverage 與 risk 的取捨稱為風險－覆蓋關係；它是一般方法框架，不能直接當作本系統的醫療安全證據。El-Yaniv 與 Wiener：Selective Classification

下圖左側使用設定檔內十筆明確標成合成的選擇分數（Selection Score），逐步接受分數較高案例。這個分數只負責排列拒答順序，讀者不能把它自動解讀為正確機率。右側則保留真實狀態：B0 與 B2 都沒有預先鎖定拒答分數，B1 則是零覆蓋。

上圖左側的 score 不是模型機率、臨床信心或 Day 26 病患輸出，只用來示範計算。coverage 從 0.4 增加到 0.5 時，新增案例剛好答對，所以 risk 從 0.25 降到 0.20；這說明曲線不必單調。

B1 的 coverage 是 0，回答數也是 0。此時 selective risk 的分母為 0，正確結果是 null，不是 0。把「完全不回答」寫成「零錯誤」，會讓沒有實用決策覆蓋的系統看起來完美。

Day 28 新增檔案與分工

本篇的設定、核心程式、執行入口、測試與圖片產生器分工如下：

results/public/ 只保存聚合指標、區間、探索性二元檢查與合成曲線。2,000 次逐次重抽樣差值會留在 gitignored 的 results/runs/day-28/<run_id>/，不把 Day 26 逐列病患相關輸出提交到公開結果。

先在自己的專案資料夾建立本篇完整檔案

接下來不會要求你前往任何程式碼網站。請在自己的電腦開啟專案資料夾，依下列順序建立檔案；每個程式碼區塊都是該檔案的完整內容，不含省略號。

本篇沿用 Day 09 的研究方案狀態、Day 26 已完成且留在本機的 6,335 列 OOF 預測，以及 Day 27 的病患消融就緒邊界。Day 28 不重新訓練模型；以下是完整統計契約、evaluation 子套件、runner、九項測試與三張技術圖生成程式。公開摘要、2,000 次重抽樣差值與 run manifest 由執行入口產生，不需要人工填入數字；Image 2 開場概念圖則是文章資產，不是統計程式的輸入。

先從專案根目錄建立需要的資料夾：

mkdir -p configs/evaluation src/triage_rag/evaluation scripts scripts/figures/day-28 tests results/public results/runs/day-28 articles/assets/day-28

如果指令沒有印出訊息是正常的。可用 test -d 資料夾路徑 && echo "資料夾已建立" 驗證單一資料夾。接著使用你熟悉的文字編輯器新增各檔案，把對應區塊完整貼入後儲存。

檔案 1：建立 configs/evaluation/day-28-safety-statistics.json

鎖定五級指標方向、record_index 分組配對 bootstrap、repeat 1 探索性二元檢查、合成 risk–coverage 教學資料、來源雜湊與解讀邊界。

請在文字編輯器建立 configs/evaluation/day-28-safety-statistics.json，貼入以下完整內容並儲存：

{
  "schema_version": 1,
  "experiment_id": "day-28-record-grouped-paired-safety-statistics",
  "scope": "exploratory_reanalysis_of_already_observed_day26_oof_predictions_not_confirmatory_or_clinical_safety_evidence",
  "sources": {
    "day09_protocol_path": "configs/experiments/day-09-research-protocol.json",
    "day09_protocol_sha256": "be9954891978cf654cf0031adaea6fa37d4e940d067b3e02293b3199a9cf0e19",
    "day26_contract_path": "configs/baselines/day-26-fair-baselines.json",
    "day26_contract_sha256": "ff6c36653da05e83276bd1f5ee671069afa5283b19a94e7103228de2a2f88dd6",
    "day26_public_path": "results/public/day-26-simple-baselines.json",
    "day26_public_sha256": "76c4a296511474b93e97a708e776e6709193bec9316665e45a9d2cb15022123a",
    "day27_registry_path": "configs/experiments/day-27-ablation-registry.json",
    "day27_registry_sha256": "40df533396bf6419ebae8c57dbe97ea04533e3dcc16b91ff125fa952294b71cd",
    "day27_public_path": "results/public/day-27-ablation-readiness.json",
    "day27_public_sha256": "53bde664581b6eb0e1d0f103fd06b6b417029d8a2f83b3b18e57edef497dd61e",
    "day26_private_result_glob": "results/runs/day-26/*/day-26-oof-baseline-results.json"
  },
  "input_contract": {
    "row_key": "record_index",
    "repeat_key": "repeat_id",
    "fold_key": "fold_id",
    "reference_field": "reference_level",
    "comparator_prediction_field": "B0_predicted_level",
    "candidate_prediction_field": "B2_predicted_level",
    "valid_levels": [1, 2, 3, 4, 5],
    "record_count": 1267,
    "repeat_count": 5,
    "fold_count": 5,
    "prediction_row_count": 6335,
    "patient_identifier_available": false,
    "independence_unit_known": false
  },
  "metric_policy": {
    "higher_is_better": [
      "accuracy",
      "balanced_accuracy",
      "macro_f1",
      "quadratic_weighted_kappa",
      "high_acuity_side_sensitivity"
    ],
    "lower_is_better": [
      "mean_absolute_error",
      "undertriage_rate",
      "overtriage_rate",
      "severe_undertriage_rate",
      "two_or_more_level_error_rate"
    ],
    "primary_descriptive_metrics": [
      "balanced_accuracy",
      "macro_f1",
      "mean_absolute_error",
      "severe_undertriage_rate"
    ],
    "level_direction": "larger_prediction_number_is_less_urgent",
    "undertriage_definition": "predicted_level_greater_than_reference_level",
    "overtriage_definition": "predicted_level_less_than_reference_level",
    "severe_undertriage_definition": "reference_in_1_2_and_prediction_in_3_4_5",
    "high_acuity_side_sensitivity_definition": "reference_in_1_2_and_prediction_in_1_2",
    "per_class_recall_required": true,
    "numerators_and_denominators_required": true
  },
  "bootstrap_policy": {
    "method": "paired_record_index_grouped_percentile_bootstrap",
    "resampling_unit": "record_index",
    "preserve_all_repeats_for_sampled_record": true,
    "preserve_method_pairing": true,
    "replications": 2000,
    "seed": 20260828,
    "confidence_level": 0.95,
    "point_estimate": "mean_of_five_repeat_specific_metrics",
    "interval": "percentile",
    "patient_level_claim_allowed": false
  },
  "paired_binary_policy": {
    "status": "exploratory_only",
    "analysis_repeat_id": 1,
    "tests": [
      {
        "id": "exact_correctness",
        "event": "prediction_equals_reference",
        "favorable_event": true
      },
      {
        "id": "two_or_more_level_error",
        "event": "absolute_error_at_least_2",
        "favorable_event": false
      },
      {
        "id": "severe_undertriage",
        "event": "reference_in_1_2_and_prediction_in_3_4_5",
        "eligible_reference_levels": [1, 2],
        "favorable_event": false
      }
    ],
    "test": "exact_two_sided_mcnemar_binomial",
    "multiplicity_adjustment": "holm_within_three_predeclared_day28_binary_checks",
    "confirmatory_family_claim_allowed": false
  },
  "selective_prediction_policy": {
    "coverage_definition": "answered_count_divided_by_all_eligible_records",
    "selective_risk_definition": "mean_zero_one_error_among_answered_records",
    "patient_risk_coverage_curve_status": "not_computable_without_prelocked_selection_score_and_nonzero_coverage",
    "b1_zero_coverage_risk": null,
    "teaching_example_is_synthetic": true,
    "teaching_example": [
      {"case_id": "S01", "reference_level": 2, "predicted_level": 2, "selection_score": 0.98},
      {"case_id": "S02", "reference_level": 3, "predicted_level": 3, "selection_score": 0.94},
      {"case_id": "S03", "reference_level": 4, "predicted_level": 4, "selection_score": 0.90},
      {"case_id": "S04", "reference_level": 3, "predicted_level": 2, "selection_score": 0.84},
      {"case_id": "S05", "reference_level": 5, "predicted_level": 5, "selection_score": 0.79},
      {"case_id": "S06", "reference_level": 2, "predicted_level": 3, "selection_score": 0.72},
      {"case_id": "S07", "reference_level": 4, "predicted_level": 4, "selection_score": 0.65},
      {"case_id": "S08", "reference_level": 1, "predicted_level": 3, "selection_score": 0.57},
      {"case_id": "S09", "reference_level": 3, "predicted_level": 3, "selection_score": 0.49},
      {"case_id": "S10", "reference_level": 5, "predicted_level": 4, "selection_score": 0.41}
    ]
  },
  "candidate_evaluation_policy": {
    "top_k_candidate_coverage_status": "not_available_in_day26_prediction_artifact",
    "straddling_error_status": "not_available_in_day26_prediction_artifact",
    "must_not_infer_candidates_from_final_prediction": true
  },
  "claim_policy": {
    "day09_status_must_remain_planning": true,
    "day09_multiplicity_was_not_locked_before_results": true,
    "day09_safety_margins_were_not_locked_before_results": true,
    "day26_results_already_observed": true,
    "confirmatory_p_values_allowed": false,
    "clinical_safety_or_deployment_claim_allowed": false,
    "confidence_interval_crossing_zero_interpretation": "compatible_with_both_effect_directions_under_this_exploratory_resampling_analysis",
    "confidence_interval_not_crossing_zero_interpretation": "descriptive_stability_signal_only_not_confirmatory_proof"
  },
  "outputs": {
    "public_summary_path": "results/public/day-28-safety-statistics.json",
    "run_output_root": "results/runs/day-28",
    "full_result_filename": "day-28-bootstrap-replicates.json"
  }
}

儲存後先確認檔名與相對路徑完全一致，再繼續建立下一個檔案。

檔案 2：建立 src/triage_rag/evaluation/__init__.py

建立 evaluation 子套件入口，公開 Day 28 指標、驗證、配對重抽樣、McNemar、Holm 與 risk–coverage 函式。

請在文字編輯器建立 src/triage_rag/evaluation/__init__.py，貼入以下完整內容並儲存：

"""Evaluation utilities for ordinal triage safety analyses."""

from .safety_statistics import (
    EvaluationContractError,
    classification_metrics,
    exact_mcnemar,
    holm_adjust,
    paired_record_bootstrap,
    risk_coverage_curve,
    validate_day28_inputs,
)

__all__ = [
    "EvaluationContractError",
    "classification_metrics",
    "exact_mcnemar",
    "holm_adjust",
    "paired_record_bootstrap",
    "risk_coverage_curve",
    "validate_day28_inputs",
]

儲存後先確認檔名與相對路徑完全一致，再繼續建立下一個檔案。

檔案 3：建立 src/triage_rag/evaluation/safety_statistics.py

實作 Accuracy、Balanced Accuracy、Macro F1、QWK、MAE、方向與安全分母，並驗證 6,335 列配對後執行 record_index 分組 bootstrap。

請在文字編輯器建立 src/triage_rag/evaluation/safety_statistics.py，貼入以下完整內容並儲存：

"""Ordinal safety metrics and record-grouped paired statistics for Day 28."""

from __future__ import annotations

import math
from collections import Counter, defaultdict
from typing import Any, Iterable, Mapping, Sequence

import numpy as np


JsonObject = dict[str, Any]


class EvaluationContractError(ValueError):
    """Raised when a Day 28 input violates the locked evaluation contract."""


def _safe_rate(numerator: int, denominator: int) -> float | None:
    return round(numerator / denominator, 6) if denominator else None


def _as_levels(values: Sequence[int], levels: Sequence[int], name: str) -> np.ndarray:
    array = np.asarray(values, dtype=int)
    if array.ndim != 1 or len(array) == 0:
        raise EvaluationContractError(f"{name} 必須是一維且非空")
    if not set(array.tolist()).issubset(set(levels)):
        raise EvaluationContractError(f"{name} 包含契約外級數")
    return array


def _quadratic_weighted_kappa(confusion: np.ndarray) -> float | None:
    count = int(confusion.sum())
    if count == 0:
        return None
    level_count = confusion.shape[0]
    indices = np.arange(level_count, dtype=float)
    weights = (indices[:, None] - indices[None, :]) ** 2
    if level_count > 1:
        weights /= float((level_count - 1) ** 2)
    expected = np.outer(confusion.sum(axis=1), confusion.sum(axis=0)) / count
    expected_disagreement = float((weights * expected).sum())
    if expected_disagreement == 0:
        return None
    observed_disagreement = float((weights * confusion).sum())
    return round(1.0 - observed_disagreement / expected_disagreement, 6)


def classification_metrics(
    reference: Sequence[int] | np.ndarray,
    predicted: Sequence[int] | np.ndarray,
    levels: Sequence[int],
) -> JsonObject:
    """Compute exact, ordinal, direction, class-recall, and safety metrics."""

    y_true = _as_levels(reference, levels, "reference")
    y_pred = _as_levels(predicted, levels, "predicted")
    if len(y_true) != len(y_pred):
        raise EvaluationContractError("reference 與 predicted 長度不同")

    level_to_index = {int(level): index for index, level in enumerate(levels)}
    confusion = np.zeros((len(levels), len(levels)), dtype=int)
    true_index = np.fromiter(
        (level_to_index[int(value)] for value in y_true), dtype=int, count=len(y_true)
    )
    pred_index = np.fromiter(
        (level_to_index[int(value)] for value in y_pred), dtype=int, count=len(y_pred)
    )
    np.add.at(confusion, (true_index, pred_index), 1)

    recalls: dict[str, float | None] = {}
    f1_values: list[float] = []
    for level in levels:
        index = level_to_index[int(level)]
        true_positive = int(confusion[index, index])
        actual_count = int(confusion[index, :].sum())
        predicted_count = int(confusion[:, index].sum())
        recalls[str(level)] = _safe_rate(true_positive, actual_count)
        f1_denominator = actual_count + predicted_count
        f1_values.append(
            (2 * true_positive / f1_denominator) if f1_denominator else 0.0
        )

    exact_count = int((y_true == y_pred).sum())
    under_count = int((y_pred > y_true).sum())
    over_count = int((y_pred < y_true).sum())
    severe_reference = np.isin(y_true, [1, 2])
    severe_under = severe_reference & np.isin(y_pred, [3, 4, 5])
    high_acuity_side_hit = severe_reference & np.isin(y_pred, [1, 2])
    absolute_error = np.abs(y_pred - y_true)
    severe_reference_count = int(severe_reference.sum())
    return {
        "answered_count": len(y_true),
        "exact_count": exact_count,
        "accuracy": _safe_rate(exact_count, len(y_true)),
        "balanced_accuracy": round(
            sum(float(value or 0.0) for value in recalls.values()) / len(levels), 6
        ),
        "macro_f1": round(float(np.mean(f1_values)), 6),
        "quadratic_weighted_kappa": _quadratic_weighted_kappa(confusion),
        "absolute_error_sum": int(absolute_error.sum()),
        "mean_absolute_error": round(float(absolute_error.mean()), 6),
        "undertriage_count": under_count,
        "undertriage_rate": _safe_rate(under_count, len(y_true)),
        "overtriage_count": over_count,
        "overtriage_rate": _safe_rate(over_count, len(y_true)),
        "severe_reference_count": severe_reference_count,
        "severe_undertriage_count": int(severe_under.sum()),
        "severe_undertriage_rate": _safe_rate(
            int(severe_under.sum()), severe_reference_count
        ),
        "high_acuity_side_hit_count": int(high_acuity_side_hit.sum()),
        "high_acuity_side_sensitivity": _safe_rate(
            int(high_acuity_side_hit.sum()), severe_reference_count
        ),
        "two_or_more_level_error_count": int((absolute_error >= 2).sum()),
        "two_or_more_level_error_rate": _safe_rate(
            int((absolute_error >= 2).sum()), len(y_true)
        ),
        "per_class_recall": recalls,
        "confusion_matrix_rows_reference_columns_prediction": confusion.tolist(),
    }


def validate_day28_inputs(
    rows: Sequence[Mapping[str, Any]],
    contract: Mapping[str, Any],
    day09: Mapping[str, Any],
    day26_public: Mapping[str, Any],
    day27_public: Mapping[str, Any],
) -> JsonObject:
    """Validate OOF completeness, pairing, reference stability, and claim boundaries."""

    data = contract["input_contract"]
    required_fields = [
        data["row_key"],
        data["repeat_key"],
        data["fold_key"],
        data["reference_field"],
        data["comparator_prediction_field"],
        data["candidate_prediction_field"],
    ]
    if len(rows) != int(data["prediction_row_count"]):
        raise EvaluationContractError("OOF 預測列數與 Day 28 契約不一致")
    for index, row in enumerate(rows):
        missing = [field for field in required_fields if field not in row]
        if missing:
            raise EvaluationContractError(f"OOF 第 {index} 列缺少欄位：{missing}")

    levels = set(int(value) for value in data["valid_levels"])
    expected_repeats = set(range(1, int(data["repeat_count"]) + 1))
    expected_folds = set(range(1, int(data["fold_count"]) + 1))
    by_record: dict[int, list[Mapping[str, Any]]] = defaultdict(list)
    pair_keys: set[tuple[int, int]] = set()
    fold_sets: dict[int, set[int]] = defaultdict(set)
    for row in rows:
        record = int(row[data["row_key"]])
        repeat = int(row[data["repeat_key"]])
        fold = int(row[data["fold_key"]])
        pair_key = (record, repeat)
        if pair_key in pair_keys:
            raise EvaluationContractError(f"record/repeat 重複：{pair_key}")
        pair_keys.add(pair_key)
        by_record[record].append(row)
        fold_sets[repeat].add(fold)
        for field in [
            data["reference_field"],
            data["comparator_prediction_field"],
            data["candidate_prediction_field"],
        ]:
            if int(row[field]) not in levels:
                raise EvaluationContractError(f"{field} 出現契約外級數")

    if len(by_record) != int(data["record_count"]):
        raise EvaluationContractError("record_index 數量與契約不一致")
    for record, record_rows in by_record.items():
        repeats = {int(row[data["repeat_key"]]) for row in record_rows}
        references = {int(row[data["reference_field"]]) for row in record_rows}
        if repeats != expected_repeats:
            raise EvaluationContractError(f"record {record} 的 repeat 不完整")
        if len(references) != 1:
            raise EvaluationContractError(f"record {record} 跨 repeat 的參考級數不一致")
    if set(fold_sets) != expected_repeats or any(
        folds != expected_folds for folds in fold_sets.values()
    ):
        raise EvaluationContractError("repeat/fold 組合不完整")

    if day09.get("status") != "planning":
        raise EvaluationContractError("Day 09 必須維持 planning，不能偽裝成已鎖定方案")
    if day09["analysis_policy"].get("multiplicity_adjustment") != "to-be-locked-before-test":
        raise EvaluationContractError("Day 09 多重比較狀態與 Day 28 邊界不一致")
    if day09["analysis_policy"].get("safety_margins") != "to-be-locked-with-clinical-review":
        raise EvaluationContractError("Day 09 安全界值狀態與 Day 28 邊界不一致")
    if day26_public["data_and_split_audit"].get("patient_identifier_available") is not False:
        raise EvaluationContractError("Day 26 病患識別碼狀態與契約不一致")
    if day27_public["contract_checks"].get("p0_patient_oof_ready") is not False:
        raise EvaluationContractError("Day 27 P0 狀態不應被 Day 28 改寫")
    return {
        "prediction_row_count": len(rows),
        "record_count": len(by_record),
        "repeat_count": len(expected_repeats),
        "fold_count": len(expected_folds),
        "one_prediction_per_record_per_repeat": True,
        "reference_stable_across_repeats": True,
        "methods_paired_on_every_record_repeat": True,
        "resampling_unit": data["row_key"],
        "patient_identifier_available": False,
        "patient_level_independence_claim_allowed": False,
        "day09_status": "planning",
        "confirmatory_analysis_allowed": False,
        "p0_patient_oof_ready": False,
    }


def _metric_average(
    reference: np.ndarray,
    predicted: np.ndarray,
    levels: Sequence[int],
    metrics: Sequence[str],
) -> dict[str, float]:
    if reference.shape != predicted.shape or reference.ndim != 2:
        raise EvaluationContractError("重抽樣矩陣必須同為 record × repeat")
    values: dict[str, list[float]] = {metric: [] for metric in metrics}
    for repeat_index in range(reference.shape[1]):
        result = classification_metrics(
            reference[:, repeat_index], predicted[:, repeat_index], levels
        )
        for metric in metrics:
            value = result.get(metric)
            if value is None:
                raise EvaluationContractError(f"{metric} 在本次重抽樣未定義")
            values[metric].append(float(value))
    return {metric: float(np.mean(items)) for metric, items in values.items()}


def paired_record_bootstrap(
    rows: Sequence[Mapping[str, Any]],
    contract: Mapping[str, Any],
    *,
    replications: int | None = None,
) -> tuple[JsonObject, JsonObject]:
    """Bootstrap record_index groups while preserving repeats and method pairing."""

    data = contract["input_contract"]
    levels = [int(value) for value in data["valid_levels"]]
    row_key = str(data["row_key"])
    repeat_key = str(data["repeat_key"])
    reference_field = str(data["reference_field"])
    comparator_field = str(data["comparator_prediction_field"])
    candidate_field = str(data["candidate_prediction_field"])
    repeat_ids = list(range(1, int(data["repeat_count"]) + 1))
    record_ids = sorted({int(row[row_key]) for row in rows})
    record_position = {record: index for index, record in enumerate(record_ids)}
    repeat_position = {repeat: index for index, repeat in enumerate(repeat_ids)}
    shape = (len(record_ids), len(repeat_ids))
    reference = np.empty(shape, dtype=int)
    comparator = np.empty(shape, dtype=int)
    candidate = np.empty(shape, dtype=int)
    for row in rows:
        position = (
            record_position[int(row[row_key])],
            repeat_position[int(row[repeat_key])],
        )
        reference[position] = int(row[reference_field])
        comparator[position] = int(row[comparator_field])
        candidate[position] = int(row[candidate_field])

    metric_policy = contract["metric_policy"]
    metrics = [
        *metric_policy["higher_is_better"],
        *metric_policy["lower_is_better"],
    ]
    metrics = list(dict.fromkeys(metrics))
    higher = set(metric_policy["higher_is_better"])
    comparator_point = _metric_average(reference, comparator, levels, metrics)
    candidate_point = _metric_average(reference, candidate, levels, metrics)

    policy = contract["bootstrap_policy"]
    replicate_count = int(replications or policy["replications"])
    if replicate_count <= 0:
        raise EvaluationContractError("bootstrap replications 必須大於零")
    rng = np.random.default_rng(int(policy["seed"]))
    raw_differences: dict[str, list[float]] = {metric: [] for metric in metrics}
    favorable_effects: dict[str, list[float]] = {metric: [] for metric in metrics}
    for _ in range(replicate_count):
        sample = rng.integers(0, len(record_ids), size=len(record_ids))
        comparator_sample = _metric_average(
            reference[sample], comparator[sample], levels, metrics
        )
        candidate_sample = _metric_average(
            reference[sample], candidate[sample], levels, metrics
        )
        for metric in metrics:
            raw = candidate_sample[metric] - comparator_sample[metric]
            favorable = raw if metric in higher else -raw
            raw_differences[metric].append(raw)
            favorable_effects[metric].append(favorable)

    alpha = 1.0 - float(policy["confidence_level"])
    public_metrics: JsonObject = {}
    for metric in metrics:
        raw_point = candidate_point[metric] - comparator_point[metric]
        favorable_point = raw_point if metric in higher else -raw_point
        raw_array = np.asarray(raw_differences[metric], dtype=float)
        favorable_array = np.asarray(favorable_effects[metric], dtype=float)
        public_metrics[metric] = {
            "higher_is_better": metric in higher,
            "B0_point_estimate": round(comparator_point[metric], 6),
            "B2_point_estimate": round(candidate_point[metric], 6),
            "difference_B2_minus_B0": round(raw_point, 6),
            "difference_B2_minus_B0_percentile_95_ci": [
                round(float(np.quantile(raw_array, alpha / 2)), 6),
                round(float(np.quantile(raw_array, 1 - alpha / 2)), 6),
            ],
            "effect_in_favorable_direction": round(favorable_point, 6),
            "effect_in_favorable_direction_percentile_95_ci": [
                round(float(np.quantile(favorable_array, alpha / 2)), 6),
                round(float(np.quantile(favorable_array, 1 - alpha / 2)), 6),
            ],
        }

    public = {
        "method": policy["method"],
        "resampling_unit": policy["resampling_unit"],
        "patient_level_claim_allowed": False,
        "preserved_repeat_count_per_sampled_record": len(repeat_ids),
        "replications": replicate_count,
        "seed": int(policy["seed"]),
        "confidence_level": float(policy["confidence_level"]),
        "interval": policy["interval"],
        "point_estimate": policy["point_estimate"],
        "metrics": public_metrics,
    }
    private = {
        "record_ids": record_ids,
        "raw_differences": {
            metric: [round(value, 8) for value in values]
            for metric, values in raw_differences.items()
        },
        "favorable_effects": {
            metric: [round(value, 8) for value in values]
            for metric, values in favorable_effects.items()
        },
    }
    return public, private


def exact_mcnemar(
    comparator_event: Sequence[bool] | np.ndarray,
    candidate_event: Sequence[bool] | np.ndarray,
) -> JsonObject:
    """Return an exact two-sided McNemar binomial test and discordant counts."""

    comparator = np.asarray(comparator_event, dtype=bool)
    candidate = np.asarray(candidate_event, dtype=bool)
    if comparator.ndim != 1 or comparator.shape != candidate.shape or len(comparator) == 0:
        raise EvaluationContractError("McNemar 需要等長且非空的一維配對事件")
    neither = int((~comparator & ~candidate).sum())
    comparator_only = int((comparator & ~candidate).sum())
    candidate_only = int((~comparator & candidate).sum())
    both = int((comparator & candidate).sum())
    discordant = comparator_only + candidate_only
    if discordant == 0:
        p_value = 1.0
    else:
        tail = sum(
            math.comb(discordant, index) for index in range(min(comparator_only, candidate_only) + 1)
        ) / (2**discordant)
        p_value = min(1.0, 2.0 * tail)
    odds_ratio = (
        round(candidate_only / comparator_only, 6) if comparator_only else None
    )
    return {
        "neither_event": neither,
        "comparator_only_event": comparator_only,
        "candidate_only_event": candidate_only,
        "both_event": both,
        "discordant_pair_count": discordant,
        "candidate_to_comparator_discordant_odds_ratio": odds_ratio,
        "exact_two_sided_p_value": p_value,
    }


def holm_adjust(p_values: Mapping[str, float]) -> dict[str, float]:
    """Adjust one declared family of p-values with Holm's step-down procedure."""

    if not p_values:
        raise EvaluationContractError("Holm 校正至少需要一個 p-value")
    for name, value in p_values.items():
        if not math.isfinite(float(value)) or not 0 <= float(value) <= 1:
            raise EvaluationContractError(f"{name} 的 p-value 超出 [0, 1]")
    ordered = sorted(p_values.items(), key=lambda item: (float(item[1]), item[0]))
    adjusted: dict[str, float] = {}
    running = 0.0
    count = len(ordered)
    for rank, (name, value) in enumerate(ordered):
        candidate = min(1.0, (count - rank) * float(value))
        running = max(running, candidate)
        adjusted[name] = running
    return adjusted


def risk_coverage_curve(
    reference: Sequence[int],
    predicted: Sequence[int],
    selection_scores: Sequence[float],
) -> list[JsonObject]:
    """Build a tie-aware empirical zero-one risk–coverage curve."""

    if not (len(reference) == len(predicted) == len(selection_scores)) or len(reference) == 0:
        raise EvaluationContractError("risk-coverage 輸入必須等長且非空")
    rows = []
    for actual, guess, score in zip(reference, predicted, selection_scores, strict=True):
        numeric_score = float(score)
        if not math.isfinite(numeric_score):
            raise EvaluationContractError("selection score 必須是有限數值")
        rows.append((numeric_score, int(actual) != int(guess)))
    rows.sort(key=lambda item: item[0], reverse=True)
    total = len(rows)
    accepted = 0
    errors = 0
    curve: list[JsonObject] = []
    for score in sorted({item[0] for item in rows}, reverse=True):
        tied = [item for item in rows if item[0] == score]
        accepted += len(tied)
        errors += sum(int(item[1]) for item in tied)
        curve.append(
            {
                "minimum_selection_score": round(score, 6),
                "answered_count": accepted,
                "error_count": errors,
                "coverage": round(accepted / total, 6),
                "selective_risk": round(errors / accepted, 6),
            }
        )
    return curve


def direction_counts(
    reference: Iterable[int], predicted: Iterable[int]
) -> dict[str, int]:
    """Small helper used by tests and reporting code."""

    counter: Counter[str] = Counter()
    for actual, guess in zip(reference, predicted, strict=True):
        if guess == actual:
            counter["exact"] += 1
        elif guess > actual:
            counter["undertriage"] += 1
        else:
            counter["overtriage"] += 1
    return {key: counter[key] for key in ["exact", "undertriage", "overtriage"]}

儲存後先確認檔名與相對路徑完全一致，再繼續建立下一個檔案。

檔案 4：建立 scripts/run_day28_safety_statistics.py

核對 Day 09、26、27 上游雜湊，找出最新相容 Day 26 OOF，產生公開聚合、本機 bootstrap replicates 與 run manifest。

請在文字編輯器建立 scripts/run_day28_safety_statistics.py，貼入以下完整內容並儲存：

#!/usr/bin/env python3
"""Run Day 28 exploratory ordinal safety and paired uncertainty analyses."""

from __future__ import annotations

import argparse
import glob
import platform
import sys
from datetime import datetime, timezone
from pathlib import Path
from typing import Any

from triage_rag.evaluation.safety_statistics import (
    classification_metrics,
    exact_mcnemar,
    holm_adjust,
    paired_record_bootstrap,
    risk_coverage_curve,
    validate_day28_inputs,
)
from triage_rag.reproducibility import (
    canonical_json_bytes,
    file_record,
    git_state,
    installed_versions,
    load_json,
    sha256_bytes,
    sha256_file,
    write_json,
)


PROJECT_ROOT = Path(__file__).resolve().parents[1]
DEFAULT_CONTRACT = "configs/evaluation/day-28-safety-statistics.json"


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="執行 Day 28 序位安全指標、record_index 分組配對 bootstrap 與探索性 McNemar。"
    )
    parser.add_argument("--contract", default=DEFAULT_CONTRACT)
    parser.add_argument(
        "--oof",
        default=None,
        help="指定 Day 26 完整 OOF JSON；未指定時使用排序後最後一份相容產物。",
    )
    parser.add_argument(
        "--bootstrap-replications",
        type=int,
        default=None,
        help="只供測試或快速重跑覆寫；正式結果使用契約的 2000 次。",
    )
    return parser.parse_args()


def _validate_source_hashes(contract: dict[str, Any]) -> list[str]:
    paths: list[str] = []
    for key, relative_path in contract["sources"].items():
        if not key.endswith("_path"):
            continue
        hash_key = f"{key[:-5]}_sha256"
        if hash_key not in contract["sources"]:
            raise ValueError(f"{key} 缺少配對 SHA-256")
        observed = sha256_file(PROJECT_ROOT / relative_path)
        expected = contract["sources"][hash_key]
        if observed != expected:
            raise ValueError(f"{relative_path} SHA-256 不符：{observed}")
        paths.append(relative_path)
    return paths


def _resolve_oof_path(contract: dict[str, Any], requested: str | None) -> Path:
    if requested:
        candidate = (PROJECT_ROOT / requested).resolve()
        if PROJECT_ROOT.resolve() not in candidate.parents:
            raise ValueError("--oof 必須位於專案資料夾內")
        if not candidate.is_file():
            raise FileNotFoundError(f"找不到 OOF：{candidate}")
        return candidate
    pattern = str(PROJECT_ROOT / contract["sources"]["day26_private_result_glob"])
    candidates = [Path(path) for path in sorted(glob.glob(pattern))]
    if not candidates:
        raise FileNotFoundError(
            "找不到 Day 26 完整 OOF；請先執行 poetry run python scripts/run_day26_baselines.py"
        )
    for candidate in reversed(candidates):
        payload = load_json(candidate)
        if (
            payload.get("schema_version") == 1
            and payload.get("experiment_id")
            == "day-26-locked-oof-simple-baselines"
        ):
            return candidate
    raise ValueError("找到 Day 26 檔案，但沒有相容的完整 OOF 產物")


def _per_repeat_metrics(
    rows: list[dict[str, Any]], contract: dict[str, Any]
) -> list[dict[str, Any]]:
    data = contract["input_contract"]
    levels = [int(value) for value in data["valid_levels"]]
    output: list[dict[str, Any]] = []
    for repeat_id in range(1, int(data["repeat_count"]) + 1):
        selected = sorted(
            (row for row in rows if int(row[data["repeat_key"]]) == repeat_id),
            key=lambda row: int(row[data["row_key"]]),
        )
        reference = [int(row[data["reference_field"]]) for row in selected]
        output.append(
            {
                "repeat_id": repeat_id,
                "record_count": len(selected),
                "B0": classification_metrics(
                    reference,
                    [int(row[data["comparator_prediction_field"]]) for row in selected],
                    levels,
                ),
                "B2": classification_metrics(
                    reference,
                    [int(row[data["candidate_prediction_field"]]) for row in selected],
                    levels,
                ),
            }
        )
    return output


def _paired_binary_results(
    rows: list[dict[str, Any]], contract: dict[str, Any]
) -> dict[str, Any]:
    data = contract["input_contract"]
    policy = contract["paired_binary_policy"]
    repeat_id = int(policy["analysis_repeat_id"])
    selected = sorted(
        (row for row in rows if int(row[data["repeat_key"]]) == repeat_id),
        key=lambda row: int(row[data["row_key"]]),
    )
    outputs: dict[str, Any] = {}
    p_values: dict[str, float] = {}
    for test in policy["tests"]:
        test_id = str(test["id"])
        eligible = selected
        if "eligible_reference_levels" in test:
            allowed = {int(value) for value in test["eligible_reference_levels"]}
            eligible = [
                row for row in selected if int(row[data["reference_field"]]) in allowed
            ]
        reference = [int(row[data["reference_field"]]) for row in eligible]
        comparator_prediction = [
            int(row[data["comparator_prediction_field"]]) for row in eligible
        ]
        candidate_prediction = [
            int(row[data["candidate_prediction_field"]]) for row in eligible
        ]
        if test_id == "exact_correctness":
            comparator_event = [a == b for a, b in zip(reference, comparator_prediction)]
            candidate_event = [a == b for a, b in zip(reference, candidate_prediction)]
        elif test_id == "two_or_more_level_error":
            comparator_event = [abs(a - b) >= 2 for a, b in zip(reference, comparator_prediction)]
            candidate_event = [abs(a - b) >= 2 for a, b in zip(reference, candidate_prediction)]
        elif test_id == "severe_undertriage":
            comparator_event = [prediction >= 3 for prediction in comparator_prediction]
            candidate_event = [prediction >= 3 for prediction in candidate_prediction]
        else:
            raise ValueError(f"未知 paired binary test：{test_id}")
        result = exact_mcnemar(comparator_event, candidate_event)
        comparator_count = sum(comparator_event)
        candidate_count = sum(candidate_event)
        result.update(
            {
                "eligible_count": len(eligible),
                "event": test["event"],
                "favorable_event": bool(test["favorable_event"]),
                "B0_event_count": comparator_count,
                "B2_event_count": candidate_count,
                "B0_event_rate": round(comparator_count / len(eligible), 6),
                "B2_event_rate": round(candidate_count / len(eligible), 6),
                "event_rate_difference_B2_minus_B0": round(
                    (candidate_count - comparator_count) / len(eligible), 6
                ),
            }
        )
        outputs[test_id] = result
        p_values[test_id] = float(result["exact_two_sided_p_value"])
    adjusted = holm_adjust(p_values)
    for test_id, value in adjusted.items():
        outputs[test_id]["holm_adjusted_p_value"] = value
    return {
        "status": policy["status"],
        "analysis_repeat_id": repeat_id,
        "test": policy["test"],
        "multiplicity_adjustment": policy["multiplicity_adjustment"],
        "confirmatory_family_claim_allowed": False,
        "tests": outputs,
    }


def _synthetic_risk_coverage(contract: dict[str, Any]) -> dict[str, Any]:
    policy = contract["selective_prediction_policy"]
    rows = policy["teaching_example"]
    curve = risk_coverage_curve(
        [int(row["reference_level"]) for row in rows],
        [int(row["predicted_level"]) for row in rows],
        [float(row["selection_score"]) for row in rows],
    )
    return {
        "is_synthetic_teaching_example": True,
        "case_count": len(rows),
        "selection_score_is_model_probability": False,
        "curve": curve,
        "interpretation": "只示範先固定 selection score，再逐步降低 coverage 並重算已回答錯誤率；不是病患結果。",
    }


def main() -> int:
    args = parse_args()
    contract_path = PROJECT_ROOT / args.contract
    contract = load_json(contract_path)
    if contract.get("schema_version") != 1:
        raise ValueError("目前只支援 Day 28 schema_version=1")
    source_paths = _validate_source_hashes(contract)
    oof_path = _resolve_oof_path(contract, args.oof)
    oof_payload = load_json(oof_path)
    rows = list(oof_payload["oof_predictions"])
    sources = contract["sources"]
    day09 = load_json(PROJECT_ROOT / sources["day09_protocol_path"])
    day26_public = load_json(PROJECT_ROOT / sources["day26_public_path"])
    day27_public = load_json(PROJECT_ROOT / sources["day27_public_path"])
    validation = validate_day28_inputs(
        rows, contract, day09, day26_public, day27_public
    )
    per_repeat = _per_repeat_metrics(rows, contract)
    bootstrap_public, bootstrap_private = paired_record_bootstrap(
        rows, contract, replications=args.bootstrap_replications
    )
    paired_binary = _paired_binary_results(rows, contract)
    synthetic_curve = _synthetic_risk_coverage(contract)

    public = {
        "schema_version": 1,
        "experiment_id": contract["experiment_id"],
        "scope": contract["scope"],
        "analysis_status": {
            "exploratory_descriptive_only": True,
            "confirmatory_analysis": False,
            "clinical_safety_or_deployment_claim_allowed": False,
            "reasons": [
                "Day 09 仍是 planning，多重校正與安全界值未在結果前鎖定。",
                "Day 26 OOF 結果已被觀察，本篇是既有結果的探索性再分析。",
                "資料只有 record_index，沒有病患識別碼，不能排除同一病患重複紀錄。",
                "完整 P0 與消融 arms 沒有病患 OOF，因此本篇只比較 B2 與 B0。",
            ],
        },
        "input_audit": validation,
        "source_oof_artifact": {
            "path": str(oof_path.relative_to(PROJECT_ROOT)),
            "sha256": sha256_file(oof_path),
            "contains_patient_rows": True,
            "published_in_public_summary": False,
        },
        "metric_definitions": contract["metric_policy"],
        "per_repeat_metrics": per_repeat,
        "paired_record_bootstrap": bootstrap_public,
        "paired_binary_checks": paired_binary,
        "candidate_evaluation": {
            **contract["candidate_evaluation_policy"],
            "top_k_candidate_coverage": None,
            "straddling_error": None,
        },
        "selective_prediction": {
            "coverage_definition": contract["selective_prediction_policy"]["coverage_definition"],
            "selective_risk_definition": contract["selective_prediction_policy"]["selective_risk_definition"],
            "patient_result": {
                "B0": {"coverage": 1.0, "risk_coverage_curve": None},
                "B1": {"coverage": 0.0, "selective_risk": None, "risk_coverage_curve": None},
                "B2": {"coverage": 1.0, "risk_coverage_curve": None},
                "status": contract["selective_prediction_policy"]["patient_risk_coverage_curve_status"],
            },
            "teaching_example": synthetic_curve,
        },
        "claim_policy": contract["claim_policy"],
        "limitations": [
            "Percentile bootstrap 描述這份資料與重抽樣設計下的變動，不修復資料代表性、標籤誤差或未識別群聚。",
            "五次重複的預測共享相同 1,267 筆紀錄；程式按 record_index 成組抽樣並保留五次預測，不把 6,335 列假裝成獨立病患。",
            "McNemar 只在契約指定的 repeat 1 做三項探索性二元檢查；Holm 校正不會把事後分析變成確認性研究。",
            "Day 26 沒有候選集合與預先鎖定的拒答分數，不能由最終級數倒推出 Top-k coverage、straddling error 或病患 risk-coverage curve。",
            "B2 的描述性指標改善不等於安全、有效、可部署，也不取代護理師或醫師的專業判斷。",
        ],
    }
    public_hash = sha256_bytes(canonical_json_bytes(public))
    public_path = PROJECT_ROOT / contract["outputs"]["public_summary_path"]
    write_json(public_path, public)

    started_at = datetime.now(timezone.utc)
    run_id = f"{started_at.strftime('%Y%m%dT%H%M%S%fZ')}-{public_hash[:8]}"
    run_directory = PROJECT_ROOT / contract["outputs"]["run_output_root"] / run_id
    run_directory.mkdir(parents=True, exist_ok=False)
    private = {
        "schema_version": 1,
        "experiment_id": contract["experiment_id"],
        "source_oof_sha256": sha256_file(oof_path),
        "bootstrap": bootstrap_private,
        "paired_binary_checks": paired_binary,
    }
    private_path = run_directory / contract["outputs"]["full_result_filename"]
    write_json(private_path, private)

    tracked_inputs = [args.contract, *source_paths, str(oof_path.relative_to(PROJECT_ROOT))]
    manifest = {
        "manifest_schema_version": 1,
        "run_id": run_id,
        "experiment_id": contract["experiment_id"],
        "started_at_utc": started_at.isoformat().replace("+00:00", "Z"),
        "command": [sys.executable, *sys.argv],
        "documented_command": [
            "poetry",
            "run",
            "python",
            "scripts/run_day28_safety_statistics.py",
        ],
        "git": git_state(PROJECT_ROOT),
        "runtime": {
            "python": platform.python_version(),
            "packages": installed_versions(["numpy"]),
        },
        "inputs": [file_record(PROJECT_ROOT, path) for path in tracked_inputs],
        "parameters": {
            "comparison": "B2 versus B0",
            "resampling_unit": "record_index",
            "bootstrap_replications": bootstrap_public["replications"],
            "seed": bootstrap_public["seed"],
            "mcnemar_repeat_id": paired_binary["analysis_repeat_id"],
        },
        "outputs": [
            file_record(PROJECT_ROOT, str(public_path.relative_to(PROJECT_ROOT))),
            file_record(PROJECT_ROOT, str(private_path.relative_to(PROJECT_ROOT))),
        ],
        "privacy": "公開檔只有聚合指標、區間與合成教學曲線；bootstrap replicates 留在 gitignored results/runs/day-28/。",
        "scope": contract["scope"],
    }
    write_json(run_directory / "run-manifest.json", manifest)

    metrics = public["paired_record_bootstrap"]["metrics"]
    print("Day 28 安全指標與配對統計：完成")
    print(f"輸入：{validation['record_count']:,} 個 record_index × {validation['repeat_count']} 次重複")
    print(f"Bootstrap：{bootstrap_public['replications']:,} 次，單位 = record_index")
    print(
        "B2 − B0 balanced accuracy："
        f"{metrics['balanced_accuracy']['difference_B2_minus_B0']:.6f} "
        f"{metrics['balanced_accuracy']['difference_B2_minus_B0_percentile_95_ci']}"
    )
    print(
        "B2 − B0 severe undertriage rate："
        f"{metrics['severe_undertriage_rate']['difference_B2_minus_B0']:.6f} "
        f"{metrics['severe_undertriage_rate']['difference_B2_minus_B0_percentile_95_ci']}"
    )
    print(f"公開摘要：{public_path.relative_to(PROJECT_ROOT)}")
    print(f"本機完整重