diff --git a/.github/workflows/backend-tests.yml b/.github/workflows/backend-tests.yml index 49426f5e..4948cee7 100644 --- a/.github/workflows/backend-tests.yml +++ b/.github/workflows/backend-tests.yml @@ -31,8 +31,6 @@ jobs: run: uv sync --frozen --python 3.13 # No warehouse.duckdb or DuckLake exists here: tests that read real data # are marked `requires_warehouse` and skip (see tests/conftest.py). - # test_base.py still imports the deleted data_collection_OLD package; drop - # the --ignore once it is ported to data_collection.base. - name: Run the backend test suite working-directory: backend - run: uv run --frozen pytest -q --ignore=tests/test_base.py + run: uv run --frozen pytest -q diff --git a/backend/data_collection/acs5.py b/backend/data_collection/acs5.py index b147273e..55040095 100644 --- a/backend/data_collection/acs5.py +++ b/backend/data_collection/acs5.py @@ -2,7 +2,8 @@ Fetch ACS 5-Year Data Profile tables (DP02-DP05) for Vermont Geographies: counties + county subdivisions + Vermont statewide + United States Years: 2009 - Latest published data -Output: one wide table and tidy table +Output: one wide table and tidy table. Each tidy row carries its published +90% margin of error as MOE (the M code for an estimate, PM for a percent). Geography selection ------------------- diff --git a/backend/data_collection/base.py b/backend/data_collection/base.py index adaad3b8..cfab3632 100644 --- a/backend/data_collection/base.py +++ b/backend/data_collection/base.py @@ -1,5 +1,12 @@ """ Shared utilities for ACS 5-Year Census B-table scrapers. + +Census annotation values (negative sentinels such as -666666666) are missing +data, so a row with any missing component is null rather than a partial sum. + +Each tidy row carries its 90% margin of error: Value_MOE for Value and +Percent_MOE (in percentage points) for Percent, derived from the published +margins with the Census Bureau's approximations for sums and proportions. """ import logging @@ -8,6 +15,7 @@ from dataclasses import dataclass from datetime import datetime +import numpy as np import pandas as pd import requests @@ -23,6 +31,10 @@ STATE_FIPS = "50" STORAGE_LOCATION = "Data/Census/ACS_5" +# The one negative margin of error that is not missing: the estimate is +# controlled to an independent count, so it has no sampling error. +CONTROLLED_MOE = -555555555 + # --------------------------------------------------------------------------- # Geography registry # --------------------------------------------------------------------------- @@ -67,11 +79,24 @@ class VarGroup: # --------------------------------------------------------------------------- +def moe_code(code: str) -> str: + """The margin-of-error variable for an estimate: B19013_001E -> B19013_001M.""" + return code[:-1] + "M" + + def fetch( - year: int, variables: list[str], for_clause: str, in_clause: str + year: int, table: str, codes: list[str], for_clause: str, in_clause: str ) -> pd.DataFrame | None: + """ + Estimates and margins of error for `codes` from one table at one geography. + + group() returns every estimate and margin in the table, which keeps + within the API's 50-variable limit (B01001 alone needs 49 estimates). + Codes the year doesn't publish come back as all-NaN columns. + """ params = { - "get": ",".join(variables) + ",NAME", + # group() already includes NAME and GEO_ID. + "get": f"group({table})", "for": for_clause, "key": API_KEY, } @@ -83,19 +108,26 @@ def fetch( params["in"] = in_clause try: - r = requests.get(BASE_URL.format(year=year), params=params, timeout=30) + r = requests.get(BASE_URL.format(year=year), params=params, timeout=60) r.raise_for_status() data = r.json() df = pd.DataFrame(data[1:], columns=data[0]) - df["year"] = year - for c in df.columns: - if c[0] == "B": - df[c] = pd.to_numeric(df[c], errors="coerce") - return df + df = df.loc[:, ~df.columns.duplicated()] except Exception as e: # noqa: BLE001 -- skip and keep the run going - logger.warning("SKIP %s / %s: %s", year, for_clause, e) + logger.warning("SKIP %s / %s / %s: %s", year, table, for_clause, e) return None + variables = [v for code in codes for v in (code, moe_code(code))] + if missing := [v for v in variables if v not in df.columns]: + logger.warning("%s / %s: no %s", year, table, ", ".join(missing)) + # Geography columns (NAME, state, county, ...), without GEO_ID or the + # table's annotation columns. + ids = [c for c in df.columns if not c.startswith(f"{table}_") and c != "GEO_ID"] + df = df.reindex(columns=ids + variables) + df[variables] = df[variables].apply(pd.to_numeric, errors="coerce") + df["year"] = year + return df + # --------------------------------------------------------------------------- # Tidy computation @@ -115,6 +147,52 @@ def pct(val: float, total: float) -> float | None: return None +def estimates(df: pd.DataFrame, codes: list[str]) -> pd.DataFrame: + """Estimates for `codes`: NaN where absent or a Census annotation value + (https://www.census.gov/data/developers/data-sets/acs-1year/notes-on-acs-estimate-and-annotation-values.html). + No real estimate is negative.""" + est = df.reindex(columns=codes) + return est.where(est >= 0) + + +def margins(df: pd.DataFrame, codes: list[str]) -> pd.DataFrame: + """Margins of error for the estimates `codes`: NaN where absent or not + computable (e.g. -222222222, too few sample cases), 0 where controlled.""" + moe = df.reindex(columns=[moe_code(c) for c in codes]) + return moe.mask(moe == CONTROLLED_MOE, 0).where(lambda m: m >= 0) + + +def proportion_moe( + p: pd.Series, num_moe: pd.Series, den: pd.Series, den_moe: pd.Series +) -> pd.Series: + """ + Margin of error of the proportion p = num / den, where num is a subset of + den: sqrt(num_moe^2 - p^2 * den_moe^2) / den, or the ratio form (+) when + the term under the root is negative, per the ACS handbook. + + The ratio form also covers a zero term, which happens when num = den + (a 100% share): the proportion form would then claim no sampling error + at all, even for a town of 50 people. + """ + under = num_moe**2 - p**2 * den_moe**2 + under = under.where(under > 0, num_moe**2 + p**2 * den_moe**2) + return np.sqrt(under) / den + + +def sum_moe(est: pd.DataFrame, moe: pd.DataFrame) -> pd.Series: + """ + Margin of error of a sum of estimates: the root sum of squared margins, + counting only the largest margin among zero estimates. A zero estimate's + margin is modeled, and summing many of them overstates the error (Census + guidance, e.g. the 2017-2021 Variance Replicate Tables documentation). + """ + sq = moe.to_numpy() ** 2 + zero = est.to_numpy() == 0 + nonzero_sq = np.where(zero, 0, sq).sum(axis=1) + zero_sq = np.where(zero, sq, 0).max(axis=1, initial=0) + return pd.Series(np.sqrt(nonzero_sq + zero_sq), index=est.index) + + def compute_tidy_generic(df: pd.DataFrame, var_groups: list[VarGroup]) -> pd.DataFrame: base = pd.DataFrame( { @@ -126,21 +204,36 @@ def compute_tidy_generic(df: pd.DataFrame, var_groups: list[VarGroup]) -> pd.Dat } ) - def row_sum(codes: list[str]) -> pd.Series: - # reindex: codes absent from df become all-NaN columns, then count as 0 - return df.reindex(columns=codes).fillna(0).sum(axis=1) + def row_sum(codes: list[str]) -> tuple[pd.Series, pd.Series]: + """The sum of `codes` and its margin of error; each null when any of + its components is missing, since an incomplete sum is not an + estimate.""" + est, moe = estimates(df, codes), margins(df, codes) + value = est.sum(axis=1, min_count=len(codes)) + value_moe = sum_moe(est, moe).where(moe.notna().all(axis=1) & value.notna()) + return value, value_moe frames = [] for var in var_groups: - value = row_sum(var.codes) + value, value_moe = row_sum(var.codes) if var.denom: - denom = row_sum(var.denom) + denom, denom_moe = row_sum(var.denom) percent = (value / denom * 100).round(1).where(denom > 0) + p = (value / denom).where(denom > 0) + p_moe = proportion_moe(p, value_moe, denom, denom_moe) + percent_moe = (100 * p_moe).round(1) else: - percent = None + percent = percent_moe = None + if len(var.codes) > 1: + value_moe = value_moe.round() # a sum of counts frames.append( base.assign( - Section=var.section, Variable=var.label, Value=value, Percent=percent + Section=var.section, + Variable=var.label, + Value=value, + Value_MOE=value_moe, + Percent=percent, + Percent_MOE=percent_moe, ) ) if not frames: @@ -179,7 +272,7 @@ def fetch_geo(geo): merged = None for table_name, codes in fetch_specs.items(): logger.debug("%s / %s / %s...", year, table_name, geo_label) - df = fetch(year, codes, for_clause, in_clause) + df = fetch(year, table_name, codes, for_clause, in_clause) if df is None: return None if merged is None: diff --git a/backend/data_collection/census.py b/backend/data_collection/census.py index 9005b068..3bfef49c 100644 --- a/backend/data_collection/census.py +++ b/backend/data_collection/census.py @@ -110,6 +110,17 @@ def merge_census_cols(name_df, data_gdf, id_vars: list | None): value_name="Value", ) + # The variable list names only estimates (E, PE), so the margins of error + # (M, PM) would drop out below. Attach each one to its estimate's row as + # MOE instead: DP03_0062M -> DP03_0062E, DP03_0062PM -> DP03_0062PE. + is_moe = df_long["Code"].str.fullmatch(r"\w+_\d+P?M") + moe = ( + df_long[is_moe] + .assign(Code=lambda d: d["Code"].str[:-1] + "E") + .rename(columns={"Value": "MOE"}) + ) + df_long = df_long[~is_moe].merge(moe, on=[*id_vars, "Code"], how="left") + # Merge to get the right names and drop the cols (keep Variable_Code - the # original census variable code - for downstream lineage) return ( diff --git a/backend/tests/test_base.py b/backend/tests/test_base.py index 665147dc..98a439b4 100644 --- a/backend/tests/test_base.py +++ b/backend/tests/test_base.py @@ -1,11 +1,15 @@ """ Tests for data_collection/base.py: - pct(), VarGroup, compute_tidy_generic() + pct(), VarGroup, compute_tidy_generic() and its margins of error, fetch() """ +import math + import pandas as pd import pytest -from data_collection_OLD.base import VarGroup, compute_tidy_generic, pct + +from data_collection import base +from data_collection.base import VarGroup, compute_tidy_generic, pct # --------------------------------------------------------------------------- # pct() @@ -115,12 +119,52 @@ def test_multi_group_produces_one_row_each(self): assert len(result) == 2 assert set(result["Variable"]) == {"Metric A", "Metric B"} - def test_missing_code_treated_as_zero(self): - """A code absent from the row (NaN/None) should contribute 0.""" - df = _make_df(B01001_002E=[None]) + def test_missing_code_makes_value_null(self): + """An absent code makes the sum null, not a partial sum.""" + df = _make_df(B01001_002E=[500.0]) groups = [VarGroup("X", "S", ["B01001_002E", "B_MISSING"])] result = compute_tidy_generic(df, groups) - assert result.iloc[0]["Value"] == 0 + assert pd.isna(result.iloc[0]["Value"]) + + def test_sentinel_estimate_is_null(self): + """-666666666 (too few sample observations) is not a value.""" + df = _make_df(B19013_001E=[-666666666.0]) + groups = [VarGroup("Median HH Income", "Income", ["B19013_001E"])] + result = compute_tidy_generic(df, groups) + assert pd.isna(result.iloc[0]["Value"]) + + def test_one_sentinel_component_makes_sum_and_percent_null(self): + df = _make_df( + B01001_001E=[1000.0], B01001_003E=[40.0], B01001_027E=[-666666666.0] + ) + groups = [ + VarGroup("Under 5", "Age", ["B01001_003E", "B01001_027E"], ["B01001_001E"]) + ] + row = compute_tidy_generic(df, groups).iloc[0] + assert pd.isna(row["Value"]) + assert pd.isna(row["Percent"]) + + def test_two_sentinel_components_do_not_sum(self): + """Two sentinels once summed to -1,333,333,332.""" + df = _make_df(B25003_002E=[-666666666.0], B25003_003E=[-666666666.0]) + groups = [VarGroup("Occupied", "Housing", ["B25003_002E", "B25003_003E"])] + assert pd.isna(compute_tidy_generic(df, groups).iloc[0]["Value"]) + + def test_sentinel_denominator_makes_percent_null(self): + df = _make_df(B23025_001E=[-999999999.0], B23025_002E=[650.0]) + groups = [VarGroup("LFP Rate", "Labor", ["B23025_002E"], ["B23025_001E"])] + row = compute_tidy_generic(df, groups).iloc[0] + assert row["Value"] == 650.0 + assert pd.isna(row["Percent"]) + + def test_real_zero_is_kept(self): + df = _make_df(B01001_001E=[1000.0], B01001_003E=[0.0], B01001_027E=[0.0]) + groups = [ + VarGroup("Under 5", "Age", ["B01001_003E", "B01001_027E"], ["B01001_001E"]) + ] + row = compute_tidy_generic(df, groups).iloc[0] + assert row["Value"] == 0 + assert row["Percent"] == 0.0 def test_passthrough_base_columns(self): df = _make_df(B01001_001E=[100.0]) @@ -151,3 +195,169 @@ def test_multiple_geographies(self): result = compute_tidy_generic(df, groups) assert len(result) == 2 assert set(result["geo_type"]) == {"county", "county_subdivision"} + + +# --------------------------------------------------------------------------- +# Margins of error +# --------------------------------------------------------------------------- + + +def _rate(num, num_moe, den, den_moe): + """One proportion row: num of den, each with its margin of error.""" + df = _make_df( + B23025_001E=[den], + B23025_001M=[den_moe], + B23025_002E=[num], + B23025_002M=[num_moe], + ) + groups = [VarGroup("LFP Rate", "Labor", ["B23025_002E"], ["B23025_001E"])] + return compute_tidy_generic(df, groups).iloc[0] + + +class TestMarginsOfError: + def test_single_code_passes_its_margin_through(self): + df = _make_df(B19013_001E=[75000.0], B19013_001M=[4120.0]) + groups = [VarGroup("Median HH Income", "Income", ["B19013_001E"])] + row = compute_tidy_generic(df, groups).iloc[0] + assert row["Value_MOE"] == 4120.0 + assert row["Percent_MOE"] is None + + def test_sum_margin_is_root_sum_of_squares(self): + df = _make_df( + B01001_002E=[500.0], + B01001_002M=[30.0], + B01001_026E=[520.0], + B01001_026M=[40.0], + ) + groups = [VarGroup("Total Pop", "Pop", ["B01001_002E", "B01001_026E"])] + assert compute_tidy_generic(df, groups).iloc[0]["Value_MOE"] == 50.0 + + def test_proportion_margin_uses_census_formula(self): + row = _rate(150.0, 50.0, 1000.0, 100.0) + # sqrt(MOE_num^2 - p^2 * MOE_den^2) / den, in percentage points + expected = 100 * math.sqrt(50**2 - 0.15**2 * 100**2) / 1000 + assert row["Percent"] == 15.0 + assert row["Percent_MOE"] == round(expected, 1) + + def test_negative_radicand_falls_back_to_ratio_formula(self): + row = _rate(900.0, 10.0, 1000.0, 100.0) + # 10^2 - 0.9^2 * 100^2 < 0, so the + form applies + expected = 100 * math.sqrt(10**2 + 0.9**2 * 100**2) / 1000 + assert row["Percent_MOE"] == round(expected, 1) + + def test_sum_counts_only_the_largest_zero_estimate_margin(self): + """Census guidance: zero estimates' margins are modeled, so a sum uses + only the largest of them instead of all of them.""" + df = _make_df( + B01001_003E=[0.0], + B01001_003M=[12.0], + B01001_004E=[0.0], + B01001_004M=[15.0], + B01001_005E=[10.0], + B01001_005M=[5.0], + ) + codes = ["B01001_003E", "B01001_004E", "B01001_005E"] + groups = [VarGroup("Under 15", "Age", codes)] + row = compute_tidy_generic(df, groups).iloc[0] + assert row["Value_MOE"] == round(math.sqrt(15**2 + 5**2)) + + def test_full_share_margin_is_not_zero(self): + """num = den makes the proportion form zero; the ratio form applies.""" + row = _rate(40.0, 12.0, 40.0, 12.0) + assert row["Percent"] == 100.0 + assert row["Percent_MOE"] == round(100 * math.sqrt(12**2 + 12**2) / 40, 1) + + def test_controlled_margin_is_zero(self): + """-555555555: the estimate is controlled and has no sampling error.""" + row = _rate(150.0, 50.0, 1000.0, -555555555.0) + assert row["Percent_MOE"] == round(100 * 50 / 1000, 1) + + def test_uncomputable_margin_is_null_but_estimate_stands(self): + """-222222222: too few sample cases to compute a margin.""" + row = _rate(0.0, -222222222.0, 1000.0, 100.0) + assert row["Value"] == 0 + assert row["Percent"] == 0.0 + assert pd.isna(row["Value_MOE"]) + assert pd.isna(row["Percent_MOE"]) + + def test_missing_margin_column_is_null(self): + df = _make_df(B19013_001E=[75000.0]) + groups = [VarGroup("Median HH Income", "Income", ["B19013_001E"])] + row = compute_tidy_generic(df, groups).iloc[0] + assert row["Value"] == 75000.0 + assert pd.isna(row["Value_MOE"]) + + def test_suppressed_estimate_nulls_value_and_margin(self): + df = _make_df(B19013_001E=[-666666666.0], B19013_001M=[-222222222.0]) + groups = [VarGroup("Median HH Income", "Income", ["B19013_001E"])] + row = compute_tidy_generic(df, groups).iloc[0] + assert pd.isna(row["Value"]) + assert pd.isna(row["Value_MOE"]) + + +# --------------------------------------------------------------------------- +# fetch() +# --------------------------------------------------------------------------- + + +class _Response: + def __init__(self, payload): + self._payload = payload + + def raise_for_status(self): + pass + + def json(self): + return self._payload + + +def test_fetch_keeps_estimates_margins_and_geography(monkeypatch): + # Like a real group() response: annotation columns, GEO_ID, then geography. + header = [ + "B19013_001E", + "B19013_001EA", + "B19013_001M", + "B19013_001MA", + "GEO_ID", + "NAME", + "state", + "county", + ] + rows = [ + ["68854", None, "4120", None, "0500000US50007", + "Chittenden County, Vermont", "50", "007"], + ["-666666666", "-", "-222222222", "**", "0500000US50009", + "Essex County, Vermont", "50", "009"], + ] # fmt: skip + calls = [] + + def fake_get(url, params, timeout): + calls.append(params) + return _Response([header, *rows]) + + monkeypatch.setattr(base, "API_KEY", "test-key") + monkeypatch.setattr(base.requests, "get", fake_get) + + df = base.fetch(2023, "B19013", ["B19013_001E"], "county:*", "state:50") + + assert calls[0]["get"] == "group(B19013)" + assert list(df.columns) == [ + "NAME", "state", "county", "B19013_001E", "B19013_001M", "year" + ] # fmt: skip + assert df["B19013_001E"].tolist() == [68854, -666666666] + assert df["B19013_001M"].tolist() == [4120, -222222222] + assert (df["year"] == 2023).all() + + +def test_fetch_fills_unpublished_codes_with_nan(monkeypatch): + header = ["B19013_001E", "B19013_001M", "GEO_ID", "NAME", "state"] + rows = [["70000", "500", "0400000US50", "Vermont", "50"]] + monkeypatch.setattr(base, "API_KEY", "test-key") + monkeypatch.setattr( + base.requests, "get", lambda *a, **k: _Response([header, *rows]) + ) + + df = base.fetch(2023, "B19013", ["B19013_001E", "B19013_002E"], "state:50", "") + + assert pd.isna(df.loc[0, "B19013_002E"]) + assert pd.isna(df.loc[0, "B19013_002M"]) diff --git a/backend/tests/test_census_tidy.py b/backend/tests/test_census_tidy.py new file mode 100644 index 00000000..4409d160 --- /dev/null +++ b/backend/tests/test_census_tidy.py @@ -0,0 +1,59 @@ +""" +Tests for data_collection/census.py: the DP profile tidy step keeps each +estimate's published margin of error as an MOE column. +""" + +import pandas as pd + +from data_collection.acs5 import ID_VARS +from data_collection.census import merge_census_cols, relabel_census_cols + +LABEL = "INCOME AND BENEFITS!!Total households!!Median household income (dollars)" + + +def _tidy(**codes) -> pd.DataFrame: + # The profile variables.json lists only E and PE codes. + labels = pd.DataFrame( + { + "Name": ["DP03_0062E", "DP03_0062PE"], + "Label": [f"Estimate!!{LABEL}", f"Percent!!{LABEL}"], + } + ) + # Like a group() response: annotation columns, GEO_ID and geography. + wide = pd.DataFrame( + { + "year": [2023], + "geo_type": ["county"], + "table": ["DP03"], + "NAME": ["Addison County, Vermont"], + "state": ["50"], + "county": ["001"], + "GEO_ID": ["0500000US50001"], + **{code: [value] for code, value in codes.items()}, + } + ) + tidy = merge_census_cols(relabel_census_cols(labels), wide, ID_VARS) + return tidy.set_index("Variable_Code") + + +def test_estimates_and_percents_carry_their_margin(): + tidy = _tidy( + DP03_0062E="85000", + DP03_0062EA=None, + DP03_0062M="3100", + DP03_0062MA=None, + DP03_0062PE="-888888888", + DP03_0062PM="-888888888", + ) + + assert sorted(tidy.index) == ["DP03_0062E", "DP03_0062PE"] + assert tidy.loc["DP03_0062E", ["Value", "MOE"]].tolist() == ["85000", "3100"] + assert tidy.loc["DP03_0062PE", "MOE"] == "-888888888" + assert tidy.loc["DP03_0062E", "Measure"] == "Estimate" + + +def test_missing_margin_is_null(): + tidy = _tidy(DP03_0062E="85000", DP03_0062PE="-888888888") + + assert tidy.loc["DP03_0062E", "Value"] == "85000" + assert pd.isna(tidy.loc["DP03_0062E", "MOE"])