# -*- coding: utf-8 -*- """ figures.py — every number the published pages print, computed here, once. Written 2026-08-23 after the round-two reviews found the THIRD recurrence of the same failure: make_pages.py held every published figure as a hardcoded string literal, so the pages were still prose with HTML around them. A5's own basis number was typed into that file on the day the rule forbidding prose numbers was written. Nothing downstream may type a number. make_pages.py imports FIGURES from here and interpolates. If a figure is not in this dict it does not go on a page. Run: python3 scripts/figures.py # prints every value python3 scripts/figures.py --json # writes web/figures.json """ import re import functools import json, math, os, sys import numpy as np from scipy import stats as sps HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, HERE) import money_analysis as M # noqa: E402 import phase2 # noqa: E402 import nominees # noqa: E402 SEED = 20260823 ROOT_DIR = os.path.normpath(os.path.join(HERE, "..")) def data_path(name): """Beside this file first, then the repo root. The published bundle is flat: scripts and data sit in one directory. Resolving only against the repo root means the bundle runs in place (falling through to the parent checkout) and dies the moment a reader copies it anywhere else — which is what a round-four reviewer found, on a bundle whose README says that exact gap is closed. build_all.py now runs the bundle from a temporary copy to keep it honest.""" here = os.path.join(HERE, name) return here if os.path.exists(here) else os.path.join(ROOT_DIR, name) # The register. Confidences live here and the page reads them from here — the round-two # review confirmed all 25 matched v3, but matching by hand is not a mechanism. CONF = dict(A1a=92, A1b=88, A2=65, A3=65, A4=70, A1ap=85, A1bp=88, A5=88, A3p=78, A4p=65, B1=80, B2=70, B3=65, B4=55, C1=78, C2=60, D1=88, E1=55, E2=57, E3=70, F1=88, F2p=60, F2pp=93, G1=62, G2=80) KNOWN_MISS = ("A1a", "A1b", "A3", "A4") # The binding spec's enumeration of what counts as a recorded legal matter for G2. G2_ENUM = {"Conviction", "Charge filed", "Settlement", "Civil judgment"} # Hits the register itself declares worthless — a free hit is as much a distortion # of a calibration score as a forced miss, and the page counted only the misses. DECLARED_FREE = ("G1", "A5", "F2p", "F2pp", "E3") CALIBRATION_SELF = ("F2p", "F2pp") # excluded from the a-hat fitting set — closes the paradox V2_CONF = [92, 88, 65, 65, 70, 80, 70, 65, 55, 78, 60, 88, 55, 57, 70, 88, 60, 62, 80] def _z(r): return math.atanh(max(min(r, 0.999999), -0.999999)) def _poisson_binomial(p): d = np.zeros(len(p) + 1) d[0] = 1.0 for x in p: d[1:] = d[1:] * (1 - x) + d[:-1] * x d[0] *= (1 - x) return d def _ci_hits(p): c = np.cumsum(_poisson_binomial(np.asarray(p, float))) return int(np.searchsorted(c, 0.025)), int(np.searchsorted(c, 0.975)) def _ahat(y, p): """MLE of the log-odds calibration shift: y ~ logistic(logit(p) + a).""" lo = np.log(p / (1 - p)) a = 0.0 for _ in range(300): q = 1 / (1 + np.exp(-(lo + a))) g = float((y - q).sum()) h = float(-(q * (1 - q)).sum()) if abs(h) < 1e-13: break s = g / h a -= s if abs(s) < 1e-13: break return a def vintage(): """The finance file's actual coverage dates, computed — the pages used to say "receipts through 30 June 2026" while 35 committees' coverage ran later (the fact-check gate caught the four sentences).""" import csv as _csv from datetime import datetime as _dt raw = list(_csv.DictReader(open(M.CSV))) dates = sorted(_dt.strptime(r["FEC_Coverage_Through"], "%m/%d/%Y") for r in raw) modal = max(set(dates), key=dates.count) return dict(n_files=len(dates), n_at_modal=dates.count(modal), modal=modal.strftime("%d %B %Y").lstrip("0"), earliest=dates[0].strftime("%d %B %Y").lstrip("0"), latest=dates[-1].strftime("%d %B %Y").lstrip("0")) def missing_finance(): """Who has no row in the frozen finance file, and why — the page used to say all three "have no federal campaign committee"; the master coding says two do.""" have = {d["name"] for d in M.load()} names = [n for n in phase2.AGE_TENURE if n not in have] return dict(n=len(names), names=", ".join(sorted(names))) def _band_ranking(ratios): """Least- and most-lopsided band by plain distance |ratio - 1|, with the distances.""" label = {"A1ap": "A1a′", "A1bp": "A1b′", "A3p": "A3′", "A4p": "A4′"} dist = {k: abs(v - 1.0) for k, v in ratios.items()} least = min(dist, key=dist.get) most = max(dist, key=dist.get) return dict(band_least_name=label[least], band_least_dist=dist[least], band_most_name=label[most], band_most_dist=dist[most]) def money(): rows = M.load() inc = [d for d in rows if d["tenure"] > 0] non = [d for d in rows if d["tenure"] == 0] x = np.array([d["tenure"] for d in inc], float) y = np.array([d["corp_pct"] for d in inc], float) xa = np.array([d["tenure"] for d in rows], float) ya = np.array([d["corp_pct"] for d in rows], float) r_inc, p_inc = sps.pearsonr(x, y) r_pool, p_pool = sps.pearsonr(xa, ya) rho_inc = sps.spearmanr(x, y).statistic rho_pool = sps.spearmanr(xa, ya).statistic slope, icept, rr, pp, se = sps.linregress(x, y) rng = np.random.default_rng(SEED) bs = [] for _ in range(10000): i = rng.integers(0, len(x), len(x)) if np.std(x[i]) and np.std(y[i]): bs.append(np.corrcoef(x[i], y[i])[0, 1]) lo, hi = np.percentile(bs, [2.5, 97.5]) RD = {"CA", "OH", "TX", "UT", "NC", "FL", "LA", "MO", "TN"} sub = [d for d in rows if d["district"][:2] not in RD] r_nord = sps.pearsonr(np.array([d["tenure"] for d in sub], float), np.array([d["corp_pct"] for d in sub], float))[0] # "Three longest careers removed." Round three found this had NO defined value: # Costa and Cuellar tie at 21.8 years, so the result depended on CSV row order # (0.355 vs 0.502 — a 0.15 swing in a statistic used in the correction story). # Rule, fixed 2026-08-23: ties break by (tenure desc, name asc), and the value # under the opposite tie-break is COMPUTED AND PUBLISHED beside it. order = sorted(inc, key=lambda d: (-d["tenure"], d["name"])) trim = order[3:] r_trim = sps.pearsonr(np.array([d["tenure"] for d in trim], float), np.array([d["corp_pct"] for d in trim], float))[0] # opposite tie-break: name DESC within equal tenure order2 = sorted(inc, key=lambda d: d["name"], reverse=True) order2 = sorted(order2, key=lambda d: -d["tenure"]) # stable: keeps name-desc within ties trim2 = order2[3:] r_trim_alt = sps.pearsonr(np.array([d["tenure"] for d in trim2], float), np.array([d["corp_pct"] for d in trim2], float))[0] trim_dropped = ", ".join(d["name"] for d in order[:3]) trim_dropped_alt = ", ".join(d["name"] for d in order2[:3]) out = dict( n_rows=len(rows), n_inc=len(inc), n_non=len(non), r_inc=r_inc, p_inc=p_inc, rho_inc=rho_inc, z_inc=_z(r_inc), r_pool=r_pool, p_pool=p_pool, rho_pool=rho_pool, z_pool=_z(r_pool), dz=_z(r_pool) - _z(r_inc), # A1a' band edges. The z form is what binds; these are its r image, and the # card first printed them rounded to 2dp, where BOTH edges fell outside the rule. r_band_lo=math.tanh(_z(r_inc) - 0.20), r_band_hi=math.tanh(_z(r_inc) + 0.20), # A band symmetric in z is NOT symmetric once converted back to r. Published # because the asymmetry runs in my favour and an undisclosed favour is a defect. r_band_down=r_inc - math.tanh(_z(r_inc) - 0.20), r_band_up=math.tanh(_z(r_inc) + 0.20) - r_inc, boot_lo=lo, boot_hi=hi, mean_inc=float(np.mean(y)), mean_non=float(np.mean([d["corp_pct"] for d in non])), step=float(np.mean(y) - np.mean([d["corp_pct"] for d in non])), slope=slope, slope_se=se, slope_lo=slope - 1.96 * se, slope_hi=slope + 1.96 * se, r2=rr ** 2, tenure_min=float(x.min()), tenure_max=float(x.max()), # Mean tenure of the served. The RAW group gap (`step` above, mean_inc - mean_non) # and the FITTED step (fit.b1_step, the jump at zero years) are different numbers # that both got called "the step"; they differ by slope x this mean tenure, and the # page now prints that reconciliation rather than leaving two steps unreconciled. tenure_mean=float(x.mean()), climb=slope * (x.max() - x.min()), scatter=float(np.std(y - (icept + slope * x), ddof=2)), r_no_redistricted=r_nord, n_no_redistricted=len(sub), r_trim3=r_trim, r_trim3_alt=r_trim_alt, trim3_dropped=trim_dropped, trim3_dropped_alt=trim_dropped_alt, trim3_tied=int(r_trim != r_trim_alt), top_name=max(inc, key=lambda d: d["tenure"])["name"], top_tenure=max(inc, key=lambda d: d["tenure"])["tenure"], top_share=max(inc, key=lambda d: d["tenure"])["corp_pct"], ) # ── Leverage and influence, published rather than asserted ────────────── # Round-four review: the page named its highest-leverage point, kept it, and # said keeping it was the discipline — without ever printing what removing it # would do. An unquantified outlier is an unanswered objection. These are # DIAGNOSTICS, published beside the estimator and never substituted for it; # the frozen specification's "no outlier removal ever" governs what is scored. # Leverage and influence are different quantities and both are reported, # because a point can be extreme in x and still move nothing. n_i = len(x) sxx = float(np.sum((x - x.mean()) ** 2)) hat = 1.0 / n_i + (x - x.mean()) ** 2 / sxx # hat-matrix diagonal drop_r = np.array([sps.pearsonr(np.delete(x, i), np.delete(y, i))[0] for i in range(n_i)]) # leave-one-out Pearson infl = np.abs(drop_r - r_inc) i_lev, i_inf = int(np.argmax(hat)), int(np.argmax(infl)) ord_lev = np.argsort(-hat) ord_inf = np.argsort(-infl) ten2 = float(np.sort(x)[-2]) # A zero on the x-axis is not a zero on the y-axis. Published because the glossary # printed n_non (never served) where it needed this, and nearly doubled the stated # point mass at zero — the same defect as every other number that lived in prose. n_zero = sum(1 for d in rows if d["corp_pct"] == 0) # The leverage-appropriate check the page owed A5: delete the point and REFIT. # Holding the full-sample slope and truncating the range is not a robustness test. ix2 = [i for i in range(n_i) if i != i_lev] xr, yr = x[ix2], y[ix2] sl_r, ic_r, r_r, _, _ = sps.linregress(xr, yr) out.update( n_zero=n_zero, pct_zero=100.0 * n_zero / len(rows), # Pearson-vs-Spearman gaps, with and without the one far-out career. The # glossary asserted the pooled gap was all structural; half of it is one man. sp_gap_inc=rho_inc - r_inc, sp_gap_pool=rho_pool - r_pool, # A4' completes the band table. Round four: the page said the two hand-set # bands were the lopsided ones; A4' is the LEAST lopsided and the only one # whose extra room runs upward. Computing three of four is how that happened. a4p_ratio=(r_inc - 0.10) / (0.55 - r_inc), # What deleting the leverage point does to the two LIVE registered bands. # Both end up outside their own top edge. Nothing here is scored that way, # and a reader is owed the number rather than the reassurance. r_pool_drop_lev=float(sps.pearsonr(np.delete(xa, int(np.argmax(xa))), np.delete(ya, int(np.argmax(xa))))[0]), sp_gap_pool_drop=float(sps.spearmanr(np.delete(xa, int(np.argmax(xa))), np.delete(ya, int(np.argmax(xa)))).statistic - sps.pearsonr(np.delete(xa, int(np.argmax(xa))), np.delete(ya, int(np.argmax(xa))))[0]), zero_overstate_pct=100.0 * (len(non) - n_zero) / n_zero, pct_rise=100.0 * (float(np.mean(y)) / float(np.mean([d['corp_pct'] for d in non])) - 1), n_inc_minus_lev=n_i - 1, # Every band on the page, and how lopsided each one is around its basis. # Round four: the page confessed the SMALLEST asymmetry (A1a', forced by the # z transform) and was silent on the two larger hand-chosen ones. Selective # disclosure reads as fuller disclosure than it is, so all of them are computed. a3p_down=r_inc - 0.10, a3p_up=0.45 - r_inc, a3p_ratio=(r_inc - 0.10) / (0.45 - r_inc), a4p_down=r_inc - 0.10, a4p_up=0.55 - r_inc, a1ap_ratio=(r_inc - math.tanh(_z(r_inc) - 0.20)) / (math.tanh(_z(r_inc) + 0.20) - r_inc), a1bp_lo=math.tanh(_z(r_pool) - 0.20), a1bp_hi=math.tanh(_z(r_pool) + 0.20), a1bp_down=r_pool - math.tanh(_z(r_pool) - 0.20), a1bp_up=math.tanh(_z(r_pool) + 0.20) - r_pool, a1bp_ratio=((r_pool - math.tanh(_z(r_pool) - 0.20)) / (math.tanh(_z(r_pool) + 0.20) - r_pool)), # THE RANKING OF LOPSIDEDNESS, COMPUTED. The glossary said "ranked by distance from # 1.00, A1a' is the least lopsided" — which is true on a log scale and false on the # plain distance the sentence names (A4' at 0.88 is nearer to 1.00 than A1a' at # 1.13). The final pre-publication reviewer caught it; this entry had already been # mis-worded twice. So the ranking is computed here on the plain |ratio - 1| the page # says it uses, and the page prints the names and the distances, not a memory. **_band_ranking(dict( A1ap=((r_inc - math.tanh(_z(r_inc) - 0.20)) / (math.tanh(_z(r_inc) + 0.20) - r_inc)), A1bp=((r_pool - math.tanh(_z(r_pool) - 0.20)) / (math.tanh(_z(r_pool) + 0.20) - r_pool)), A3p=(r_inc - 0.10) / (0.45 - r_inc), A4p=(r_inc - 0.10) / (0.55 - r_inc))), n_zero_inc=sum(1 for d in inc if d["corp_pct"] == 0), slope_refit=float(sl_r), climb_refit=float(sl_r * (xr.max() - xr.min())), lev_name=inc[i_lev]["name"], lev_h=float(hat[i_lev]), lev_h_next=float(hat[ord_lev[1]]), lev_next_name=inc[int(ord_lev[1])]["name"], infl_name=inc[i_inf]["name"], infl_dr=float(infl[i_inf]), infl_next_name=inc[int(ord_inf[1])]["name"], infl_next_dr=float(infl[ord_inf[1]]), sp_gap_inc_drop=float(sps.spearmanr(np.delete(x, i_lev), np.delete(y, i_lev)).statistic - drop_r[i_lev]), r_inc_drop_lev=float(drop_r[i_lev]), # r with the leverage point gone r_inc_drop_infl=float(drop_r[i_inf]), lev_is_infl=int(i_lev == i_inf), # How much of the "the slope moves further than the step" claim rests on # ground only one person stands on. tenure_2nd=ten2, n_inc_above_2nd=int(np.sum(x > ten2)), n_inc_over_20=int(np.sum(x > 20.0)), climb_to_2nd=float(slope * (ten2 - x.min())), fitted_at_top=float(icept + slope * x.max()), resid_top=float(y[int(np.argmax(x))] - (icept + slope * x.max())), ) for pty in ("DEM", "REP"): s = [d for d in rows if d["party"] == pty] si = [d for d in inc if d["party"] == pty] out[f"r_pool_{pty}"] = sps.pearsonr(np.array([d["tenure"] for d in s], float), np.array([d["corp_pct"] for d in s], float))[0] out[f"n_pool_{pty}"] = len(s) out[f"r_inc_{pty}"] = sps.pearsonr(np.array([d["tenure"] for d in si], float), np.array([d["corp_pct"] for d in si], float))[0] out[f"mean_inc_{pty}"] = float(np.mean([d["corp_pct"] for d in si])) # Challenger/open-seat, counted WITHIN the tenure==0 population — the definition the # spec binds. Round three: counting FEC codes over all 112 rows gave "64 ... of whom # 50 and 16" (=66) under a heading reading "both are correct", because Gillen and # Luria are FEC-coded C while carrying prior House service (tenure > 0). for code, key in (("C", "chal"), ("O", "open")): g = [d for d in non if d["inc_code"] == code] out[f"n_{key}"] = len(g) out[f"mean_{key}"] = float(np.mean([d["corp_pct"] for d in g])) # the FEC-code groupings over all rows, kept under their own honest names — these are # what the August record's 0.3% / 1.3% means actually were computed over for code, key in (("C", "chal_feccode"), ("O", "open_feccode")): g = [d for d in rows if d["inc_code"] == code] out[f"n_{key}"] = len(g) out[f"mean_{key}"] = float(np.mean([d["corp_pct"] for d in g])) assert out["n_chal"] + out["n_open"] == len(non), "challenger+open must equal never-served" return out def roster(): from tier_gate import TIER_UNVERIFIED P = phase2.PHASE2 party = {n: p for lst in nominees.NOMINEES.values() for _, p, n in lst} at = phase2.AGE_TENURE mil = {k: sum(1 for v in P.values() if v["mil"] == k) for k in ("Yes", "No", "Unclear")} def fisher(pred, field, drop=(), drop_names=()): t = {} for n, v in P.items(): if party[n] not in ("DEM", "REP") or v[field] in drop: continue # A row whose sourcing the index has declared NOT ESTABLISHED is not evidence. # It stays in the published record, with its full detail; it does not stay in # a statistic. Found by an independent media lawyer, who pointed out that a # claim published beside its own notice of insufficient sourcing was also an # INPUT to a party-asymmetry figure this project bets on. if n in drop_names and pred(v[field]): continue t[(party[n], pred(v[field]))] = t.get((party[n], pred(v[field])), 0) + 1 R = [t.get(("REP", True), 0), t.get(("REP", False), 0)] D = [t.get(("DEM", True), 0), t.get(("DEM", False), 0)] orr, p = sps.fisher_exact([R, D]) # The odds ratio is what the register bets on; it is NOT what a reader hears. # Round-four review: a reader told "odds ratio 2.11" concludes "twice as likely", # and the plain-English figure is the risk ratio. Both are published so the # threshold and the sentence a reader can repeat are never the same number by # accident. Shares are printed too, because a ratio of two small shares is the # easiest number on this page to mislead with. pR = R[0] / sum(R) if sum(R) else float("nan") pD = D[0] / sum(D) if sum(D) else float("nan") return dict(R_yes=R[0], R_n=sum(R), D_yes=D[0], D_n=sum(D), OR=orr, p=p, R_pct=100 * pR, D_pct=100 * pD, RR=(pR / pD) if pD else float("nan"), or_over_rr_pct=100 * (orr / (pR / pD) - 1) if pD and pR else float("nan"), pp_diff=100 * (pR - pD)) return dict( n_roster=len(P), n_districts=len(nominees.NOMINEES), # Tier claims whose cited source is tertiary. Counted so the page can print # the number rather than leaving the ladder to be spot-checked by a reviewer. n_tier_unverified=len(__import__('tier_gate').TIER_UNVERIFIED), # How many of them claim TIER 1 -- a court or government record -- on nothing but # a wiki page. That is the strongest claim in the ladder resting on the weakest # citation, and it was not separated out until a red-teamer asked for it. n_tier1_unverified=sum( 1 for n in __import__('tier_gate').TIER_UNVERIFIED if (P.get(n, {}).get("tier") or "").startswith("1 - ")), n_tier_unverified_was=30, # the count the page carried before Phase A closed it # PHASE A, 2026-08-30 — how the 30 wiki-only rows were resolved, computed from the # tier strings so the account of the fix cannot itself go stale. n_tier_resourced=sum(1 for v in P.values() if "re-sourced 2026-08-30" in (v.get("status") or "")), n_tier_gov=sum(1 for v in P.values() if (v.get("tier") or "").startswith("1 - Government record")), # Of those government-record rows, the ones resting on the member's OWN official # House biography (5) versus a court record newly surfaced by the stricter gate # (Gillen, 1). Split so the page arithmetic sums to the 30 originally flagged: # 2 verified + 5 House-bio + 23 downgraded = 30, and Gillen is the +1 the tighter # rule caught, not part of the original count. n_tier_govbio=sum(1 for v in P.values() if "official House biography" in (v.get("tier") or "")), n_tier_gov_new=sum(1 for v in P.values() if (v.get("tier") or "").startswith("1 - Government record") and "court" in (v.get("tier") or "").lower()), # 2026-09-01: rows relabelled by the final pre-publication review from "court/govt # record" to the member's own House.gov page, because their citations were the # member's own page plus outlets and no court record. Counted separately so the # government-record total on the page still sums. n_tier1_all=sum(1 for v in P.values() if (v.get("tier") or "").startswith("1 - ")), n_tier_court=sum(1 for v in P.values() if (v.get("tier") or "").startswith("1 - Court")), n_tier_gov_ownpage=sum(1 for v in P.values() if "House.gov page" in (v.get("tier") or "")), # Rows whose free-text notes had wiki-only reported matters struck under the # source-it-or-drop-it rule on 2026-09-02, counted from the rows so the page's # statement of how many cannot go stale (it did once, at "two"). # The legal column's floor. A row's "None found" means nothing surfaced in a # general-press search; most rows say in their own detail that no court-docket search # was performed. Counted from the rows so the page can state the floor as a number. n_no_docket_search=sum(1 for v in P.values() if __import__("re").search( r"no (PACER|docket|state docket)[^.]*search|general-press search only", v.get("legal_d") or "", __import__("re").I)), n_notes_struck=sum(1 for v in P.values() if "struck 2026-09-02" in (v.get("notes") or "").lower()), n_tier3_wiki=sum(1 for v in P.values() if (v.get("tier") or "").startswith("3 - Tertiary (Wikipedia only")), # Two districts hold no candidate: their primaries postdate the roster freeze. # The pages said 'all 115 candidates in the 59 districts', implying coverage # the index does not have. Computed so the disclosure cannot go stale. n_districts_empty=sum(1 for v in nominees.NOMINEES.values() if not v), n_districts_filled=sum(1 for v in nominees.NOMINEES.values() if v), n_dem=sum(1 for v in party.values() if v == "DEM"), n_rep=sum(1 for v in party.values() if v == "REP"), n_ind=sum(1 for v in party.values() if v == "IND"), n_age_blank=sum(1 for v in at.values() if v[0] == ""), n_age_marker=sum(1 for v in at.values() if isinstance(v[0], str) and v[0].startswith("~")), n_age_star=sum(1 for v in at.values() if isinstance(v[0], str) and v[0].endswith("*")), mil_yes=mil["Yes"], mil_no=mil["No"], mil_unclear=mil["Unclear"], # G1 under the spec's own rule: "Unclear" EXCLUDED. The published basis counted them in. G1=fisher(lambda v: v == "Yes", "mil", drop=("Unclear",)), G1_as_published=fisher(lambda v: v == "Yes", "mil"), # G2 counts EXACTLY the binding spec's enumeration (ANALYSIS_SPEC_v3 §4 G2: "the # legal column carries one of the enumerated adverse values — Conviction, Charge # filed, Settlement, or Civil judgment. The positive enumeration governs"). This used # to count anything that was not "None found", which swept in one row coded "Civil # suit (pending)" — a pending suit is not in the enumeration — and so printed a basis # of 7 Democrats where the spec's own rule yields 6. Found by the final pre-publication # reviewer, 2026-09-02. The verdict is invariant (p = 1.0 either way); the printed basis # was not the binding rule's, and the spec's own printed basis had the same defect. G2=fisher(lambda v: v in G2_ENUM, "legal", drop_names=TIER_UNVERIFIED), # As it stood before the not-established rows were excluded. Published beside the # scored basis so nobody can say a universe was picked after seeing which flattered. G2_with_unverified=fisher(lambda v: v in G2_ENUM, "legal"), # And as it stood under the broad "anything but None found" count, so the change # above is itself checkable: the basis before the enumeration was applied. G2_broad=fisher(lambda v: v != "None found", "legal", drop_names=TIER_UNVERIFIED), ) # Dependence clusters. The published interval convolves 25 INDEPENDENT Bernoullis, # which is the one assumption this register spends five cards refuting: A5 is implied # by A1a' in 20,000 of 20,000 resamples, A2 is nested, the four dead bets all trace to # one roster-completion event, B1-B4 share a single scoreability floor, C1 and C2 come # from one uncollected sweep, and F2'/F2" are deterministic functions of the rest. # Round-six review found the caveat the page volunteers (sample size) is the weaker of # the two problems. Predictions sharing a cluster move together; the clustered interval # is what a reader should actually read. CLUSTERS = { "money": ("A1a", "A1b", "A2", "A3", "A4", "A1ap", "A1bp", "A3p", "A4p", "A5"), "writing": ("B1", "B2", "B3", "B4"), "collection": ("C1", "C2"), "roster": ("D1",), "ages": ("E1", "E2", "E3"), "self": ("F1",), "calibration":("F2p", "F2pp"), "coding": ("G1", "G2"), } def scoring(): keys = list(CONF) p = np.array([CONF[k] for k in keys], float) / 100 free = [i for i, k in enumerate(keys) if k not in KNOWN_MISS] fit = [i for i, k in enumerate(keys) if k not in CALIBRATION_SELF] n = len(keys) lo, hi = _ci_hits(p) lo_f, hi_f = _ci_hits(p[free]) n_free = len(free) # Clustered interval: each cluster resolves as one draw at its mean confidence, # carrying all of its members with it. This is the pessimistic reading and the # honest one; the truth sits between it and the independent version above. # Conditioned exactly like the 44-76% interval above it: the four known misses are # zeroed, the rest of each cluster resolves together. Round-two review found the # first version left them live, so its top read 100% against a stated ceiling of # 84% two paragraphs earlier — a contradiction, not a wider interval. csizes, cprobs, cdead = [], [], [] for members in CLUSTERS.values(): live = [k for k in members if k not in KNOWN_MISS] idx = [keys.index(k) for k in live] csizes.append(len(idx)) cdead.append(len(members) - len(live)) cprobs.append(float(np.mean(p[idx])) if idx else 0.0) rngc = np.random.default_rng(SEED) draws = (rngc.random((40000, len(csizes))) < np.array(cprobs)) @ np.array(csizes, float) c_lo, c_hi = np.percentile(draws, [2.5, 97.5]) n_eff = len(CLUSTERS) rng = np.random.default_rng(SEED) km = [keys.index(k) for k in KNOWN_MISS] ia, ib = keys.index("F2p"), keys.index("F2pp") bad_v3 = bad_v4 = neg = 0 N = 20000 for _ in range(N): y = (rng.random(n) < p).astype(float) y[km] = 0.0 # v4: a-hat fitted WITHOUT the two self-referential predictions. Test it the same # way as v3 — vary their outcomes and check the estimator moves. It cannot, by # construction, but "by construction" is an argument and this is a measurement. a4 = _ahat(y[fit], p[fit]) neg += a4 < 0 for guess in (1.0, 0.0): yy = y.copy(); yy[ia] = guess; yy[ib] = guess if abs(_ahat(yy[fit], p[fit]) - a4) > 1e-12: bad_v4 += 1 break # v3: a-hat fitted over the whole set, including the two predictions it scores ok = False for guess in (1.0, 0.0): yy = y.copy(); yy[ia] = guess; yy[ib] = guess if (1.0 if _ahat(yy, p) < 0 else 0.0) == guess: ok = True break if not ok: bad_v3 += 1 # what freezing the four dead bets costs, versus deleting them brier_all = float(np.mean((p - np.where(np.isin(np.arange(n), km), 0.0, 1.0)) ** 2)) brier_del = float(np.mean((p[free] - 1.0) ** 2)) return dict( n_scored=n, mean_conf=float(p.mean() * 100), v2_n=len(V2_CONF), v2_mean=float(np.mean(V2_CONF)), n_free=n_free, n_free_hits=len(DECLARED_FREE), free_hit_conf=float(np.mean([CONF[k] for k in DECLARED_FREE])), ci_clustered_lo_pct=100.0*c_lo/n, ci_clustered_hi_pct=100.0*c_hi/n, n_clusters=n_eff, ci_lo_hits=lo, ci_hi_hits=hi, ci_lo_pct=100 * lo / n, ci_hi_pct=100 * hi / n, ci_free_lo_pct=100 * lo_f / n, ci_free_hi_pct=100 * hi_f / n, max_attainable_pct=100 * len(free) / n, p_hitrate_le_47=float(np.cumsum(_poisson_binomial(p))[int(0.47 * n)]), stale_lo=float(np.mean(V2_CONF)) - 24, stale_hi=float(np.mean(V2_CONF)) + 16, paradox_v3_pct=100 * bad_v3 / N, paradox_v4_pct=100 * bad_v4 / N, p_ahat_negative=neg / N, n_fit=len(fit), n_known_miss=len(KNOWN_MISS), brier_frozen=brier_all, brier_if_deleted=brier_del, brier_cost_of_freezing=brier_all - brier_del, ) def _flat(prefix, d, out): """Flatten one level of nesting so page placeholders can address roster.G1.p directly.""" for k, v in d.items(): if isinstance(v, dict): _flat(f"{prefix}{k}.", v, out) else: out[f"{prefix}{k}"] = v return out # ── Phase 2 additions, 2026-08-24 ──────────────────────────────────────────── # The ten districts coded in batches 9 and 10, named in 00_HANDOFF_BRIEF.md. # Restricting to everything else reconstructs the batches 1-8 coding universe the # August headline figures were computed on — the reconstruction that refutes the # story's "I cannot distinguish between two explanations." LATE_DISTRICTS = {"IA-02", "NY-01", "OH-15", "TX-34", "TX-28", "OH-09", "MT-01", "NC-01", "AL-02", "AZ-02"} # The three districts that violate the spec's printed selection rule # (2024 House margin <= 10 OR presidential margin <= 10). Their margins were TYPED here, # duplicating districts_59.csv, and then typed a third time into the page prose -- three # copies of a number with nothing holding them together. Read from the frozen file now. def _off_rule(): import csv as _csv out = {} with open(data_path(os.path.join("public_data", "districts_59.csv")) if os.path.exists(data_path(os.path.join("public_data", "districts_59.csv"))) else data_path("districts_59.csv"), encoding="utf-8") as fh: for r in _csv.DictReader(fh): if "does not meet" in (r.get("selection_rule_status") or "").lower() \ or "violates" in (r.get("selection_rule_status") or "").lower() \ or "exception" in (r.get("selection_rule_status") or "").lower(): out[r["district"]] = float(r["pres_margin_2024"]) if not out: raise SystemExit("figures: districts_59.csv marks no off-rule district. The three " "named exceptions are load-bearing for the dual-universe scoring.") return out OFF_RULE = _off_rule() def august(): """The August universe, reconstructed. August used FEC incumbency codes, not tenure > 0 — that is how its n_inc was 38 and its D-incumbent mean 13.2. HISTORICAL CODING PINNED: the August record carried Kevin Kiley as REP (the party-status flag was known but the column followed the Phase 1 source). The roster has since been corrected to IND; this reconstruction keeps the August coding, because it reproduces the August record, not today's.""" rows = M.load() rows = [dict(d, party="REP") if d["name"] == "Kevin Kiley" else d for d in rows] old = [d for d in rows if d["district"] not in LATE_DISTRICTS] inc = [d for d in old if d["inc_code"] == "I"] def pr(g): return float(sps.pearsonr(np.array([d["tenure"] for d in g], float), np.array([d["corp_pct"] for d in g], float))[0]) def rho(g): return float(sps.spearmanr([d["tenure"] for d in g], [d["corp_pct"] for d in g]).statistic) trim = sorted(inc, key=lambda d: (-d["tenure"], d["name"]))[3:] out = dict(n=len(old), n_inc=len(inc), r_pool=pr(old), rho_pool=rho(old), r_inc=pr(inc), rho_inc=rho(inc), r_trim3=pr(trim)) for pty in ("DEM", "REP"): out[f"r_inc_{pty}"] = pr([d for d in inc if d["party"] == pty]) # the means August published beside the correlations, under BOTH definitions, # on the FULL current data — this is where the vintages mix (R10's finding) full_inc_fec = [d for d in rows if d["inc_code"] == "I"] for pty in ("DEM", "REP"): out[f"mean_inc_feccode_{pty}"] = float(np.mean( [d["corp_pct"] for d in full_inc_fec if d["party"] == pty])) return out def rule_compliant(): """The universe under the spec's printed rule, i.e. WITHOUT the three named off-rule districts. Reported beside the as-published figures, never scored — the dual-universe ruling registered before October can see the answer.""" rows = [d for d in M.load() if d["district"] not in OFF_RULE] inc = [d for d in rows if d["tenure"] > 0] def pr(g): return float(sps.pearsonr(np.array([d["tenure"] for d in g], float), np.array([d["corp_pct"] for d in g], float))[0]) r_inc = pr(inc) # Round-six review: v5's ruling says the universe choice "decides two of the # verdicts" and names A1b and A1a'. A1b' flips too, and was named nowhere — # the flip table was published one row short, in the direction that favours the # universe actually chosen. All three are computed here so it cannot recur. r_pool_rc = pr(rows) all_rows0 = M.load() r_pool_pub = pr(all_rows0) a1bp_lo, a1bp_hi = (math.tanh(_z(r_pool_pub) - 0.20), math.tanh(_z(r_pool_pub) + 0.20)) # referent computed, never quoted: the as-published incumbents-only r all_rows = M.load() pub_inc = [d for d in all_rows if d["tenure"] > 0] r_pub = pr(pub_inc) return dict(n=len(rows), n_inc=len(inc), n_off_rule=len(OFF_RULE), r_pool=r_pool_rc, r_inc=r_inc, z_inc=_z(r_inc), dz_vs_published=abs(_z(r_inc) - _z(r_pub)), # A1b': is the rule-compliant pooled r outside its registered band? a1bp_lo=a1bp_lo, a1bp_hi=a1bp_hi, a1bp_flips=int(not (a1bp_lo <= r_pool_rc <= a1bp_hi)), n_flips=1 + 1 + int(not (a1bp_lo <= r_pool_rc <= a1bp_hi))) def murphy(y, p, bins=None): """Murphy decomposition of the Brier score: BS = reliability - resolution + uncertainty. Reinstated 2026-08-30. ANALYSIS_SPEC_v1 committed to publishing this; v2 dropped it six days after the register froze, without a changelog. It is the only listed statistic that separates whether my confidences are UNBIASED from whether they carry any INFORMATION, and with four dead bets at high confidence it is the one most likely to look bad. Implemented now, and self-tested, so October is a run and not a promise. y: outcomes (0/1). p: stated confidences in [0,1]. """ y = np.asarray(y, float) p = np.asarray(p, float) n = len(y) obar = float(y.mean()) groups = {} for pi, yi in zip(p, y): groups.setdefault(round(float(pi), 6), []).append(yi) rel = res = 0.0 for pk, ys in groups.items(): nk = len(ys) ok = float(np.mean(ys)) rel += nk * (pk - ok) ** 2 res += nk * (ok - obar) ** 2 rel /= n res /= n unc = obar * (1 - obar) bs = float(np.mean((p - y) ** 2)) return dict(brier=bs, reliability=rel, resolution=res, uncertainty=unc, identity_residual=abs(bs - (rel - res + unc))) def _murphy_selftest(): """The decomposition identity must hold, or the number published in October is junk.""" rng = np.random.default_rng(SEED) for _ in range(200): p = rng.choice([0.55, 0.62, 0.70, 0.78, 0.88, 0.92], size=25) y = (rng.random(25) < p).astype(float) d = murphy(y, p) assert d["identity_residual"] < 1e-12, d # a perfectly calibrated constant forecaster has zero resolution d = murphy(np.array([1., 0., 1., 0.]), np.array([0.5] * 4)) assert abs(d["resolution"]) < 1e-12 and abs(d["reliability"]) < 1e-12, d return True def style_b4(): """B4's candidate conjunct under the shipped pooling and the corrected one. style_analysis.py pools a cell's samples with a bare newline join; 72 samples (66 of them social) end without terminal punctuation, so sentences merge across sample boundaries, post-cell FK inflates, and the quote-to-post gap SHRINKS — in exactly the direction that makes B4's candidate-median conjunct hold. Declared a defect in the frozen instrument per ANALYSIS_SPEC §5.1; B4 is scored on the corrected pipeline with both values published.""" import csv as _csv import re as _re import statistics as _st path = data_path("style_corpus_all.csv") rows = list(_csv.DictReader(open(path))) def fk(text): words = _re.findall(r"[\w'@#-]+", text) n = len(words) if n == 0: return None sents = [x for x in _re.split(r"[.!?]+", text) if x.strip()] lower = [w.lower().strip("'-") for w in words if _re.match(r"[A-Za-z]", w)] syll = sum(max(1, len(_re.findall(r"[aeiouy]+", w.lower()))) for w in lower) return 0.39 * (n / max(1, len(sents))) + 11.8 * (syll / max(1, len(lower))) - 15.59 TERM = _re.compile(r"[.!?][\"\')\]]*\s*$") cells = {} for r in rows: sub = ({"quote": "quote", "spoken": "spoken"}.get(r["subtype"], "release") if r["register"] == "official" else "post") cells.setdefault((r["candidate"], sub), []).append(r["text"]) def gap(cand, fixed): vals = {} for reg in ("quote", "post"): ts = cells.get((cand, reg)) if not ts: return None if fixed: ts = [t if TERM.search(t) else t.rstrip() + "." for t in ts] vals[reg] = fk("\n".join(ts)) return vals["quote"] - vals["post"] cands = sorted({c for c, reg in cells if reg in ("quote",)}) shipped = [g for c in cands if (g := gap(c, False)) is not None] fixed = [g for c in cands if (g := gap(c, True)) is not None] # SECOND instrument defect, found 2026-08-30 by an independent reviewer. The frozen # protocol sets a per-register minimum and says a cell below it is reported as # "insufficient corpus", never as a number. The medians above pool every candidate # with both cells present, minimum or not — so cells the protocol forbids reporting # are bounding a comparison the protocol scores. # # THE MINIMUM ITSELF WAS THEN GOT WRONG, found 2026-08-30 by an independent # fact-checker: this used 300 words for BOTH registers. STYLE_PROTOCOL.md corpus # rule 4 sets "8 samples and 600 total words per register"; the v0.2 changelog # recalibrates the SOCIAL minimum only, 600 -> 300, because "eight posts of ~40 # words cannot reach 600", and says the sample minimum stays 8. Applying 300 to the # official register waives half the frozen protocol, and it waives it in the # direction that makes my own disclosure smaller: 7 candidates clear the loose bar, # 3 clear the protocol as written. Both are published. MIN_SAMPLES = 8 MIN_WORDS = {"quote": 600, "post": 300} # official / social, per the protocol MIN_WORDS_LOOSE = {"quote": 300, "post": 300} # what this function used until 08-30 def compliant(cand, minw): for reg in ("quote", "post"): ts = cells.get((cand, reg)) or [] if len(ts) < MIN_SAMPLES or sum(len(t.split()) for t in ts) < minw[reg]: return False return True ok = [c for c in cands if compliant(c, MIN_WORDS) and gap(c, True) is not None] loose = [c for c in cands if compliant(c, MIN_WORDS_LOOSE) and gap(c, True) is not None] fixed_ok = [gap(c, True) for c in ok] loose_ok = [gap(c, True) for c in loose] return dict( n_gap_compliant=len(ok), n_gap_compliant_loose=len(loose), n_gap_below_minimum=len(fixed) - len(ok), cand_median_compliant=float(_st.median(fixed_ok)) if fixed_ok else float("nan"), cand_median_loose=float(_st.median(loose_ok)) if loose_ok else float("nan"), n_samples=len(rows), n_unterminated=sum(1 for r in rows if not TERM.search(r["text"])), n_unterminated_social=sum(1 for r in rows if r["register"] != "official" and not TERM.search(r["text"])), n_gap_candidates=len(shipped), cand_median_shipped=float(_st.median(shipped)), cand_median_fixed=float(_st.median(fixed)), shift=float(_st.median(fixed) - _st.median(shipped)), b4_conjunct_shipped=int(_st.median(shipped) <= 2.5), b4_conjunct_fixed=int(_st.median(fixed) <= 2.5), ) def care(): """E4's care-professions split, recomputed rather than retyped. The number shipped on the page as a bare literal — "10 Democrats to 3 Republicans", "odds ratio 3.83", "p = 0.074" — typed out of VALIDATION_RESULTS.md and never recomputed. That is this project's exact recurring defect: a figure living in prose, going stale when its inputs move. Its input DID move. The frozen cross-tab was built when the roster read 57 DEM / 57 REP. On 2026-08-24 Kevin Kiley was corrected to IND (he left the GOP on 2026-03-09 and filed CA-06 as an independent; the roster column had fossilised the Phase 1 party), so the live roster is 57 / 56 / 2 IND and one Republican leaves the table. Under H1 his most-recent career is legislator and attorney, so he sat in REP-other, not REP-care: the care cells are untouched and only the REP denominator falls by one. The 2x3 table below is the FROZEN artifact — it is the cross-tab recorded in VALIDATION_RESULTS.md, which is covered by timestamp/COMMITMENT-2.txt. It is carried here verbatim rather than re-coded. Re-coding the 113 candidates today would be MY coding, after seeing the result, which is the precise failure the blind freeze exists to prevent (and which E4's uniformed half already demonstrated). So: frozen coding, scripted arithmetic, declared correction. Both the as-frozen and the corrected figures are published. """ # DEM, REP as frozen in VALIDATION_RESULTS.md frozen = dict(care=(10, 3), uniformed=(6, 9), other=(41, 45)) party = {n: p for lst in nominees.NOMINEES.values() for _, p, n in lst} live = {k: sum(1 for v in party.values() if v == k) for k in ("DEM", "REP", "IND")} fD = sum(v[0] for v in frozen.values()) fR = sum(v[1] for v in frozen.values()) if (fD, fR) != (57, 57): raise SystemExit("care(): the frozen cross-tab no longer sums to the 57/57 " "it was recorded under — it is frozen, so this means it was edited") # The correction is exactly one Republican, out of REP-other. corr = dict(frozen, other=(frozen["other"][0], frozen["other"][1] - 1)) cD = sum(v[0] for v in corr.values()) cR = sum(v[1] for v in corr.values()) if (cD, cR) != (live["DEM"], live["REP"]): raise SystemExit("care(): the corrected cross-tab (%d/%d) does not match the live " "roster (%d DEM / %d REP). The roster moved and this table did not." % (cD, cR, live["DEM"], live["REP"])) def stats(t): D, R = t["care"] dn, rn = (sum(v[0] for v in t.values()) - D), (sum(v[1] for v in t.values()) - R) orr, p = sps.fisher_exact([[D, dn], [R, rn]]) chi = sps.chi2_contingency([[t[k][0] for k in ("uniformed", "care", "other")], [t[k][1] for k in ("uniformed", "care", "other")]])[1] return dict(D=D, R=R, D_n=D + dn, R_n=R + rn, OR=orr, p=p, chi2_p=chi, D_pct=100 * D / (D + dn), R_pct=100 * R / (R + rn)) out = {} for k, v in stats(frozen).items(): out["frozen_" + k] = v out.update(stats(corr)) # The uniformed half, printed because the refutation is the more useful half and # a page that prints only the surviving result is doing the thing it warns about. # Oriented Republican-first, the way the register's claim was stated ("an R # uniformed-service cluster"), so the odds ratio here is the same number the # frozen record prints rather than its reciprocal. uD, uR = corr["uniformed"] fuD, fuR = frozen["uniformed"] out["u_OR"], out["u_p"] = sps.fisher_exact([[uR, cR - uR], [uD, cD - uD]]) out["u_frozen_OR"], out["u_frozen_p"] = sps.fisher_exact( [[fuR, fR - fuR], [fuD, fD - fuD]]) out["u_D"], out["u_R"] = uD, uR out["kiley_shift_OR"] = out["frozen_OR"] - out["OR"] return out def dependence(): """A5 and A1a-prime are not independent bets. How not-independent, measured. The page has been asserting "A5 holds in 20,000 of 20,000 resamples in which A1a-prime holds" as a bare literal since it was written -- a figure typed in prose with no line of code behind it, on a page whose entire argument is that numbers must be recomputed rather than retyped. Found by an independent reviewer. It is now computed here, and the count that ships is whatever this function returns. WHAT THIS IS, stated narrowly so it cannot be read as more: candidates are resampled with replacement from the frozen file, both statistics are recomputed on each resample, and the two conditions are evaluated. It measures how tightly the two conditions move together in THIS data. It is not a forecast of the October refresh, which is the same people with one more quarter of receipts, not a fresh draw. The dependence is the point: two cards that rise and fall together are close to one bet counted twice, and a calibration score that treats them as two is flattered by it. """ rows = M.load() inc = [d for d in rows if d["tenure"] > 0] x = np.array([d["tenure"] for d in inc], float) y = np.array([d["corp_pct"] for d in inc], float) xa = np.array([d["tenure"] for d in rows], float) ya = np.array([d["corp_pct"] for d in rows], float) z_inc0, z_pool0 = _z(float(np.corrcoef(x, y)[0, 1])), _z(float(np.corrcoef(xa, ya)[0, 1])) rng = np.random.default_rng(SEED) N = 20000 n_a, n_both, n_a5 = 0, 0, 0 ii = np.array([k for k, d in enumerate(rows) if d["tenure"] > 0]) for _ in range(N): j = rng.integers(0, len(rows), len(rows)) xs, ys = xa[j], ya[j] m = xs > 0 if m.sum() < 3 or np.std(xs) == 0 or np.std(ys) == 0 or np.std(xs[m]) == 0 \ or np.std(ys[m]) == 0: continue zi = _z(float(np.corrcoef(xs[m], ys[m])[0, 1])) zp = _z(float(np.corrcoef(xs, ys)[0, 1])) a1 = abs(zi - z_inc0) <= 0.20 # A1a-prime: incumbents-only r stays put a5 = (zp - zi) >= 0.20 # A5: pooled z exceeds incumbents-only by 0.20 n_a += a1 n_a5 += a5 n_both += (a1 and a5) return dict(N=N, n_a1ap=int(n_a), n_a5=int(n_a5), n_both=int(n_both), pct_given=100.0 * n_both / n_a if n_a else float("nan"), pct_a5=100.0 * n_a5 / N) def validation(): """The inter-coder replication, READ from the frozen validation file. "46 of 50 (92%)" and "alpha 0.91" were typed onto the page three times. A red-teamer changed one to "49 of 50 (98%)" and the build passed, shipping a page that contradicted the timestamped file beside it. VALIDATION_RESULTS.md is covered by timestamp/COMMITMENT-2.txt, so reading it here means the page cannot disagree with the record it cites without failing the commitment gate first.""" for base in (ROOT_DIR, HERE, os.path.join(ROOT_DIR, "public_data")): c = os.path.join(base, "VALIDATION_RESULTS.md") if os.path.exists(c): break else: raise SystemExit("figures: VALIDATION_RESULTS.md not found — the page cites it") txt = open(c, encoding="utf-8").read() m = re.search(r"(\d+)\s*(?:of|/)\s*(\d+)\s+sampled", txt) or \ re.search(r"\b(\d+)\s*(?:of|/)\s*(\d+)\b[^.\n]{0,60}(?:decision|agree)", txt, re.I) a = re.search(r"(?:alpha|α)\s*=?\s*([01]\.\d+)", txt) if not m or not a: raise SystemExit("figures: could not read the replication result out of " "VALIDATION_RESULTS.md — the page must not state it from memory") agree, n = int(m.group(1)), int(m.group(2)) return dict(agree=agree, n=n, pct=100.0 * agree / n, alpha=float(a.group(1))) def laplace(): """The Laplace rule on the review rounds run so far, computed rather than typed.""" rounds = 6 # blind review rounds completed before F1 was priced return dict(rounds=rounds, pct=100.0 * (rounds + 1) / (rounds + 2)) @functools.lru_cache(maxsize=1) def step_and_slope(): """share ~ b0 + b1*served + b2*tenure, OLS, all rows with a finance row. `served` is 1 for any prior congressional service, 0 otherwise; `tenure` is years, which is 0 for everyone with served=0. So b1 is the STEP -- the jump on arriving -- and b2 is the SLOPE, the annual gain among those who have served. The two coefficients answer the question the page has been arguing about with correlations, which cannot distinguish them: r is one number and this is two. NO CAUSAL CLAIM. This is a cross-section of 112 candidates with no controls for committee position, majority status, seat safety, district industry mix or total receipts. It describes the shape of the association. It does not explain it, and nothing here says which way any arrow points. """ rows = M.load() x_served = np.array([1.0 if d["tenure"] > 0 else 0.0 for d in rows]) x_ten = np.array([float(d["tenure"]) for d in rows]) y = np.array([float(d["corp_pct"]) for d in rows]) X = np.column_stack([np.ones(len(y)), x_served, x_ten]) beta, *_ = np.linalg.lstsq(X, y, rcond=None) resid = y - X @ beta n, k = X.shape dof = n - k sigma2 = float(resid @ resid) / dof cov = sigma2 * np.linalg.inv(X.T @ X) se = np.sqrt(np.diag(cov)) t = beta / se pvals = [float(2 * sps.t.sf(abs(v), dof)) for v in t] # HETEROSKEDASTICITY-ROBUST inference. The dependent variable is a corporate-PAC SHARE # with a large point mass at ~0 (most challengers take almost none), so the residuals # are heteroskedastic by construction and the classical SE above understates the # coefficients' true uncertainty. The step is large enough that it survives either way; # the SLOPE's significance does not, and reporting only the classical p on the slope # overstates it. HC3 (the most conservative common sandwich) and a pairs bootstrap are # computed here so the page can state the slope honestly. Added 2026-08-31 after a # review found the page called the per-year slope "measurable, p = 0.0011" when it is # not significant under robust inference. XtXinv = np.linalg.inv(X.T @ X) hat = np.einsum("ij,jk,ik->i", X, XtXinv, X) # leverages h_ii w_hc3 = (resid ** 2) / (1.0 - hat) ** 2 cov_hc3 = XtXinv @ (X.T @ (X * w_hc3[:, None])) @ XtXinv se_hc3 = np.sqrt(np.diag(cov_hc3)) p_hc3 = [float(2 * sps.t.sf(abs(beta[i] / se_hc3[i]), dof)) for i in range(k)] rng = np.random.default_rng(20260831) # fixed: the figure must reproduce B = 8000 boot = np.empty(B) for i in range(B): idx = rng.integers(0, n, n) bb, *_ = np.linalg.lstsq(X[idx], y[idx], rcond=None) boot[i] = bb[2] slope_boot_lo = float(np.percentile(boot, 2.5)) slope_boot_hi = float(np.percentile(boot, 97.5)) slope_boot_frac_pos = float((boot > 0).mean()) # What the slope is worth across the observed range of service, so the two # coefficients are comparable in the units a reader thinks in: percentage points. span = float(x_ten.max()) out = dict(n=n, dof=dof, b0=float(beta[0]), b1_step=float(beta[1]), b2_slope=float(beta[2]), se0=float(se[0]), se1=float(se[1]), se2=float(se[2]), t1=float(t[1]), t2=float(t[2]), p1=pvals[1], p2=pvals[2], step_lo=float(beta[1] - 1.96 * se[1]), step_hi=float(beta[1] + 1.96 * se[1]), slope_lo=float(beta[2] - 1.96 * se[2]), slope_hi=float(beta[2] + 1.96 * se[2]), span=span, slope_worth=float(beta[2]) * span, slope_worth_lo=float(beta[2] - 1.96 * se[2]) * span, slope_worth_hi=float(beta[2] + 1.96 * se[2]) * span, r2=float(1 - (resid @ resid) / (((y - y.mean()) ** 2).sum())), sigma=float(np.sqrt(sigma2)), # Years of service at which the accumulated slope equals the step. The # single most useful number the fit produces: below it the jump on # arriving dominates, above it seniority does, and it says in one figure # which of the two the phrase "the pattern" should mean. crossover=float(beta[1] / beta[2]), n_above_crossover=int((x_ten > float(beta[1] / beta[2])).sum()), # robust inference (HC3 + pairs bootstrap) se1_hc3=float(se_hc3[1]), se2_hc3=float(se_hc3[2]), p1_hc3=float(p_hc3[1]), p2_hc3=float(p_hc3[2]), step_hc3_lo=float(beta[1] - 1.96 * se_hc3[1]), step_hc3_hi=float(beta[1] + 1.96 * se_hc3[1]), slope_hc3_lo=float(beta[2] - 1.96 * se_hc3[2]), slope_hc3_hi=float(beta[2] + 1.96 * se_hc3[2]), slope_boot_lo=slope_boot_lo, slope_boot_hi=slope_boot_hi, slope_boot_frac_pos=slope_boot_frac_pos) # And the same fit with the highest-leverage point removed, because this page names # that point and keeps it everywhere else and must here too. top = int(np.argmax(x_ten)) keep = [i for i in range(n) if i != top] Xk, yk = X[keep], y[keep] bk, *_ = np.linalg.lstsq(Xk, yk, rcond=None) rk = yk - Xk @ bk ck = (float(rk @ rk) / (len(yk) - k)) * np.linalg.inv(Xk.T @ Xk) # WHERE THE SLOPE'S ROBUST UNCERTAINTY COMES FROM. An earlier build blamed the # "floor of zeros" (the never-served) for the slope's wide HC3 interval. A reviewer # decomposed it and that was wrong: never-served rows have zero tenure, so they carry # NO information about a per-year rate and contribute exactly nothing to the slope's # variance. The HC3 variance of the slope is a sum of per-observation terms # a_i^2 * e_i^2 / (1-h_i)^2 (a = the slope row of (X'X)^-1 X'), and one observation — # the longest-serving member, high leverage AND a large residual — is most of it. # Printed on the page, with the same fit refitted without her under HC3, so the reader # can see that the slope's significance hinges on one 44-year career, not on the zeros. A2 = (XtXinv @ X.T)[2] contrib = (A2 ** 2) * w_hc3 share_top = float(contrib[top] / contrib.sum()) share_never = float(contrib[x_served == 0].sum() / contrib.sum()) Xi = np.linalg.inv(Xk.T @ Xk) hk = np.einsum("ij,jk,ik->i", Xk, Xi, Xk) wk = (rk ** 2) / (1.0 - hk) ** 2 ck_hc3 = Xi @ (Xk.T @ (Xk * wk[:, None])) @ Xi se2_drop_hc3 = float(np.sqrt(np.diag(ck_hc3))[2]) dofk = len(yk) - k out.update(b1_step_drop=float(bk[1]), b2_slope_drop=float(bk[2]), se2_drop=float(np.sqrt(np.diag(ck))[2]), drop_name=rows[top]["name"], slope_hc3_share_top=share_top, slope_hc3_share_never=share_never, se2_drop_hc3=se2_drop_hc3, p2_drop_hc3=float(2 * sps.t.sf(abs(bk[2] / se2_drop_hc3), dofk)), slope_drop_hc3_lo=float(bk[2] - 1.96 * se2_drop_hc3), slope_drop_hc3_hi=float(bk[2] + 1.96 * se2_drop_hc3)) return out def build(): f = {} _flat("money.", money(), f) _flat("roster.", roster(), f) _flat("score.", scoring(), f) _flat("vintage.", vintage(), f) _flat("missing.", missing_finance(), f) _flat("august.", august(), f) _flat("rule.", rule_compliant(), f) _flat("style.", style_b4(), f) _flat("care.", care(), f) _flat("dep.", dependence(), f) _flat("val.", validation(), f) _flat("laplace.", laplace(), f) _flat("fit.", step_and_slope(), f) _flat("offrule.", {k.replace("-", "_"): v for k, v in OFF_RULE.items()}, f) _flat("conf.", CONF, f) return f FIGURES = None def get(): global FIGURES if FIGURES is None: FIGURES = build() return FIGURES if __name__ == "__main__": F = get() if "--json" in sys.argv: # Beside this file if there is no sibling web/ — the published bundle is flat, # and a reader who unpacks it anywhere else got a FileNotFoundError from a # command RUN_ME.md tells them to run. make_figures.py already guarded this; # this one did not. web = os.path.normpath(os.path.join(HERE, "..", "web")) out = (os.path.join(web, "figures.json") if os.path.isdir(web) else os.path.join(HERE, "figures.json")) json.dump(F, open(out, "w"), indent=1, default=float, sort_keys=True) print("wrote", out) for k in sorted(F): v = F[k] if isinstance(v, dict): print(f"{k:34s} " + " ".join(f"{a}={b:.4f}" if isinstance(b, float) else f"{a}={b}" for a, b in v.items())) elif isinstance(v, float): print(f"{k:34s} {v:.6f}") else: print(f"{k:34s} {v}")