"""Does the published RR carry a composition term? The shuffle test. Written 2026-09-19 (session 70), after the rain-before-seven study found that a pooled difference of two conditional rates reads a confident effect when the exposure labels are shuffled inside their blocks and carry no information at all. That defect is a property of the ESTIMATOR, not of the data, and this pack's headline — RR 1.0231 for a follow-up within 168 h after a Friday release — is a pooled ratio of two conditional rates over a very heterogeneous population of projects. So it has the shape. This asks whether it has the defect. python research/friday-releases-replication/audit-shuffle.py The test: permute each release's weekday label WITHIN ITS OWN PROJECT. Every project keeps its exact weekday multiset and its exact set of outcomes; only the pairing between them is destroyed. Under a world where the weekday means nothing, the pooled RR must come back at 1.0000. Whatever it comes back as is the part of the published number that is composition rather than effect. Reported whatever it says. """ import collections import datetime as dt import random import sys import analyse import rule REPS = 200 SEED = 20260919 FRI, MONTHU = {4}, {0, 1, 2, 3} def rows_with_project(rows, window_h): """analyse.follow_ups, but keeping the project each release belongs to.""" by_project = collections.defaultdict(list) for project, t, v, yanked in rows: by_project[project].append((t, v, yanked)) horizon = analyse.END - dt.timedelta(hours=rule.HORIZON_H) out = [] for project, times in by_project.items(): times.sort() for i, (t, v, yanked) in enumerate(times): if t > horizon: continue nxt = times[i + 1] if i + 1 < len(times) else None inside = nxt is not None and (nxt[0] - t) <= dt.timedelta(hours=window_h) out.append((project, t.weekday(), inside, yanked)) return out def pooled_rr(recs, day_field=1, out_field=2): k1 = n1 = k2 = n2 = 0 for r in recs: d = r[day_field] if d in FRI: n1 += 1 k1 += bool(r[out_field]) elif d in MONTHU: n2 += 1 k2 += bool(r[out_field]) return rule.rr_ci(k1, n1, k2, n2), (k1, n1, k2, n2) def shuffle_within_project(recs, rng): by_project = collections.defaultdict(list) for i, r in enumerate(recs): by_project[r[0]].append(i) days = [r[1] for r in recs] for idxs in by_project.values(): pool = [days[i] for i in idxs] rng.shuffle(pool) for i, d in zip(idxs, pool): days[i] = d return [(r[0], days[i], r[2], r[3]) for i, r in enumerate(recs)] def run(recs, label, out_field): (rr, lo, hi), (k1, n1, k2, n2) = pooled_rr(recs, out_field=out_field) print(f"\n{label}") print(f" observed RR = {rr:.4f} CI [{lo:.4f}, {hi:.4f}] " f"Friday {k1:,}/{n1:,}, Mon-Thu {k2:,}/{n2:,}") rng = random.Random(SEED) vals = [] for _ in range(REPS): s = shuffle_within_project(recs, rng) vals.append(pooled_rr(s, out_field=out_field)[0][0]) vals.sort() mean = sum(vals) / len(vals) p2, p97 = vals[int(0.025 * len(vals))], vals[int(0.975 * len(vals)) - 1] print(f" shuffled RR = {mean:.4f} 2.5-97.5 % of {REPS} shuffles " f"[{p2:.4f}, {p97:.4f}]") off = mean - 1.0 print(f" composition term = {off:+.4f} " f"({off / (rr - 1.0) * 100:+.1f} % of the published effect)" if abs(rr - 1) > 1e-9 else "") verdict = ("CLEAN — the estimator reads 1.0000 on information-free labels" if p2 <= 1.0 <= p97 else "CONFOUNDED — the estimator reads an effect on information-free labels") print(f" {verdict}") return p2 <= 1.0 <= p97 def mh_rr(recs, out_field, projects=None): """Mantel-Haenszel risk ratio, strata = projects. RR_MH = sum_i (k1_i * n2_i / n_i) / sum_i (k2_i * n1_i / n_i), over strata i. A project that releases only on Fridays contributes nothing to either sum, which is exactly the composition the pooled ratio was reading. """ num = den = 0.0 per = collections.defaultdict(lambda: [0, 0, 0, 0]) # k1 n1 k2 n2 for r in recs: c = per[r[0]] if r[1] in FRI: c[1] += 1 c[0] += bool(r[out_field]) elif r[1] in MONTHU: c[3] += 1 c[2] += bool(r[out_field]) for p in (projects if projects is not None else per): k1, n1, k2, n2 = per[p] n = n1 + n2 if n == 0 or n1 == 0 or n2 == 0: continue num += k1 * n2 / n den += k2 * n1 / n return float("nan") if den == 0 else num / den def mh_ci(recs, out_field, reps=500, seed=SEED): """Percentile CI by resampling PROJECTS with replacement.""" rng = random.Random(seed) by_project = collections.defaultdict(list) for r in recs: by_project[r[0]].append(r) names = list(by_project) vals = [] for _ in range(reps): pick = [names[rng.randrange(len(names))] for _ in names] sub, seen = [], collections.Counter() for p in pick: seen[p] += 1 tag = f"{p}#{seen[p]}" # a project drawn twice is two strata sub.extend((tag, r[1], r[2], r[3]) for r in by_project[p]) v = mh_rr(sub, out_field) if v == v: vals.append(v) vals.sort() return vals[int(0.025 * len(vals))], vals[int(0.975 * len(vals)) - 1] def corrected(recs, label, out_field, published): rr = mh_rr(recs, out_field) lo, hi = mh_ci(recs, out_field) print(f"\n{label}") print(f" published (pooled) RR = {published}") print(f" corrected (within project) RR = {rr:.4f} CI [{lo:.4f}, {hi:.4f}]") print(" the corrected interval " + ("EXCLUDES 1 — the effect survives stratification" if lo > 1 or hi < 1 else "CONTAINS 1 — what was left after stratification is not distinguishable " "from nothing")) return rr, lo, hi def main(): rows, drops = analyse.load() print(f"{drops['kept']:,} releases kept of {drops['rows read']:,} rows") recs = rows_with_project(rows, 168) print(f"{len(recs):,} releases scored (inside the 168 h horizon)") print(f"{len({r[0] for r in recs}):,} projects — the blocks") ok1 = run(recs, "THE HEADLINE: any follow-up within 168 h, Friday vs Mon-Thu", out_field=2) ok2 = run(recs, "THE YANK RATE: a release later withdrawn, Friday vs Mon-Thu", out_field=3) print("\n" + ("Both estimators are clean under the shuffle." if ok1 and ok2 else "At least one estimator moves on information-free labels.")) print("\n" + "=" * 70) print("THE CORRECTED NUMBERS, stratified by project") corrected(recs, "any follow-up within 168 h, Friday vs Mon-Thu", 2, "1.0231 CI [1.0113, 1.0350]") corrected(recs, "a release later withdrawn, Friday vs Mon-Thu", 3, "1.2133 CI [1.0081, 1.4603] (reported, never sealed)") return 0 if __name__ == "__main__": sys.exit(main())