From 594c044f9fb927bde30c106b3afcd397bf3a70db Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Thu, 15 Oct 2026 08:02:54 -0400 Subject: [PATCH 1/5] Measure each build's HBAI distribution and single ages against its base Calibration fits grouped targets; these checks measure what the groups hide, on the built dataset. Each PR build and its base build report person-weighted equivalised income medians, means and Gini, and relative low income for everyone, children, working-age adults and pensioners, against DWP's HBAI FYE 2025 estimates and 95% intervals (tables 8a and 8b, read from the published ODS files), at period 2024. They also report people by sex and single year of age against the ONS 2024-based projection at the calibration year. The comparison job adds them to the PR comment and flags a measure that leaves its interval; it fails only on a broken or mismatched report until the same-code spread is measured (#598). On uk-data 1.58.1, 5 of 14 HBAI measures are outside their intervals (median AHC, both Gini coefficients, child low income BHC and AHC), and single ages 0 to 79 are off by 19.5% on average. Co-Authored-By: Claude Opus 5.5 --- .github/workflows/pull_request.yaml | 38 ++ changelog.d/distribution-checks.added.md | 1 + .../targets/reference_data/hbai_fye_2025.yaml | 27 ++ .../tests/test_distribution_validation.py | 162 ++++++++ policyengine_uk_data/tests/test_fit_report.py | 6 +- .../utils/distribution_validation.py | 377 ++++++++++++++++++ 6 files changed, 609 insertions(+), 2 deletions(-) create mode 100644 changelog.d/distribution-checks.added.md create mode 100644 policyengine_uk_data/targets/reference_data/hbai_fye_2025.yaml create mode 100644 policyengine_uk_data/tests/test_distribution_validation.py create mode 100644 policyengine_uk_data/utils/distribution_validation.py diff --git a/.github/workflows/pull_request.yaml b/.github/workflows/pull_request.yaml index 4b035dc3..c1a63b81 100644 --- a/.github/workflows/pull_request.yaml +++ b/.github/workflows/pull_request.yaml @@ -86,6 +86,13 @@ jobs: run: uv run --frozen make data env: TESTING: "1" + - name: Measure the income distribution and single ages + run: uv run --frozen python -m policyengine_uk_data.utils.distribution_validation measure policyengine_uk_data/storage/enhanced_frs_2024_25.h5 --out distribution_report.json + - name: Save the distribution report + uses: actions/upload-artifact@v7 + with: + name: distribution_report.json + path: distribution_report.json - name: Save calibration log (constituencies) uses: actions/upload-artifact@v7 with: @@ -150,6 +157,19 @@ jobs: with: name: base_constituency_calibration_log.csv path: base/constituency_calibration_log.csv + - name: Measure the base's income distribution and single ages + working-directory: base + # The head's measurement code, on the base's build. + run: | + mkdir -p policyengine_uk_data/targets/reference_data + cp ../head/policyengine_uk_data/utils/distribution_validation.py policyengine_uk_data/utils/ + cp ../head/policyengine_uk_data/targets/reference_data/hbai_fye_2025.yaml policyengine_uk_data/targets/reference_data/ + uv run --frozen python -m policyengine_uk_data.utils.distribution_validation measure policyengine_uk_data/storage/enhanced_frs_2024_25.h5 --out ../base_distribution_report.json + - name: Save the base distribution report + uses: actions/upload-artifact@v7 + with: + name: base_distribution_report.json + path: base_distribution_report.json fit-comparison: name: Calibration fit against the base needs: [test, base-build] @@ -183,6 +203,13 @@ jobs: with: name: base_constituency_calibration_log.csv path: logs/base + - name: Download the distribution reports + uses: actions/download-artifact@v8 + continue-on-error: true + with: + pattern: "*distribution_report.json" + path: reports + merge-multiple: true - name: Compare the final fits id: compare env: @@ -206,6 +233,17 @@ jobs: set -e { echo "$marker"; cat body.md 2>/dev/null || echo '## Calibration fit against the base build: FAIL'; echo; echo "[Workflow run]($RUN_URL)"; } > fit.md fi + if [ -f reports/base_distribution_report.json ] && [ -f reports/distribution_report.json ]; then + set +e + uv run --frozen python -m policyengine_uk_data.utils.distribution_validation compare \ + reports/base_distribution_report.json reports/distribution_report.json --markdown distribution.md + distribution_status=$? + set -e + { echo; cat distribution.md; } >> fit.md + if [ "$distribution_status" != 0 ]; then echo "status=1" >> "$GITHUB_OUTPUT"; fi + else + { echo; echo "The income distribution and single-age reports are missing for one build."; } >> fit.md + fi cat fit.md >> "$GITHUB_STEP_SUMMARY" - name: Post or update the pull request comment if: ${{ !cancelled() && github.event.pull_request.head.repo.full_name == github.repository }} diff --git a/changelog.d/distribution-checks.added.md b/changelog.d/distribution-checks.added.md new file mode 100644 index 00000000..5077bf91 --- /dev/null +++ b/changelog.d/distribution-checks.added.md @@ -0,0 +1 @@ +Measure each pull request build's HBAI income distribution (medians, means, Gini and relative low income by group against DWP's HBAI FYE 2025 estimates and intervals) and single-age population (against the ONS 2024-based projection), and report them against the base build in the pull request comment. They flag but do not fail yet. diff --git a/policyengine_uk_data/targets/reference_data/hbai_fye_2025.yaml b/policyengine_uk_data/targets/reference_data/hbai_fye_2025.yaml new file mode 100644 index 00000000..c209ab76 --- /dev/null +++ b/policyengine_uk_data/targets/reference_data/hbai_fye_2025.yaml @@ -0,0 +1,27 @@ +# Households below average income, FYE 2025 (DWP, 26 March 2026), from the +# published data tables: +# https://assets.publishing.service.gov.uk/media/69c2902a7e02b81c0d1c7762/hbai-data_tables-2025.zip +# uncertainty-hbai-2024-25-tables.ods, tables 8a (medians, means, Gini) and +# 8b (low income), read on 15 October 2026. Estimates and 95% confidence +# intervals for the United Kingdom, FRS 2024/25, which policyengine-uk models +# as period 2024. Incomes are equivalised, in £ a week at 2024/25 prices; +# rates and Gini coefficients are per cent. +schema_version: 1 +source: DWP Households below average income FYE 2025 +period: 2024 +measures: + median_bhc: {estimate: 719.481062, low: 703.567322, high: 735.575892, table: 8a} + median_ahc: {estimate: 623.144059, low: 605.209749, high: 637.071794, table: 8a} + mean_bhc: {estimate: 854.979494, low: 836.367495, high: 870.351587, table: 8a} + mean_ahc: {estimate: 749.82155, low: 732.329379, high: 764.496483, table: 8a} + gini_bhc: {estimate: 32.736632, low: 31.889881, high: 33.978308, table: 8a} + gini_ahc: {estimate: 36.780956, low: 35.897706, high: 38.293054, table: 8a} + # Percentage of individuals below 60% of the relative median. + relative_low_income_bhc_all: {estimate: 15.911572, low: 14.446526, high: 17.637657, table: 8b} + relative_low_income_bhc_children: {estimate: 20.519181, low: 17.356404, high: 23.955585, table: 8b} + relative_low_income_bhc_working_age: {estimate: 14.281508, low: 12.908347, high: 16.018814, table: 8b} + relative_low_income_bhc_pensioners: {estimate: 15.895071, low: 13.690455, high: 17.991754, table: 8b} + relative_low_income_ahc_all: {estimate: 19.640844, low: 17.9939, high: 21.271606, table: 8b} + relative_low_income_ahc_children: {estimate: 27.41076, low: 23.49369, high: 30.644085, table: 8b} + relative_low_income_ahc_working_age: {estimate: 18.579734, low: 17.012601, high: 20.343162, table: 8b} + relative_low_income_ahc_pensioners: {estimate: 13.866371, low: 11.618244, high: 15.898207, table: 8b} diff --git a/policyengine_uk_data/tests/test_distribution_validation.py b/policyengine_uk_data/tests/test_distribution_validation.py new file mode 100644 index 00000000..3c86db33 --- /dev/null +++ b/policyengine_uk_data/tests/test_distribution_validation.py @@ -0,0 +1,162 @@ +"""Single-age and HBAI distribution checks (utils/distribution_validation.py).""" + +import numpy as np +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from policyengine_uk_data.utils.distribution_validation import ( + SCHEMA, + TOP_AGE_BAND, + compare, + hbai_measures, + load_hbai_reference, + render, + single_ages, + summarise_single_ages, + weighted_gini, + weighted_median, +) + +samples = st.lists( + st.tuples(st.floats(-1e4, 1e5, allow_nan=False), st.integers(1, 5)), + min_size=1, + max_size=40, +) + + +def _expand(pairs): + return np.repeat([v for v, _ in pairs], [w for _, w in pairs]).astype(float) + + +@settings(max_examples=300, deadline=None) +@given(pairs=samples) +def test_weighted_median_equals_the_median_of_the_expanded_sample(pairs): + """Differential: integer weights are repeated records.""" + values = np.array([v for v, _ in pairs]) + weights = np.array([w for _, w in pairs], float) + expanded = np.sort(_expand(pairs)) + # The lower median of the expanded sample. + lower = expanded[(len(expanded) - 1) // 2] + upper = expanded[len(expanded) // 2] + assert lower <= weighted_median(values, weights) <= upper + # Scaling every weight leaves it unchanged. + assert weighted_median(values, weights * 7.5) == weighted_median(values, weights) + + +@settings(max_examples=300, deadline=None) +@given( + pairs=st.lists( + st.tuples(st.floats(0.01, 1e5), st.integers(1, 5)), min_size=2, max_size=40 + ) +) +def test_weighted_gini_equals_the_mean_difference_gini_of_the_expanded_sample(pairs): + expanded = _expand(pairs) + n = len(expanded) + mean_difference = np.abs(expanded[:, None] - expanded[None, :]).sum() / (n * n) + reference = mean_difference / (2 * expanded.mean()) + values = np.array([v for v, _ in pairs]) + weights = np.array([w for _, w in pairs], float) + assert weighted_gini(values, weights) == pytest.approx(reference, abs=1e-9) + + +def test_gini_of_equal_incomes_is_zero(): + assert weighted_gini([5.0, 5.0, 5.0], [1, 2, 3]) == pytest.approx(0) + + +def test_hbai_measures_on_a_known_population(): + income = np.array([100.0, 200.0, 300.0, 400.0, 1000.0]) + weights = np.ones(5) + child = np.array([1, 0, 0, 0, 0], bool) + pensioner = np.array([0, 0, 0, 0, 1], bool) + m = hbai_measures(income, income - 50, weights, child, pensioner) + assert m["median_bhc"] == 300.0 + assert m["mean_bhc"] == 400.0 + # Below 60% of 300: 100 only, a child. + assert m["relative_low_income_bhc_all"] == pytest.approx(20.0) + assert m["relative_low_income_bhc_children"] == pytest.approx(100.0) + assert m["relative_low_income_bhc_working_age"] == pytest.approx(0.0) + assert m["relative_low_income_bhc_pensioners"] == pytest.approx(0.0) + + +def test_the_hbai_reference_is_the_published_table(): + """Spot values from DWP HBAI FYE 2025, tables 8a and 8b.""" + ref = load_hbai_reference() + assert ref["period"] == 2024 + m = ref["measures"] + assert m["median_bhc"]["estimate"] == pytest.approx(719.481062) + assert m["median_ahc"]["estimate"] == pytest.approx(623.144059) + assert m["relative_low_income_ahc_children"]["estimate"] == pytest.approx(27.41076) + assert len(m) == 14 + + +def test_single_ages_band_the_oldest_and_split_by_sex(): + age = np.array([0, 0, 79, 80, 85, 90]) + male = np.array([1, 0, 1, 1, 0, 1], bool) + weights = np.array([2.0, 3.0, 1.0, 1.0, 4.0, 5.0]) + reference = { + (s, a): 1.0 for s in ("male", "female") for a in range(TOP_AGE_BAND + 1) + } + rows = single_ages(age, male, weights, reference) + assert rows["male"]["0"]["value"] == 2.0 and rows["female"]["0"]["value"] == 3.0 + assert rows["male"]["80"]["value"] == 6.0 and rows["female"]["80"]["value"] == 4.0 + summary = summarise_single_ages(rows) + assert summary["error_80_plus"] == {"male": 5.0, "female": 3.0} + + +def _report(**hbai_values): + ref = load_hbai_reference()["measures"] + hbai = {} + for name, r in ref.items(): + value = hbai_values.get(name, r["estimate"]) + hbai[name] = { + "value": value, + "reference": r["estimate"], + "low": r["low"], + "high": r["high"], + "inside": r["low"] <= value <= r["high"], + } + single = { + "mean_abs_error_0_79": 0.2, + "max_abs_error_0_79": 0.9, + "max_abs_error_15_19": 0.6, + "error_80_plus": {"male": 0.0, "female": 0.0}, + } + return { + "schema": SCHEMA, + "periods": {"hbai": 2024, "single_age": 2025}, + "hbai": hbai, + "summary": {"single_age": single}, + } + + +def test_compare_flags_a_measure_that_leaves_its_interval_but_does_not_fail(): + result = compare(_report(), _report(median_bhc=760.0)) + assert result["status"] == "pass" + assert result["flags"] == [ + "HBAI median_bhc leaves its 95% interval: 719.48 -> 760.00 (HBAI 719.48, 703.57 to 735.58)." + ] + text = render(result) + assert "760.00 (outside)" in text and "do not fail yet" in text + + +def test_compare_with_itself_passes_with_no_flags(): + result = compare(_report(), _report()) + assert result["status"] == "pass" and result["flags"] == [] + assert all(r["further_from_reference"] == 0 for r in result["rows"]) + + +@pytest.mark.parametrize( + "change, message", + [ + (lambda r: r.update(schema="other"), "malformed"), + (lambda r: r["periods"].update(hbai=2025), "different periods"), + (lambda r: r["hbai"]["median_bhc"].update(value=float("nan")), "not finite"), + ], +) +def test_compare_fails_on_a_broken_or_mismatched_report(change, message): + head = _report() + change(head) + result = compare(_report(), head) + assert result["status"] == "fail" + assert any(message in f for f in result["failures"]) diff --git a/policyengine_uk_data/tests/test_fit_report.py b/policyengine_uk_data/tests/test_fit_report.py index 276540e2..21ccf262 100644 --- a/policyengine_uk_data/tests/test_fit_report.py +++ b/policyengine_uk_data/tests/test_fit_report.py @@ -83,7 +83,7 @@ def test_builds_of_different_epochs_cannot_be_compared(): assert "saves the final epoch" in result["notes"][0] -def test_the_pull_request_workflow_uploads_only_calibration_logs(): +def test_the_pull_request_workflow_uploads_only_aggregate_reports(): """The datasets are licensed microdata: no job may upload or cache one.""" from pathlib import Path @@ -98,7 +98,9 @@ def test_the_pull_request_workflow_uploads_only_calibration_logs(): assert not uses.startswith("actions/cache"), step if uses.startswith("actions/upload-artifact"): uploads.append(step["with"]["path"]) - assert uploads and all(p.endswith("calibration_log.csv") for p in uploads), uploads + assert uploads and all( + p.endswith(("calibration_log.csv", "distribution_report.json")) for p in uploads + ), uploads assert jobs["fit-comparison"]["needs"] == ["test", "base-build"] diff --git a/policyengine_uk_data/utils/distribution_validation.py b/policyengine_uk_data/utils/distribution_validation.py new file mode 100644 index 00000000..f6f1dae8 --- /dev/null +++ b/policyengine_uk_data/utils/distribution_validation.py @@ -0,0 +1,377 @@ +"""Single-age population and HBAI income distribution of a built dataset. + +Calibration fits grouped targets. These checks measure what the groups can +hide, on the built dataset: + +- **HBAI.** Person-weighted medians, means and Gini coefficients of + equivalised household income, before and after housing costs. Also the + share of individuals below 60% of the median, for everyone, HBAI + dependent children, working-age adults and those at State Pension age. + Each is compared with DWP's HBAI FYE 2025 estimate and 95% confidence + interval (``targets/reference_data/hbai_fye_2025.yaml``), at period 2024. +- **Single ages.** Weighted people by sex and single year of age, at the + calibration year, against the ONS 2024-based projection the calibration + targets come from. Ages from 80 up are one band, because the dataset's + oldest records are in the 80-84 band (#593). + +``compare`` sets a pull request's build against its base's. While the +same-code spread of these measures is unmeasured (#598), it reports moves +and flags a measure that leaves the HBAI confidence interval, but fails +only on a broken or mismatched report. + + python -m policyengine_uk_data.utils.distribution_validation measure DATASET.h5 --out report.json + python -m policyengine_uk_data.utils.distribution_validation compare BASE.json HEAD.json [--markdown out.md] + +Reports hold aggregates only: no household's values. +""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import math +import sys +from pathlib import Path +from typing import Mapping, Optional + +import numpy as np +import yaml + +SCHEMA = "policyengine-uk-data-distribution/1" +HBAI_REFERENCE = ( + Path(__file__).parents[1] / "targets" / "reference_data" / "hbai_fye_2025.yaml" +) +WEEKS_PER_YEAR = 52 +TOP_AGE_BAND = 80 +SEXES = {"male": ("MALE", "Males"), "female": ("FEMALE", "Females")} + + +def weighted_median(values, weights) -> float: + values, weights = np.asarray(values, float), np.asarray(weights, float) + order = np.argsort(values, kind="mergesort") + cumulative = np.cumsum(weights[order]) + return float(values[order][np.searchsorted(cumulative, cumulative[-1] / 2)]) + + +def weighted_gini(values, weights) -> float: + """Gini coefficient of a weighted distribution (0 to 1).""" + values, weights = np.asarray(values, float), np.asarray(weights, float) + order = np.argsort(values, kind="mergesort") + x, w = values[order], weights[order] + cw = np.cumsum(w) + cxw = np.cumsum(x * w) + total_w, total_xw = cw[-1], cxw[-1] + # Area under the Lorenz curve by the trapezium rule over weighted points. + lorenz = cxw / total_xw + previous = np.concatenate([[0.0], lorenz[:-1]]) + area = np.sum(w / total_w * (lorenz + previous) / 2) + return float(1 - 2 * area) + + +def hbai_measures(income_bhc, income_ahc, weights, child, pensioner) -> dict: + """The HBAI reference measures from person-level arrays. + + ``income_*`` are each person's household equivalised income in £ a week; + ``weights`` are person weights; ``child`` and ``pensioner`` are the HBAI + groups (working-age adults are everyone else). + """ + child, pensioner = np.asarray(child, bool), np.asarray(pensioner, bool) + groups = { + "all": np.ones_like(child), + "children": child, + "working_age": ~child & ~pensioner, + "pensioners": pensioner, + } + out = {} + for suffix, income in (("bhc", income_bhc), ("ahc", income_ahc)): + income = np.asarray(income, float) + median = weighted_median(income, weights) + out[f"median_{suffix}"] = median + out[f"mean_{suffix}"] = float(np.average(income, weights=weights)) + out[f"gini_{suffix}"] = 100 * weighted_gini(income, weights) + low = income < 0.6 * median + for group, mask in groups.items(): + w = weights * mask + out[f"relative_low_income_{suffix}_{group}"] = float( + 100 * np.sum(w * low) / np.sum(w) + ) + return out + + +def load_hbai_reference(path=None) -> dict: + with open(path or HBAI_REFERENCE) as f: + raw = yaml.safe_load(f) + if raw.get("schema_version") != 1: + raise ValueError("HBAI reference: schema_version must be 1") + for name, m in raw["measures"].items(): + if not m["low"] <= m["estimate"] <= m["high"]: + raise ValueError(f"HBAI reference {name}: estimate outside its interval") + return raw + + +def single_ages(age, is_male, weights, reference: Mapping) -> dict: + """People by sex and single age, ages from ``TOP_AGE_BAND`` up as one band. + + ``reference`` maps (sex, age) to the ONS count, with ages from + ``TOP_AGE_BAND`` up summed by the caller. + """ + age = np.minimum(np.floor(np.asarray(age, float)), TOP_AGE_BAND).astype(int) + out = {} + for sex, mask in (("male", is_male), ("female", ~np.asarray(is_male, bool))): + counts = np.bincount( + age[mask], weights=np.asarray(weights)[mask], minlength=TOP_AGE_BAND + 1 + ) + out[sex] = { + str(a): {"value": float(counts[a]), "reference": float(reference[(sex, a)])} + for a in range(TOP_AGE_BAND + 1) + } + return out + + +def ons_single_ages(year: int) -> dict: + """The ONS 2024-based projection by sex and age, with 80+ as one band.""" + import pandas as pd + + from policyengine_uk_data.targets.sources.ons_demographics import ( + _download_uk_projection, + _year_column, + ) + + df = _download_uk_projection() + column = _year_column(df, year) + if column is None: + raise ValueError(f"the ONS projection has no {year}") + ages = pd.to_numeric( + df["Age"].astype(str).str.extract(r"(\d+)")[0], errors="coerce" + ) + out = {} + for sex, (_, ons_sex) in SEXES.items(): + rows = df["Sex"] == ons_sex + for a in range(TOP_AGE_BAND): + out[(sex, a)] = float(df.loc[rows & (ages == a), column].sum()) + out[(sex, TOP_AGE_BAND)] = float( + df.loc[rows & (ages >= TOP_AGE_BAND), column].sum() + ) + return out + + +def summarise_single_ages(rows: Mapping) -> dict: + errors = { + (sex, int(a)): cell["value"] / cell["reference"] - 1 + for sex, by_age in rows.items() + for a, cell in by_age.items() + if cell["reference"] > 0 + } + under = [abs(e) for (_, a), e in errors.items() if a < TOP_AGE_BAND] + teens = [abs(e) for (_, a), e in errors.items() if 15 <= a <= 19] + return { + "mean_abs_error_0_79": float(np.mean(under)), + "max_abs_error_0_79": float(np.max(under)), + "max_abs_error_15_19": float(np.max(teens)), + "error_80_plus": { + sex: float( + rows[sex][str(TOP_AGE_BAND)]["value"] + / rows[sex][str(TOP_AGE_BAND)]["reference"] + - 1 + ) + for sex in rows + }, + } + + +def measure( + dataset_path, hbai_period: Optional[int] = None, age_year: int = 2025 +) -> dict: + """The distribution report of a built dataset (aggregates only).""" + from policyengine_uk import Microsimulation + + reference = load_hbai_reference() + period = hbai_period or int(reference["period"]) + sim = Microsimulation(dataset=str(dataset_path)) + + def person(variable, year): + return np.asarray(sim.calculate(variable, year, map_to="person").values, float) + + weights = person("household_weight", period) + measures = hbai_measures( + person("equiv_hbai_household_net_income", period) / WEEKS_PER_YEAR, + person("equiv_hbai_household_net_income_ahc", period) / WEEKS_PER_YEAR, + weights, + person("is_hbai_dependent_child", period) > 0, + person("is_SP_age", period) > 0, + ) + hbai = {} + for name, ref in reference["measures"].items(): + value = measures[name] + hbai[name] = { + "value": value, + "reference": ref["estimate"], + "low": ref["low"], + "high": ref["high"], + "inside": bool(ref["low"] <= value <= ref["high"]), + } + gender = np.asarray(sim.calculate("gender", age_year).values).astype(str) + ages = single_ages( + person("age", age_year), + gender == "MALE", + person("household_weight", age_year), + ons_single_ages(age_year), + ) + data = Path(dataset_path).read_bytes() + return { + "schema": SCHEMA, + "dataset": Path(dataset_path).name, + "sha256": hashlib.sha256(data).hexdigest(), + "periods": {"hbai": period, "single_age": age_year}, + "hbai": hbai, + "single_age": ages, + "summary": { + "hbai_outside_interval": sorted( + n for n, m in hbai.items() if not m["inside"] + ), + "single_age": summarise_single_ages(ages), + }, + } + + +def _finite(x) -> bool: + return isinstance(x, (int, float)) and math.isfinite(x) + + +def compare(base: Mapping, head: Mapping) -> dict: + """A head build's distribution against its base's. + + Fails only on a broken or mismatched report; flags a measure that + leaves its HBAI interval, or moves away from its reference. + """ + failures, flags = [], [] + for side, report in (("base", base), ("head", head)): + if report.get("schema") != SCHEMA: + failures.append(f"The {side} distribution report is missing or malformed.") + if failures: + return { + "status": "fail", + "failures": failures, + "flags": [], + "rows": [], + "single_age": {}, + } + if base["periods"] != head["periods"]: + failures.append( + f"The reports measured different periods ({base['periods']} and {head['periods']})." + ) + rows = [] + for name, h in head["hbai"].items(): + b = base["hbai"].get(name) + if b is None or not (_finite(h["value"]) and _finite(b["value"])): + failures.append(f"HBAI {name} is missing or not finite on one side.") + continue + worse = abs(h["value"] - h["reference"]) - abs(b["value"] - b["reference"]) + row = { + "measure": name, + "base": b["value"], + "head": h["value"], + "reference": h["reference"], + "low": h["low"], + "high": h["high"], + "base_inside": b["inside"], + "head_inside": h["inside"], + "further_from_reference": worse, + } + rows.append(row) + if b["inside"] and not h["inside"]: + flags.append( + f"HBAI {name} leaves its 95% interval: {b['value']:.2f} -> {h['value']:.2f} " + f"(HBAI {h['reference']:.2f}, {h['low']:.2f} to {h['high']:.2f})." + ) + sb, sh = base["summary"]["single_age"], head["summary"]["single_age"] + single = { + k: {"base": sb[k], "head": sh[k]} + for k in ("mean_abs_error_0_79", "max_abs_error_0_79", "max_abs_error_15_19") + } + return { + "status": "fail" if failures else "pass", + "failures": failures, + "flags": flags, + "rows": rows, + "single_age": single, + } + + +def render(result: Mapping) -> str: + lines = [ + f"## Income distribution and single ages against the base build: {result['status']}", + "", + ] + if result["failures"]: + lines += ["### Failures", ""] + [f"- {f}" for f in result["failures"]] + [""] + if result["rows"]: + lines += [ + "HBAI FYE 2025 (DWP), period 2024: equivalised income in £ a week, rates and Gini in per cent. " + "A value outside HBAI's 95% interval is marked (outside).", + "", + "| Measure | Base | Head | HBAI | 95% interval |", + "|---|---:|---:|---:|---|", + ] + for r in result["rows"]: + mark = lambda inside: "" if inside else " (outside)" # noqa: E731 + lines.append( + f"| {r['measure']} | {r['base']:.2f}{mark(r['base_inside'])} | {r['head']:.2f}{mark(r['head_inside'])} " + f"| {r['reference']:.2f} | {r['low']:.2f} to {r['high']:.2f} |" + ) + lines.append("") + if result["single_age"]: + lines += [ + "| Single ages against the ONS projection | Base | Head |", + "|---|---:|---:|", + ] + labels = { + "mean_abs_error_0_79": "Mean error, ages 0 to 79", + "max_abs_error_0_79": "Largest error, ages 0 to 79", + "max_abs_error_15_19": "Largest error, ages 15 to 19", + } + for key, label in labels.items(): + cell = result["single_age"][key] + lines.append( + f"| {label} | {100 * cell['base']:.1f}% | {100 * cell['head']:.1f}% |" + ) + lines.append("") + lines += ["### Flags", ""] + ([f"- {f}" for f in result["flags"]] or ["None."]) + lines += [ + "", + "These measures report and flag but do not fail yet: their same-code spread is unmeasured (#598).", + ] + return "\n".join(lines) + "\n" + + +def main(argv=None) -> int: + parser = argparse.ArgumentParser( + prog="python -m policyengine_uk_data.utils.distribution_validation" + ) + sub = parser.add_subparsers(dest="command", required=True) + m = sub.add_parser("measure", help="measure a built dataset") + m.add_argument("dataset") + m.add_argument("--out", required=True) + c = sub.add_parser("compare", help="compare a head build's report with its base's") + c.add_argument("base") + c.add_argument("head") + c.add_argument("--markdown") + args = parser.parse_args(argv) + if args.command == "measure": + report = measure(args.dataset) + Path(args.out).write_text(json.dumps(report, indent=1)) + print(json.dumps(report["summary"], indent=1)) + return 0 + base = json.loads(Path(args.base).read_text()) + head = json.loads(Path(args.head).read_text()) + result = compare(base, head) + text = render(result) + if args.markdown: + Path(args.markdown).write_text(text) + print(text) + return 1 if result["status"] == "fail" else 0 + + +if __name__ == "__main__": + sys.exit(main()) From b6c83c4ee492859d437707b5a66c45a95d3cc4c9 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Thu, 15 Oct 2026 08:33:37 -0400 Subject: [PATCH 2/5] Measure HBAI incomes as DWP does; fail on incomplete reports Sol review of #600: - Negative incomes before housing costs are reset to zero, and income after housing costs is derived from the adjusted figure (DWP HBAI methodology FYE 2025, "Negative incomes"). Weekly amounts use 365.25 / 7, as the FRS loader annualises them. - A report missing any reference measure, or with a missing or non-finite single-age summary, fails the comparison. - A missing distribution report fails the job. On uk-data 1.58.1, four of 14 measures are now outside their intervals (both Gini coefficients, child low income BHC and AHC); median AHC moves inside. Co-Authored-By: Claude Opus 5.5 --- .github/workflows/pull_request.yaml | 1 + .../tests/test_distribution_validation.py | 45 ++++++++++++++++ .../utils/distribution_validation.py | 51 +++++++++++++++++-- 3 files changed, 94 insertions(+), 3 deletions(-) diff --git a/.github/workflows/pull_request.yaml b/.github/workflows/pull_request.yaml index c1a63b81..23d28aeb 100644 --- a/.github/workflows/pull_request.yaml +++ b/.github/workflows/pull_request.yaml @@ -243,6 +243,7 @@ jobs: if [ "$distribution_status" != 0 ]; then echo "status=1" >> "$GITHUB_OUTPUT"; fi else { echo; echo "The income distribution and single-age reports are missing for one build."; } >> fit.md + echo "status=1" >> "$GITHUB_OUTPUT" fi cat fit.md >> "$GITHUB_STEP_SUMMARY" - name: Post or update the pull request comment diff --git a/policyengine_uk_data/tests/test_distribution_validation.py b/policyengine_uk_data/tests/test_distribution_validation.py index 3c86db33..ae8103a8 100644 --- a/policyengine_uk_data/tests/test_distribution_validation.py +++ b/policyengine_uk_data/tests/test_distribution_validation.py @@ -160,3 +160,48 @@ def test_compare_fails_on_a_broken_or_mismatched_report(change, message): result = compare(_report(), head) assert result["status"] == "fail" assert any(message in f for f in result["failures"]) + + +def test_negative_incomes_before_housing_costs_are_reset_to_zero_as_hbai_does(): + """Sol review of #600: HBAI resets negative BHC to zero and derives AHC from + the adjusted figure, so AHC can still be negative.""" + from policyengine_uk_data.utils.distribution_validation import ( + WEEKS_PER_YEAR, + hbai_incomes, + ) + + weekly_bhc = np.array([-100.0, 100.0, 200.0, 300.0, 400.0]) + annual_bhc = weekly_bhc * WEEKS_PER_YEAR + annual_ahc = (weekly_bhc - 50) * WEEKS_PER_YEAR + bhc, ahc = hbai_incomes(annual_bhc, annual_ahc, np.ones(5), np.ones(5)) + np.testing.assert_allclose(bhc, [0, 100, 200, 300, 400]) + np.testing.assert_allclose(ahc, [-50, 50, 150, 250, 350]) + m = hbai_measures(bhc, ahc, np.ones(5), np.zeros(5, bool), np.zeros(5, bool)) + assert m["mean_bhc"] == pytest.approx(200) and m["mean_ahc"] == pytest.approx(150) + assert m["gini_bhc"] == pytest.approx(40.0) + assert m["gini_ahc"] == pytest.approx(53.333333, rel=1e-6) + + +def test_an_incomplete_or_non_finite_report_fails(): + """Sol review of #600: every reference measure and age summary must be there.""" + head = _report() + del head["hbai"]["median_bhc"] + result = compare(_report(), head) + assert result["status"] == "fail" + assert result["failures"] == ["The head report has no median_bhc."] + head = _report() + head["summary"]["single_age"]["mean_abs_error_0_79"] = float("nan") + result = compare(_report(), head) + assert result["status"] == "fail" and "not finite" in result["failures"][0] + + +def test_the_workflow_fails_when_a_distribution_report_is_missing(): + from pathlib import Path + + text = ( + Path(__file__).parents[2] / ".github" / "workflows" / "pull_request.yaml" + ).read_text() + block = text[ + text.index("The income distribution and single-age reports are missing") : + ] + assert 'echo "status=1" >> "$GITHUB_OUTPUT"' in block.split("\n fi", 1)[0] diff --git a/policyengine_uk_data/utils/distribution_validation.py b/policyengine_uk_data/utils/distribution_validation.py index f6f1dae8..bc7d87a0 100644 --- a/policyengine_uk_data/utils/distribution_validation.py +++ b/policyengine_uk_data/utils/distribution_validation.py @@ -42,7 +42,8 @@ HBAI_REFERENCE = ( Path(__file__).parents[1] / "targets" / "reference_data" / "hbai_fye_2025.yaml" ) -WEEKS_PER_YEAR = 52 +# The FRS loader annualises weekly amounts with 365.25 / 7 (datasets/frs.py). +WEEKS_PER_YEAR = 365.25 / 7 TOP_AGE_BAND = 80 SEXES = {"male": ("MALE", "Males"), "female": ("FEMALE", "Females")} @@ -69,6 +70,25 @@ def weighted_gini(values, weights) -> float: return float(1 - 2 * area) +def hbai_incomes(bhc, ahc, equivalisation_bhc, equivalisation_ahc): + """Equivalised weekly incomes as HBAI defines them. + + HBAI resets negative incomes before housing costs to zero, and derives + the income after housing costs from the adjusted figure, so it can still + be negative (DWP, HBAI background information and methodology, FYE 2025, + "Negative incomes"). ``bhc`` and ``ahc`` are annual household incomes. + """ + bhc, ahc = np.asarray(bhc, float), np.asarray(ahc, float) + housing_costs = bhc - ahc + adjusted = np.maximum(bhc, 0.0) + return ( + adjusted / np.asarray(equivalisation_bhc, float) / WEEKS_PER_YEAR, + (adjusted - housing_costs) + / np.asarray(equivalisation_ahc, float) + / WEEKS_PER_YEAR, + ) + + def hbai_measures(income_bhc, income_ahc, weights, child, pensioner) -> dict: """The HBAI reference measures from person-level arrays. @@ -194,9 +214,15 @@ def person(variable, year): return np.asarray(sim.calculate(variable, year, map_to="person").values, float) weights = person("household_weight", period) + income_bhc, income_ahc = hbai_incomes( + person("hbai_household_net_income", period), + person("hbai_household_net_income_ahc", period), + person("household_equivalisation_bhc", period), + person("household_equivalisation_ahc", period), + ) measures = hbai_measures( - person("equiv_hbai_household_net_income", period) / WEEKS_PER_YEAR, - person("equiv_hbai_household_net_income_ahc", period) / WEEKS_PER_YEAR, + income_bhc, + income_ahc, weights, person("is_hbai_dependent_child", period) > 0, person("is_SP_age", period) > 0, @@ -257,6 +283,25 @@ def compare(base: Mapping, head: Mapping) -> dict: "rows": [], "single_age": {}, } + required = set(load_hbai_reference()["measures"]) + for side, report in (("base", base), ("head", head)): + missing = sorted(required - set(report.get("hbai") or {})) + if missing: + failures.append(f"The {side} report has no {', '.join(missing)}.") + summary = (report.get("summary") or {}).get("single_age") or {} + for key in ("mean_abs_error_0_79", "max_abs_error_0_79", "max_abs_error_15_19"): + if not _finite(summary.get(key)): + failures.append( + f"The {side} report's single-age {key} is missing or not finite." + ) + if failures: + return { + "status": "fail", + "failures": failures, + "flags": [], + "rows": [], + "single_age": {}, + } if base["periods"] != head["periods"]: failures.append( f"The reports measured different periods ({base['periods']} and {head['periods']})." From 48f88ad60c7319fb1c89e32f331db73d59ff11e5 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Thu, 15 Oct 2026 09:26:36 -0400 Subject: [PATCH 3/5] An empty HBAI group has no rate, not a NaN hbai_measures divided by a zero group weight, which warned and gave NaN. It now returns None for an empty group, the report marks it outside its interval, and compare fails on it as missing or not finite. Co-Authored-By: Claude Opus 5.5 --- .../tests/test_distribution_validation.py | 11 +++++++++++ policyengine_uk_data/utils/distribution_validation.py | 8 +++++--- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/policyengine_uk_data/tests/test_distribution_validation.py b/policyengine_uk_data/tests/test_distribution_validation.py index ae8103a8..05ffc884 100644 --- a/policyengine_uk_data/tests/test_distribution_validation.py +++ b/policyengine_uk_data/tests/test_distribution_validation.py @@ -180,6 +180,11 @@ def test_negative_incomes_before_housing_costs_are_reset_to_zero_as_hbai_does(): assert m["mean_bhc"] == pytest.approx(200) and m["mean_ahc"] == pytest.approx(150) assert m["gini_bhc"] == pytest.approx(40.0) assert m["gini_ahc"] == pytest.approx(53.333333, rel=1e-6) + # No children or pensioners: those groups have no rate, not a NaN. + assert m["relative_low_income_bhc_children"] is None + assert m["relative_low_income_ahc_pensioners"] is None + # Median 200, line 120: incomes 0 and 100 are below it. + assert m["relative_low_income_bhc_working_age"] == pytest.approx(40.0) def test_an_incomplete_or_non_finite_report_fails(): @@ -193,6 +198,12 @@ def test_an_incomplete_or_non_finite_report_fails(): head["summary"]["single_age"]["mean_abs_error_0_79"] = float("nan") result = compare(_report(), head) assert result["status"] == "fail" and "not finite" in result["failures"][0] + head = _report() + name = next(iter(head["hbai"])) + head["hbai"][name]["value"] = None # an empty group + result = compare(_report(), head) + assert result["status"] == "fail" + assert f"HBAI {name} is missing or not finite on one side." in result["failures"] def test_the_workflow_fails_when_a_distribution_report_is_missing(): diff --git a/policyengine_uk_data/utils/distribution_validation.py b/policyengine_uk_data/utils/distribution_validation.py index bc7d87a0..67ff43b1 100644 --- a/policyengine_uk_data/utils/distribution_validation.py +++ b/policyengine_uk_data/utils/distribution_validation.py @@ -113,8 +113,10 @@ def hbai_measures(income_bhc, income_ahc, weights, child, pensioner) -> dict: low = income < 0.6 * median for group, mask in groups.items(): w = weights * mask - out[f"relative_low_income_{suffix}_{group}"] = float( - 100 * np.sum(w * low) / np.sum(w) + total = np.sum(w) + # An empty group has no rate; compare fails on it. + out[f"relative_low_income_{suffix}_{group}"] = ( + float(100 * np.sum(w * low) / total) if total > 0 else None ) return out @@ -235,7 +237,7 @@ def person(variable, year): "reference": ref["estimate"], "low": ref["low"], "high": ref["high"], - "inside": bool(ref["low"] <= value <= ref["high"]), + "inside": value is not None and bool(ref["low"] <= value <= ref["high"]), } gender = np.asarray(sim.calculate("gender", age_year).values).astype(str) ages = single_ages( From 512e85e305c1882db86833c568f3b206447647ad Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Thu, 15 Oct 2026 15:36:58 -0400 Subject: [PATCH 4/5] Measure the distribution after the tests and report a crashed comparison (delta review of #600) The measurement downloads the ONS projection, so it now runs after the test suite with if: !cancelled(): an ONS outage no longer skips the tests and the fit summary. A crashed distribution comparison writes a FAIL line into the comment, and the final error names both checks. Co-Authored-By: Claude Opus 5.5 --- .github/workflows/pull_request.yaml | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/.github/workflows/pull_request.yaml b/.github/workflows/pull_request.yaml index 897908b7..96e6e99a 100644 --- a/.github/workflows/pull_request.yaml +++ b/.github/workflows/pull_request.yaml @@ -86,13 +86,6 @@ jobs: run: uv run --frozen make data env: TESTING: "1" - - name: Measure the income distribution and single ages - run: uv run --frozen python -m policyengine_uk_data.utils.distribution_validation measure policyengine_uk_data/storage/enhanced_frs_2024_25.h5 --out distribution_report.json - - name: Save the distribution report - uses: actions/upload-artifact@v7 - with: - name: distribution_report.json - path: distribution_report.json # The raw calibration logs stay on the runner: a target supported by a # few households, divided by a published count, could give one # household's value back. Only the summary the fit comparison needs is @@ -111,6 +104,17 @@ jobs: # Relaxes built-dataset outcome tolerances to match the reduced # fidelity of the TESTING build above (see test_cgt_band_donors). TESTING: "1" + # After the tests, and whatever they did: the measurement downloads the + # ONS projection, and an ONS outage must not skip the test suite. + - name: Measure the income distribution and single ages + if: ${{ !cancelled() }} + run: uv run --frozen python -m policyengine_uk_data.utils.distribution_validation measure policyengine_uk_data/storage/enhanced_frs_2024_25.h5 --out distribution_report.json + - name: Save the distribution report + if: ${{ !cancelled() }} + uses: actions/upload-artifact@v7 + with: + name: distribution_report.json + path: distribution_report.json # The no-worsening check (W3). The base is built the same way as the head: # TESTING, one clone and the head's lockfile, so both builds run one model # and environment and the comparison measures the data change. As in the @@ -238,7 +242,7 @@ jobs: reports/base_distribution_report.json reports/distribution_report.json --markdown distribution.md distribution_status=$? set -e - { echo; cat distribution.md; } >> fit.md + { echo; cat distribution.md 2>/dev/null || echo '## Income distribution and single ages: FAIL (the comparison did not finish)'; } >> fit.md if [ "$distribution_status" != 0 ]; then echo "status=1" >> "$GITHUB_OUTPUT"; fi else { echo; echo "The income distribution and single-age reports are missing for one build."; } >> fit.md @@ -272,5 +276,5 @@ jobs: - name: Fail on a fit regression if: steps.compare.outputs.status != '0' run: | - echo "::error::The head build's calibration fit is worse than the base's; see the job summary." + echo "::error::The calibration fit or distribution comparison failed; see the job summary." exit 1 From 45d474284aed3527680a0b2e9439aa444bba88c2 Mon Sep 17 00:00:00 2001 From: Max Ghenis Date: Thu, 15 Oct 2026 23:17:22 -0400 Subject: [PATCH 5/5] Merge #599; test the fit summary's and the base's step order (delta review of #600) Co-Authored-By: Claude Opus 5.5 --- .../tests/test_distribution_validation.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/policyengine_uk_data/tests/test_distribution_validation.py b/policyengine_uk_data/tests/test_distribution_validation.py index 7d61f0a7..7738583d 100644 --- a/policyengine_uk_data/tests/test_distribution_validation.py +++ b/policyengine_uk_data/tests/test_distribution_validation.py @@ -232,3 +232,14 @@ def test_the_measurement_runs_after_the_tests_and_cannot_skip_them(): assert names.index("Save the fit summary") < measure assert steps[measure]["if"] == "${{ !cancelled() }}" assert steps[measure]["continue-on-error"] is True + # The fit summary runs after the tests too, whatever they did. + summarise = names.index("Summarise the calibration fit") + assert names.index("Run tests") < summarise + assert steps[summarise]["if"] == "${{ !cancelled() }}" + assert steps[names.index("Save the fit summary")]["if"] == "${{ !cancelled() }}" + # On the base, the distribution measurement runs after the fit summary. + base = yaml.safe_load(workflow.read_text())["jobs"]["base-build"]["steps"] + base_names = [step.get("name") for step in base] + assert base_names.index("Save the base fit summary") < base_names.index( + "Measure the base's income distribution and single ages" + )