plot_pvalues.py

"""Phase 1 p-values against the empirical null — the picture of the project's thesis.

The null here is NOT the textbook uniform. It is 2000 label permutations pushed through
the identical per-signal test battery, so it carries SECOM's real correlation structure
(112 of 590 columns are near-duplicates at |r| > 0.99). That structure does not shift the
null's centre but roughly doubles its spread, which is exactly what makes a naive
"how many hits did I get" reading unsafe.

  (a) where the p-values sit, against the null's own 95% envelope
  (b) how many hits at p<0.05 a TRUE null produces — the number the naive screen omits

Needs phase1_permnull_p.npy, written by phase1.py (gitignored: 7.6MB, regenerable).
"""
import numpy as np, pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from secom_common import load

OBS, NUL, INK, MUTED = "#4a3aa7", "#8a8a85", "#0b0b0b", "#52514e"
SURF = "#fcfcfb"
NB = 20

D = pd.read_csv("phase1_univariate.csv", index_col=0)
P = np.load("phase1_permnull_p.npy")            # B x m, null p-values
p = D.p.values
m, B = len(p), P.shape[0]
pi0 = min(1., 2*np.mean(p > .5))
edges = np.linspace(0, 1, NB+1)
exp_per_bin = m/NB

obs_h, _ = np.histogram(p, bins=edges)
null_h = np.array([np.histogram(P[b], bins=edges)[0] for b in range(B)])
lo, hi = np.percentile(null_h, [2.5, 97.5], axis=0)
mid = null_h.mean(0)

hits = (P < .05).sum(1)
obs_hits = int((p < .05).sum())

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11.5, 4.6))
fig.patch.set_facecolor(SURF)

# --- (a) p-value histogram vs null envelope ------------------------------
ctr = (edges[:-1]+edges[1:])/2
# append the closing edge, else step="post" leaves the last bin unfilled
ax1.fill_between(np.r_[edges[:-1], 1.], np.r_[lo, lo[-1]], np.r_[hi, hi[-1]],
                 step="post", color=NUL, alpha=.35, lw=0)
ax1.step(np.r_[edges[:-1], 1], np.r_[mid, mid[-1]], where="post", color=NUL, lw=1.6)
ax1.bar(ctr, obs_h, width=.045, color=OBS, alpha=.85, zorder=3)
ax1.axhline(exp_per_bin, color=INK, ls="--", lw=1.1, zorder=4)
ax1.axhline(pi0*exp_per_bin, color=INK, ls=":", lw=1.1, zorder=4)
# reference lines described in the empty upper-right, not on top of the bars
ax1.annotate(f"- - -  uniform expectation, {exp_per_bin:.1f} per bin\n"
             f"·····  Storey null component, pi0 = {pi0:.3f}\n"
             f"         -> at most ~{m*(1-pi0):.0f} of {m} signals are non-null",
             xy=(.97, .93), xycoords="axes fraction", ha="right", va="top",
             fontsize=8, color=INK, linespacing=1.6)
ax1.set_xlabel("univariate p-value", color=MUTED)
ax1.set_ylabel(f"signals per bin (of {m})", color=MUTED)
ax1.set_xlim(0, 1)
ax1.set_title("The signal is one spike at the left; the rest is flat", fontsize=11,
              color=INK, loc="left", pad=26)
ax1.annotate(f"bars = observed   ·   grey band = 95% of {B} permuted nulls   ·   "
             f"first bin {obs_h[0]} vs null {mid[0]:.0f}",
             xy=(0, 1.012), xycoords="axes fraction", ha="left", va="bottom",
             fontsize=8.5, color=MUTED)

# --- (b) hits at p<0.05: observed vs null --------------------------------
bins2 = np.arange(hits.min()-1, max(hits.max(), obs_hits)+4, 2)
ax2.hist(hits, bins=bins2, color=NUL, alpha=.55, lw=0)
ax2.hist(hits, bins=bins2, histtype="step", color=NUL, lw=1.6)
ax2.axvline(obs_hits, color=OBS, lw=2.4)
ax2.set_xlim(bins2[0], obs_hits+8)
ax2.annotate(f"observed\n{obs_hits} hits", xy=(obs_hits, ax2.get_ylim()[1]*.80),
             xytext=(-7, 0), textcoords="offset points", ha="right", va="top",
             fontsize=9, color=OBS, linespacing=1.4)
ax2.axvline(hits.max(), color=INK, lw=1.1, ls=":")
ax2.annotate(f"a true null reached\n{hits.max()} hits once in {B}", xy=(hits.max(), ax2.get_ylim()[1]*.45),
             xytext=(-6, 0), textcoords="offset points", ha="right", va="top",
             fontsize=8, color=INK, linespacing=1.4)
ax2.set_xlabel("signals with p < 0.05", color=MUTED)
ax2.set_ylabel(f"permutations (of {B})", color=MUTED)
ax2.set_title("What noise alone produces at p < 0.05", fontsize=11,
              color=INK, loc="left", pad=26)
ax2.annotate(f"null mean {hits.mean():.1f}, sd {hits.std():.1f}   ·   "
             f"permutation p < {1/B:.4f}",
             xy=(0, 1.012), xycoords="axes fraction", ha="left", va="bottom",
             fontsize=8.5, color=MUTED)

for ax in (ax1, ax2):
    ax.set_facecolor(SURF)
    ax.grid(axis="y", 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)

fig.suptitle("Phase 1 p-values against a permuted null — why 86 'significant' signals is not 86 findings",
             fontsize=13, color=INK, x=.008, ha="left", y=.978)
fig.text(.008, .012,
         f"{m} testable signals, {B} label permutations through the identical per-signal test battery "
         f"(Welch t or Mann-Whitney, chosen by shape).\n"
         f"The null keeps the real correlation structure, so its sd is {hits.std():.1f} where "
         f"Binomial({m}, 0.05) would give {np.sqrt(m*.05*.95):.1f} — the near-duplicate columns vote in blocks.\n"
         f"BH q=0.05 cuts the {obs_hits} to 28. See REPORT.md.",
         fontsize=7.5, color=MUTED, ha="left", linespacing=1.5)
fig.tight_layout(rect=[0, .09, 1, .935])
fig.savefig("figures/pvalue_null.png", dpi=200, facecolor=fig.get_facecolor())
print(f"observed hits {obs_hits}   null mean {hits.mean():.2f} sd {hits.std():.2f} "
      f"max {hits.max()}   pi0 {pi0:.4f}")
print(f"first bin: observed {obs_h[0]}, null mean {mid[0]:.1f}, null 95% [{lo[0]:.0f}, {hi[0]:.0f}]")
print(f"draws matching or beating observed: {(hits >= obs_hits).sum()}/{B}")
print("  -> figures/pvalue_null.png")