"""exercise_0a.py Exercise 0-A. Spatial dependence before anything else. Runs the whole exercise and writes both the report and the figures: 1. Moran's I on real data under five weights schemes 2. The same measure on three synthetic fields with known answers 3. Getis-Ord Gi* with and without a false discovery correction 4. Local Moran cluster map 5. The last mile handoff package for a human cartographer python exercise_0a.py """ from __future__ import annotations import json import sys from pathlib import Path import matplotlib import numpy as np import pandas as pd matplotlib.use("Agg") import matplotlib.pyplot as plt # noqa: E402 from matplotlib.collections import PolyCollection # noqa: E402 from scipy import ndimage # noqa: E402 sys.path.insert(0, str(Path(__file__).resolve().parent)) from spatial_stats import (WEIGHT_SCHEMES, build_weights, fdr_mask, # noqa getis_ord, local_morans, morans_i, read_polygons) ROOT = Path(__file__).resolve().parent.parent SNAP, FIG, OUT = ROOT / "snapshots", ROOT / "figures", ROOT / "handoff" GEO = SNAP / "ms_tracts.geojson" SEED = 20260824 ALPHA = 0.05 QCOLOR = {"HH": "#b2182b", "LL": "#2166ac", "HL": "#ef8a62", "LH": "#67a9cf", "ns": "#e8e8e8"} def load_values(): df = pd.read_csv(SNAP / "acs2023_mississippi_asi.csv", dtype={"code": str, "sumlev": str}) tr = df[df.sumlev == "140"].copy() tr["pct_pov"] = 100 * (tr.C17002_E002 + tr.C17002_E003) / tr.C17002_E001 return tr.set_index("code")["pct_pov"] def draw(ax, rings, facecolors, title, subtitle="", source=""): polys, fc = [], [] for r, c in zip(rings, facecolors): for ring in r: polys.append(ring) fc.append(c) ax.add_collection(PolyCollection(polys, facecolors=fc, edgecolors="white", linewidths=0.2)) ax.autoscale_view() ax.set_aspect("equal") ax.axis("off") ax.set_title(title, fontsize=11, loc="left", fontweight="bold") if subtitle: ax.text(0, 1.005, subtitle, transform=ax.transAxes, fontsize=8, va="bottom", color="#444444") if source: ax.text(0, -0.03, source, transform=ax.transAxes, fontsize=7, va="top", color="#666666") def synthetic_block(n=30): rng = np.random.default_rng(SEED) cb = np.indices((n, n)).sum(0) % 2 fields = { "random": rng.normal(size=(n, n)), "dispersed": cb * 2.0 + rng.normal(scale=0.3, size=(n, n)), "clustered": ndimage.gaussian_filter( rng.normal(size=(n, n)), sigma=3) * 8, } idx = np.arange(n * n).reshape(n, n) Wr = np.zeros((n * n, n * n)) Wq = np.zeros((n * n, n * n)) for dr in (-1, 0, 1): for dc in (-1, 0, 1): if dr == dc == 0: continue r0, r1 = max(0, dr), min(n, n + dr) c0, c1 = max(0, dc), min(n, n + dc) a = idx[r0:r1, c0:c1].ravel() b = idx[r0 - dr:r1 - dr, c0 - dc:c1 - dc].ravel() Wq[a, b] = 1.0 if dr == 0 or dc == 0: Wr[a, b] = 1.0 return fields, Wr, Wq def main(): for d in (FIG, OUT): d.mkdir(exist_ok=True) ids, rings = read_polygons(GEO) vals = load_values() x = np.array([vals.get(i, np.nan) for i in ids], dtype=float) keep = np.isfinite(x) ids_k = [i for i, k in zip(ids, keep) if k] rings_k = [r for r, k in zip(rings, keep) if k] xk = x[keep] import hashlib for snap in (GEO, SNAP / "acs2023_mississippi_asi.csv"): print(f"snapshot {snap.name} sha256 " f"{hashlib.sha256(snap.read_bytes()).hexdigest()}") print(f"Mississippi tracts with a poverty estimate: {len(xk)}\n") # ---- 1. weights matrix sensitivity on REAL data -------------------- print("[1] GLOBAL MORAN'S I UNDER FIVE WEIGHTS SCHEMES") print(" Same data. The analyst picks the scheme; the tool rarely asks.") print(f" {'scheme':<12}{'neighbors':>11}{'I':>9}{'z':>9}" f"{'p':>9} verdict") Ws, rows = {}, [] for s in WEIGHT_SCHEMES: W = build_weights(rings_k, s) Ws[s] = W m = morans_i(xk, W, permutations=499, seed=SEED) avg = float((W > 0).sum(1).mean()) v = ("clustered" if m["p_sim"] < ALPHA and m["I"] > m["EI"] else "dispersed" if m["p_sim"] < ALPHA else "random") rows.append({"scheme": s, "neighbors": avg, **m, "verdict": v}) print(f" {s:<12}{avg:>11.2f}{m['I']:>+9.4f}{m['z_sim']:>+9.2f}" f"{m['p_sim']:>9.3f} {v}") # ---- 2. synthetic fields with known answers ------------------------ print("\n[2] THE SAME MEASURES ON DATA WHOSE ANSWER IS KNOWN") fields, Wr, Wq = synthetic_block() print(f" {'field':<11}{'rook I':>9}{'queen I':>10}" f"{'Gi* sig':>9}{'after FDR':>11}") syn = {} for name, f in fields.items(): v = f.ravel().astype(float) ir = morans_i(v, Wr, permutations=199, seed=SEED) iq = morans_i(v, Wq, permutations=199, seed=SEED) g = getis_ord(v, Wq) sig = int((g["p"] < ALPHA).sum()) kept = int(fdr_mask(g["p"], ALPHA).sum()) syn[name] = (ir["I"], iq["I"], sig, kept) print(f" {name:<11}{ir['I']:>+9.4f}{iq['I']:>+10.4f}" f"{sig:>9d}{kept:>11d}") print(f" 900 units; chance alone puts about {ALPHA*900:.0f} outside " f"the critical values at alpha {ALPHA}") # ---- 3. Gi* on real data, corrected and uncorrected ---------------- print("\n[3] GETIS-ORD Gi* ON THE REAL DATA, QUEEN CONTIGUITY") g = getis_ord(xk, Ws["queen"]) raw_hot = g["z"] > 1.96 raw_cold = g["z"] < -1.96 keptm = fdr_mask(g["p"], ALPHA) print(f" uncorrected hot {int(raw_hot.sum())}, " f"cold {int(raw_cold.sum())}, total " f"{int(raw_hot.sum() + raw_cold.sum())}") print(f" surviving Benjamini-Hochberg at q={ALPHA}: " f"{int(keptm.sum())}") print(f" lost to the correction: " f"{int(raw_hot.sum() + raw_cold.sum() - keptm.sum())}") # ---- 4. LISA ------------------------------------------------------- lm = local_morans(xk, Ws["queen"], permutations=499, seed=SEED) sig = lm["p"] < ALPHA print("\n[4] LOCAL MORAN CLUSTER TYPES, p < 0.05") for q in ("HH", "LL", "HL", "LH"): print(f" {q} {int(((lm['quadrant'] == q) & sig).sum()):>4d}") print(f" not significant {int((~sig).sum()):>4d}") # ---- figures ------------------------------------------------------- fig, axes = plt.subplots(1, 3, figsize=(18, 6.5)) src = ("Source: U.S. Census Bureau, ACS 2023 5-year. " f"{len(xk)} tracts. EPSG:5070.") draw(axes[0], rings_k, [QCOLOR[q] if s else QCOLOR["ns"] for q, s in zip(lm["quadrant"], sig)], "Local Moran cluster map", "Queen contiguity, 499 permutations, p < 0.05", src) handles = [plt.Rectangle((0, 0), 1, 1, fc=QCOLOR[q]) for q in ("HH", "LL", "HL", "LH", "ns")] axes[0].legend(handles, ["high near high", "low near low", "high near low", "low near high", "not significant"], loc="upper left", bbox_to_anchor=(1.01, 1.0), fontsize=8, frameon=False) draw(axes[1], rings_k, ["#b2182b" if h else "#2166ac" if c else "#e8e8e8" for h, c in zip(raw_hot, raw_cold)], "Gi* hot spots, uncorrected", f"{int(raw_hot.sum()+raw_cold.sum())} tracts significant at 0.05", src) draw(axes[2], rings_k, ["#b2182b" if (h and k) else "#2166ac" if (c and k) else "#e8e8e8" for h, c, k in zip(raw_hot, raw_cold, keptm)], "Gi* hot spots, false discovery corrected", f"{int(keptm.sum())} tracts survive Benjamini-Hochberg", src) fig.suptitle("Cluster detection depends on choices the tool does not ask " "you to make", fontsize=14, fontweight="bold") fig.tight_layout() fig.savefig(FIG / "exercise_0a_clusters.png", dpi=140, bbox_inches="tight") plt.close(fig) # ---- 5. the last mile handoff -------------------------------------- hand = pd.DataFrame({ "geoid": ids_k, "pct_pov": xk, "gi_z": g["z"], "gi_p": g["p"], "gi_significant_raw": raw_hot | raw_cold, "gi_significant_fdr": keptm, "lisa_quadrant": lm["quadrant"], "lisa_p": lm["p"], "lisa_significant": sig, }) hand.to_csv(OUT / "exercise_0a_for_review.csv", index=False) (OUT / "exercise_0a_decisions.json").write_text(json.dumps({ "decisions_the_machine_made_without_being_asked": { "weights_scheme": "queen contiguity", "row_standardized": True, "permutations": 499, "alpha": ALPHA, "multiple_comparison_correction": "Benjamini-Hochberg", "class_count_on_maps": 5, "color_ramp": "diverging red and blue, unjustified", "projection": "EPSG:5070", }, "what_a_human_must_check_before_this_map_is_published": [ "Is queen contiguity right for this phenomenon, or does the " "process travel along roads and rivers instead of across " "shared borders?", "Do the tracts with no estimate bias the cluster pattern, and " "should they be shown as absent rather than as background?", "Does the diverging ramp imply a midpoint that means anything?", "Are the surviving clusters explicable on the ground, and does " "local knowledge contradict any of them?", "Should the uncorrected map be published at all?", ], "global_moran_by_scheme": {r["scheme"]: round(r["I"], 4) for r in rows}, }, indent=2)) print(f"\n[5] LAST MILE HANDOFF written to {OUT}") print(" exercise_0a_for_review.csv every unit, every statistic") print(" exercise_0a_decisions.json what the machine chose, and " "what a human must check") print(f"\nwrote {FIG / 'exercise_0a_clusters.png'}") if __name__ == "__main__": main()