phase0c.py

"""PHASE 0c — exact replacement for the CMH test, which is invalid at ~3 missing rows/col.
Permute the label WITHIN time octiles: kills the time confound, keeps every cell count exact.
"""
import numpy as np, pandas as pd
from secom_common import load

rng = np.random.default_rng(0)
B = 20000
X, y, ts = load()
o = np.argsort(ts.values)
X, y = X.iloc[o].reset_index(drop=True), y.iloc[o].values.astype(float)
n = len(X)
strat = np.repeat(np.arange(8), n // 8 + 1)[:n]          # 8 chronological octiles

miss = X.isna()
cols = X.columns[(miss.mean() > 0) & (miss.mean() < 1)]
M = miss[cols].values.astype(float)                       # n x k

# B permutations of y, shuffled within each stratum
Y = np.tile(y, (B, 1))
for s in range(8):
    idx = np.nonzero(strat == s)[0]
    for r in range(B):
        Y[r, idx] = rng.permutation(Y[r, idx])

obs = y @ M                                               # fails among missing rows, per col
null = Y @ M                                              # B x k
dev_obs = np.abs(obs - null.mean(0))
p = ((np.abs(null - null.mean(0)) >= dev_obs - 1e-9).sum(0) + 1) / (B + 1)

R = pd.DataFrame({"miss_rate": miss[cols].mean().values, "p_perm": p}, index=cols)


def bh(pv, q):
    pv = np.asarray(pv, float); o = np.argsort(pv)
    ok = pv[o] <= q * np.arange(1, len(pv) + 1) / len(pv)
    k = np.nonzero(ok)[0].max() + 1 if ok.any() else 0
    r = np.zeros(len(pv), bool); r[o[:k]] = True; return r


R["bh05"] = bh(R.p_perm, .05)
print("[7b] MISSING vs LABEL, time-stratified PERMUTATION test (B=20000, labels shuffled within time octile)")
print(f"     tested columns            : {len(R)}")
print(f"     raw p<.05                 : {int((R.p_perm<.05).sum())}   [null expectation {.05*len(R):.0f}]")
print(f"     BH q=0.05 survivors       : {int(R.bh05.sum())}")
print(f"     KS test of p-values vs U[0,1]: p={__import__('scipy.stats', fromlist=['x']).kstest(R.p_perm,'uniform').pvalue:.3f}"
      "   (large p = p-values uniform = no label information in the NA pattern)")
print("\n     smallest p, time-adjusted:")
print(R.sort_values("p_perm").head(6).to_string(float_format=lambda v: f"{v:.4g}"))
R.to_csv("phase0c_perm.csv")