plot_missingness.py

"""Missingness over time — why Phase 0 concluded MAR-on-time rather than MNAR.

The two populations are drawn separately on purpose. In one raster of all 590 signals
the 28 heavily-missing columns are 5% of the height and the block structure — the whole
point — is invisible. Split, each population is legible at its own scale:

  (a) the 28 signals above 50% missing: long contiguous outages, excluded from modeling
  (b) the other 562: scattered single-lot dropouts
  (c) missing values per lot over time — the -0.417 Spearman drift

The rasters are evidence for "instrument off for a stretch" and against "random
dropout". What they CANNOT show is any relation to the label — that needed the
time-stratified permutation test in phase0c.py, which found none.
"""
import numpy as np, pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from matplotlib.colors import ListedColormap
from scipy import stats
from secom_common import load

INK, MUTED, MARK, SURF = "#0b0b0b", "#52514e", "#2f2e2b", "#fcfcfb"
CMAP = ListedColormap([SURF, MARK])

X, y, ts = load()
o = np.argsort(ts.values)
X, ts = X.iloc[o].reset_index(drop=True), ts.iloc[o].reset_index(drop=True)
n, p = X.shape
M = X.isna().values
mrate = M.mean(0)
xnum = mdates.date2num(ts)
cut = int(.7*n)


def longest_run(col):
    best = cur = 0
    for v in col:
        cur = cur + 1 if v else 0
        best = max(best, cur)
    return best


def blockiness(cols):
    """Median share of a column's gaps that sit in its single longest run."""
    v = [longest_run(M[:, j])/M[:, j].sum() for j in cols if M[:, j].any()]
    return float(np.median(v)) if v else np.nan


hi = np.nonzero(mrate > .5)[0]
lo = np.nonzero(mrate <= .5)[0]
hi = hi[np.lexsort((np.array([M[:, j].argmax() for j in hi]), -mrate[hi]))]
lo = lo[np.argsort(-mrate[lo])]
b_hi, b_lo = blockiness(hi), blockiness(lo)
nz = int((mrate == 0).sum())
rho, rp = stats.spearmanr(M.sum(1), np.arange(n))

fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(11, 8.4), sharex=True,
                                    gridspec_kw={"height_ratios": [1.0, 1.5, .95]})
fig.patch.set_facecolor(SURF)

for ax, idx, ttl, sub in [
    (ax1, hi, f"The {len(hi)} signals above 50% missing: long contiguous outages",
     f"excluded from all modeling, never imputed   ·   median {b_hi:.0%} of each "
     f"signal's gaps fall in its single longest run"),
    (ax2, lo, f"The other {len(lo)} signals: scattered single-lot dropouts",
     f"median {b_lo:.0%} of gaps in the longest run   ·   {nz} of these are never "
     f"missing at all (the clean band at the foot)")]:
    ax.imshow(M[:, idx].T, aspect="auto", interpolation="nearest", cmap=CMAP,
              vmin=0, vmax=1, extent=[xnum[0], xnum[-1], len(idx), 0])
    ax.set_ylim(len(idx), 0)
    ax.set_ylabel(f"{len(idx)} signals\n(most missing at top)", color=MUTED, fontsize=9)
    ax.set_title(ttl, fontsize=11, color=INK, loc="left", pad=26)
    ax.annotate(sub, xy=(0, 1.012), xycoords="axes fraction", ha="left", va="bottom",
                fontsize=8.5, color=MUTED)

cnt = M.sum(1)
ax3.plot(ts, cnt, color=MARK, lw=.6, alpha=.30)
ax3.plot(ts, pd.Series(cnt).rolling(50, center=True).median(), color=MARK, lw=2)
ax3.set_ylabel("missing values\nper lot", color=MUTED, fontsize=9)
ax3.set_xlabel("lot timestamp", color=MUTED)
ax3.set_title("Fewer values go missing as the record goes on", fontsize=11,
              color=INK, loc="left", pad=26)
ax3.annotate(f"Spearman(missing count, time) rho = {rho:+.3f}, p = {rp:.1e}   ·   "
             f"thin = per lot, thick = 50-lot rolling median",
             xy=(0, 1.012), xycoords="axes fraction", ha="left", va="bottom",
             fontsize=8.5, color=MUTED)
ax3.grid(axis="y", color="#e6e5e1", lw=.8)
ax3.set_axisbelow(True)
ax3.annotate("Phase 4 temporal split", xy=(xnum[cut], ax3.get_ylim()[1]*.95),
             xytext=(-8, 0), textcoords="offset points", ha="right", va="top",
             fontsize=8.5, color=INK)

for ax in (ax1, ax2, ax3):
    ax.axvline(xnum[cut], color=INK, lw=1.3, ls=":")
    ax.set_facecolor(SURF)
    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)
ax3.xaxis.set_major_locator(mdates.WeekdayLocator(byweekday=mdates.MO, interval=2))
ax3.xaxis.set_major_formatter(mdates.DateFormatter("%b %d"))

fig.suptitle("SECOM missingness over time — instrument availability, not lost data",
             fontsize=13, color=INK, x=.008, ha="left", y=.987)
fig.text(.008, .012,
         f"{n} lots x {p} signals, chronological; each dark mark is one absent value "
         f"({int(M.sum()):,} of {n*p:,} cells, {M.mean():.2%}).\n"
         f"Conditioned on time, 0 of 538 partially-missing columns show any missing-vs-label "
         f"association (stratified permutation test, B=20,000).\n"
         f"That is what licenses imputation — not the block structure alone. See REPORT.md.",
         fontsize=7.5, color=MUTED, ha="left", linespacing=1.5)
fig.tight_layout(rect=[0, .072, 1, .95])
fig.savefig("figures/missingness_over_time.png", dpi=200, facecolor=fig.get_facecolor())
print(f"cells {M.sum():,}/{n*p:,} ({M.mean():.4%})  >50%-missing {len(hi)} (blockiness "
      f"{b_hi:.3f})  rest {len(lo)} (blockiness {b_lo:.3f})  never-missing {nz}")
print(f"spearman rho={rho:+.4f} p={rp:.3e}")
print("  -> figures/missingness_over_time.png")