1+ #!/usr/bin/env python3
2+ """
3+ Statistical validation of classification performance.
4+ Computes:
5+ 1. Fold-wise summary: mean ± SD, 95% CI (t-dist) for F1/precision/recall
6+ 2. Bootstrap 95% CI for test set macro F1
7+ 3. DeLong 95% CI for test set AUC (requires probability scores)
8+ 4. Wilcoxon signed-rank + Holm correction for model pairs (requires all models)
9+
10+ Usage:
11+ python extra_stat_validation.py \
12+ \
13+ \
14+ \
15+
16+ """
17+ import argparse
18+ import json
19+ import numpy as np
20+ import pandas as pd
21+ from pathlib import Path
22+ from scipy import stats
23+ from scipy .stats import wilcoxon
24+ from sklearn .metrics import f1_score , roc_auc_score
25+ from itertools import combinations
26+
27+
28+ # ─────────────────────────────────────────────
29+ # 1. Fold-wise summary
30+ # ─────────────────────────────────────────────
31+ def fold_summary (fold_results_path : str , label : str ) -> list :
32+ with open (fold_results_path ) as f :
33+ data = json .load (f )
34+
35+ # Handle both formats
36+ if isinstance (data , list ):
37+ # fold_results.json — plain list of fold dicts
38+ folds = data
39+ elif 'fold_results' in data :
40+ # cv_evaluation_results.json — nested under fold_results key
41+ folds = data ['fold_results' ]
42+ else :
43+ raise ValueError (f"Unrecognised format in { fold_results_path } )
44+
45+ # Extract per-fold metrics — handle nested 'deterministic' key if present
46+ def get_metric (fold , metric ):
47+ if 'deterministic' in fold :
48+ return fold ['deterministic' ][metric ]
49+ return fold [metric ]
50+
51+ metrics = ['f1' , 'precision' , 'recall' ]
52+ rows = []
53+ for metric in metrics :
54+ values = np .array ([get_metric (f , metric ) for f in folds ])
55+ n = len (values )
56+ mean , std = values .mean (), values .std (ddof = 1 )
57+ se = std / np .sqrt (n )
58+ ci_lo , ci_hi = stats .t .interval (0.95 , df = n - 1 , loc = mean , scale = se )
59+ rows .append ({
60+ 'model' : label ,
61+ 'metric' : metric ,
62+ 'mean' : round (mean , 4 ),
63+ 'std' : round (std , 4 ),
64+ 'ci_lo' : round (ci_lo , 4 ),
65+ 'ci_hi' : round (ci_hi , 4 ),
66+ 'values' : values .tolist (),
67+ })
68+ return rows
69+
70+
71+ # ─────────────────────────────────────────────
72+ # 2. Bootstrap CI for test macro F1
73+ # ─────────────────────────────────────────────
74+ def bootstrap_f1_ci (y_true , y_pred , n_bootstrap = 10000 , seed = 42 ):
75+ """Paired bootstrap 95% CI for macro F1."""
76+ rng = np .random .default_rng (seed )
77+ n = len (y_true )
78+ scores = []
79+ for _ in range (n_bootstrap ):
80+ idx = rng .integers (0 , n , size = n )
81+ scores .append (f1_score (y_true [idx ], y_pred [idx ],
82+ average = 'macro' , zero_division = 0 ))
83+ scores = np .array (scores )
84+ return {
85+ 'f1_observed' : round (f1_score (y_true , y_pred ,
86+ average = 'macro' , zero_division = 0 ), 4 ),
87+ 'ci_lo' : round (np .percentile (scores , 2.5 ), 4 ),
88+ 'ci_hi' : round (np .percentile (scores , 97.5 ), 4 ),
89+ 'n_bootstrap' : n_bootstrap ,
90+ }
91+
92+
93+ # ─────────────────────────────────────────────
94+ # 3. DeLong AUC CI
95+ # ─────────────────────────────────────────────
96+ def delong_auc_ci (y_true , y_score , alpha = 0.05 ):
97+ """
98+ DeLong et al. (1988) method for AUC confidence interval.
99+ y_score: predicted probability of positive class.
100+ """
101+ y_true = np .array (y_true )
102+ y_score = np .array (y_score )
103+
104+ pos = y_score [y_true == 1 ]
105+ neg = y_score [y_true == 0 ]
106+ n_pos , n_neg = len (pos ), len (neg )
107+
108+ # Placement values
109+ psi_pos = np .array ([np .mean (p > neg ) + 0.5 * np .mean (p == neg ) for p in pos ])
110+ psi_neg = np .array ([np .mean (n < pos ) + 0.5 * np .mean (n == pos ) for n in neg ])
111+
112+ auc = psi_pos .mean ()
113+
114+ # Variance via structural components
115+ var_pos = np .var (psi_pos , ddof = 1 ) / n_pos
116+ var_neg = np .var (psi_neg , ddof = 1 ) / n_neg
117+ se = np .sqrt (var_pos + var_neg )
118+
119+ z = stats .norm .ppf (1 - alpha / 2 )
120+ return {
121+ 'auc' : round (auc , 4 ),
122+ 'se' : round (se , 6 ),
123+ 'ci_lo' : round (max (0 , auc - z * se ), 4 ),
124+ 'ci_hi' : round (min (1 , auc + z * se ), 4 ),
125+ }
126+
127+
128+ # ─────────────────────────────────────────────
129+ # 4. Wilcoxon signed-rank + Holm correction
130+ # ─────────────────────────────────────────────
131+ def pairwise_wilcoxon (fold_data : dict , metric = 'f1' ):
132+ """
133+ fold_data: {label: [fold_0_score, fold_1_score, ...]}
134+ Returns DataFrame with pairwise Wilcoxon results + Holm correction.
135+ """
136+ labels = list (fold_data .keys ())
137+ rows = []
138+ for a , b in combinations (labels , 2 ):
139+ xa = np .array (fold_data [a ])
140+ xb = np .array (fold_data [b ])
141+ diff = xa - xb
142+ if np .all (diff == 0 ):
143+ stat , p = np .nan , 1.0
144+ else :
145+ stat , p = wilcoxon (xa , xb , alternative = 'two-sided' )
146+ rows .append ({
147+ 'model_a' : a ,
148+ 'model_b' : b ,
149+ 'mean_diff' : round ((xa - xb ).mean (), 4 ),
150+ 'statistic' : stat ,
151+ 'p_value' : p ,
152+ })
153+
154+ df = pd .DataFrame (rows ).sort_values ('p_value' )
155+
156+ # Holm correction
157+ n = len (df )
158+ df ['p_holm' ] = [min (1.0 , p * (n - i ))
159+ for i , p in enumerate (df ['p_value' ])]
160+ df ['significant' ] = df ['p_holm' ] < 0.05
161+ df ['p_holm' ] = df ['p_holm' ].round (4 )
162+ return df
163+
164+
165+ # ─────────────────────────────────────────────
166+ # Main
167+ # ─────────────────────────────────────────────
168+ def main ():
169+ parser = argparse .ArgumentParser ()
170+ parser .add_argument ('--fold_results' , nargs = '+' , required = True ,
171+ help = 'fold_results.json files, one per model' )
172+ parser .add_argument ('--test_preds' , nargs = '+' , required = True ,
173+ help = 'test_predictions.csv files, one per model' )
174+ parser .add_argument ('--labels' , nargs = '+' , required = True ,
175+ help = 'Model labels (same order as above)' )
176+ parser .add_argument ('--prob_col' , default = 'mean_confidence' ,
177+ help = 'Column with predicted probability for DeLong' )
178+ parser .add_argument ('--output_dir' , default = 'stat_results' )
179+ args = parser .parse_args ()
180+
181+ output_dir = Path (args .output_dir )
182+ output_dir .mkdir (parents = True , exist_ok = True )
183+
184+ assert len (args .fold_results ) == len (args .test_preds ) == len (args .labels ), \
185+ "Must provide equal numbers of fold_results, test_preds, and labels"
186+
187+ # ── 1 & 2 & 3: per-model analysis ──
188+ all_fold_summary = []
189+ all_bootstrap = []
190+ all_delong = []
191+ fold_f1_data = {} # for Wilcoxon
192+
193+ for fold_path , pred_path , label in zip (
194+ args .fold_results , args .test_preds , args .labels ):
195+
196+ # Fold summary
197+ with open (fold_path ) as f :
198+ fold_results = json .load (f )
199+ rows = fold_summary (fold_path , label )
200+ all_fold_summary .extend (rows )
201+ fold_f1_data [label ] = [r ['values' ] for r in rows
202+ if r ['metric' ] == 'f1' ][0 ]
203+
204+ # Test set
205+ df = pd .read_csv (pred_path )
206+ y_true = (df ['true_label' ] == 'pc' ).astype (int ).values
207+ y_pred = (df ['consensus_prediction' ] == 'pc' ).astype (int ).values
208+
209+ # Bootstrap F1
210+ boot = bootstrap_f1_ci (y_true , y_pred )
211+ boot ['model' ] = label
212+ all_bootstrap .append (boot )
213+
214+ # DeLong AUC — only if probability column is meaningful
215+ if args .prob_col in df .columns :
216+ # mean_confidence is max(P(lnc), P(pc)) — need to check direction
217+ # If consensus_prediction == 'pc', confidence is P(pc); else P(lnc)
218+ y_score = np .where (
219+ df ['consensus_prediction' ] == 'pc' ,
220+ df [args .prob_col ],
221+ 1 - df [args .prob_col ]
222+ )
223+ dl = delong_auc_ci (y_true , y_score )
224+ dl ['model' ] = label
225+ all_delong .append (dl )
226+ else :
227+ print (f" WARNING: { args .prob_col } { pred_path }
228+ f"skipping DeLong for { label } )
229+
230+ # ── 4: Wilcoxon pairwise ──
231+ wilcoxon_df = pairwise_wilcoxon (fold_f1_data , metric = 'f1' )
232+
233+ # ── Save outputs ──
234+ fold_df = pd .DataFrame (all_fold_summary ).drop (columns = 'values' )
235+ fold_df .to_csv (output_dir / 'fold_summary.csv' , index = False )
236+
237+ boot_df = pd .DataFrame (all_bootstrap )
238+ boot_df .to_csv (output_dir / 'bootstrap_f1_ci.csv' , index = False )
239+
240+ if all_delong :
241+ delong_df = pd .DataFrame (all_delong )
242+ delong_df .to_csv (output_dir / 'delong_auc_ci.csv' , index = False )
243+
244+ wilcoxon_df .to_csv (output_dir / 'wilcoxon_pairwise.csv' , index = False )
245+
246+ # ── Print summary ──
247+ print ("\n ── Fold-wise summary ──" )
248+ print (fold_df .to_string (index = False ))
249+
250+ print ("\n ── Bootstrap F1 CI (test set) ──" )
251+ print (boot_df .to_string (index = False ))
252+
253+ if all_delong :
254+ print ("\n ── DeLong AUC CI (test set) ──" )
255+ print (delong_df .to_string (index = False ))
256+
257+ print ("\n ── Wilcoxon pairwise (fold F1, Holm-corrected) ──" )
258+ print (wilcoxon_df .to_string (index = False ))
259+
260+
261+ if __name__ == '__main__' :
262+ main ()
0 commit comments