R7/analysis: RQ1/RQ2 confirmatory decision layer (frozen §6 gates in code)

Encodes the four lead-paper confirmatory tests exactly per the frozen prereg §6,
calibrated on synthetic ground truth only (never re-fit to cell data):

- RQ1-P1 leak: correlation AUC + BCa CI; gate = CI excludes 0.5 (D2); 0.60 CI
  lower bound = separate materiality label (material / weak-but-real), not the gate.
- RQ1-P2 padding: paired ΔAUC = AUC(no-pad) − AUC(pad) over circuits; effective
  iff CI > 0.
- RQ2-P1 anonymity set: ΔH = H(federated) − H(single, matched N) with Miller-Madow
  per-circuit entropy; two-sided grow / honest-shrink / inconclusive by CI sign.
- RQ2-P3 mechanism: Spearman ρ(top-3 bridge concentration, per-circuit H) + CI.

apply_holm corrects the reported RQ1/RQ2 subset against the full frozen family of
7 (family_size default 7) — never re-optimised to the 4 reported. Bootstrap
p-values order the Holm step-down only; every decision is a CI gate, never a bare p.

stats.py: bootstrap CIs can now return their resample distribution so the Holm
ordering p-value comes from the same resamples as the CI.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

cmd_chat/sor/analysis/stats.py

@@ -128,6 +128,19 @@ def _bca_endpoints(
     return q_lo, q_hi, "bca"
 
 
+def two_sided_bootstrap_p(thetas: Sequence[float], null: float) -> float:
+    """A bootstrap two-sided p-value for H0: θ = ``null`` from a bootstrap
+    distribution ``thetas`` — 2·min(P(θ* ≤ null), P(θ* ≥ null)), capped at 1.0.
+    Used only to *order* the Holm family; the pre-registered decision gates are
+    the CIs, never a p-value alone (§6)."""
+    b = len(thetas)
+    if b == 0:
+        return 1.0
+    le = sum(1 for t in thetas if t <= null) / b
+    ge = sum(1 for t in thetas if t >= null) / b
+    return min(1.0, 2.0 * min(le, ge))
+
+
 def bootstrap_ci(
     units: Sequence[T],
     statistic: Callable[[Sequence[T]], float],
@@ -136,16 +149,19 @@ def bootstrap_ci(
     alpha: float = DEFAULT_ALPHA,
     seed: int = 0,
     method: str = "bca",
-) -> CIResult:
+    return_dist: bool = False,
+):
     """One-sample bootstrap CI of ``statistic`` over ``units`` (the unit of
     analysis — a circuit-pair for RQ1, a circuit for RQ2). Resamples ``units`` with
     replacement ``n_resamples`` times. ``method="bca"`` applies bias-correction +
     acceleration (falling back to percentile if degenerate); ``"percentile"``
-    forces the plain interval."""
+    forces the plain interval. With ``return_dist=True`` also returns the sorted
+    bootstrap distribution (so a p-value can be derived from the same resamples)."""
     n = len(units)
     if n == 0:
-        return CIResult(float("nan"), float("nan"), float("nan"),
-                        alpha, n_resamples, "empty", seed)
+        res = CIResult(float("nan"), float("nan"), float("nan"),
+                       alpha, n_resamples, "empty", seed)
+        return (res, []) if return_dist else res
     theta_hat = float(statistic(units))
     rng = random.Random(seed)
     thetas: List[float] = []
@@ -160,8 +176,9 @@ def bootstrap_ci(
     else:
         q_lo, q_hi, used = alpha / 2.0, 1.0 - alpha / 2.0, "percentile"
 
-    return CIResult(theta_hat, _percentile(thetas, q_lo), _percentile(thetas, q_hi),
-                    alpha, n_resamples, used, seed)
+    res = CIResult(theta_hat, _percentile(thetas, q_lo), _percentile(thetas, q_hi),
+                   alpha, n_resamples, used, seed)
+    return (res, thetas) if return_dist else res
 
 
 def two_sample_diff_ci(
@@ -173,15 +190,18 @@ def two_sample_diff_ci(
     alpha: float = DEFAULT_ALPHA,
     seed: int = 0,
     method: str = "bca",
-) -> CIResult:
+    return_dist: bool = False,
+):
     """Bootstrap CI for the difference ``statistic(A) - statistic(B)`` of two
     independent arms (RQ2-P1: ΔH = H(federated) − H(single-house, matched N)).
     Each arm is resampled independently. BCa uses a combined leave-one-out
-    jackknife across both arms (each point dropped from its own arm)."""
+    jackknife across both arms (each point dropped from its own arm). With
+    ``return_dist=True`` also returns the sorted bootstrap distribution."""
     na, nb = len(units_a), len(units_b)
     if na == 0 or nb == 0:
-        return CIResult(float("nan"), float("nan"), float("nan"),
-                        alpha, n_resamples, "empty", seed)
+        res = CIResult(float("nan"), float("nan"), float("nan"),
+                       alpha, n_resamples, "empty", seed)
+        return (res, []) if return_dist else res
     theta_hat = float(statistic(units_a)) - float(statistic(units_b))
     rng = random.Random(seed)
     thetas: List[float] = []
@@ -203,8 +223,9 @@ def two_sample_diff_ci(
     else:
         q_lo, q_hi, used = alpha / 2.0, 1.0 - alpha / 2.0, "percentile"
 
-    return CIResult(theta_hat, _percentile(thetas, q_lo), _percentile(thetas, q_hi),
-                    alpha, n_resamples, used, seed)
+    res = CIResult(theta_hat, _percentile(thetas, q_lo), _percentile(thetas, q_hi),
+                   alpha, n_resamples, used, seed)
+    return (res, thetas) if return_dist else res
 
 
 def mean(xs: Sequence[float]) -> float: