plot_stability.py

"""Stability selection frequencies against the permuted null, and per training fold.

  (a) every signal's selection frequency, ranked, against the range 40 permuted nulls
      reach at the same rank. The 60% threshold the brief specifies sits INSIDE that
      range — noise cleared it 5 times in one permutation against 6 observed — so the
      panel draws the noise ceiling rather than just the threshold.
  (b) the same signals re-run inside each chronological training fold, penalty
      calibrated to equal sparsity so folds are comparable.

Reads phase3_stability.csv, phase3_nullfreq.csv, phase3_folds_calibrated.csv (phase3.py,
phase3_folds_calib.py).
"""
import numpy as np, pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

OBS, NUL, INK, MUTED, SURF = "#4a3aa7", "#8a8a85", "#0b0b0b", "#52514e", "#fcfcfb"
RAMP = ["#c9c3ea", "#9c92d4", "#6f61bd", "#4a3aa7"]     # sequential: fold is ordered
PI, NSHOW = .60, 60

F = pd.read_csv("phase3_stability.csv", index_col=0)
NF = pd.read_csv("phase3_nullfreq.csv").values
FD = pd.read_csv("phase3_folds_calibrated.csv", index_col=0)
folds = [c for c in FD.columns if c.endswith("%")]

obs = F.freq.sort_values(ascending=False)
null_sorted = -np.sort(-NF, axis=1)                       # each permutation, ranked
n_all, B = len(obs), NF.shape[0]
rank = np.arange(1, n_all+1)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11.5, 4.9),
                               gridspec_kw={"width_ratios": [1.35, 1]})
fig.patch.set_facecolor(SURF)

# --- (a) ranked frequencies vs the null's own ranked frequencies ---------
ax1.fill_between(rank, null_sorted.min(0), null_sorted.max(0), color=NUL, alpha=.30, lw=0)
ax1.plot(rank, null_sorted.mean(0), color=NUL, lw=1.6)
ax1.plot(rank, obs.values, color=OBS, lw=2)
ax1.plot(rank[:9], obs.values[:9], "o", ms=5, color=OBS,
         markeredgecolor=SURF, markeredgewidth=1)
ax1.axhline(PI, color=INK, ls="--", lw=1.1)
ax1.axhline(NF.max(), color=INK, ls=":", lw=1.1)
ax1.annotate(f"60% threshold", xy=(NSHOW, PI), xytext=(-4, 4), textcoords="offset points",
             ha="right", fontsize=8, color=INK)
ax1.annotate(f"highest frequency noise ever reached ({NF.max():.2f})", xy=(NSHOW, NF.max()),
             xytext=(-4, 4), textcoords="offset points", ha="right", fontsize=8, color=INK)
# ranks 4 and 5 differ by 0.01, so a fixed offset stacks their labels: alternate
for i, (s, v) in enumerate(obs.head(6).items()):
    ax1.annotate(s, xy=(i+1, v), xytext=(7, 7 if i % 2 == 0 else -13),
                 textcoords="offset points", fontsize=8.5, color=OBS)
ax1.set_xlim(0, NSHOW)
ax1.set_ylim(0, 1.04)
ax1.set_xlabel(f"signal rank (of {n_all} in the pool)", color=MUTED)
ax1.set_ylabel("selection frequency over 100 subsamples", color=MUTED)
ax1.yaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax1.set_title("Only S59 clears what noise never reached", fontsize=11,
              color=INK, loc="left", pad=26)
ax1.annotate(f"line = observed   ·   grey = range over {B} permuted nulls at the same rank",
             xy=(0, 1.012), xycoords="axes fraction", ha="left", va="bottom",
             fontsize=8.5, color=MUTED)

# --- (b) per training fold ----------------------------------------------
order = FD.sort_values(["folds>60%"] + folds, ascending=False).index.tolist()[::-1]
for i, sig in enumerate(order):
    v = FD.loc[sig, folds].astype(float).values
    ax2.plot([v.min(), v.max()], [i, i], color=MUTED, lw=1.2, zorder=1)
    for k, (f, c) in enumerate(zip(folds, RAMP)):
        ax2.plot(v[k], i, "o", ms=7, color=c, markeredgecolor=SURF,
                 markeredgewidth=1, zorder=3, label=f"train {f}" if i == 0 else None)
    ax2.annotate(f"{int(FD.loc[sig,'folds>60%'])}/4", xy=(1.03, i),
                 xycoords=("axes fraction", "data"), ha="left", va="center",
                 fontsize=8.5, color=INK if FD.loc[sig, "folds>60%"] >= 3 else MUTED)
ax2.axvline(PI, color=INK, ls="--", lw=1.1)
ax2.set_yticks(range(len(order)))
ax2.set_yticklabels(order, fontsize=9.5)
ax2.set_xlim(0, 1.02)
ax2.set_ylim(-.6, len(order)-.4)
ax2.xaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax2.set_xlabel("selection frequency within the fold", color=MUTED)
ax2.set_title("Reproduced inside each training fold", fontsize=11,
              color=INK, loc="left", pad=26)
ax2.annotate("penalty calibrated per fold to equal sparsity (q≈20)",
             xy=(0, 1.012), xycoords="axes fraction", ha="left", va="bottom",
             fontsize=8.5, color=MUTED)
ax2.annotate("folds >60%", xy=(1.03, 1.012), xycoords="axes fraction",
             ha="left", va="bottom", fontsize=8, color=MUTED)
ax2.legend(frameon=False, fontsize=8, loc="lower right", handletextpad=.2,
           borderaxespad=.2, labelspacing=.25)

for ax in (ax1, ax2):
    ax.set_facecolor(SURF)
    ax.grid(axis="y" if ax is ax1 else "x", color="#e6e5e1", lw=.8)
    ax.set_axisbelow(True)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    for s in ("left", "bottom"):
        ax.spines[s].set_color("#d5d4cf")
    ax.tick_params(colors=MUTED, labelsize=9)
ax2.spines["left"].set_visible(False)

fig.suptitle("Stability selection — frequencies against noise, and whether they reproduce",
             fontsize=13, color=INK, x=.008, ha="left", y=.975)
fig.text(.008, .012,
         f"L1 logistic on 100 subsamples at 50%, stratified; imputation and scaling refit inside every "
         f"subsample. {int((F.freq==0).sum())} of {n_all} signals were never selected once.\n"
         f"The empirical FDR at the 60% line is 0.29 — it is a soft threshold, not a test. S64 reproduces "
         f"in 4/4 folds but is univariately null (q_BH = 0.50): a suppressor, not a finding. See REPORT.md.",
         fontsize=7.5, color=MUTED, ha="left", linespacing=1.5)
fig.tight_layout(rect=[0, .075, .965, .935])
fig.savefig("figures/stability_frequencies.png", dpi=200, facecolor=fig.get_facecolor())
print(obs.head(8).round(3).to_string())
print(f"null: max {NF.max():.2f}, mean top-rank {null_sorted[:,0].mean():.2f}, "
      f"signals >60% under null max {int((NF>PI).sum(1).max())}")
print(f"never selected: {int((F.freq==0).sum())}/{n_all}")
print("  -> figures/stability_frequencies.png")