"""asi_toolkit.py Reference implementation of the Aggregation Sensitivity Index that Chapter 4 of GeoAI Risks suggests (section 4.4) and Chapter 5 proposes as a family of measures (section 5.6.5). The toolkit runs three tests the book names in sections 4.4, 5.6 and 6.7, a rank-order stability check drawn from the index components section 4.4 lists, and an evidence resolution check drawn from section 4.5, then closes with the index itself: 1. Multiscale validation (4.4, 5.6.1, 6.7.1) 2. Boundary perturbation testing (4.4, 5.6.2, 6.7.2) 3. Cross-zonal transfer testing (4.4, 5.6.3, 6.7.3) 4. Rank-order stability (4.4, an index component) 5. Evidence resolution (4.5 recommendation 3, 6.6.3) then the Aggregation Sensitivity Index (4.4, 5.6.5) Requires pandas and numpy. Nothing else. No GIS package, no shapefile, no projection library, no API key. Every input is a plain text table. Usage ----- python asi_toolkit.py --pull download and freeze the snapshots python asi_toolkit.py --report run the five diagnostics and print Data source ----------- American Community Survey 5-year estimates, table-based Summary File, https://www2.census.gov/programs-surveys/acs/summary_file/ Public, no registration, estimate and margin of error in adjacent columns. """ from __future__ import annotations import argparse import hashlib import io import json from pathlib import Path import urllib.request import numpy as np import pandas as pd HERE = Path(__file__).resolve().parent SNAP = HERE.parent / "snapshots" YEAR = 2023 STATE = "28" # Mississippi STATE_NAME = "Mississippi" ZIP_LO, ZIP_HI = 38601, 39776 # Mississippi ZCTA range SEED = 20260823 # Weights for the four index components. Equal weighting is a choice # Exercise II-A asks students to defend or change in Step 8. ASI_WEIGHTS = {"scale_spread": 1.0, "zoning_spread": 1.0, "transfer_gap": 1.0, "rank_instability": 1.0} BASE = (f"https://www2.census.gov/programs-surveys/acs/summary_file/{YEAR}" f"/table-based-SF/data/5YRData/acsdt5y{YEAR}-{{table}}.dat") SUMLEV = {"050": "county", "140": "tract", "150": "block group", "860": "ZCTA"} # Only the columns each table contributes. Everything else is discarded at # download time so the frozen snapshot stays small enough to ship. TABLES = { "b19013": ["B19013_E001", "B19013_M001"], "b25044": ["B25044_E001", "B25044_E003", "B25044_E010"], "c17002": ["C17002_E001", "C17002_E002", "C17002_E003"], "b01001": (["B01001_E001"] + [f"B01001_E{n:03d}" for n in range(20, 26)] + [f"B01001_E{n:03d}" for n in range(44, 50)]), } # -------------------------------------------------------------------------- # Stage 0. Freeze the data. # -------------------------------------------------------------------------- def _get(url: str, timeout: int) -> bytes: """Fetch a URL with the standard library, so no extra package is needed.""" with urllib.request.urlopen(url, timeout=timeout) as resp: return resp.read() def pull() -> None: SNAP.mkdir(parents=True, exist_ok=True) frames = [] for table, cols in TABLES.items(): url = BASE.format(table=table) print(f" downloading {table} ...", flush=True) raw = _get(url, timeout=900) df = pd.read_csv(io.BytesIO(raw), sep="|", dtype={"GEO_ID": str}, low_memory=False) keep = ["GEO_ID"] + [c for c in cols if c in df.columns] df = df[keep].copy() df["sumlev"] = df["GEO_ID"].str[:3] df = df[df["sumlev"].isin(SUMLEV)] df["code"] = df["GEO_ID"].str.split("US").str[-1] is_zcta = df["sumlev"] == "860" zc = pd.to_numeric(df["code"], errors="coerce") # A county FIPS and a ZCTA can share the same five digits, so the # state filter must never reach the ZCTA rows and the merge key must # carry the summary level alongside the code. keep_state = ~is_zcta & df["code"].str.startswith(STATE) keep_zcta = is_zcta & zc.between(ZIP_LO, ZIP_HI) df = df[keep_state | keep_zcta].drop(columns=["GEO_ID"]) frames.append(df) m = frames[0] for f in frames[1:]: m = m.merge(f, on=["sumlev", "code"], how="outer") out = SNAP / f"acs{YEAR}_{STATE_NAME.lower()}_asi.csv" m.to_csv(out, index=False) digest = hashlib.sha256(out.read_bytes()).hexdigest() (SNAP / "manifest.json").write_text(json.dumps({ "file": out.name, "rows": int(len(m)), "cols": int(m.shape[1]), "sha256": digest, "source": BASE.format(table=""), "vintage": f"ACS {YEAR} 5-year", "state": STATE_NAME, }, indent=2)) print(f"\nfrozen: {out}") print(f"rows {len(m)} cols {m.shape[1]} sha256 {digest}") def load() -> pd.DataFrame: path = SNAP / f"acs{YEAR}_{STATE_NAME.lower()}_asi.csv" if not path.exists(): raise SystemExit("No snapshot. Run: python asi_toolkit.py --pull") m = pd.read_csv(path, dtype={"code": str, "sumlev": str}) age = [c for c in m.columns if c.startswith("B01001_E") and c != "B01001_E001"] m["level"] = m["sumlev"].map(SUMLEV) m["pop"] = m["B01001_E001"] m["pop65"] = m[age].sum(axis=1, min_count=1) m["hh"] = m["B25044_E001"] m["hh_nv"] = m["B25044_E003"] + m["B25044_E010"] m["pov_univ"] = m["C17002_E001"] m["pov_below"] = m["C17002_E002"] + m["C17002_E003"] m["pct65"] = 100 * m["pop65"] / m["pop"] m["pct_nv"] = 100 * m["hh_nv"] / m["hh"] m["pct_pov"] = 100 * m["pov_below"] / m["pov_univ"] m["medinc"] = m["B19013_E001"].where(m["B19013_E001"] > 0) m["medinc_moe"] = m["B19013_M001"].where(m["B19013_M001"] > 0) return m # -------------------------------------------------------------------------- # Helpers. # -------------------------------------------------------------------------- def corr(a: pd.Series, b: pd.Series) -> tuple[float, int]: ok = a.notna() & b.notna() & np.isfinite(a) & np.isfinite(b) if ok.sum() < 3: return float("nan"), int(ok.sum()) return float(np.corrcoef(a[ok], b[ok])[0, 1]), int(ok.sum()) def need_index(d: pd.DataFrame, method: str = "rank") -> pd.Series: """Equal weight composite of three need variables.""" cols = ["pct65", "pct_pov", "pct_nv"] if method == "rank": parts = [d[c].rank(pct=True) for c in cols] elif method == "z": parts = [(d[c] - d[c].mean()) / d[c].std() for c in cols] elif method == "minmax": parts = [(d[c] - d[c].min()) / (d[c].max() - d[c].min()) for c in cols] else: raise ValueError(method) return sum(parts) def jaccard(a, b) -> float: a, b = set(a), set(b) return len(a & b) / len(a | b) if (a | b) else float("nan") # -------------------------------------------------------------------------- # The five diagnostics. # -------------------------------------------------------------------------- def multiscale(m, xvar="pct_pov", yvar="pct_nv"): rows = [] for lev in ("block group", "tract", "county"): d = m[m.level == lev] r, n = corr(d[xvar], d[yvar]) rows.append({"level": lev, "n": n, "r": r, "r2": r * r}) return pd.DataFrame(rows) def boundary_perturbation(m, n_iter=500, seed=SEED): bg = m[m.level == "block group"].dropna( subset=["pov_univ", "pov_below", "hh", "hh_nv"]).copy() n_zones = int((m.level == "tract").sum()) rng = np.random.default_rng(seed) cols = ["pov_univ", "pov_below", "hh", "hh_nv"] vals = [] for _ in range(n_iter): g = bg.groupby(rng.integers(0, n_zones, len(bg)))[cols].sum() r, _ = corr(100 * g.pov_below / g.pov_univ, 100 * g.hh_nv / g.hh) vals.append(r) return np.array(vals), n_zones def cross_zonal(m, xvar="pct_pov", yvar="pct_nv"): rows = [] for lev in ("tract", "county", "ZCTA"): d = m[m.level == lev] r, n = corr(d[xvar], d[yvar]) rows.append({"partition": lev, "n": n, "r": r}) return pd.DataFrame(rows) def rank_stability(m): bg = m[m.level == "block group"].dropna( subset=["pct65", "pct_pov", "pct_nv"]).copy() tr = m[m.level == "tract"].dropna( subset=["pct65", "pct_pov", "pct_nv"]).copy() bg["idx"] = need_index(bg) bg["parent_tract"] = bg["code"].str[:11] out = {"geography": jaccard(bg.nlargest(10, "idx")["parent_tract"], tr.assign(idx=need_index(tr)) .nlargest(10, "idx")["code"])} tops = {} for meth in ("rank", "z", "minmax"): tops[meth] = tr.assign(i=need_index(tr, meth)).nlargest(10, "i")["code"] out["rank_vs_z"] = jaccard(tops["rank"], tops["z"]) out["rank_vs_minmax"] = jaccard(tops["rank"], tops["minmax"]) out["z_vs_minmax"] = jaccard(tops["z"], tops["minmax"]) return out def reliability(m): rows = [] for lev in ("block group", "tract", "county"): d = m[m.level == lev] dd = d.dropna(subset=["medinc", "medinc_moe"]) cv = 100 * (dd.medinc_moe / 1.645) / dd.medinc rows.append({"level": lev, "n_estimates": len(dd), "n_suppressed": int(d.medinc.isna().sum()), "pct_cv_over_30": 100 * float((cv > 30).mean()), "median_cv": float(cv.median())}) return pd.DataFrame(rows) def asi(ms, pert, xz, rk): """Aggregation Sensitivity Index, higher meaning more sensitive. Four components, combined with the weights in ASI_WEIGHTS: scale spread range of |r| across resolutions, over the largest |r| zoning spread interquartile range of r across random partitions, over the absolute median r transfer gap absolute difference between the tract and ZCTA coefficients, over the largest absolute coefficient among the tract, county and ZCTA partitions rank instability 1 minus the mean Jaccard across one geography comparison and two standardization comparisons Three components fall between 0 and 1 by construction. The zoning spread can exceed 1 when the median correlation sits near zero. """ r = ms["r"].abs() scale_spread = float((r.max() - r.min()) / r.max()) zoning_spread = float(np.subtract(*np.percentile(pert, [75, 25])) / abs(np.median(pert))) tr = xz.set_index("partition")["r"].abs() transfer_gap = float(abs(tr["tract"] - tr["ZCTA"]) / tr.max()) rank_inst = 1 - float(np.mean([rk["geography"], rk["rank_vs_z"], rk["rank_vs_minmax"]])) comps = {"scale_spread": scale_spread, "zoning_spread": abs(zoning_spread), "transfer_gap": transfer_gap, "rank_instability": rank_inst} comps["ASI"] = float(np.average([comps[k] for k in ASI_WEIGHTS], weights=list(ASI_WEIGHTS.values()))) return comps def report(): m = load() man = json.loads((SNAP / "manifest.json").read_text()) print("=" * 74) print("AGGREGATION SENSITIVITY REPORT") print(f" source {man['vintage']}, {man['state']}") print(f" snapshot {man['file']} rows {man['rows']}") print(f" sha256 {man['sha256']}") print(f" seed {SEED}") print("=" * 74) print("\n[1] MULTISCALE VALIDATION (book 4.4, 5.6.1, 6.7.1)") print(" Does the result change when geographic resolution changes?") ms = multiscale(m) for _, x in ms.iterrows(): print(f" {x.level:<12s} n={x.n:>5d} r={x.r:+.4f} r2={x.r2:.4f}") print("\n[2] BOUNDARY PERTURBATION (book 4.4, 5.6.2, 6.7.2)") print(" Does the result change when boundaries move?") pert, nz = boundary_perturbation(m) obs = float(ms.loc[ms.level == "tract", "r"].iloc[0]) pct = float((pert < obs).mean() * 100) print(f" {len(pert)} random partitions into {nz} zones") print(f" r p5={np.percentile(pert,5):+.4f} " f"median={np.median(pert):+.4f} p95={np.percentile(pert,95):+.4f}") if pct >= 100.0: print(f" observed tract r = {obs:+.4f}, above all {len(pert)} " f"random zonings, the largest of which reaches " f"{pert.max():+.4f}") else: print(f" observed tract r = {obs:+.4f}, which is the " f"{pct:.1f}th percentile of random zonings") print("\n[3] CROSS-ZONAL TRANSFER (book 4.4, 5.6.3, 6.7.3)") print(" Does the relationship survive a different partition?") xz = cross_zonal(m) for _, x in xz.iterrows(): print(f" {x.partition:<12s} n={x.n:>5d} r={x.r:+.4f}") print("\n[4] RANK-ORDER STABILITY (book 4.4, ASI component)") print(" Do the same places stay in the top ten?") rk = rank_stability(m) print(f" block group vs tract Jaccard {rk['geography']:.3f}") print(f" percentile rank vs z score Jaccard {rk['rank_vs_z']:.3f}") print(f" percentile rank vs min-max Jaccard {rk['rank_vs_minmax']:.3f}") print(f" z score vs min-max Jaccard {rk['z_vs_minmax']:.3f}") print("\n[5] EVIDENCE RESOLUTION (book 4.5 recommendation 3)") print(" Is the evidence as fine as the display?") for _, x in reliability(m).iterrows(): print(f" {x.level:<12s} estimates={x.n_estimates:>5d} " f"suppressed={x.n_suppressed:>4d} " f"CV>30%={x.pct_cv_over_30:5.1f}% median CV={x.median_cv:5.1f}%") print("\n" + "-" * 74) a = asi(ms, pert, xz, rk) for k, v in a.items(): tag = " <== AGGREGATION SENSITIVITY INDEX" if k == "ASI" else "" print(f" {k:<18s} {v:.3f}{tag}") print("-" * 74) return a if __name__ == "__main__": ap = argparse.ArgumentParser(description=__doc__) ap.add_argument("--pull", action="store_true", help="freeze the snapshot") ap.add_argument("--report", action="store_true", help="run diagnostics") args = ap.parse_args() if args.pull: pull() if args.report or not args.pull: report()