plot_failrate.py

"""Fail rate over calendar time — the confounder that shapes every other phase.

Phase 0 [6] found the fail rate varies 7x across the record (chi2 p = 5.7e-10). That is
why the temporal split is mandatory, why missingness had to be tested conditional on
time, and why S129's association exists only in the later half.

Two panels:
  (a) equal-count bins (~196 lots each) drawn at their true calendar extent, with Wilson
      95% intervals — bin WIDTH shows how unevenly lots arrive, bin HEIGHT the rate.
      Wilson rather than normal: at ~13 failures per bin the normal interval is wrong
      and can run below zero.
  (b) a 150-lot rolling rate, so the structure is visible without any bin choice at all.
"""
import numpy as np, pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from scipy import stats
from secom_common import load

FAIL_C, INK, MUTED = "#eb6834", "#0b0b0b", "#52514e"
NBIN, ROLL = 8, 150

X, y, ts = load()
o = np.argsort(ts.values)
y, ts = y.iloc[o].reset_index(drop=True), ts.iloc[o].reset_index(drop=True)
yv = y.values
n = len(y)
base = yv.mean()


def wilson(k, m, z=1.96):
    """Wilson score interval. Correct at small counts; normal approx is not."""
    if m == 0:
        return 0., 0.
    p = k/m
    d = 1 + z**2/m
    c = (p + z**2/(2*m))/d
    h = z*np.sqrt(p*(1-p)/m + z**2/(4*m**2))/d
    return c-h, c+h


# --- equal-count bins ----------------------------------------------------
# qcut, matching phase0b.py exactly: a different edge rule shifts chi2 in the 3rd
# digit and the figure would then disagree with the report over the same statistic.
oct_ = pd.qcut(np.arange(n), NBIN, labels=False)
rows = []
for i in range(NBIN):
    idx = np.nonzero(oct_ == i)[0]
    s, e = idx[0], idx[-1]+1
    k, m = int(yv[s:e].sum()), e-s
    lo, hi = wilson(k, m)
    rows.append((ts[s], ts[e-1], k, m, k/m, lo, hi))
B = pd.DataFrame(rows, columns=["t0", "t1", "fails", "lots", "rate", "lo", "hi"])
ct = np.array([[r.fails, r.lots-r.fails] for r in B.itertuples()])
chi2_p = stats.chi2_contingency(ct)[1]
rho, rho_p = stats.spearmanr(np.arange(n), yv)
cut = int(.7*n)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(11, 6.6), sharex=True,
                               gridspec_kw={"height_ratios": [1.25, 1]})
fig.patch.set_facecolor("#fcfcfb")

# (a) binned rate with Wilson CI
for r in B.itertuples():
    w = mdates.date2num(r.t1) - mdates.date2num(r.t0)
    ax1.bar(mdates.date2num(r.t0), r.rate, width=w, align="edge",
            color=FAIL_C, alpha=.45, edgecolor=FAIL_C, lw=1.4)
    mid = mdates.date2num(r.t0) + w/2
    ax1.plot([mid, mid], [r.lo, r.hi], color=FAIL_C, lw=1.6, solid_capstyle="butt")
    ax1.plot([mid], [r.rate], "o", ms=4, color=FAIL_C,
             markeredgecolor="#fcfcfb", markeredgewidth=1.2)
ax1.axhline(base, color=MUTED, ls="--", lw=1.2)
ax1.annotate(f"overall {base:.2%}", xy=(ts.iloc[-1], base), xytext=(-4, 5),
             textcoords="offset points", ha="right", fontsize=8.5, color=MUTED)
ax1.set_ylabel("fail rate (Wilson 95% CI)", color=MUTED)
ax1.set_ylim(0, max(B.hi)*1.18)
ax1.yaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax1.set_title(f"Fail rate is not stationary: {B.rate.min():.1%} to {B.rate.max():.1%} "
              f"across {NBIN} equal-count bins", fontsize=11, color=INK, loc="left", pad=30)
ax1.annotate(f"chi-square homogeneity p = {chi2_p:.1e}   ·   "
             f"Spearman(time, fail) rho = {rho:+.3f}, p = {rho_p:.1e}   ·   "
             f"bin width = calendar span of ~{n//NBIN} lots",
             xy=(0, 1.012), xycoords="axes fraction", ha="left", va="bottom",
             fontsize=8.5, color=MUTED)

# (b) rolling rate — no bin choice
roll = pd.Series(yv).rolling(ROLL, center=True).mean()
ax2.plot(ts, roll, color=FAIL_C, lw=2)
ax2.axhline(base, color=MUTED, ls="--", lw=1.2)
ax2.set_ylabel(f"{ROLL}-lot rolling fail rate", color=MUTED)
ax2.yaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax2.set_ylim(0, np.nanmax(roll)*1.15)
ax2.set_xlabel("lot timestamp", color=MUTED)
ax2.set_title("Same signal without binning", fontsize=11, color=INK, loc="left", pad=10)

# the Phase 4 temporal split, on both panels
for ax in (ax1, ax2):
    ax.axvline(ts[cut], color=INK, lw=1.3, ls=":")
ax1.annotate(f"Phase 4 temporal split (70%)\ntrain {yv[:cut].mean():.2%}  →  "
             f"test {yv[cut:].mean():.2%}",
             xy=(ts[cut], ax1.get_ylim()[1]*.97), xytext=(-8, 0), textcoords="offset points",
             ha="right", va="top", fontsize=8.5, color=INK, linespacing=1.4)

for ax in (ax1, ax2):
    ax.set_facecolor("#fcfcfb")
    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)
ax2.xaxis.set_major_locator(mdates.WeekdayLocator(byweekday=mdates.MO, interval=2))
ax2.xaxis.set_major_formatter(mdates.DateFormatter("%b %d"))

fig.suptitle("SECOM fail rate over time — the confounder behind every other result",
             fontsize=13, color=INK, x=.008, ha="left", y=.985)
fig.text(.008, .012,
         f"{n} lots, {int(yv.sum())} failures, {ts.min():%Y-%m-%d} to {ts.max():%Y-%m-%d}. "
         f"Bins hold ~{n//NBIN} lots each, so unequal widths show uneven lot arrival, not missing data.\n"
         f"A model trained before the dotted line is tested on a period with a "
         f"{yv[:cut].mean()/yv[cut:].mean():.1f}x lower fail rate — some of the temporal AUC drop is this "
         f"shift, not leakage.\n"
         f"The rolling line is centred, so it stops {ROLL//2} lots short of each end. See REPORT.md.",
         fontsize=7.5, color=MUTED, ha="left", linespacing=1.5)
fig.tight_layout(rect=[0, .075, 1, .945])
fig.savefig("figures/failrate_over_time.png", dpi=200, facecolor=fig.get_facecolor())
print(B.assign(t0=B.t0.dt.strftime("%m-%d"), t1=B.t1.dt.strftime("%m-%d")).to_string(
    float_format=lambda v: f"{v:.4f}"))
print(f"chi2 p={chi2_p:.3e}  spearman rho={rho:+.4f} p={rho_p:.3e}  "
      f"train={yv[:cut].mean():.4f} test={yv[cut:].mean():.4f}")
print("  -> figures/failrate_over_time.png")