Add reviewer variability statistics reporting

Author: LMBooth

analysis_pipeline/stage6_build_publication_report.py

@@ -3,12 +3,18 @@
 import argparse
 import glob
 import json
+import math
 from datetime import datetime, timezone
 from pathlib import Path
 from typing import Any
 
 import pandas as pd
 
+try:
+    from scipy.stats import t as student_t
+except ImportError:  # pragma: no cover - report still works with normal approximation
+    student_t = None
+
 
 def _analysis_root() -> Path:
     return Path(__file__).resolve().parent
@@ -179,6 +185,57 @@ def _as_float(value: Any) -> float | None:
         return None
 
 
+def _finite_series(values: pd.Series) -> pd.Series:
+    numeric = pd.to_numeric(values, errors="coerce")
+    return numeric.dropna()
+
+
+def _ci_multiplier(n: int) -> float:
+    if n <= 1:
+        return 0.0
+    if student_t is not None:
+        return float(student_t.ppf(0.975, df=n - 1))
+    return 1.96
+
+
+def _describe_metric(values: pd.Series) -> dict[str, float | int | None]:
+    clean = _finite_series(values)
+    n = int(clean.shape[0])
+    if n == 0:
+        return {
+            "n": 0,
+            "mean": None,
+            "sd": None,
+            "se": None,
+            "ci95_low": None,
+            "ci95_high": None,
+            "min": None,
+            "max": None,
+        }
+    mean_value = float(clean.mean())
+    if n == 1:
+        sd_value = 0.0
+        se_value = 0.0
+        ci_low = mean_value
+        ci_high = mean_value
+    else:
+        sd_value = float(clean.std(ddof=1))
+        se_value = float(sd_value / math.sqrt(n))
+        margin = _ci_multiplier(n) * se_value
+        ci_low = mean_value - margin
+        ci_high = mean_value + margin
+    return {
+        "n": n,
+        "mean": mean_value,
+        "sd": sd_value,
+        "se": se_value,
+        "ci95_low": float(ci_low),
+        "ci95_high": float(ci_high),
+        "min": float(clean.min()),
+        "max": float(clean.max()),
+    }
+
+
 def _collect_stage6_result_paths(explicit: list[str] | None, glob_pattern: str | None) -> list[Path]:
     paths: list[Path] = []
     for item in explicit or []:
@@ -303,6 +360,42 @@ def _build_global_best_table(best_df: pd.DataFrame) -> pd.DataFrame:
     return pd.DataFrame(rows).sort_values(["dataset", "protocol"]).reset_index(drop=True)
 
 
+def _build_modality_protocol_variability_table(best_df: pd.DataFrame) -> pd.DataFrame:
+    if best_df.empty:
+        return pd.DataFrame()
+    required = {"dataset", "protocol", "scenario", "balanced_accuracy_mean", "macro_f1_mean"}
+    if not required.issubset(set(best_df.columns)):
+        return pd.DataFrame()
+
+    rows: list[dict[str, Any]] = []
+    grouped = best_df.groupby(["dataset", "protocol"], sort=True)
+    for (dataset, protocol), group in grouped:
+        ba = _describe_metric(group["balanced_accuracy_mean"])
+        f1 = _describe_metric(group["macro_f1_mean"])
+        scenario_names = sorted(str(item) for item in group["scenario"].dropna().unique())
+        rows.append(
+            {
+                "dataset": str(dataset),
+                "protocol": str(protocol),
+                "n_scenarios": int(ba["n"]),
+                "mean_balanced_accuracy": ba["mean"],
+                "sd_balanced_accuracy": ba["sd"],
+                "se_balanced_accuracy": ba["se"],
+                "ci95_low_balanced_accuracy": ba["ci95_low"],
+                "ci95_high_balanced_accuracy": ba["ci95_high"],
+                "min_balanced_accuracy": ba["min"],
+                "max_balanced_accuracy": ba["max"],
+                "mean_macro_f1": f1["mean"],
+                "sd_macro_f1": f1["sd"],
+                "se_macro_f1": f1["se"],
+                "ci95_low_macro_f1": f1["ci95_low"],
+                "ci95_high_macro_f1": f1["ci95_high"],
+                "scenarios": "; ".join(scenario_names),
+            }
+        )
+    return pd.DataFrame(rows).sort_values(["dataset", "protocol"]).reset_index(drop=True)
+
+
 def _participants_with_drops(
     epoch_payload: dict[str, Any] | None,
 ) -> tuple[list[str], list[str], list[str], pd.DataFrame]:
@@ -426,6 +519,7 @@ def _markdown_from_components(
     stage6_scenarios: pd.DataFrame,
     stage6_best: pd.DataFrame,
     stage6_global_best: pd.DataFrame,
+    stage6_variability: pd.DataFrame,
     dropped_windows_table: pd.DataFrame,
 ) -> str:
     lines: list[str] = []
@@ -640,6 +734,12 @@ def _markdown_from_components(
     if not stage6_global_best.empty:
         lines.append("### Global Best Pipeline per Dataset and Protocol")
         lines.append(_to_markdown_table(stage6_global_best))
+    if not stage6_variability.empty:
+        lines.append("### Modality/Protocol Variability Across Class Scenarios")
+        lines.append(
+            "_Values summarise the best-pipeline balanced-accuracy and macro-F1 means across loaded class-scenario result files._"
+        )
+        lines.append(_to_markdown_table(stage6_variability, decimals=4))
     lines.append("")
 
     stage1_strict_state = _stage1_strict_mode(run_manifest_payload)
@@ -743,6 +843,7 @@ def main() -> None:
 
     stage6_scenarios, stage6_best, _stage6_agg = _build_stage6_frames(result_paths)
     stage6_global_best = _build_global_best_table(stage6_best)
+    stage6_variability = _build_modality_protocol_variability_table(stage6_best)
 
     eeg_drop_subjects, ecg_drop_subjects, pupil_drop_subjects, dropped_windows_table = _participants_with_drops(
         epoch_payload
@@ -765,6 +866,7 @@ def main() -> None:
         stage6_scenarios=stage6_scenarios,
         stage6_best=stage6_best,
         stage6_global_best=stage6_global_best,
+        stage6_variability=stage6_variability,
         dropped_windows_table=dropped_windows_table,
     )
 
@@ -786,6 +888,7 @@ def main() -> None:
             "scenario_rows": stage6_scenarios.to_dict(orient="records"),
             "best_rows": stage6_best.to_dict(orient="records"),
             "global_best_by_dataset_protocol": stage6_global_best.to_dict(orient="records"),
+            "modality_protocol_variability": stage6_variability.to_dict(orient="records"),
         },
         "segmentation_dropped_windows": {
             "subjects_with_eeg_drops": eeg_drop_subjects,