"""Independent re-computation of the Vellora metropolitan–nonmetropolitan study. Recomputes every statistic in public/data/research/rural-metro-results.v1.json from the raw verified input files with SciPy and NumPy (a separate implementation from src/lib/research/rural-study.ts) and exits non-zero on any mismatch. Run from the repository root, in a throwaway environment: uvx --with scipy --with numpy python scripts/verify-rural-study.py """ import json import math import sys import numpy as np import scipy from scipy import stats load = lambda p: json.load(open(p, encoding="utf-8")) cardio = load("public/data/places-tx-cardio.v2025.1.json") places = load("public/data/places-tx.v2025.1.json") heart = load("public/data/cdc-heart-mortality-tx.2022-2024.v1.json") cancer = load("public/data/scp-cancer-incidence-tx.2019-2023.v1.json") facts = load("public/data/research/cardio-disparities-inputs.v1.json") res = load("public/data/research/rural-metro-results.v1.json") order = facts["countyOrder"] assert all(f["countyOrder"] == order for f in (cardio, places, heart, cancer)) P25 = "places-2025|all" cell = lambda f, i, per=P25: np.array([np.nan if v is None else v for v in f["observations"][i][per]["values"]], dtype=float) AA = { "HDMORT": cell(heart, "hd-mort", "hd-2022-2024|all"), "LUNG": cell(cancer, "cancer-inc-lung", "scp-2019-2023|all"), "CHD": cell(cardio, "chd-aa"), "BPHIGH": cell(cardio, "bphigh-aa"), "DIABETES": cell(cardio, "diabetes-aa"), "OBESITY": cell(cardio, "obesity-aa"), "CSMOKING": cell(cardio, "smoking-aa"), "LPA": cell(cardio, "lpa-aa"), } CRUDE = {"CHD": cell(places, "chd"), "BPHIGH": cell(places, "bphigh"), "DIABETES": cell(places, "diabetes"), "OBESITY": cell(places, "obesity"), "CSMOKING": cell(cardio, "smoking"), "LPA": cell(cardio, "lpa")} rucc = np.array([c["rucc"] for c in facts["counties"]]) pop = np.array([c["pop2020"] for c in facts["counties"]]) metro = rucc <= 3 failures = [] def check(label, got, want, tol=5e-5): if got is None or want is None or not math.isclose(got, want, rel_tol=0, abs_tol=tol): failures.append(f"{label}: results file {got} vs recomputed {want}") def pcheck(label, got, want): check(label, got, want, max(1e-12, 1e-6 * want)) def hl(a, b): d = np.sort(np.subtract.outer(a, b).ravel()) n1, n2 = len(a), len(b) c = math.floor(n1 * n2 / 2 - stats.norm.ppf(0.975) * math.sqrt(n1 * n2 * (n1 + n2 + 1) / 12)) return np.median(d), (d[c - 1], d[len(d) - c]) def comparison(label, a, b, r): mw = stats.mannwhitneyu(a, b, alternative="two-sided", method="asymptotic", use_continuity=True) est, ci = hl(a, b) check(f"{label} n nonmetro", r["nNonmetro"], len(a), 0) check(f"{label} n metro", r["nMetro"], len(b), 0) check(f"{label} median difference", r["medianDifference"], np.median(a) - np.median(b)) check(f"{label} HL", r["hl"], est) check(f"{label} HL lo", r["hlCI"][0], ci[0]) check(f"{label} HL hi", r["hlCI"][1], ci[1]) check(f"{label} U", r["U"], mw.statistic, 1e-9) check(f"{label} P", r["probSuperiority"], mw.statistic / (len(a) * len(b))) pcheck(f"{label} p", r["p"], mw.pvalue) check("metro counties", res["counties"]["metro"], int(metro.sum()), 0) check("nonmetro counties", res["counties"]["nonmetro"], int((~metro).sum()), 0) raw = {} for r in res["indicators"]: v, i = AA[r["id"]], r["id"] ok = ~np.isnan(v) a, b = v[ok & ~metro], v[ok & metro] for g, x, mask in (("nonmetro", a, ~metro), ("metro", b, metro)): s = r[g] check(f"{i} {g} n", s["n"], len(x), 0) check(f"{i} {g} suppressed", s["suppressed"], int((mask & ~ok).sum()), 0) check(f"{i} {g} median", s["median"], np.median(x)) check(f"{i} {g} q1", s["q1"], np.quantile(x, 0.25)) check(f"{i} {g} q3", s["q3"], np.quantile(x, 0.75)) check(f"{i} {g} mean", s["mean"], x.mean()) check(f"{i} {g} sd", s["sd"], x.std(ddof=1)) check(f"{i} {g} min", s["min"], x.min(), 0) check(f"{i} {g} max", s["max"], x.max(), 0) comparison(i, a, b, r["comparison"]) raw[i] = r["comparison"]["p"] m = r["missingness"] if (~ok).any(): sn, sm = int((~ok & ~metro).sum()), int((~ok & metro).sum()) pcheck(f"{i} Fisher p", m["fisherP"], stats.fisher_exact([[sn, int((~metro).sum()) - sn], [sm, int(metro.sum()) - sm]]).pvalue) check(f"{i} median pop suppressed", m["medianPop2020Suppressed"], np.median(pop[~ok]), 1e-9) check(f"{i} median pop published", m["medianPop2020Published"], np.median(pop[ok]), 1e-9) elif m is not None: failures.append(f"{i}: missingness reported without suppression") for fam in (["HDMORT", "LUNG"], ["CHD", "BPHIGH", "DIABETES", "OBESITY", "CSMOKING", "LPA"]): adj = stats.false_discovery_control([raw[i] for i in fam], method="bh") for i, p in zip(fam, adj): pcheck(f"{i} BH-adjusted p", next(r for r in res["indicators"] if r["id"] == i)["comparison"]["pAdj"], p) e = res["exploratory"] for r in e["crude"]: v = CRUDE[r["id"]] ok = ~np.isnan(v) comparison(f"crude {r['id']}", v[ok & ~metro], v[ok & metro], r) for r in e["gradient"]: v = AA[r["id"]] ok = ~np.isnan(v) rho, p = stats.spearmanr(rucc[ok], v[ok]) check(f"gradient {r['id']} rho", r["rho"], rho) pcheck(f"gradient {r['id']} p", r["p"], p) for r in e["adjacency"]: v = AA[r["id"]] ok = ~np.isnan(v) & ~metro adj, non = v[ok & np.isin(rucc, [4, 6, 8])], v[ok & np.isin(rucc, [5, 7, 9])] est, ci = hl(non, adj) check(f"adjacency {r['id']} n", r["nAdjacent"] + r["nNonadjacent"], len(adj) + len(non), 0) check(f"adjacency {r['id']} HL", r["hl"], est) check(f"adjacency {r['id']} HL lo", r["hlCI"][0], ci[0]) pcheck(f"adjacency {r['id']} p", r["p"], stats.mannwhitneyu(non, adj, alternative="two-sided", method="asymptotic", use_continuity=True).pvalue) if failures: print(f"RURAL STUDY CHECK FAILED ({len(failures)}):\n" + "\n".join(failures[:40])) sys.exit(1) print(f"OK: every statistic in the rural study results file reproduced with SciPy {scipy.__version__} and NumPy {np.__version__}: " f"{res['counties']['metro']} metropolitan and {res['counties']['nonmetro']} nonmetropolitan counties; {len(res['indicators'])} indicators " f"(descriptives, Hodges–Lehmann shifts and CIs, Mann–Whitney, BH within families, Fisher missingness tests), crude, gradient and adjacency analyses.")