"""exercise_4a.py Exercise IV-A. The Uncertainty Cascade Assessment. Propagates the published margins of error of three survey numerators, with the denominators and tract populations held at their published estimates, through an index, a threshold and a population count, and reports what arrives at the far end beside the single number the point estimate produced. python exercise_4a.py """ from __future__ import annotations 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 sys.path.insert(0, str(Path(__file__).resolve().parent)) import common as C # noqa: E402 from mapping import load_geometry # noqa: E402 from matplotlib.collections import PolyCollection # noqa: E402 GEO = C.SNAP / "ms_tracts.geojson" N_SIM = 1000 CUTOFF = 0.80 # tracts above this percentile qualify MOE_TO_SD = 1.645 # the Census publishes a 90 percent margin WEIGHTINGS = {"equal": (1, 1, 1), "poverty heavy": (1, 2, 1), "age heavy": (2, 1, 1)} def load(): base = C.load_acs() tr = base[base.level == "tract"].copy() moe = pd.read_csv(C.SNAP / "acs2023_mississippi_moe.csv", dtype={"code": str, "sumlev": str}) moe = moe[moe.sumlev == "140"] age_e = [f"B01001_E{n:03d}" for n in list(range(20, 26)) + list(range(44, 50))] age_m = [f"B01001_M{n:03d}" for n in list(range(20, 26)) + list(range(44, 50))] # Margins on a derived sum combine in quadrature, which is the Census # Bureau's own published rule for aggregating counts. moe["pop65_moe"] = np.sqrt((moe[age_m] ** 2).sum(axis=1)) moe["nv_moe"] = np.sqrt(moe["B25044_M003"] ** 2 + moe["B25044_M010"] ** 2) moe["pov_moe"] = np.sqrt(moe["C17002_M002"] ** 2 + moe["C17002_M003"] ** 2) keep = ["code", "pop65_moe", "nv_moe", "pov_moe", "B01001_M001", "B25044_M001", "C17002_M001"] d = tr.merge(moe[keep], on="code", how="inner") need = ["pop", "pop65", "hh", "hh_nv", "pov_univ", "pov_below"] d = d.dropna(subset=need + ["pop65_moe", "nv_moe", "pov_moe"]) # Water and zero-population tracts, whose GEOID ends 990000, carry a zero # denominator. Dividing by it yields NaN, one NaN poisons the percentile # rank, and every downstream comparison then returns False. Keep only # tracts whose three shares are finite. d = d[(d["pop"] > 0) & (d["hh"] > 0) & (d["pov_univ"] > 0)] return d.reset_index(drop=True) def index_from(p65, ppov, pnv, w=(1, 1, 1)): r = [pd.Series(v).rank(pct=True).to_numpy() for v in (p65, ppov, pnv)] return w[0] * r[0] + w[1] * r[1] + w[2] * r[2] def main(): C.ensure_dirs() d = load() C.banner("exercise IV-A: the uncertainty cascade assessment", "ACS 2023 5-year, Mississippi, estimates and margins of error", "acs2023_mississippi_asi.csv", {"tracts": len(d), "draws": N_SIM}) print(f" moe file acs2023_mississippi_moe.csv sha256 " f"{C.sha256(C.SNAP / 'acs2023_mississippi_moe.csv')[:32]}") out = {"tracts": int(len(d)), "draws": N_SIM, "cutoff_pct": CUTOFF * 100} # ---- 1. the point estimate ------------------------------------------ p65 = 100 * d.pop65 / d["pop"] ppov = 100 * d.pov_below / d.pov_univ pnv = 100 * d.hh_nv / d.hh idx0 = index_from(p65, ppov, pnv) bar0 = np.quantile(idx0, CUTOFF) qual0 = idx0 >= bar0 pop0 = int(d["pop"][qual0].sum()) first0 = int(np.argmax(idx0)) print("\n[1] THE POINT ESTIMATE") print(f" qualifying tracts {int(qual0.sum())}") print(f" population in them {pop0:,}") print(f" highest ranked tract {d.code.iloc[first0]}") out["point"] = {"qualifying": int(qual0.sum()), "population": pop0, "top_tract": str(d.code.iloc[first0])} # ---- 2. the cascade -------------------------------------------------- rng = np.random.default_rng(C.SEED) n = len(d) tot = np.empty(N_SIM) first_hits = 0 ever = np.zeros(n, bool) for s in range(N_SIM): pop65 = np.maximum(rng.normal(d.pop65, d.pop65_moe / MOE_TO_SD), 0) below = np.maximum(rng.normal(d.pov_below, d.pov_moe / MOE_TO_SD), 0) nov = np.maximum(rng.normal(d.hh_nv, d.nv_moe / MOE_TO_SD), 0) i = index_from(100 * pop65 / d["pop"], 100 * below / d.pov_univ, 100 * nov / d.hh) q = i >= np.quantile(i, CUTOFF) tot[s] = d["pop"][q].sum() ever |= q if int(np.argmax(i)) == first0: first_hits += 1 p5, p50, p95 = np.percentile(tot, [5, 50, 95]) print("\n[2] ONE THOUSAND REALIZATIONS") print(f" 5th percentile {p5:>12,.0f}") print(f" median {p50:>12,.0f}") print(f" 95th percentile {p95:>12,.0f}") print(f" interval width {p95-p5:>12,.0f} " f"({100*(p95-p5)/pop0:.1f}% of the point estimate)") print(f" top rank retained {100*first_hits/N_SIM:.1f}% of draws") print(f" tracts qualifying at least once {int(ever.sum())} " f"against {int(qual0.sum())} under the point estimate") out["cascade"] = { "p5": int(p5), "p50": int(p50), "p95": int(p95), "width": int(p95 - p5), "width_pct_of_point": round(100 * (p95 - p5) / pop0, 1), "top_rank_retention_pct": round(100 * first_hits / N_SIM, 1), "ever_qualifying": int(ever.sum())} # ---- 3. reproducibility --------------------------------------------- def one(seed): r = np.random.default_rng(seed) pop65 = np.maximum(r.normal(d.pop65, d.pop65_moe / MOE_TO_SD), 0) below = np.maximum(r.normal(d.pov_below, d.pov_moe / MOE_TO_SD), 0) nov = np.maximum(r.normal(d.hh_nv, d.nv_moe / MOE_TO_SD), 0) i = index_from(100 * pop65 / d["pop"], 100 * below / d.pov_univ, 100 * nov / d.hh) return float(d["pop"][i >= np.quantile(i, CUTOFF)].sum()) a, b, c = one(C.SEED), one(C.SEED), one(1) print("\n[3] SEED DISCIPLINE") print(f" same seed twice {a:,.0f} and {b:,.0f} " f"{'identical' if a == b else 'DIFFERENT'}") print(f" a different seed {c:,.0f}") out["seed_check"] = {"same": a == b, "value": int(a), "other_seed": int(c)} # ---- 4. the ensemble -------------------------------------------------- print("\n[4] ENSEMBLE OVER WEIGHTINGS") sets = {} for name, w in WEIGHTINGS.items(): i = index_from(p65, ppov, pnv, w) q = set(d.code[i >= np.quantile(i, CUTOFF)]) sets[name] = q print(f" {name:<14} qualifying {len(q):>4d}") inter = set.intersection(*sets.values()) union = set.union(*sets.values()) print(f" qualify under every weighting {len(inter)}") print(f" qualify under at least one {len(union)}") out["ensemble"] = {"intersection": len(inter), "union": len(union), "by_weighting": {k: len(v) for k, v in sets.items()}} # ---- figures ---------------------------------------------------------- fig, ax = plt.subplots(figsize=(8.5, 5.5)) ax.hist(tot, bins=40, color="#c6dbef", edgecolor="#4292c6") ax.axvline(pop0, color="#b2182b", lw=2, label=f"point estimate {pop0:,}") ax.axvspan(p5, p95, color="#fdd0a2", alpha=0.35, label=f"5th to 95th, {p5:,.0f} to {p95:,.0f}") ax.set_xlabel("population in qualifying tracts") ax.set_ylabel("realizations") ax.set_title("One number, and the interval the survey actually supports", loc="left", fontweight="bold") ax.legend(frameon=False, fontsize=9) ax.text(0, -0.13, "Source: U.S. Census Bureau, ACS 2023 5-year summary " f"file, estimates and margins of error. {len(d)} " f"tracts, {N_SIM} draws, seed {C.SEED}.", transform=ax.transAxes, fontsize=7, color="#666") fig.savefig(C.FIG / "exercise_4a_cascade.png", dpi=150, bbox_inches="tight") plt.close(fig) tbl = d[["code", "pop"]].copy() tbl["pct65"], tbl["pct_pov"], tbl["pct_nv"] = p65, ppov, pnv tbl["index_point"] = idx0 tbl["qualifies_point"] = qual0 tbl["qualifies_in_any_draw"] = ever C.write_handoff("exercise_4a", tbl, { "draws": N_SIM, "seed": C.SEED, "cutoff": f"{CUTOFF:.0%} percentile, recomputed per draw", "error_model": "normal, sd = published 90 percent margin / 1.645, " "truncated at zero", "correlated_errors": "ignored; inputs drawn independently", "denominators": "treated as exact, their own margins not propagated", "weighting": "equal across three inputs", }, [ "Errors were drawn independently. ACS margins on variables from the " "same survey are correlated, so the true interval is wider or " "narrower than reported here. Which, and can you bound it?", "Denominators were held exact while numerators varied. Is that " "defensible for tracts with small populations?", "The top ranked tract holds its rank in only part of the draws. What " "would you write beneath a published ranking?", "How would you state the final figure to someone who asked for one " "number and will not read a footnote?", ]) C.save_checkvalues("exercise_4a", out) print(f"\nwrote {C.FIG / 'exercise_4a_cascade.png'}") if __name__ == "__main__": main()