diff --git a/changelog.d/moop-spm-experiment.added.md b/changelog.d/moop-spm-experiment.added.md new file mode 100644 index 00000000..36ddb91f --- /dev/null +++ b/changelog.d/moop-spm-experiment.added.md @@ -0,0 +1 @@ +Add an experiment that measures how the choice of method for imputing medical out-of-pocket spending (MOOP) changes Supplemental Poverty Measure poverty on held-out 2026 CPS ASEC SPM units (`experiments/moop_spm`). diff --git a/experiments/moop_spm/NOTE.md b/experiments/moop_spm/NOTE.md new file mode 100644 index 00000000..f9309ce5 --- /dev/null +++ b/experiments/moop_spm/NOTE.md @@ -0,0 +1,45 @@ +# How the MOOP imputation method moves SPM poverty + +The American Community Survey (ACS) does not ask about medical out-of-pocket spending (MOOP), so researchers who build a Supplemental Poverty Measure (SPM) on the ACS impute it. For 2014, Census researchers imputed it with predictive mean matching by insurance status, and subtracting their imputed MOOP raised the ACS SPM rate by 4.57 percentage points. The Current Population Survey's Annual Social and Economic Supplement (CPS ASEC) asks about MOOP, and for the same year, subtracting reported MOOP raised its SPM rate by 3.48 points. Yet the ACS MOOP total of $561 billion trailed the CPS ASEC's $567 billion by 1% ([Fox, Pacas and Glassman, 2017](https://www.census.gov/content/dam/Census/library/working-papers/2017/demo/Fox-ACSSPM2014.pdf), slides 10, 11 and 26). The surveys also differ in income and other resources, so we test the imputation alone. We impute MOOP within the CPS ASEC and check each method against the file's own MOOP. + +## Data and design + +We use the 134,729 people, 55,401 SPM units and 53,344 households in the 2026 CPS ASEC public-use file, which covers calendar-year 2025. We impute each unit's SPM medical expenses (`SPM_MEDXPNS`), its members' MOOP plus their Medicare premiums. Census subtracts them from SPM resources dollar for dollar, so we swap each unit's reported value for an imputed one and compare the new resources with an unchanged threshold. With reported MOOP and the ASEC person weight, we compute an SPM rate of 13.11% with a 90% margin of error of 0.31 points. Census published 13.11% and 0.31 in [P60-290](https://www.census.gov/library/publications/2026/demo/p60-290.html). Subtracting MOOP raises the rate by 2.25 points (SE 0.08), and Census reports 2.28. + +We limit every method to predictors an ACS-style file carries: members' ages, sex, disability and current health coverage, the unit's cash income on its own and relative to its SPM threshold, family type, tenure and state. We assign whole households to 5 folds and impute each fold from models we fit on the rest, without survey weights. We average 10 draws for each stochastic method and compute standard errors from the 160 replicate weights, treating each fitted model as fixed. + +- In method (a), we fit a logistic model of whether a unit has any MOOP. Each unit with a predicted probability of 0.5 or more gets the mean MOOP of donor spenders in its cell, and the rest get zero (logit plus cell means). In the OLS variant, predicted spenders get an OLS prediction among donor spenders instead. We form cells from coverage group, head's age, cash income relative to the threshold and unit size. +- In method (b), logit draws with donors, we draw participation at random from the same model's probabilities, and each spender takes a random donor spender's MOOP from its cell. +- In method (c), we fit quantile regression forests from PolicyEngine's microimpute package, alone (QRF) and inside its zero-inflated wrapper (ZI-QRF). We give each unit the value at a random quantile of its predicted distribution. In PolicyEngine's microimpute [paper](../../paper/main.pdf), María Juaristi found that quantile regression forests had the lowest average quantile loss of five methods for imputing Survey of Consumer Finances (SCF) net worth onto the CPS. +- We also apply predictive mean matching (PMM) within five coverage groups, where the Census researchers used eight. In each group, we regress log(1 + MOOP) on the predictors across all donors, including those with zero spending. Each recipient then takes the MOOP of a random one of the five donors whose predictions lie nearest its own. + +## Results + +| | Reported | (a) Logit + cell mean | (a) Logit + OLS mean | (b) Logit draw + donor | PMM | (c) QRF | (c) ZI-QRF | +|---|---:|---:|---:|---:|---:|---:|---:| +| SPM rate, all people (%) | 13.11 | 13.42 | 13.46 | 13.50 | 13.43 | 13.04 | 13.07 | +| Change, all people (pp, SE) | | +0.30 (0.08) | +0.34 (0.08) | +0.39 (0.07) | +0.32 (0.07) | -0.07 (0.07) | -0.05 (0.07) | +| Change, 65 and older (pp) | | +0.98 | +1.15 | +0.75 | +0.60 | -0.10 | -0.13 | +| Change, with a disability (pp) | | +1.20 | +1.40 | +0.59 | +0.41 | -0.21 | -0.22 | +| Change, Medicare (pp) | | +1.20 | +1.23 | +0.85 | +0.59 | -0.08 | -0.11 | +| Change, resources before MOOP 1.0 to 1.25 x threshold (pp) | | +14.69 | +14.37 | +4.88 | +4.19 | +0.55 | +0.61 | +| Mean MOOP, resources before MOOP 1.0 to 1.25 x threshold | $2,385 | $3,189 | $3,108 | $2,949 | $2,799 | $2,352 | $2,372 | +| Units with zero MOOP (%) | 6.4 | 0.5 | 4.3 | 6.5 | 6.2 | 6.2 | 6.1 | +| Median MOOP | $3,700 | $6,132 | $5,891 | $3,678 | $3,624 | $3,623 | $3,610 | +| 99th percentile MOOP | $33,640 | $15,585 | $15,769 | $33,147 | $32,666 | $31,356 | $31,424 | +| Total MOOP ($ billions) | 859 | 880 | 876 | 855 | 846 | 828 | 828 | +| People newly poor / newly not poor (%) | | 1.37 / 1.07 | 1.40 / 1.06 | 1.74 / 1.36 | 1.63 / 1.31 | 1.33 / 1.40 | 1.37 / 1.41 | + +We compute rates and changes over people and MOOP statistics over SPM units, and show standard errors in parentheses. + +We estimate that logit plus cell means raises the SPM rate by 0.30 percentage points (SE 0.08), or 13% of MOOP's own 2.25-point effect. The OLS variant adds 0.34 points, logit draws with donors 0.39 and PMM 0.32. QRF changes the rate by -0.07 points (SE 0.07) and ZI-QRF by -0.05 (SE 0.07), though their draws fall 3.6% and 3.7% short of the reported MOOP total. For people 65 and older, the rate changes by +0.98 points under logit plus cell means, +0.60 under PMM and -0.10 under QRF. For people with a disability, it changes by +1.20, +0.41 and -0.21. + +Logit draws with donors track the reported distribution, with 6.5% of units at zero against 6.4% and a 99th percentile of $33,147 against $33,640, yet raise the rate by 0.39 points. A unit falls below the line when its MOOP exceeds the margin by which its other resources clear the threshold, and reported MOOP rises with resources. Units whose resources before MOOP run 1.0 to 1.25 times their threshold spend $2,385 on average, and units at 2 to 3 times spend $6,086. In that first band above the line, cell means assign $3,189, donor draws $2,949, PMM $2,799 and QRF $2,352. The band's 24.9% poverty rate rises by 14.69, 4.88, 4.19 and 0.55 points under those four methods. + +Under the 0.5 cutoff, 0.5% of units get no MOOP, against 6.4% in the file. Drawing participation instead lowers the national change from +0.30 to +0.10 points with cell means, and from +0.59 to +0.39 with donor draws (Table 3 in `TABLE.md`). When we draw participation and assign cell means, the rate still rises 11.50 points for people at 1.0 to 1.25 times the threshold and falls 4.80 points at 1.25 to 1.5. Cell means also cut the 99th percentile from $33,640 to $15,585. Under the same cutoff, coverage-group means alone, without the age, income and size cells, raise the national rate by 2.07 points, 92% of MOOP's own effect. + +Every method moves people across the line in both directions. QRF moves 1.33% of people into poverty and 1.40% out of it, and logit plus cell means moves 1.37% in and 1.07% out. Census researchers turned to the ACS for state and substate SPM estimates, where these flows need not offset. + +Census imputes MOOP for 43% of the people in the file and the Medicare premium for 13%, so part of our benchmark comes from its imputation. In the 44% of SPM units where Census imputes neither for any member, methods (a), (b) and PMM raise the rate by 0.26 to 0.52 points. QRF and ZI-QRF change it by +0.01 and +0.07 points (SE 0.09 and 0.10). Because our donors and recipients share one survey's questions, income measures and population, we do not test the differences between surveys that an ACS imputation also faces. + +For each method's change by age, disability, coverage type and resources relative to the threshold, see Table 1 in `TABLE.md`; for the code, see `README.md`. diff --git a/experiments/moop_spm/NOTE.template.md b/experiments/moop_spm/NOTE.template.md new file mode 100644 index 00000000..6abab906 --- /dev/null +++ b/experiments/moop_spm/NOTE.template.md @@ -0,0 +1,32 @@ +# How the MOOP imputation method moves SPM poverty + +The American Community Survey (ACS) does not ask about medical out-of-pocket spending (MOOP), so researchers who build a Supplemental Poverty Measure (SPM) on the ACS impute it. For 2014, Census researchers imputed it with predictive mean matching by insurance status, and subtracting their imputed MOOP raised the ACS SPM rate by 4.57 percentage points. The Current Population Survey's Annual Social and Economic Supplement (CPS ASEC) asks about MOOP, and for the same year, subtracting reported MOOP raised its SPM rate by 3.48 points. Yet the ACS MOOP total of $561 billion trailed the CPS ASEC's $567 billion by 1% ([Fox, Pacas and Glassman, 2017](https://www.census.gov/content/dam/Census/library/working-papers/2017/demo/Fox-ACSSPM2014.pdf), slides 10, 11 and 26). The surveys also differ in income and other resources, so we test the imputation alone. We impute MOOP within the CPS ASEC and check each method against the file's own MOOP. + +## Data and design + +We use the {n_persons} people, {n_units} SPM units and {n_households} households in the 2026 CPS ASEC public-use file, which covers calendar-year 2025. We impute each unit's SPM medical expenses (`SPM_MEDXPNS`), its members' MOOP plus their Medicare premiums. Census subtracts them from SPM resources dollar for dollar, so we swap each unit's reported value for an imputed one and compare the new resources with an unchanged threshold. With reported MOOP and the ASEC person weight, we compute an SPM rate of {rate2_all_reported}% with a 90% margin of error of {reported_moe90} points. Census published {published_rate}% and {published_moe90} in [P60-290](https://www.census.gov/library/publications/2026/demo/p60-290.html). Subtracting MOOP raises the rate by {moop_effect} points (SE {moop_effect_se}), and Census reports {published_moop_effect}. + +We limit every method to predictors an ACS-style file carries: members' ages, sex, disability and current health coverage, the unit's cash income on its own and relative to its SPM threshold, family type, tenure and state. We assign whole households to {folds} folds and impute each fold from models we fit on the rest, without survey weights. We average {draws} draws for each stochastic method and compute standard errors from the 160 replicate weights, treating each fitted model as fixed. + +- In method (a), we fit a logistic model of whether a unit has any MOOP. Each unit with a predicted probability of 0.5 or more gets the mean MOOP of donor spenders in its cell, and the rest get zero (logit plus cell means). In the OLS variant, predicted spenders get an OLS prediction among donor spenders instead. We form cells from coverage group, head's age, cash income relative to the threshold and unit size. +- In method (b), logit draws with donors, we draw participation at random from the same model's probabilities, and each spender takes a random donor spender's MOOP from its cell. +- In method (c), we fit quantile regression forests from PolicyEngine's microimpute package, alone (QRF) and inside its zero-inflated wrapper (ZI-QRF). We give each unit the value at a random quantile of its predicted distribution. In PolicyEngine's microimpute [paper](../../paper/main.pdf), María Juaristi found that quantile regression forests had the lowest average quantile loss of five methods for imputing Survey of Consumer Finances (SCF) net worth onto the CPS. +- We also apply predictive mean matching (PMM) within five coverage groups, where the Census researchers used eight. In each group, we regress log(1 + MOOP) on the predictors across all donors, including those with zero spending. Each recipient then takes the MOOP of a random one of the five donors whose predictions lie nearest its own. + +## Results + +{summary_table} + +We compute rates and changes over people and MOOP statistics over SPM units, and show standard errors in parentheses. + +We estimate that logit plus cell means raises the SPM rate by {pp_all_logit_cell_mean} percentage points (SE {se_all_logit_cell_mean}), or {share_of_effect_logit_cell_mean}% of MOOP's own {moop_effect}-point effect. The OLS variant adds {pp_all_logit_reg_mean} points, logit draws with donors {pp_all_logit_draw_hotdeck} and PMM {pp_all_pmm}. QRF changes the rate by {diff_all_qrf} points (SE {se_all_qrf}) and ZI-QRF by {diff_all_zi_qrf} (SE {se_all_zi_qrf}), though their draws fall {agg_gap_qrf}% and {agg_gap_zi_qrf}% short of the reported MOOP total. For people 65 and older, the rate changes by {diff_senior_logit_cell_mean} points under logit plus cell means, {diff_senior_pmm} under PMM and {diff_senior_qrf} under QRF. For people with a disability, it changes by {diff_disabled_logit_cell_mean}, {diff_disabled_pmm} and {diff_disabled_qrf}. + +Logit draws with donors track the reported distribution, with {zero_logit_draw_hotdeck}% of units at zero against {zero_reported}% and a 99th percentile of {p99_logit_draw_hotdeck} against {p99_reported}, yet raise the rate by {pp_all_logit_draw_hotdeck} points. A unit falls below the line when its MOOP exceeds the margin by which its other resources clear the threshold, and reported MOOP rises with resources. Units whose resources before MOOP run 1.0 to 1.25 times their threshold spend {bandmean1_reported} on average, and units at 2 to 3 times spend {bandmean4_reported}. In that first band above the line, cell means assign {bandmean1_logit_cell_mean}, donor draws {bandmean1_logit_draw_hotdeck}, PMM {bandmean1_pmm} and QRF {bandmean1_qrf}. The band's {rate_band1_reported}% poverty rate rises by {pp_band1_logit_cell_mean}, {pp_band1_logit_draw_hotdeck}, {pp_band1_pmm} and {pp_band1_qrf} points under those four methods. + +Under the 0.5 cutoff, {zero_logit_cell_mean}% of units get no MOOP, against {zero_reported}% in the file. Drawing participation instead lowers the national change from {diff_all_logit_cell_mean} to {diff_all_logit_draw_cell_mean} points with cell means, and from {diff_all_logit_threshold_hotdeck} to {diff_all_logit_draw_hotdeck} with donor draws (Table 3 in `TABLE.md`). When we draw participation and assign cell means, the rate still rises {pp_band1_logit_draw_cell_mean} points for people at 1.0 to 1.25 times the threshold and falls {pp_band2_logit_draw_cell_mean} points at 1.25 to 1.5. Cell means also cut the 99th percentile from {p99_reported} to {p99_logit_cell_mean}. Under the same cutoff, coverage-group means alone, without the age, income and size cells, raise the national rate by {pp_all_logit_group_mean} points, {share_of_effect_logit_group_mean}% of MOOP's own effect. + +Every method moves people across the line in both directions. QRF moves {newly_poor_qrf}% of people into poverty and {newly_nonpoor_qrf}% out of it, and logit plus cell means moves {newly_poor_logit_cell_mean}% in and {newly_nonpoor_logit_cell_mean}% out. Census researchers turned to the ACS for state and substate SPM estimates, where these flows need not offset. + +Census imputes MOOP for {share_persons_moop_allocated}% of the people in the file and the Medicare premium for {share_persons_premium_allocated}%, so part of our benchmark comes from its imputation. In the {share_units_no_allocation}% of SPM units where Census imputes neither for any member, methods (a), (b) and PMM raise the rate by {subset_pp_low_ab_pmm} to {subset_pp_high_ab_pmm} points. QRF and ZI-QRF change it by {subset_diff_qrf} and {subset_diff_zi_qrf} points (SE {subset_se_qrf} and {subset_se_zi_qrf}). Because our donors and recipients share one survey's questions, income measures and population, we do not test the differences between surveys that an ACS imputation also faces. + +For each method's change by age, disability, coverage type and resources relative to the threshold, see Table 1 in `TABLE.md`; for the code, see `README.md`. diff --git a/experiments/moop_spm/README.md b/experiments/moop_spm/README.md new file mode 100644 index 00000000..2ddac855 --- /dev/null +++ b/experiments/moop_spm/README.md @@ -0,0 +1,129 @@ +# MOOP imputation and SPM poverty + +How much does the choice of method for imputing medical out-of-pocket +spending (MOOP) change Supplemental Poverty Measure (SPM) poverty? We +measure it on the public-use file of the 2026 Current Population Survey +Annual Social and Economic Supplement (CPS ASEC), which records SPM-unit +MOOP for calendar year 2025. + +The American Community Survey (ACS) does not ask about medical spending, +so researchers who build an SPM on the ACS impute MOOP. For 2014, +Census researchers imputed it with predictive mean matching by insurance +status ([Fox, Pacas and Glassman, 2017](https://www.census.gov/content/dam/Census/library/working-papers/2017/demo/Fox-ACSSPM2014.pdf)). +We impute MOOP within the CPS ASEC and check each method against the +file's own MOOP. + +We summarize the findings in `NOTE.md` and give the full tables in +`TABLE.md`. `render.py` writes both from `results/results.json`. + +## Data + +- We use Census's CPS ASEC 2026 public-use CSV release + (`asecpub26csv.zip`; `cps_asec.py` pins its SHA-256) and the 160 person + replicate weights it ships. +- We impute `SPM_MEDXPNS` for each SPM unit. Census labels the variable + "MOOP and Medicare Part B subsidy" in its data dictionary. On every + record, the variable equals the unit sum of person `MOOP` plus the unit + sum of `PEMCPREM` (edited Medicare premium). +- `check_spm_identity` confirms that, on every record, `SPM_RESOURCES` + equals cash income plus SNAP, housing, school lunch, energy and WIC + subsidies, less taxes, FICA, capped work and child care expenses, + `SPM_MEDXPNS` and child support paid. So an imputed MOOP `m'` changes + a unit's resources to `SPM_RESOURCES + SPM_MEDXPNS - m'`. It also + confirms that `SPM_POOR` equals `SPM_RESOURCES < SPM_POVTHRESHOLD`, so + we recompute poverty from the new resources and the unchanged + threshold. +- With reported MOOP and the ASEC person weight, we reproduce Census's + published 2025 SPM rates (P60-290 Table 5) and their margins of error, + for all people and each age group, to within 0.01 percentage point. + Subtracting MOOP raises our rates by 0.02 to 0.03 points less than it + raises Census's. `run.py` writes the published figures beside ours + into `results.json`. +- Census imputes MOOP for 43% of the file's people (`I_MOOP` of 1, 2 or + 3) and the Medicare premium for 13% (`I_MCPREM` of 2). + +## Predictors + +For each SPM unit, we build 73 predictor columns +(`cps_asec.design_matrix`) from information an ACS-style recipient file +carries: + +- head's age and sex; counts of members aged 0-17, 18-44, 45-64, 65-74, + 75+; adult women; +- members with a disability (`PRDISFLG`); +- members with each type of current coverage (`NOW_*`: employer, direct + purchase, Medicare, Medicaid or other means-tested, TRICARE/VA/CHAMPVA, + uninsured); +- unit cash income (signed log) and cash income relative to the SPM + threshold; +- SPM family type, tenure and state (indicators). + +## Methods + +We fit every method on donor units, without survey weights, and impute +held-out recipient units (`methods.py`). + +| Name | Participation | Amount | +|---|---|---| +| (a) `logit_cell_mean` | logistic P(MOOP > 0) >= 0.5 | mean of donor spenders in the unit's cell | +| (a) `logit_reg_mean` | same | max(0, OLS prediction among donor spenders) | +| (b) `logit_draw_hotdeck` | Bernoulli draw from the logistic P | random donor spender in the cell | +| `pmm` | implicit | predictive mean matching on log(1 + MOOP) within five coverage groups, random one of 5 nearest donors | +| (c) `qrf` | implicit | microimpute `QRF`, one random quantile per unit | +| (c) `zi_qrf` | microimpute gate | microimpute `ZeroInflatedImputer` around `QRF` | + +We form cells from the unit's coverage group (Medicare, employer, direct +purchase, public, uninsured), the head's age band, cash income relative +to the threshold (under 1, 1-2, 2-4, 4+) and unit size (1, 2, 3+). Each +recipient uses the finest cell with at least 30 donor spenders. Where its +full cell has fewer, we drop unit size, then income, then age and then +coverage group. + +For PMM, we match within five coverage groups, where the Census +researchers used eight groups by insurance status. + +To decompose method (a), we cross the participation rule (0.5 cutoff or +a draw) with the amount rule (cell mean or a donor draw). In two more +variants, we give every unit `p x cell mean`, or give each unit at or +above the 0.5 cutoff its coverage group's mean (no age, income or size +cells). + +## Design + +- We assign whole households to 5 folds (`run.assign_folds`) and impute + each SPM unit from models we fit on the other 4. +- We take 10 draws for each stochastic method and average poverty rates + over them. +- We weight rates with the ASEC person weight (`MARSUPWT / 100`) and + compute standard errors from the 160 replicate weights, + `SE^2 = (4 / 160) * sum_r (theta_r - theta_0)^2`, treating the fitted + imputation models as fixed. +- We break the results down by age, disability, current coverage, and + the unit's resources before MOOP relative to its threshold. +- We repeat the comparison among people in SPM units where Census + imputes no member's MOOP or Medicare premium (`I_MOOP` of 0 or -1, and + `I_MCPREM` of 0 or -1). + +## Running + +Install microimpute, then run from the repository root: + +```bash +python -m experiments.moop_spm.run --data-dir ~/data/cps_asec_2026 +python -m experiments.moop_spm.render +``` + +The first command downloads the 139 MB Census zip, checks its SHA-256, +caches each fold's imputations in `results/cache/`, which git ignores, +and writes `results/results.json`. +`run.py` refits every fold when any of these change: + +- the settings or the experiment's source files; +- any microimpute source file; +- the version of Python, numpy, pandas, scipy, scikit-learn, + quantile-forest or microimpute; +- the model inputs or their row keys. + +`run.py` also records each method's run time in `results.json`. The +tests in `tests/test_moop_spm_experiment.py` check the experiment's +invariants on synthetic data. diff --git a/experiments/moop_spm/TABLE.md b/experiments/moop_spm/TABLE.md new file mode 100644 index 00000000..ce282bfd --- /dev/null +++ b/experiments/moop_spm/TABLE.md @@ -0,0 +1,63 @@ +# MOOP imputation and SPM poverty: tables + +Source: CPS ASEC 2026 public-use file (calendar-year 2025 income), 134,729 people in 55,401 SPM units and 53,344 households. We impute each unit's MOOP from models we fit on the other 4 of 5 folds of whole households and average 10 draws for each stochastic method. We weight rates by the ASEC person weight and compute the standard errors in parentheses from the 160 replicate weights, treating the fitted imputation models as fixed. + +## Table 1. SPM poverty rate with reported MOOP (%), and the change under each imputation (percentage points, SE) + +| Group | People (n) | Reported | (a) Logit + cell mean | (a) Logit + OLS mean | (b) Logit draw + donor | PMM | (c) QRF | (c) ZI-QRF | +|---|---:|---:|---:|---:|---:|---:|---:|---:| +| All people | 134,729 | 13.1 | +0.30 (0.08) | +0.34 (0.08) | +0.39 (0.07) | +0.32 (0.07) | -0.07 (0.07) | -0.05 (0.07) | +| **Age** | | | | | | | | | +| Under 18 | 31,444 | 13.4 | +0.17 (0.13) | +0.26 (0.14) | +0.42 (0.11) | +0.41 (0.11) | +0.00 (0.11) | +0.06 (0.11) | +| 18 to 64 | 77,083 | 12.3 | +0.13 (0.09) | +0.11 (0.09) | +0.26 (0.07) | +0.19 (0.07) | -0.08 (0.07) | -0.06 (0.07) | +| 65 and older | 26,202 | 15.4 | +0.98 (0.17) | +1.15 (0.16) | +0.75 (0.13) | +0.60 (0.15) | -0.10 (0.14) | -0.13 (0.13) | +| **Disability** | | | | | | | | | +| With a disability | 14,477 | 22.7 | +1.20 (0.22) | +1.40 (0.21) | +0.59 (0.18) | +0.41 (0.20) | -0.21 (0.18) | -0.22 (0.18) | +| No disability | 94,282 | 11.7 | +0.21 (0.09) | +0.20 (0.09) | +0.35 (0.07) | +0.28 (0.07) | -0.07 (0.07) | -0.04 (0.07) | +| Not asked (mostly under 15) | 25,970 | 13.2 | +0.15 (0.14) | +0.28 (0.15) | +0.42 (0.11) | +0.41 (0.12) | +0.02 (0.12) | +0.05 (0.12) | +| **Current health insurance** | | | | | | | | | +| Medicare | 27,328 | 16.6 | +1.20 (0.17) | +1.23 (0.16) | +0.85 (0.14) | +0.59 (0.15) | -0.08 (0.15) | -0.11 (0.13) | +| Employer | 66,520 | 3.9 | -0.21 (0.08) | -0.12 (0.08) | +0.14 (0.07) | +0.23 (0.07) | -0.02 (0.07) | +0.02 (0.08) | +| Direct purchase | 8,341 | 18.3 | +0.06 (0.39) | +2.28 (0.48) | +0.20 (0.32) | +0.71 (0.36) | -0.22 (0.36) | -0.04 (0.36) | +| Medicaid or other means-tested | 18,134 | 32.1 | +1.18 (0.23) | +0.32 (0.22) | +0.99 (0.18) | +0.33 (0.20) | +0.01 (0.17) | -0.03 (0.18) | +| TRICARE, VA or CHAMPVA | 3,174 | 7.2 | -0.17 (0.21) | -0.59 (0.36) | -0.29 (0.29) | -0.15 (0.29) | -0.01 (0.26) | -0.09 (0.30) | +| Uninsured | 11,232 | 25.8 | +0.07 (0.27) | -0.22 (0.28) | +0.06 (0.23) | -0.05 (0.22) | -0.33 (0.23) | -0.28 (0.24) | +| **Resources before MOOP / threshold** | | | | | | | | | +| Under 1.0 | 13,803 | 100.0 | +0.00 (0.00) | +0.00 (0.00) | +0.00 (0.00) | +0.00 (0.00) | +0.00 (0.00) | +0.00 (0.00) | +| 1.0 to 1.25 | 8,441 | 24.9 | +14.69 (1.01) | +14.37 (1.05) | +4.88 (0.74) | +4.19 (0.75) | +0.55 (0.79) | +0.61 (0.81) | +| 1.25 to 1.5 | 9,233 | 5.8 | -4.75 (0.50) | -3.92 (0.49) | +0.50 (0.49) | +0.23 (0.51) | -0.76 (0.48) | -0.57 (0.50) | +| 1.5 to 2.0 | 17,495 | 1.6 | -1.56 (0.15) | -1.54 (0.15) | +0.19 (0.17) | +0.17 (0.18) | -0.29 (0.17) | -0.26 (0.16) | +| 2.0 to 3.0 | 30,010 | 0.3 | -0.32 (0.07) | -0.32 (0.07) | +0.05 (0.07) | +0.03 (0.08) | -0.10 (0.07) | -0.06 (0.07) | +| 3.0 or more | 55,747 | 0.1 | -0.05 (0.02) | -0.05 (0.02) | +0.03 (0.02) | +0.02 (0.02) | +0.02 (0.02) | +0.00 (0.02) | + +Every method imputes MOOP of zero or more, so people whose unit's resources before MOOP fall below its threshold stay poor, and the 'Under 1.0' row cannot change. + +## Table 2. MOOP distribution across SPM units (unit weights) and gross flows of people (person weights) + +| Statistic | Reported | (a) Logit + cell mean | (a) Logit + OLS mean | (b) Logit draw + donor | PMM | (c) QRF | (c) ZI-QRF | +|---|---:|---:|---:|---:|---:|---:|---:| +| Mean | $6,011 | $6,156 | $6,133 | $5,984 | $5,922 | $5,797 | $5,791 | +| Median | $3,700 | $6,132 | $5,891 | $3,678 | $3,624 | $3,623 | $3,610 | +| Share with zero MOOP | 6.4% | 0.5% | 4.3% | 6.5% | 6.2% | 6.2% | 6.1% | +| 75th percentile | $8,160 | $8,625 | $8,983 | $8,154 | $8,091 | $7,949 | $7,891 | +| 90th percentile | $14,300 | $11,436 | $11,507 | $14,276 | $14,158 | $13,820 | $13,765 | +| 99th percentile | $33,640 | $15,585 | $15,769 | $33,147 | $32,666 | $31,356 | $31,424 | +| Total ($ billions) | 859 | 880 | 876 | 855 | 846 | 828 | 828 | +| Mean, resources before MOOP 1.0 to 1.25 x threshold | $2,385 | $3,189 | $3,108 | $2,949 | $2,799 | $2,352 | $2,372 | +| Mean, resources before MOOP 1.25 to 1.5 x threshold | $3,381 | $3,679 | $3,667 | $3,446 | $3,452 | $3,282 | $3,272 | +| Mean, resources before MOOP 1.5 to 2.0 x threshold | $4,510 | $4,788 | $4,612 | $4,546 | $4,527 | $4,335 | $4,389 | +| Mean, resources before MOOP 2.0 to 3.0 x threshold | $6,086 | $6,116 | $5,872 | $5,953 | $5,783 | $5,809 | $5,843 | +| People newly poor (%) | 0.0 | 1.37 | 1.40 | 1.74 | 1.63 | 1.33 | 1.37 | +| People newly not poor (%) | 0.0 | 1.07 | 1.06 | 1.36 | 1.31 | 1.40 | 1.41 | +| Change in rate among people in units where Census imputes no MOOP or Medicare premium (pp) | 0.00 | +0.36 (0.11) | +0.26 (0.11) | +0.52 (0.09) | +0.48 (0.10) | +0.01 (0.09) | +0.07 (0.10) | + +## Table 3. How the participation and amount rules of the logit approach move the rate (all people) + +| Participation rule | Amount rule | Rate (%) | Change (pp) | SE | +|---|---|---:|---:|---:| +| p >= 0.5 | Cell mean | 13.42 | +0.30 | 0.08 | +| p >= 0.5 | Random donor in cell | 13.70 | +0.59 | 0.07 | +| Bernoulli(p) | Cell mean | 13.21 | +0.10 | 0.08 | +| Bernoulli(p) | Random donor in cell | 13.50 | +0.39 | 0.07 | +| None (all units) | p x cell mean | 13.19 | +0.07 | 0.08 | +| p >= 0.5 | Coverage-group mean (5 cells) | 15.18 | +2.07 | 0.10 | diff --git a/experiments/moop_spm/__init__.py b/experiments/moop_spm/__init__.py new file mode 100644 index 00000000..c1f4b6cc --- /dev/null +++ b/experiments/moop_spm/__init__.py @@ -0,0 +1,4 @@ +"""How MOOP imputation methods change SPM poverty (CPS ASEC 2026). + +Run from the repository root with ``python -m experiments.moop_spm.run``. +""" diff --git a/experiments/moop_spm/cps_asec.py b/experiments/moop_spm/cps_asec.py new file mode 100644 index 00000000..2155a5f9 --- /dev/null +++ b/experiments/moop_spm/cps_asec.py @@ -0,0 +1,441 @@ +"""Load the 2026 CPS ASEC public-use file into person and SPM-unit frames. + +The 2026 CPS ASEC (fielded February to April 2026) reports calendar-year +2025 income and carries Census's Supplemental Poverty Measure (SPM) +variables for every person. This module downloads the public-use CSV +release, checks its SHA-256, and builds: + +- a person frame with SPM resources, thresholds, poverty status, the + person weight and Census's 160 replicate weights; and +- an SPM-unit frame with reported medical out-of-pocket spending (MOOP, + ``SPM_MEDXPNS``) and the predictors an ACS-style recipient file carries. + +Variable names and codes follow the Census data dictionary +``asec2026_ddl_pub_full.pdf`` shipped with the release. +""" + +from __future__ import annotations + +import hashlib +import urllib.request +import zipfile +import zlib +from pathlib import Path + +import numpy as np +import pandas as pd + +SURVEY_YEAR = 2026 +REFERENCE_YEAR = 2025 +ZIP_URL = ( + "https://www2.census.gov/programs-surveys/cps/datasets/2026/march/asecpub26csv.zip" +) +ZIP_SHA256 = "fe819d0d2fc4470c282e76c247b8aa8b6054811619730307e9f2fe6186707d93" +PERSON_CSV = "pppub26.csv" +HOUSEHOLD_CSV = "hhpub26.csv" +REPWGT_CSV = "asec_csv_repwgt_2026.csv" +N_REPLICATES = 160 + +PERSON_COLUMNS = [ + "PH_SEQ", + "PPPOS", + "A_AGE", + "A_SEX", + "PRDISFLG", + "MARSUPWT", + "I_MOOP", + "I_MCPREM", + "MOOP", + "PEMCPREM", + "NOW_COV", + "NOW_GRP", + "NOW_DIR", + "NOW_MCARE", + "NOW_MCAID", + "NOW_MIL", + "NOW_VACARE", + "NOW_CHAMPVA", + "SPM_ID", + "SPM_HEAD", + "SPM_WEIGHT", + "SPM_RESOURCES", + "SPM_POVTHRESHOLD", + "SPM_POOR", + "SPM_MEDXPNS", + "SPM_TOTVAL", + "SPM_SNAPSUB", + "SPM_CAPHOUSESUB", + "SPM_SCHLUNCH", + "SPM_ENGVAL", + "SPM_WICVAL", + "SPM_FEDTAX", + "SPM_STTAX", + "SPM_FICA", + "SPM_CAPWKCCXPNS", + "SPM_CHILDSUPPD", + "SPM_NUMPER", + "SPM_FAMTYPE", + "SPM_TENMORTSTATUS", +] + +# Current coverage at interview (NOW_*), the closest analogue of the +# ACS coverage question. Codes: 1 = yes, 2 = no. +COVERAGE_ORDER = [ + "medicare", + "employer", + "direct", + "medicaid", + "military", + "uninsured", +] +COVERAGE_LABELS = { + "medicare": "Medicare", + "employer": "Employer", + "direct": "Direct purchase", + "medicaid": "Medicaid or other means-tested", + "military": "TRICARE, VA or CHAMPVA", + "uninsured": "Uninsured", +} +# Coverage groups for SPM units, used to define imputation cells and +# PMM groups. TRICARE/VA units are folded into the public group. +UNIT_GROUPS = ["medicare", "employer", "direct", "public", "uninsured"] + +AGE_BANDS = [(0, 17), (18, 44), (45, 64), (65, 74), (75, 200)] + + +def sha256_of(path: Path) -> str: + """Return the hex SHA-256 digest of a file.""" + digest = hashlib.sha256() + with open(path, "rb") as handle: + for chunk in iter(lambda: handle.read(1 << 20), b""): + digest.update(chunk) + return digest.hexdigest() + + +def crc32_of(path: Path) -> int: + """Return the CRC-32 of a file, as a zip archive records it.""" + crc = 0 + with open(path, "rb") as handle: + for chunk in iter(lambda: handle.read(1 << 20), b""): + crc = zlib.crc32(chunk, crc) + return crc + + +def extract_verified(zip_path: Path, data_dir: Path, names: list[str]) -> None: + """Extract ``names`` from a zip unless matching copies already exist. + + A copy matches when its size and CRC-32 equal the archive member's; + any other copy is replaced from the archive. + """ + with zipfile.ZipFile(zip_path) as archive: + for name in names: + info = archive.getinfo(name) + target = data_dir / name + if ( + target.exists() + and target.stat().st_size == info.file_size + and crc32_of(target) == info.CRC + ): + continue + archive.extract(name, data_dir) + + +def ensure_raw(data_dir: Path) -> Path: + """Download and unpack the ASEC CSV release into ``data_dir``. + + The zip is fetched once and checked against the pinned SHA-256. The + three CSVs this experiment reads are extracted next to it, and copies + that differ from the zip's members are extracted again. + + Args: + data_dir: Directory that holds (or will hold) the raw files. + + Returns: + The data directory. + + Raises: + ValueError: If the zip's digest does not match ``ZIP_SHA256``. + """ + data_dir = Path(data_dir).expanduser() + data_dir.mkdir(parents=True, exist_ok=True) + zip_path = data_dir / "asecpub26csv.zip" + if not zip_path.exists(): + urllib.request.urlretrieve(ZIP_URL, zip_path) + digest = sha256_of(zip_path) + if digest != ZIP_SHA256: + raise ValueError( + f"{zip_path} has SHA-256 {digest}, expected {ZIP_SHA256}. " + "Census may have replaced the file; check before rerunning." + ) + extract_verified(zip_path, data_dir, [PERSON_CSV, HOUSEHOLD_CSV, REPWGT_CSV]) + return data_dir + + +def coverage_category(persons: pd.DataFrame) -> pd.Series: + """Assign each person one current-coverage category. + + The hierarchy is Medicare, employer, direct purchase, Medicaid or + other means-tested, TRICARE/VA/CHAMPVA, then uninsured. A person with + Medicare and Medicaid counts as Medicare. + """ + military = ( + (persons["NOW_MIL"] == 1) + | (persons["NOW_VACARE"] == 1) + | (persons["NOW_CHAMPVA"] == 1) + ) + conditions = [ + persons["NOW_MCARE"] == 1, + persons["NOW_GRP"] == 1, + persons["NOW_DIR"] == 1, + persons["NOW_MCAID"] == 1, + military, + persons["NOW_COV"] == 2, + ] + category = np.select(conditions, COVERAGE_ORDER, default="unclassified") + return pd.Series(category, index=persons.index, name="coverage") + + +def load_persons(data_dir: Path) -> pd.DataFrame: + """Read person records, attach state and replicate weights. + + Returns one row per ASEC person record with the columns in + ``PERSON_COLUMNS`` plus ``state`` (GESTFIPS), ``weight`` (MARSUPWT + divided by 100, the person weight in people), ``coverage`` and + ``repwgt_1`` to ``repwgt_160``. + """ + data_dir = ensure_raw(data_dir) + persons = pd.read_csv(data_dir / PERSON_CSV, usecols=PERSON_COLUMNS) + households = pd.read_csv( + data_dir / HOUSEHOLD_CSV, usecols=["H_SEQ", "GESTFIPS"] + ).rename(columns={"H_SEQ": "PH_SEQ", "GESTFIPS": "state"}) + persons = persons.merge(households, on="PH_SEQ", how="left", validate="many_to_one") + if persons["state"].isna().any(): + raise ValueError("Some persons have no household record.") + + replicate_columns = [f"pwwgt{i}" for i in range(1, N_REPLICATES + 1)] + replicates = pd.read_csv( + data_dir / REPWGT_CSV, + usecols=["PPPOS", "h_seq", *replicate_columns], + ).rename(columns={"h_seq": "PH_SEQ"}) + replicates = replicates.rename( + columns={f"pwwgt{i}": f"repwgt_{i}" for i in range(1, N_REPLICATES + 1)} + ) + persons = persons.merge( + replicates, on=["PH_SEQ", "PPPOS"], how="left", validate="one_to_one" + ).copy() + if persons["repwgt_1"].isna().any(): + raise ValueError("Some persons have no replicate weights.") + + persons["weight"] = persons["MARSUPWT"] / 100.0 + persons["coverage"] = coverage_category(persons) + if (persons["coverage"] == "unclassified").any(): + raise ValueError("Coverage hierarchy left persons unclassified.") + return persons.reset_index(drop=True) + + +def check_spm_identity(persons: pd.DataFrame) -> dict[str, float]: + """Check how SPM resources and MOOP are assembled in the file. + + Returns three checks that must come out zero for the imputation + exercise to be valid: + + - ``resources_gap``: SPM_RESOURCES against cash income plus noncash + benefits less taxes, capped work and child care expenses, MOOP and + child support paid. + - ``poor_mismatch``: share of persons whose SPM_POOR differs from + ``SPM_RESOURCES < SPM_POVTHRESHOLD``. + - ``moop_gap``: SPM_MEDXPNS against the unit sum of person MOOP plus + the unit sum of PEMCPREM (edited Medicare premium). + """ + rebuilt = ( + persons["SPM_TOTVAL"] + + persons["SPM_SNAPSUB"] + + persons["SPM_CAPHOUSESUB"] + + persons["SPM_SCHLUNCH"] + + persons["SPM_ENGVAL"] + + persons["SPM_WICVAL"] + - persons["SPM_FEDTAX"] + - persons["SPM_STTAX"] + - persons["SPM_FICA"] + - persons["SPM_CAPWKCCXPNS"] + - persons["SPM_MEDXPNS"] + - persons["SPM_CHILDSUPPD"] + ) + poor = (persons["SPM_RESOURCES"] < persons["SPM_POVTHRESHOLD"]).astype(int) + by_unit = persons.groupby("SPM_ID").agg( + moop=("MOOP", "sum"), + premium=("PEMCPREM", "sum"), + reported=("SPM_MEDXPNS", "first"), + ) + return { + "resources_gap": float((rebuilt - persons["SPM_RESOURCES"]).abs().max()), + "poor_mismatch": float((poor != persons["SPM_POOR"]).mean()), + "moop_gap": float( + (by_unit["reported"] - by_unit["moop"] - by_unit["premium"]).abs().max() + ), + } + + +def allocation_flags(persons: pd.DataFrame) -> pd.DataFrame: + """Flag persons whose MOOP or Medicare premium Census imputed. + + ``SPM_MEDXPNS`` sums members' ``MOOP`` and ``PEMCPREM``. ``I_MOOP`` + codes 1 to 3 mark hot-deck, logical and whole-unit imputation of + ``MOOP``; ``I_MCPREM`` code 2 marks logical imputation of + ``PEMCPREM``. + """ + return pd.DataFrame( + { + "allocated_moop": persons["I_MOOP"].isin([1, 2, 3]).astype(int), + "allocated_premium": (persons["I_MCPREM"] == 2).astype(int), + }, + index=persons.index, + ) + + +def _signed_log1p(values: pd.Series) -> pd.Series: + return np.sign(values) * np.log1p(values.abs()) + + +def build_units(persons: pd.DataFrame) -> pd.DataFrame: + """Collapse persons to one row per SPM unit. + + The frame holds the imputation target ``moop`` (SPM_MEDXPNS), the + unit's SPM resources and threshold, its weight (``SPM_WEIGHT`` / 100, + which equals the head's person weight), and ACS-style predictors: + age, sex, disability, current coverage, cash income, family size and + type, tenure and state. + """ + persons = persons.loc[:, ~persons.columns.str.startswith("repwgt_")].copy() + persons = persons.join(allocation_flags(persons)) + persons["disabled"] = (persons["PRDISFLG"] == 1).astype(int) + persons["female_adult"] = ( + (persons["A_SEX"] == 2) & (persons["A_AGE"] >= 18) + ).astype(int) + for low, high in AGE_BANDS: + name = f"n_age_{low}_{high}" if high < 200 else f"n_age_{low}_plus" + persons[name] = persons["A_AGE"].between(low, high).astype(int) + military = ( + (persons["NOW_MIL"] == 1) + | (persons["NOW_VACARE"] == 1) + | (persons["NOW_CHAMPVA"] == 1) + ) + # Coverage counts are not mutually exclusive: a person with two + # coverage types counts toward both. + any_coverage = { + "n_any_medicare": persons["NOW_MCARE"] == 1, + "n_any_employer": persons["NOW_GRP"] == 1, + "n_any_direct": persons["NOW_DIR"] == 1, + "n_any_medicaid": persons["NOW_MCAID"] == 1, + "n_any_military": military, + "n_uninsured": persons["NOW_COV"] == 2, + } + for name, flag in any_coverage.items(): + persons[name] = flag.astype(int) + + count_columns = [ + c for c in persons.columns if c.startswith("n_age_") or c in any_coverage + ] + count_columns += [ + "disabled", + "female_adult", + "allocated_moop", + "allocated_premium", + ] + sums = persons.groupby("SPM_ID")[count_columns].sum() + sums = sums.rename( + columns={ + "disabled": "n_disabled", + "female_adult": "n_female_adult", + "allocated_moop": "n_allocated_moop", + "allocated_premium": "n_allocated_premium", + } + ) + heads = persons.loc[persons["SPM_HEAD"] == 1].set_index("SPM_ID") + if not heads.index.is_unique or len(heads) != persons["SPM_ID"].nunique(): + raise ValueError("Expected exactly one SPM head per SPM unit.") + units = pd.DataFrame( + { + "household": heads["PH_SEQ"], + "unit_weight": heads["SPM_WEIGHT"] / 100.0, + "moop": heads["SPM_MEDXPNS"].astype(float), + "resources": heads["SPM_RESOURCES"].astype(float), + "threshold": heads["SPM_POVTHRESHOLD"].astype(float), + "cash_income": heads["SPM_TOTVAL"].astype(float), + "n_persons": heads["SPM_NUMPER"].astype(int), + "famtype": heads["SPM_FAMTYPE"].astype(int), + "tenure": heads["SPM_TENMORTSTATUS"].astype(int), + "state": heads["state"].astype(int), + "head_age": heads["A_AGE"].astype(float), + "head_female": (heads["A_SEX"] == 2).astype(int), + } + ).join(sums) + units["pre_moop_ratio"] = (units["resources"] + units["moop"]) / units["threshold"] + units["log_cash_income"] = _signed_log1p(units["cash_income"]) + units["cash_income_ratio"] = (units["cash_income"] / units["threshold"]).clip( + -1.0, 10.0 + ) + units["group"] = unit_coverage_group(units) + # Units where Census imputed any part of SPM_MEDXPNS for any member. + units["any_allocated"] = ( + units["n_allocated_moop"] + units["n_allocated_premium"] > 0 + ).astype(int) + return units.sort_index() + + +def unit_coverage_group(units: pd.DataFrame) -> pd.Series: + """Classify SPM units by the coverage of their members. + + Any member with Medicare makes the unit ``medicare``; otherwise any + employer coverage makes it ``employer``; then direct purchase; then + Medicaid, other means-tested, TRICARE or VA (``public``); units with + none of these are ``uninsured``. + """ + conditions = [ + units["n_any_medicare"] > 0, + units["n_any_employer"] > 0, + units["n_any_direct"] > 0, + (units["n_any_medicaid"] > 0) | (units["n_any_military"] > 0), + ] + group = np.select(conditions, UNIT_GROUPS[:4], default="uninsured") + return pd.Series(group, index=units.index, name="group") + + +NUMERIC_FEATURES: list[str] = [ + "head_age", + "head_female", + "n_age_0_17", + "n_age_18_44", + "n_age_45_64", + "n_age_65_74", + "n_age_75_plus", + "n_female_adult", + "n_disabled", + "n_any_medicare", + "n_any_employer", + "n_any_direct", + "n_any_medicaid", + "n_any_military", + "n_uninsured", + "log_cash_income", + "cash_income_ratio", +] +CATEGORICAL_FEATURES: list[str] = ["famtype", "tenure", "state"] + + +def design_matrix(units: pd.DataFrame) -> pd.DataFrame: + """Return the predictor matrix shared by every imputation method. + + Numeric predictors enter as-is; family type, tenure and state enter + as 0/1 indicators with the lowest code dropped. Every method sees + exactly these columns. + """ + frames = [units[NUMERIC_FEATURES].astype(float)] + for column in CATEGORICAL_FEATURES: + dummies = pd.get_dummies( + units[column].astype(int), prefix=column, drop_first=True + ).astype(float) + frames.append(dummies) + matrix = pd.concat(frames, axis=1) + matrix.columns = [str(c) for c in matrix.columns] + return matrix diff --git a/experiments/moop_spm/evaluate.py b/experiments/moop_spm/evaluate.py new file mode 100644 index 00000000..b3241334 --- /dev/null +++ b/experiments/moop_spm/evaluate.py @@ -0,0 +1,281 @@ +"""SPM poverty and MOOP distribution metrics with replicate-weight SEs. + +Replacing a unit's reported MOOP ``m`` with an imputed value ``m'`` +changes its SPM resources to ``SPM_RESOURCES + m - m'``; a person is in +SPM poverty when the unit's resources fall below ``SPM_POVTHRESHOLD``. +``check_spm_identity`` in ``cps_asec`` confirms both relations hold +exactly in the public-use file. + +Standard errors use the 160 ASEC replicate weights with Census's +formula SE^2 = (4 / 160) * sum_r (theta_r - theta_0)^2. They treat the +fitted imputation models as fixed: they capture sampling variation in +the evaluation sample, not refitting variation. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from .cps_asec import COVERAGE_LABELS, COVERAGE_ORDER, N_REPLICATES + +REPLICATE_FACTOR = 4.0 / N_REPLICATES +RATIO_BANDS: list[tuple[float, float, str]] = [ + (-np.inf, 1.0, "Under 1.0"), + (1.0, 1.25, "1.0 to 1.25"), + (1.25, 1.5, "1.25 to 1.5"), + (1.5, 2.0, "1.5 to 2.0"), + (2.0, 3.0, "2.0 to 3.0"), + (3.0, np.inf, "3.0 or more"), +] +QUANTILES = {"median": 0.5, "p75": 0.75, "p90": 0.9, "p95": 0.95, "p99": 0.99} + + +def weight_matrix(persons: pd.DataFrame) -> np.ndarray: + """Columns: full-sample person weight, then the 160 replicates.""" + columns = ["weight"] + [f"repwgt_{i}" for i in range(1, N_REPLICATES + 1)] + return persons[columns].to_numpy(dtype=float) + + +def replicate_se(estimates: np.ndarray) -> np.ndarray: + """SE from an array whose last axis is (full sample, 160 replicates).""" + deviations = estimates[..., 1:] - estimates[..., :1] + return np.sqrt(REPLICATE_FACTOR * (deviations**2).sum(axis=-1)) + + +def person_groups(persons: pd.DataFrame, units: pd.DataFrame) -> pd.DataFrame: + """Boolean membership of each person in each reporting group. + + Groups: everyone; age under 18, 18 to 64, 65 or older; disability + (PRDISFLG, asked of civilians 15 and older); current coverage + category; and the unit's resources before MOOP relative to its + threshold. + """ + age = persons["A_AGE"] + ratio = persons["SPM_ID"].map(units["pre_moop_ratio"]) + groups = { + ("all", "All people"): np.ones(len(persons), dtype=bool), + ("age", "Under 18"): age < 18, + ("age", "18 to 64"): (age >= 18) & (age < 65), + ("age", "65 and older"): age >= 65, + ("disability", "With a disability"): persons["PRDISFLG"] == 1, + ("disability", "No disability"): persons["PRDISFLG"] == 2, + ("disability", "Not asked (mostly under 15)"): persons["PRDISFLG"] == -1, + } + for category in COVERAGE_ORDER: + groups[("coverage", COVERAGE_LABELS[category])] = ( + persons["coverage"] == category + ) + for low, high, label in RATIO_BANDS: + groups[("pre_moop_ratio", label)] = (ratio >= low) & (ratio < high) + frame = pd.DataFrame( + {f"{dim}::{label}": np.asarray(mask) for (dim, label), mask in groups.items()} + ) + return frame + + +def poverty_flags( + persons: pd.DataFrame, units: pd.DataFrame, unit_moop: np.ndarray +) -> np.ndarray: + """Person-level SPM poverty under imputed unit MOOP. + + Args: + persons: Person frame with SPM_ID, SPM_RESOURCES, SPM_MEDXPNS and + SPM_POVTHRESHOLD. + units: Unit frame indexed by SPM_ID, in the order of ``unit_moop``'s + last axis. + unit_moop: Imputed MOOP, shape ``(n_draws, n_units)`` or + ``(n_units,)``. + + Returns: + Boolean poverty flags of shape ``(n_draws, n_persons)``. + """ + unit_moop = np.atleast_2d(unit_moop) + position = pd.Index(units.index).get_indexer(persons["SPM_ID"]) + if (position < 0).any(): + raise ValueError("Persons reference SPM units missing from units.") + pre_moop = (persons["SPM_RESOURCES"] + persons["SPM_MEDXPNS"]).to_numpy(dtype=float) + threshold = persons["SPM_POVTHRESHOLD"].to_numpy(dtype=float) + resources = pre_moop[None, :] - unit_moop[:, position] + return resources < threshold[None, :] + + +def group_rates( + flags: np.ndarray, groups: np.ndarray, weights: np.ndarray +) -> np.ndarray: + """Weighted share flagged within each group, for every weight column. + + Args: + flags: Person values, shape ``(n_persons,)`` (may be fractional, + e.g. the share of draws in which a person is poor). + groups: Membership, shape ``(n_persons, n_groups)``. + weights: Weights, shape ``(n_persons, n_weights)``. + + Returns: + Array of shape ``(n_groups, n_weights)``. + """ + member = groups.astype(float) + numerator = (member * flags[:, None]).T @ weights + denominator = member.T @ weights + # An empty group has no rate; return NaN without a warning. + with np.errstate(invalid="ignore", divide="ignore"): + return numerator / denominator + + +def weighted_quantile(values: np.ndarray, weights: np.ndarray, q: float) -> float: + """Smallest value whose weighted cumulative share reaches ``q``.""" + order = np.argsort(values, kind="stable") + cumulative = np.cumsum(weights[order]) + index = np.searchsorted(cumulative, q * cumulative[-1], side="left") + return float(values[order][min(index, len(values) - 1)]) + + +def _mean_or_nan(values: np.ndarray, weights: np.ndarray) -> float: + """Weighted mean, or NaN for an empty selection.""" + if len(values) == 0: + return float("nan") + return float(np.average(values, weights=weights)) + + +def ratio_band_index(ratio: np.ndarray) -> np.ndarray: + """Index into ``RATIO_BANDS`` for each resources-to-threshold ratio.""" + edges = np.array([high for _, high, _ in RATIO_BANDS[:-1]]) + return np.searchsorted(edges, ratio, side="right") + + +def moop_distribution( + moop: np.ndarray, + weights: np.ndarray, + near_line: np.ndarray, + band: np.ndarray, +) -> dict[str, float]: + """Weighted SPM-unit MOOP statistics for one set of values. + + ``band_means`` gives mean MOOP within each band of resources before + MOOP relative to the threshold (``RATIO_BANDS`` order). + """ + stats = { + "band_means": [ + _mean_or_nan(moop[band == b], weights[band == b]) + for b in range(len(RATIO_BANDS)) + ], + "mean": float(np.average(moop, weights=weights)), + "share_zero": float(np.average(moop == 0, weights=weights)), + "aggregate_billions": float((moop * weights).sum() / 1e9), + "near_line_mean": _mean_or_nan(moop[near_line], weights[near_line]), + } + for name, q in QUANTILES.items(): + stats[name] = weighted_quantile(moop, weights, q) + return stats + + +def average_distribution( + draws: np.ndarray, + weights: np.ndarray, + near_line: np.ndarray, + band: np.ndarray, +) -> dict[str, float]: + """Average the distribution statistics over draws.""" + per_draw = [moop_distribution(d, weights, near_line, band) for d in draws] + out: dict = { + k: float(np.mean([s[k] for s in per_draw])) + for k in per_draw[0] + if k != "band_means" + } + out["band_means"] = np.mean([s["band_means"] for s in per_draw], axis=0).tolist() + if len(per_draw) > 1: + out["mean_draw_sd"] = float(np.std([s["mean"] for s in per_draw], ddof=1)) + return out + + +def evaluate( + persons: pd.DataFrame, + units: pd.DataFrame, + imputations: dict[str, np.ndarray], +) -> dict: + """Compute every reported metric for reported and imputed MOOP. + + Args: + persons: Person frame from ``load_persons``. + units: Unit frame from ``build_units``. + imputations: Method name to array ``(n_draws, n_units)``. + + Returns: + A JSON-serializable dict of results. + """ + weights = weight_matrix(persons) + groups = person_groups(persons, units) + group_matrix = groups.to_numpy() + reported = poverty_flags(persons, units, units["moop"].to_numpy())[0] + published_check = (reported == (persons["SPM_POOR"] == 1).to_numpy()).all() + if not published_check: + raise ValueError("Recomputed poverty does not match SPM_POOR.") + reported_f = reported.astype(float) + reported_rates = group_rates(reported_f, group_matrix, weights) + + unit_weight = units["unit_weight"].to_numpy() + near_line = (units["pre_moop_ratio"] >= 1.0) & (units["pre_moop_ratio"] < 2.0) + near_line = near_line.to_numpy() + band = ratio_band_index(units["pre_moop_ratio"].to_numpy()) + no_allocation_person = persons["SPM_ID"].map(units["any_allocated"]).to_numpy() == 0 + subset_groups = np.column_stack( + [np.ones(len(persons), dtype=bool), no_allocation_person] + ) + reported_subset = group_rates(reported_f, subset_groups, weights) + + results: dict = { + "groups": [c.split("::") for c in groups.columns], + "group_person_counts": groups.sum().astype(int).tolist(), + "methods": {}, + } + + def pack(rates: np.ndarray, diff: np.ndarray) -> dict[str, list[float]]: + return { + "rate": rates[:, 0].tolist(), + "rate_se": replicate_se(rates).tolist(), + "diff": diff[:, 0].tolist(), + "diff_se": replicate_se(diff).tolist(), + } + + results["methods"]["reported"] = { + **pack(reported_rates, np.zeros_like(reported_rates)), + "distribution": moop_distribution( + units["moop"].to_numpy(), unit_weight, near_line, band + ), + "n_draws": 1, + } + w0 = weights[:, 0] + for method, draws in imputations.items(): + flags = poverty_flags(persons, units, draws) + mean_flag = flags.mean(axis=0) + rates = group_rates(mean_flag, group_matrix, weights) + entry = pack(rates, rates - reported_rates) + # Gross flows relative to reported MOOP, averaged over draws. + newly_poor = (flags & ~reported[None, :]).mean(axis=0) + newly_nonpoor = (~flags & reported[None, :]).mean(axis=0) + entry["newly_poor_share"] = float(np.average(newly_poor, weights=w0)) + entry["newly_nonpoor_share"] = float(np.average(newly_nonpoor, weights=w0)) + per_draw_overall = [ + float(np.average(f.astype(float), weights=w0)) for f in flags + ] + entry["overall_draw_sd"] = ( + float(np.std(per_draw_overall, ddof=1)) if len(flags) > 1 else 0.0 + ) + entry["distribution"] = average_distribution( + draws, unit_weight, near_line, band + ) + subset_rates = group_rates(mean_flag, subset_groups, weights) + subset_diff = subset_rates - reported_subset + entry["no_allocation_subset"] = { + "rate": float(subset_rates[1, 0]), + "diff": float(subset_diff[1, 0]), + "diff_se": float(replicate_se(subset_diff)[1]), + } + entry["n_draws"] = len(draws) + results["methods"][method] = entry + results["methods"]["reported"]["no_allocation_subset"] = { + "rate": float(reported_subset[1, 0]), + "diff": 0.0, + "diff_se": 0.0, + } + return results diff --git a/experiments/moop_spm/methods.py b/experiments/moop_spm/methods.py new file mode 100644 index 00000000..90af7811 --- /dev/null +++ b/experiments/moop_spm/methods.py @@ -0,0 +1,398 @@ +"""MOOP imputation methods compared on held-out CPS ASEC SPM units. + +Every method takes the same donor and recipient predictor matrices and +returns an array of shape ``(n_draws, n_recipients)`` of imputed SPM-unit +MOOP. Deterministic methods repeat one row. + +Methods: + +- ``logit_cell_mean``: a logistic regression for any MOOP; recipients + with predicted probability of at least 0.5 are "predicted spenders" + and get the mean MOOP of donor spenders in their cell. +- ``logit_reg_mean``: the same participation rule; predicted spenders + get the OLS prediction (in dollars) from a regression on donor + spenders. +- ``logit_draw_hotdeck``: participation drawn from the logistic + probability; spenders get the MOOP of a random donor spender in their + cell. +- ``pmm``: predictive mean matching within coverage groups: an OLS of + log(1 + MOOP) on the predictors, then each recipient takes the + reported MOOP of one of its five nearest donors by predicted value. +- ``qrf``: microimpute's ``QRF`` fit on all donors, one random quantile + draw per recipient. +- ``zi_qrf``: microimpute's ``ZeroInflatedImputer`` around ``QRF``. + +Decomposition variants cross the participation rule (threshold or +draw) with the amount rule (cell mean or donor draw), and +``logit_expected_value`` assigns every recipient the probability times +the cell mean, and ``logit_group_mean`` uses the 0.5 rule with means by +coverage group only (five cells). +""" + +from __future__ import annotations + +from collections.abc import Callable, Sequence +from dataclasses import dataclass + +import numpy as np +import pandas as pd +from sklearn.linear_model import LinearRegression, LogisticRegression +from sklearn.pipeline import Pipeline, make_pipeline +from sklearn.preprocessing import StandardScaler + +from microimpute.models.qrf import QRF +from microimpute.models.zero_inflated import ZeroInflatedImputer + +MIN_CELL_DONORS = 30 +PMM_NEIGHBORS = 5 +PARTICIPATION_CUTOFF = 0.5 + +# Imputation cells, finest first. A recipient uses the finest level +# whose cell holds at least MIN_CELL_DONORS donor spenders. +CELL_LEVELS: list[tuple[str, ...]] = [ + ("group", "age_band", "income_band", "size_band"), + ("group", "age_band", "income_band"), + ("group", "age_band"), + ("group",), + (), +] +# Coarse cells for the sensitivity check: coverage group only. +COARSE_CELL_LEVELS: list[tuple[str, ...]] = [("group",), ()] + + +def cell_keys(units: pd.DataFrame) -> pd.DataFrame: + """Return the cell variables for each SPM unit. + + - ``group``: coverage group of the unit (see ``unit_coverage_group``). + - ``age_band``: head's age under 45, 45 to 64, 65 or older. + - ``income_band``: cash income relative to the SPM threshold, under 1, + 1 to 2, 2 to 4, 4 or more. + - ``size_band``: unit size 1, 2, or 3 or more. + """ + age_band = pd.cut( + units["head_age"], [-1, 44, 64, 200], labels=["u45", "45_64", "65p"] + ) + income_band = pd.cut( + units["cash_income"] / units["threshold"], + [-np.inf, 1.0, 2.0, 4.0, np.inf], + right=False, + labels=["lt1", "1_2", "2_4", "4p"], + ) + size_band = pd.cut(units["n_persons"], [0, 1, 2, np.inf], labels=["1", "2", "3p"]) + return pd.DataFrame( + { + "group": units["group"].astype(str), + "age_band": age_band.astype(str), + "income_band": income_band.astype(str), + "size_band": size_band.astype(str), + }, + index=units.index, + ) + + +def _key(frame: pd.DataFrame, columns: Sequence[str]) -> np.ndarray: + if not columns: + return np.full(len(frame), "all", dtype=object) + return frame[list(columns)].astype(str).agg("|".join, axis=1).to_numpy() + + +@dataclass +class CellAssignment: + """Donor spender amounts available to each recipient.""" + + pools: list[np.ndarray] + levels: np.ndarray + + def means(self) -> np.ndarray: + return np.array([pool.mean() for pool in self.pools]) + + def draw(self, rng: np.random.Generator) -> np.ndarray: + return np.array([pool[rng.integers(len(pool))] for pool in self.pools]) + + +def assign_cells( + donor_cells: pd.DataFrame, + donor_moop: np.ndarray, + recipient_cells: pd.DataFrame, + min_donors: int = MIN_CELL_DONORS, + cell_levels: Sequence[tuple[str, ...]] = tuple(CELL_LEVELS), +) -> CellAssignment: + """Match each recipient to the finest cell with enough donor spenders. + + Args: + donor_cells: Cell variables of donors. + donor_moop: Reported MOOP of donors. + recipient_cells: Cell variables of recipients. + min_donors: Minimum number of donor spenders a cell needs. + cell_levels: Cell definitions to try, finest first. + + Returns: + For each recipient, the array of positive donor MOOP values in its + cell and the index into ``cell_levels`` used (0 = finest). + """ + spender = donor_moop > 0 + pools: list[np.ndarray] = [None] * len(recipient_cells) # type: ignore + levels = np.full(len(recipient_cells), -1) + unresolved = np.ones(len(recipient_cells), dtype=bool) + for level, columns in enumerate(cell_levels): + donor_key = _key(donor_cells.loc[spender], columns) + amounts = donor_moop[spender] + grouped: dict[str, np.ndarray] = { + key: amounts[donor_key == key] for key in np.unique(donor_key) + } + recipient_key = _key(recipient_cells, columns) + for i in np.flatnonzero(unresolved): + pool = grouped.get(recipient_key[i]) + if pool is not None and len(pool) >= min_donors: + pools[i] = pool + levels[i] = level + unresolved[i] = False + if not unresolved.any(): + break + if unresolved.any(): + raise ValueError("Some recipients found no donor cell.") + return CellAssignment(pools=pools, levels=levels) + + +def fit_participation(X: pd.DataFrame, moop: np.ndarray) -> Pipeline: + """Fit an (effectively unpenalized) logistic regression for MOOP > 0.""" + model = make_pipeline( + StandardScaler(), + LogisticRegression(C=1e6, max_iter=10_000), + ) + model.fit(X.to_numpy(), (moop > 0).astype(int)) + return model + + +def participation_probability(model: Pipeline, X: pd.DataFrame) -> np.ndarray: + return model.predict_proba(X.to_numpy())[:, 1] + + +@dataclass +class FoldData: + """Donor and recipient data for one cross-fitting fold.""" + + X_donor: pd.DataFrame + moop_donor: np.ndarray + cells_donor: pd.DataFrame + group_donor: np.ndarray + X_recipient: pd.DataFrame + cells_recipient: pd.DataFrame + group_recipient: np.ndarray + + +def _repeat(values: np.ndarray, n_draws: int) -> np.ndarray: + return np.tile(values, (n_draws, 1)) + + +def logit_family( + fold: FoldData, n_draws: int, rng: np.random.Generator +) -> dict[str, np.ndarray]: + """Run every logistic-participation method on one fold. + + Sharing the logistic fit and the cell assignment across the variants + keeps their differences down to the participation and amount rules. + """ + model = fit_participation(fold.X_donor, fold.moop_donor) + prob = participation_probability(model, fold.X_recipient) + cells = assign_cells(fold.cells_donor, fold.moop_donor, fold.cells_recipient) + cell_mean = cells.means() + group_mean = assign_cells( + fold.cells_donor, + fold.moop_donor, + fold.cells_recipient, + cell_levels=COARSE_CELL_LEVELS, + ).means() + + spenders = fold.moop_donor > 0 + regression = LinearRegression().fit( + fold.X_donor.to_numpy()[spenders], fold.moop_donor[spenders] + ) + reg_mean = np.clip(regression.predict(fold.X_recipient.to_numpy()), 0.0, None) + + predicted_spender = prob >= PARTICIPATION_CUTOFF + out = { + "logit_cell_mean": _repeat( + np.where(predicted_spender, cell_mean, 0.0), n_draws + ), + "logit_reg_mean": _repeat(np.where(predicted_spender, reg_mean, 0.0), n_draws), + "logit_expected_value": _repeat(prob * cell_mean, n_draws), + "logit_group_mean": _repeat( + np.where(predicted_spender, group_mean, 0.0), n_draws + ), + } + draw_hotdeck = np.empty((n_draws, len(prob))) + threshold_hotdeck = np.empty((n_draws, len(prob))) + draw_cell_mean = np.empty((n_draws, len(prob))) + for r in range(n_draws): + participates = rng.random(len(prob)) < prob + donor_amount = cells.draw(rng) + draw_hotdeck[r] = np.where(participates, donor_amount, 0.0) + threshold_hotdeck[r] = np.where(predicted_spender, donor_amount, 0.0) + draw_cell_mean[r] = np.where(participates, cell_mean, 0.0) + out["logit_draw_hotdeck"] = draw_hotdeck + out["logit_threshold_hotdeck"] = threshold_hotdeck + out["logit_draw_cell_mean"] = draw_cell_mean + return out + + +def nearest_donors( + donor_pred: np.ndarray, recipient_pred: np.ndarray, k: int +) -> np.ndarray: + """Indices of each recipient's ``k`` nearest donors by predicted value. + + Args: + donor_pred: Predicted values of donors. + recipient_pred: Predicted values of recipients. + k: Number of neighbors; capped at the number of donors. + + Returns: + Integer array of shape ``(n_recipients, min(k, n_donors))`` with + positions into ``donor_pred``. + """ + k = min(k, len(donor_pred)) + order = np.argsort(donor_pred, kind="stable") + sorted_pred = donor_pred[order] + position = np.searchsorted(sorted_pred, recipient_pred) + # The k nearest donors sit in the window of 2k sorted donors around the + # insertion point. + offsets = np.arange(-k, k) + candidates = np.clip(position[:, None] + offsets[None, :], 0, len(sorted_pred) - 1) + distance = np.abs(sorted_pred[candidates] - recipient_pred[:, None]) + # Clipping at either end repeats an index; count each donor once. + repeated = np.concatenate( + [ + np.zeros((len(candidates), 1), dtype=bool), + np.diff(candidates, axis=1) == 0, + ], + axis=1, + ) + distance = np.where(repeated, np.inf, distance) + closest = np.argsort(distance, axis=1, kind="stable")[:, :k] + return order[np.take_along_axis(candidates, closest, axis=1)] + + +def pmm( + fold: FoldData, + n_draws: int, + rng: np.random.Generator, + neighbors: int = PMM_NEIGHBORS, +) -> np.ndarray: + """Predictive mean matching within coverage groups. + + Within each coverage group, an OLS regression of log(1 + MOOP) on the + predictors (all donors, zeros included) gives a predicted value for + donors and recipients. Each recipient takes the reported MOOP of one + donor drawn at random from its ``neighbors`` nearest donors by + predicted value, so imputed values are always reported values. + """ + draws = np.empty((n_draws, len(fold.X_recipient))) + for group in np.unique(fold.group_recipient): + donor_mask = fold.group_donor == group + recipient_mask = fold.group_recipient == group + X_d = fold.X_donor.to_numpy()[donor_mask] + y_d = fold.moop_donor[donor_mask] + model = LinearRegression().fit(X_d, np.log1p(y_d)) + nearest = nearest_donors( + model.predict(X_d), + model.predict(fold.X_recipient.to_numpy()[recipient_mask]), + neighbors, + ) + rows = np.arange(len(nearest)) + for r in range(n_draws): + pick = nearest[rows, rng.integers(nearest.shape[1], size=len(rows))] + draws[r, recipient_mask] = y_d[pick] + return draws + + +def _as_frame(prediction) -> pd.DataFrame: + if isinstance(prediction, dict): + return next(iter(prediction.values())) + return prediction + + +def qrf( + fold: FoldData, + n_draws: int, + seed: int, + n_jobs: int = 4, +) -> np.ndarray: + """microimpute QRF on all donors; each draw samples a random quantile.""" + train = fold.X_donor.copy() + train["moop"] = fold.moop_donor + predictors = list(fold.X_donor.columns) + fitted = QRF(seed=seed, log_level="ERROR").fit( + train, + predictors=predictors, + imputed_variables=["moop"], + not_numeric_categorical=predictors, + n_jobs=n_jobs, + ) + draws = np.empty((n_draws, len(fold.X_recipient))) + for r in range(n_draws): + draws[r] = _as_frame(fitted.predict(fold.X_recipient))["moop"].to_numpy() + return draws + + +def zi_qrf(fold: FoldData, n_draws: int, seed: int) -> np.ndarray: + """microimpute ZeroInflatedImputer (gate + QRF on positive donors).""" + train = fold.X_donor.copy() + train["moop"] = fold.moop_donor + predictors = list(fold.X_donor.columns) + fitted = ZeroInflatedImputer( + base_imputer_class=QRF, seed=seed, log_level="ERROR" + ).fit( + train, + predictors=predictors, + imputed_variables=["moop"], + not_numeric_categorical=predictors, + ) + draws = np.empty((n_draws, len(fold.X_recipient))) + for r in range(n_draws): + draws[r] = _as_frame(fitted.predict(fold.X_recipient))["moop"].to_numpy() + return draws + + +MAIN_METHODS = [ + "logit_cell_mean", + "logit_reg_mean", + "logit_draw_hotdeck", + "pmm", + "qrf", + "zi_qrf", +] +DECOMPOSITION_METHODS = [ + "logit_threshold_hotdeck", + "logit_draw_cell_mean", + "logit_expected_value", + "logit_group_mean", +] +DETERMINISTIC = { + "logit_cell_mean", + "logit_reg_mean", + "logit_expected_value", + "logit_group_mean", +} +METHOD_LABELS = { + "reported": "Reported (CPS ASEC)", + "logit_cell_mean": "(a) Logit, p >= 0.5, cell mean", + "logit_reg_mean": "(a) Logit, p >= 0.5, regression mean", + "logit_draw_hotdeck": "(b) Logit draw, donor amount in cell", + "pmm": "PMM within coverage groups", + "qrf": "(c) QRF (microimpute)", + "zi_qrf": "(c) Zero-inflated QRF (microimpute)", + "logit_threshold_hotdeck": "Logit, p >= 0.5, donor amount in cell", + "logit_draw_cell_mean": "Logit draw, cell mean", + "logit_expected_value": "p x cell mean for every unit", + "logit_group_mean": "Logit, p >= 0.5, coverage-group mean", +} +STOCHASTIC = { + "logit_draw_hotdeck", + "pmm", + "qrf", + "zi_qrf", + "logit_threshold_hotdeck", + "logit_draw_cell_mean", +} + +ImputeFn = Callable[[FoldData, int, np.random.Generator], np.ndarray] diff --git a/experiments/moop_spm/published_2025.json b/experiments/moop_spm/published_2025.json new file mode 100644 index 00000000..bf4918a7 --- /dev/null +++ b/experiments/moop_spm/published_2025.json @@ -0,0 +1,36 @@ +{ + "source": "U.S. Census Bureau, Poverty in the United States: 2025 (P60-290), September 2026", + "report_url": "https://www.census.gov/library/publications/2026/demo/p60-290.html", + "tables": { + "table_4": "https://www2.census.gov/programs-surveys/demo/tables/p60/290/table_4_spm_person.xlsx", + "table_5": "https://www2.census.gov/programs-surveys/demo/tables/p60/290/table_5_spm_program_effect_rates.xlsx", + "table_8": "https://www2.census.gov/programs-surveys/demo/tables/p60/290/table_8_spm_hist.xlsx" + }, + "accessed": "2026-11-05", + "spm_rate_percent": { + "All people": 13.11, + "Under 18": 13.4, + "18 to 64": 12.28, + "65 and older": 15.37 + }, + "spm_rate_moe90": { + "All people": 0.31, + "Under 18": 0.56, + "18 to 64": 0.34, + "65 and older": 0.53 + }, + "medical_expenses_effect_pp": { + "All people": 2.28, + "Under 18": 1.99, + "18 to 64": 1.88, + "65 and older": 3.83 + }, + "medical_expenses_effect_moe90": { + "All people": 0.13, + "Under 18": 0.21, + "18 to 64": 0.14, + "65 and older": 0.26 + }, + "number_in_poverty_thousands": 44390, + "population_thousands": 338600 +} diff --git a/experiments/moop_spm/render.py b/experiments/moop_spm/render.py new file mode 100644 index 00000000..04831c1b --- /dev/null +++ b/experiments/moop_spm/render.py @@ -0,0 +1,410 @@ +"""Render the results table and note from ``results/results.json``. + +Every number in ``TABLE.md`` and ``NOTE.md`` comes from the results file +written by ``run.py``. The note's prose lives in ``NOTE.template.md`` +with ``{placeholders}``; rendering fails on any placeholder the results +cannot fill. + +Usage (from the repository root): + + python -m experiments.moop_spm.render +""" + +from __future__ import annotations + +import json +import string +from pathlib import Path + +HERE = Path(__file__).resolve().parent +MAIN = [ + "logit_cell_mean", + "logit_reg_mean", + "logit_draw_hotdeck", + "pmm", + "qrf", + "zi_qrf", +] +SHORT = { + "reported": "Reported", + "logit_cell_mean": "(a) Logit + cell mean", + "logit_reg_mean": "(a) Logit + OLS mean", + "logit_draw_hotdeck": "(b) Logit draw + donor", + "pmm": "PMM", + "qrf": "(c) QRF", + "zi_qrf": "(c) ZI-QRF", +} +DECOMPOSITION = [ + "logit_cell_mean", + "logit_threshold_hotdeck", + "logit_draw_cell_mean", + "logit_draw_hotdeck", + "logit_expected_value", + "logit_group_mean", +] +# Bands of resources before MOOP shown in Table 2 (index into RATIO_BANDS). +BAND_ROWS = [ + (1, "1.0 to 1.25"), + (2, "1.25 to 1.5"), + (3, "1.5 to 2.0"), + (4, "2.0 to 3.0"), +] +DIMENSION_TITLES = { + "all": None, + "age": "Age", + "disability": "Disability", + "coverage": "Current health insurance", + "pre_moop_ratio": "Resources before MOOP / threshold", +} + + +def pct(x: float, digits: int = 1) -> str: + return f"{100 * x:.{digits}f}" + + +def signed(x: float, digits: int = 2) -> str: + value = 100 * x + if round(value, digits) == 0: + value = 0.0 + return f"{value:+.{digits}f}" + + +def dollars(x: float) -> str: + return f"${x:,.0f}" + + +def main_table(results: dict) -> list[str]: + methods = results["methods"] + groups = results["groups"] + counts = results["group_person_counts"] + header = ["Group", "People (n)", SHORT["reported"]] + [SHORT[m] for m in MAIN] + lines = [ + "| " + " | ".join(header) + " |", + "|" + "|".join(["---"] + ["---:"] * (len(header) - 1)) + "|", + ] + last_dim = None + for i, (dim, label) in enumerate(groups): + if dim != last_dim and DIMENSION_TITLES.get(dim): + lines.append(f"| **{DIMENSION_TITLES[dim]}** |" + " |" * (len(header) - 1)) + last_dim = dim + row = [label, f"{counts[i]:,}", pct(methods["reported"]["rate"][i])] + for method in MAIN: + entry = methods[method] + row.append(f"{signed(entry['diff'][i])} ({100 * entry['diff_se'][i]:.2f})") + lines.append("| " + " | ".join(row) + " |") + return lines + + +def distribution_table(results: dict) -> list[str]: + methods = results["methods"] + columns = ["reported"] + MAIN + header = ["Statistic"] + [SHORT[m] for m in columns] + lines = [ + "| " + " | ".join(header) + " |", + "|" + "|".join(["---"] + ["---:"] * len(columns)) + "|", + ] + rows = [ + ("Mean", lambda d: dollars(d["mean"])), + ("Median", lambda d: dollars(d["median"])), + ("Share with zero MOOP", lambda d: f"{100 * d['share_zero']:.1f}%"), + ("75th percentile", lambda d: dollars(d["p75"])), + ("90th percentile", lambda d: dollars(d["p90"])), + ("99th percentile", lambda d: dollars(d["p99"])), + ("Total ($ billions)", lambda d: f"{d['aggregate_billions']:,.0f}"), + ] + for b, label in BAND_ROWS: + rows.append( + ( + f"Mean, resources before MOOP {label} x threshold", + lambda d, b=b: dollars(d["band_means"][b]), + ) + ) + for name, fmt in rows: + lines.append( + "| " + + " | ".join([name] + [fmt(methods[m]["distribution"]) for m in columns]) + + " |" + ) + flows = [ + ("People newly poor (%)", "newly_poor_share"), + ("People newly not poor (%)", "newly_nonpoor_share"), + ] + for name, key in flows: + lines.append( + "| " + + " | ".join([name, "0.0"] + [f"{100 * methods[m][key]:.2f}" for m in MAIN]) + + " |" + ) + subset = [ + ( + "Change in rate among people in units where Census imputes no " + "MOOP or Medicare premium (pp)" + ), + "0.00", + ] + [ + f"{signed(methods[m]['no_allocation_subset']['diff'])} " + f"({100 * methods[m]['no_allocation_subset']['diff_se']:.2f})" + for m in MAIN + ] + lines.append("| " + " | ".join(subset) + " |") + return lines + + +def decomposition_table(results: dict) -> list[str]: + methods = results["methods"] + lines = [ + "| Participation rule | Amount rule | Rate (%) | Change (pp) | SE |", + "|---|---|---:|---:|---:|", + ] + rules = { + "logit_cell_mean": ("p >= 0.5", "Cell mean"), + "logit_threshold_hotdeck": ("p >= 0.5", "Random donor in cell"), + "logit_draw_cell_mean": ("Bernoulli(p)", "Cell mean"), + "logit_draw_hotdeck": ("Bernoulli(p)", "Random donor in cell"), + "logit_expected_value": ("None (all units)", "p x cell mean"), + "logit_group_mean": ("p >= 0.5", "Coverage-group mean (5 cells)"), + } + for method in DECOMPOSITION: + entry = methods[method] + participation, amount = rules[method] + lines.append( + f"| {participation} | {amount} | {pct(entry['rate'][0], 2)} | " + f"{signed(entry['diff'][0])} | {100 * entry['diff_se'][0]:.2f} |" + ) + return lines + + +def summary_table(results: dict) -> list[str]: + """Compact table for the note: headline rates, key groups, shape.""" + methods = results["methods"] + keyed = {"::".join(g): i for i, g in enumerate(results["groups"])} + columns = ["reported"] + MAIN + header = [""] + [SHORT[m] for m in columns] + lines = [ + "| " + " | ".join(header) + " |", + "|" + "|".join(["---"] + ["---:"] * len(columns)) + "|", + ] + + def change(group: str, with_se: bool = False): + i = keyed[group] + + def fmt(m: str) -> str: + if m == "reported": + return "" + entry = methods[m] + text = signed(entry["diff"][i]) + if with_se: + text += f" ({100 * entry['diff_se'][i]:.2f})" + return text + + return fmt + + def dist(key: str, fmt): + return lambda m: fmt(methods[m]["distribution"][key]) + + def flows(m: str) -> str: + if m == "reported": + return "" + entry = methods[m] + return ( + f"{100 * entry['newly_poor_share']:.2f} / " + f"{100 * entry['newly_nonpoor_share']:.2f}" + ) + + rows = [ + ( + "SPM rate, all people (%)", + lambda m: pct(methods[m]["rate"][keyed["all::All people"]], 2), + ), + ("Change, all people (pp, SE)", change("all::All people", with_se=True)), + ("Change, 65 and older (pp)", change("age::65 and older")), + ("Change, with a disability (pp)", change("disability::With a disability")), + ("Change, Medicare (pp)", change("coverage::Medicare")), + ( + "Change, resources before MOOP 1.0 to 1.25 x threshold (pp)", + change("pre_moop_ratio::1.0 to 1.25"), + ), + ( + "Mean MOOP, resources before MOOP 1.0 to 1.25 x threshold", + lambda m: dollars(methods[m]["distribution"]["band_means"][1]), + ), + ("Units with zero MOOP (%)", dist("share_zero", lambda x: f"{100 * x:.1f}")), + ("Median MOOP", dist("median", dollars)), + ("99th percentile MOOP", dist("p99", dollars)), + ("Total MOOP ($ billions)", dist("aggregate_billions", lambda x: f"{x:,.0f}")), + ("People newly poor / newly not poor (%)", flows), + ] + for name, fmt in rows: + lines.append("| " + " | ".join([name] + [fmt(m) for m in columns]) + " |") + return lines + + +def placeholder_values(results: dict) -> dict[str, str]: + """Flat name-to-string map the note template can reference.""" + methods = results["methods"] + groups = results["groups"] + data = results["data"] + published = results["published"] + values: dict[str, str] = { + "n_persons": f"{data['n_persons']:,}", + "n_units": f"{data['n_units']:,}", + "n_households": f"{data['n_households']:,}", + "population_millions": f"{data['population_millions']:.1f}", + "folds": str(results["config"]["folds"]), + "draws": str(results["config"]["draws"]), + "share_units_any_allocated": f"{100 * data['share_units_any_allocated']:.0f}", + "share_persons_moop_allocated": f"{100 * data['share_persons_moop_allocated']:.0f}", + "share_persons_premium_allocated": f"{100 * data['share_persons_premium_allocated']:.0f}", + "share_persons_any_allocated": f"{100 * data['share_persons_any_allocated']:.0f}", + "share_units_reported_zero": f"{100 * data['share_units_reported_zero']:.1f}", + "published_rate": f"{published['spm_rate_percent']['All people']:.2f}", + "published_moe90": f"{published['spm_rate_moe90']['All people']:.2f}", + "published_moop_effect": f"{published['medical_expenses_effect_pp']['All people']:.2f}", + "published_poor_thousands": f"{published['number_in_poverty_thousands']:,}", + } + keyed = {"::".join(g): i for i, g in enumerate(groups)} + slug = { + "all::All people": "all", + "age::Under 18": "child", + "age::18 to 64": "adult", + "age::65 and older": "senior", + "disability::With a disability": "disabled", + "disability::No disability": "nondisabled", + "coverage::Medicare": "medicare", + "coverage::Employer": "employer", + "coverage::Direct purchase": "direct", + "coverage::Medicaid or other means-tested": "medicaid", + "coverage::Uninsured": "uninsured", + "pre_moop_ratio::1.0 to 1.25": "band1", + "pre_moop_ratio::1.25 to 1.5": "band2", + "pre_moop_ratio::1.5 to 2.0": "band3", + "pre_moop_ratio::2.0 to 3.0": "band4", + "pre_moop_ratio::3.0 or more": "band5", + } + for key, name in slug.items(): + i = keyed[key] + for method, entry in methods.items(): + values[f"rate_{name}_{method}"] = pct(entry["rate"][i]) + values[f"rate2_{name}_{method}"] = pct(entry["rate"][i], 2) + values[f"diff_{name}_{method}"] = signed(entry["diff"][i]) + values[f"diff1_{name}_{method}"] = signed(entry["diff"][i], 1) + # Unsigned magnitudes for prose such as "raises the rate by X". + values[f"pp_{name}_{method}"] = f"{100 * abs(entry['diff'][i]):.2f}" + values[f"pp1_{name}_{method}"] = f"{100 * abs(entry['diff'][i]):.1f}" + values[f"se_{name}_{method}"] = f"{100 * entry['diff_se'][i]:.2f}" + rep = methods["reported"] + values["reported_rate_se"] = f"{100 * rep['rate_se'][0]:.2f}" + values["reported_moe90"] = f"{1.645 * 100 * rep['rate_se'][0]:.2f}" + values["moop_effect"] = f"{-100 * methods['no_moop']['diff'][0]:.2f}" + values["moop_effect_se"] = f"{100 * methods['no_moop']['diff_se'][0]:.2f}" + values["reported_poor_thousands"] = ( + f"{rep['rate'][0] * data['population_millions'] * 1000:,.0f}" + ) + values["summary_table"] = "\n".join(summary_table(results)) + reported_total = rep["distribution"]["aggregate_billions"] + moop_effect = -methods["no_moop"]["diff"][0] + for method, entry in methods.items(): + dist = entry["distribution"] + values[f"share_of_effect_{method}"] = ( + f"{100 * entry['diff'][0] / moop_effect:.0f}" + ) + values[f"agg_{method}"] = f"{dist['aggregate_billions']:,.0f}" + values[f"agg_gap_{method}"] = ( + f"{100 * (1 - dist['aggregate_billions'] / reported_total):.1f}" + ) + values[f"mean_{method}"] = dollars(dist["mean"]) + values[f"median_{method}"] = dollars(dist["median"]) + values[f"p99_{method}"] = dollars(dist["p99"]) + values[f"p90_{method}"] = dollars(dist["p90"]) + values[f"zero_{method}"] = f"{100 * dist['share_zero']:.1f}" + values[f"nearline_{method}"] = dollars(dist["near_line_mean"]) + for b, _ in BAND_ROWS: + values[f"bandmean{b}_{method}"] = dollars(dist["band_means"][b]) + if method != "reported": + values[f"newly_poor_{method}"] = f"{100 * entry['newly_poor_share']:.2f}" + values[f"newly_nonpoor_{method}"] = ( + f"{100 * entry['newly_nonpoor_share']:.2f}" + ) + values[f"drawsd_{method}"] = f"{100 * entry['overall_draw_sd']:.3f}" + sub = entry["no_allocation_subset"] + values[f"subset_diff_{method}"] = signed(sub["diff"]) + values[f"subset_se_{method}"] = f"{100 * sub['diff_se']:.2f}" + values["share_units_no_allocation"] = ( + f"{100 * (1 - data['share_units_any_allocated']):.0f}" + ) + # The note says methods (a), (b) and PMM "raise" the subset's rate. + ab_pmm = ["logit_cell_mean", "logit_reg_mean", "logit_draw_hotdeck", "pmm"] + subset_ab_pmm = [methods[m]["no_allocation_subset"]["diff"] for m in ab_pmm] + if min(subset_ab_pmm) <= 0: + raise ValueError("A method (a), (b) or PMM no longer raises the subset rate.") + values["subset_pp_low_ab_pmm"] = f"{100 * min(subset_ab_pmm):.2f}" + values["subset_pp_high_ab_pmm"] = f"{100 * max(subset_ab_pmm):.2f}" + return values + + +def main() -> None: + results = json.loads((HERE / "results" / "results.json").read_text()) + data = results["data"] + folds = results["config"]["folds"] + source = ( + f"Source: {data['survey']} (calendar-year {data['reference_year']} " + f"income), {data['n_persons']:,} people in {data['n_units']:,} SPM " + f"units and {data['n_households']:,} households. We impute each " + f"unit's MOOP from models we fit on the other {folds - 1} of {folds} " + "folds of whole households and average " + f"{results['config']['draws']} draws for each stochastic method. We " + "weight rates by the ASEC person weight and compute the standard " + "errors in parentheses from the 160 replicate weights, treating the " + "fitted imputation models as fixed." + ) + table_1 = ( + "## Table 1. SPM poverty rate with reported MOOP (%), and the change " + "under each imputation (percentage points, SE)" + ) + table_1_note = ( + "Every method imputes MOOP of zero or more, so people whose unit's " + "resources before MOOP fall below its threshold stay poor, and the " + "'Under 1.0' row cannot change." + ) + table_2 = ( + "## Table 2. MOOP distribution across SPM units (unit weights) and " + "gross flows of people (person weights)" + ) + table_3 = ( + "## Table 3. How the participation and amount rules of the logit " + "approach move the rate (all people)" + ) + table = [ + "# MOOP imputation and SPM poverty: tables", + "", + source, + "", + table_1, + "", + *main_table(results), + "", + table_1_note, + "", + table_2, + "", + *distribution_table(results), + "", + table_3, + "", + *decomposition_table(results), + "", + ] + (HERE / "TABLE.md").write_text("\n".join(table)) + + template = (HERE / "NOTE.template.md").read_text() + values = placeholder_values(results) + fields = {name for _, name, _, _ in string.Formatter().parse(template) if name} + missing = sorted(fields - set(values)) + if missing: + raise KeyError(f"Template placeholders without values: {missing}") + (HERE / "NOTE.md").write_text(template.format(**values)) + print("wrote TABLE.md and NOTE.md") + + +if __name__ == "__main__": + main() diff --git a/experiments/moop_spm/results/.gitignore b/experiments/moop_spm/results/.gitignore new file mode 100644 index 00000000..e934adfd --- /dev/null +++ b/experiments/moop_spm/results/.gitignore @@ -0,0 +1 @@ +cache/ diff --git a/experiments/moop_spm/results/results.json b/experiments/moop_spm/results/results.json new file mode 100644 index 00000000..0015c2c2 --- /dev/null +++ b/experiments/moop_spm/results/results.json @@ -0,0 +1,1685 @@ +{ + "groups": [ + [ + "all", + "All people" + ], + [ + "age", + "Under 18" + ], + [ + "age", + "18 to 64" + ], + [ + "age", + "65 and older" + ], + [ + "disability", + "With a disability" + ], + [ + "disability", + "No disability" + ], + [ + "disability", + "Not asked (mostly under 15)" + ], + [ + "coverage", + "Medicare" + ], + [ + "coverage", + "Employer" + ], + [ + "coverage", + "Direct purchase" + ], + [ + "coverage", + "Medicaid or other means-tested" + ], + [ + "coverage", + "TRICARE, VA or CHAMPVA" + ], + [ + "coverage", + "Uninsured" + ], + [ + "pre_moop_ratio", + "Under 1.0" + ], + [ + "pre_moop_ratio", + "1.0 to 1.25" + ], + [ + "pre_moop_ratio", + "1.25 to 1.5" + ], + [ + "pre_moop_ratio", + "1.5 to 2.0" + ], + [ + "pre_moop_ratio", + "2.0 to 3.0" + ], + [ + "pre_moop_ratio", + "3.0 or more" + ] + ], + "group_person_counts": [ + 134729, + 31444, + 77083, + 26202, + 14477, + 94282, + 25970, + 27328, + 66520, + 8341, + 18134, + 3174, + 11232, + 13803, + 8441, + 9233, + 17495, + 30010, + 55747 + ], + "methods": { + "reported": { + "rate": [ + 0.1311167481463543, + 0.13391314768360826, + 0.1227728935055599, + 0.15377903597654416, + 0.22650265077441, + 0.11662416320335221, + 0.13179831972283498, + 0.1655511081179155, + 0.03927344743889666, + 0.1830291783452912, + 0.3211404432041963, + 0.07158687578614623, + 0.2575160613078942, + 1.0, + 0.24904795316391506, + 0.0580351973091584, + 0.015613799700366192, + 0.003205602288123372, + 0.0005050823529117521 + ], + "rate_se": [ + 0.0018588778777483958, + 0.0034248493266703637, + 0.002047004606950538, + 0.0031913248726654795, + 0.004720027401982268, + 0.0018577453350814479, + 0.0036695832505496302, + 0.003346339647614366, + 0.0013512909560803852, + 0.008464100416272931, + 0.00694969820172331, + 0.009599933498527842, + 0.00690480253979412, + 0.0, + 0.008908376186531767, + 0.004758269857449709, + 0.0015372929958834547, + 0.0006807169472628365, + 0.00016539128557088874 + ], + "diff": [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ], + "diff_se": [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ], + "distribution": { + "band_means": [ + 1550.2326498624302, + 2384.907329611911, + 3381.045403045494, + 4510.285658463513, + 6086.405622898044, + 8758.680056741825 + ], + "mean": 6010.649674747584, + "share_zero": 0.06397480538010987, + "aggregate_billions": 858.9108113234299, + "near_line_mean": 3662.534944614901, + "median": 3700.0, + "p75": 8160.0, + "p90": 14300.0, + "p95": 19363.0, + "p99": 33640.0 + }, + "n_draws": 1, + "no_allocation_subset": { + "rate": 0.11888757471294231, + "diff": 0.0, + "diff_se": 0.0 + } + }, + "logit_cell_mean": { + "rate": [ + 0.13415628217653797, + 0.13562351436516815, + 0.12408642421216211, + 0.1636152292389468, + 0.23851159983907055, + 0.11870002015577995, + 0.1332762113328937, + 0.17757514293977786, + 0.03717757342208421, + 0.1836198550055982, + 0.3329829809921909, + 0.06983741107438445, + 0.2581705915154295, + 1.0, + 0.3959271137135803, + 0.01049172765052451, + 0.0, + 0.0, + 0.0 + ], + "rate_se": [ + 0.0019468154070516092, + 0.0035186951222816853, + 0.002119999471433626, + 0.0036362454305232336, + 0.0049436709905289495, + 0.0019408721939548593, + 0.003726839656372921, + 0.003704759827623985, + 0.0013549184712824693, + 0.00781083735948231, + 0.007200169953566473, + 0.009276990524667015, + 0.007398082313775466, + 0.0, + 0.01058725137953196, + 0.0017747313212909811, + 0.0, + 0.0, + 0.0 + ], + "diff": [ + 0.0030395340301836715, + 0.0017103666815598828, + 0.0013135307066022145, + 0.009836193262402648, + 0.012008949064660546, + 0.002075856952427735, + 0.001477891610058707, + 0.012024034821862345, + -0.002095874016812445, + 0.000590676660306988, + 0.0118425377879946, + -0.0017494647117617806, + 0.0006545302075353088, + 0.0, + 0.14687916054966524, + -0.04754346965863389, + -0.015613799700366192, + -0.003205602288123372, + -0.0005050823529117521 + ], + "diff_se": [ + 0.0008180059629516925, + 0.0013035249101429728, + 0.0008513585174365269, + 0.001700760175690585, + 0.0022166555193685843, + 0.0008534195314263675, + 0.0013639264474922183, + 0.0017400452539473661, + 0.0007995634521281728, + 0.003912246604921555, + 0.002264939995008925, + 0.0021096058855819425, + 0.0026591946094825444, + 0.0, + 0.010079266639459208, + 0.005009044521643856, + 0.0015372929958834547, + 0.0006807169472628365, + 0.00016539128557088874 + ], + "newly_poor_share": 0.013729688297802767, + "newly_nonpoor_share": 0.010690154267619214, + "overall_draw_sd": 0.0, + "distribution": { + "mean": 6155.6312349156, + "share_zero": 0.005279168755883774, + "aggregate_billions": 879.6284102868192, + "near_line_mean": 4081.0472789812484, + "median": 6131.531531531532, + "p75": 8625.363829787235, + "p90": 11435.809585492229, + "p95": 11843.143099787685, + "p99": 15585.390476190476, + "band_means": [ + 1870.5061230776732, + 3189.307895967847, + 3678.8893039791515, + 4788.33309364516, + 6115.830696118508, + 8743.743631067751 + ] + }, + "no_allocation_subset": { + "rate": 0.12248282022546712, + "diff": 0.003595245512524811, + "diff_se": 0.0010590374111963492 + }, + "n_draws": 1 + }, + "logit_reg_mean": { + "rate": [ + 0.13456367230036045, + 0.1365381597759004, + 0.12391407168349722, + 0.1652512730300891, + 0.2405044491594176, + 0.11864655276566594, + 0.13459096502966644, + 0.1778728633731612, + 0.03802978929505018, + 0.20580694960070725, + 0.3243529512179724, + 0.06567192702516342, + 0.25535292027178713, + 1.0, + 0.39272830542182846, + 0.018800882761431173, + 0.00025718107874482494, + 0.0, + 0.0 + ], + "rate_se": [ + 0.0019297947742283497, + 0.00355469058130079, + 0.002133128126743632, + 0.0033258321521071033, + 0.004961029840896307, + 0.0019341852344244394, + 0.0037624407372463054, + 0.003396460333037017, + 0.0013933921165993413, + 0.008902797249971825, + 0.006976568959238853, + 0.00886977606713273, + 0.007321829076126535, + 0.0, + 0.011268433729019143, + 0.002408981277020418, + 0.00016852617912994163, + 0.0, + 0.0 + ], + "diff": [ + 0.003446924154006159, + 0.0026250120922921283, + 0.0011411781779373242, + 0.011472237053544942, + 0.014001798385007591, + 0.002022389562313723, + 0.0027926453068314594, + 0.012321755255245675, + -0.001243658143846478, + 0.02277777125541605, + 0.003212508013776094, + -0.005914948760982808, + -0.0021631410361070413, + 0.0, + 0.1436803522579134, + -0.039234314547727224, + -0.015356618621621366, + -0.003205602288123372, + -0.0005050823529117521 + ], + "diff_se": [ + 0.0008360480955614244, + 0.0013638438599782386, + 0.0009299518217500681, + 0.001608095393742314, + 0.00214726813111919, + 0.000873479869299981, + 0.0014862955043137044, + 0.0016425867760879685, + 0.0008410904704522769, + 0.004762662420736158, + 0.0021883155293683674, + 0.003551474809015206, + 0.0027622365662070465, + 0.0, + 0.010534575496027043, + 0.004859843827809802, + 0.0015157676026204886, + 0.0006807169472628365, + 0.00016539128557088874 + ], + "newly_poor_share": 0.014004473042537253, + "newly_nonpoor_share": 0.010557548888531124, + "overall_draw_sd": 0.0, + "distribution": { + "mean": 6133.498352364934, + "share_zero": 0.04271336726641513, + "aggregate_billions": 876.4656619755368, + "near_line_mean": 3973.2215435922753, + "median": 5891.305254140832, + "p75": 8982.978379356733, + "p90": 11507.143411865734, + "p95": 12902.1198417475, + "p99": 15769.244207504948, + "band_means": [ + 2139.2481601516083, + 3108.0399971311067, + 3667.4699536830517, + 4611.998286184863, + 5871.831558908602, + 8795.034112455074 + ] + }, + "no_allocation_subset": { + "rate": 0.12144512309537356, + "diff": 0.002557548382431246, + "diff_se": 0.0010741198580138027 + }, + "n_draws": 1 + }, + "logit_expected_value": { + "rate": [ + 0.13185490814571796, + 0.1340100566600272, + 0.12231537098228347, + 0.1589165608945135, + 0.23249253432766445, + 0.11673157231191805, + 0.13189301657242325, + 0.17201108386170375, + 0.03668481478116816, + 0.1812791557198023, + 0.32896777520192655, + 0.06777401072948212, + 0.2558179312672401, + 1.0, + 0.36123834196046917, + 0.0086775137220621, + 0.0, + 0.0, + 0.0 + ], + "rate_se": [ + 0.0019495264562355505, + 0.003529941001211207, + 0.0021269012424618566, + 0.003577058161495216, + 0.004870658946944573, + 0.0019489126467759999, + 0.0037170153300361067, + 0.0035972427715047148, + 0.0013511746806648649, + 0.007710759224980344, + 0.007159537944493091, + 0.00920916717000466, + 0.0073054250535077615, + 0.0, + 0.010462499115272446, + 0.0016501149655657615, + 0.0, + 0.0, + 0.0 + ], + "diff": [ + 0.0007381599993636634, + 9.69089764189246e-05, + -0.00045752252327642406, + 0.005137524917969344, + 0.005989883553254444, + 0.00010740910856583963, + 9.469684958826607e-05, + 0.006459975743788232, + -0.002588632657728497, + -0.0017500226254889006, + 0.007827331997730225, + -0.0038128650566641026, + -0.0016981300406540556, + 0.0, + 0.1121903887965541, + -0.049357683587096295, + -0.015613799700366192, + -0.003205602288123372, + -0.0005050823529117521 + ], + "diff_se": [ + 0.0007778415172662687, + 0.0012991654857934859, + 0.0008238767003684714, + 0.0016431087363203355, + 0.0020354030361025804, + 0.0008159018976255727, + 0.0013812427602133454, + 0.0016624905852716915, + 0.0007827821044075652, + 0.0038670575181979287, + 0.002228537265668587, + 0.00243842213940639, + 0.002573172904484096, + 0.0, + 0.009694225268483923, + 0.0050625812983993304, + 0.0015372929958834547, + 0.0006807169472628365, + 0.00016539128557088874 + ], + "newly_poor_share": 0.011867941998969457, + "newly_nonpoor_share": 0.011129781999605948, + "overall_draw_sd": 0.0, + "distribution": { + "mean": 5989.323876012299, + "share_zero": 0.0, + "aggregate_billions": 855.8633938086812, + "near_line_mean": 3836.8891068483226, + "median": 5858.598852547809, + "p75": 8517.372864312221, + "p90": 11408.300447091582, + "p95": 11775.329137318375, + "p99": 15538.421207496805, + "band_means": [ + 1620.767903258226, + 2915.0225616475104, + 3425.4666682227853, + 4565.641459478447, + 5942.080219538605, + 8653.654440016455 + ] + }, + "no_allocation_subset": { + "rate": 0.12023174015284101, + "diff": 0.0013441654398986952, + "diff_se": 0.0009865229670288846 + }, + "n_draws": 1 + }, + "logit_group_mean": { + "rate": [ + 0.1518298935292951, + 0.14601269934809355, + 0.13926499088189884, + 0.19704230707608433, + 0.286219075942204, + 0.13357081365245865, + 0.14399342555107872, + 0.2160072843241851, + 0.046957056508792255, + 0.208330282937395, + 0.3508961995965275, + 0.07296066665490294, + 0.2703396947181051, + 1.0, + 0.5771262029116778, + 0.10059036457500722, + 0.0005277538124302876, + 0.0, + 0.0 + ], + "rate_se": [ + 0.001916084285788687, + 0.003547638555171311, + 0.002115717557843004, + 0.003627462361800205, + 0.004745597110354895, + 0.0019809167963349226, + 0.003798596637989846, + 0.0036815264688785784, + 0.0014671792570957265, + 0.008086433613936772, + 0.007056145438873903, + 0.009442087264323983, + 0.0072075135199208195, + 0.0, + 0.011365787820070965, + 0.004306975496823684, + 0.00015386128427577844, + 0.0, + 0.0 + ], + "diff": [ + 0.020713145382940806, + 0.012099551664485286, + 0.016492097376338943, + 0.043263271099540174, + 0.05971642516779402, + 0.01694665044910644, + 0.012195105828243735, + 0.050456176206269576, + 0.007683609069895597, + 0.0253011045921038, + 0.02975575639233119, + 0.001373790868756719, + 0.012823633410210922, + 0.0, + 0.3280782497477628, + 0.04255516726584883, + -0.015086045887935905, + -0.003205602288123372, + -0.0005050823529117521 + ], + "diff_se": [ + 0.00095603451557515, + 0.0015988107137469827, + 0.0009946392894943218, + 0.0020006928290158283, + 0.002907794154664922, + 0.0009726709067019774, + 0.0017208996151363, + 0.0021354647707454303, + 0.0008988577873044474, + 0.004127679611779838, + 0.0027895248385098673, + 0.0030548004960280002, + 0.0027954392571739235, + 0.0, + 0.01146497111396994, + 0.006173757445266712, + 0.0015156002644129511, + 0.0006807169472628365, + 0.00016539128557088874 + ], + "newly_poor_share": 0.028362078213072756, + "newly_nonpoor_share": 0.007648932830132147, + "overall_draw_sd": 0.0, + "distribution": { + "mean": 6451.750259697265, + "share_zero": 0.005279168755883774, + "aggregate_billions": 921.9432756651952, + "near_line_mean": 6062.9886365521415, + "median": 7178.6743306764765, + "p75": 7205.237434308045, + "p90": 7243.382978723404, + "p95": 7964.515263644774, + "p99": 8464.614849187936, + "band_means": [ + 5027.304500409313, + 5750.231114286323, + 5914.754000891628, + 6315.091795006906, + 6729.75087732737, + 6984.767444836202 + ] + }, + "no_allocation_subset": { + "rate": 0.14047762113486845, + "diff": 0.021590046421926135, + "diff_se": 0.0013852092871485232 + }, + "n_draws": 1 + }, + "logit_draw_hotdeck": { + "rate": [ + 0.13497467589728124, + 0.13814477963016772, + 0.12532665451021344, + 0.16124920271937723, + 0.23242608040649704, + 0.12008007406032241, + 0.1360302538519213, + 0.17403654147552666, + 0.040687702850948436, + 0.18504352852387337, + 0.331033217033513, + 0.06873632655570744, + 0.258162272860673, + 1.0, + 0.2978935666823896, + 0.06299810859034971, + 0.0175340572414454, + 0.003662284961536977, + 0.000767416416662028 + ], + "rate_se": [ + 0.0017805675412109463, + 0.0032243233506604567, + 0.0019673887165264654, + 0.0032126776482786328, + 0.004500292768659439, + 0.0017899470206137963, + 0.003422163214020956, + 0.003278953470888837, + 0.0011815555906528653, + 0.007481880847773399, + 0.006837136384722952, + 0.008812471093231241, + 0.006940516519463729, + 0.0, + 0.005641736013713324, + 0.0020386174199169336, + 0.0006332141982323176, + 0.0002129556930161621, + 6.433945459216442e-05 + ], + "diff": [ + 0.0038579277509269427, + 0.004231631946559455, + 0.002553761004653543, + 0.007470166742833073, + 0.0059234296320870294, + 0.0034559108569701957, + 0.004231934129086318, + 0.00848543335761115, + 0.0014142554120517778, + 0.002014350178582175, + 0.009892773829316692, + -0.002850549230438787, + 0.0006462115527788304, + 0.0, + 0.04884561351847452, + 0.004962911281191314, + 0.001920257541079209, + 0.00045668267341360475, + 0.0002623340637502759 + ], + "diff_se": [ + 0.0006741552158614499, + 0.0010599015483495466, + 0.0007187656148785475, + 0.001333375023197983, + 0.001823368133438438, + 0.000694479860762997, + 0.0011260203638576377, + 0.001371428107696201, + 0.0007142118639856538, + 0.0032482593705501616, + 0.0017843270113983, + 0.002882017251427382, + 0.002314944995923872, + 0.0, + 0.007372116522066939, + 0.004912543629542214, + 0.0016638442120374415, + 0.0007259623274631244, + 0.00017522472361162703 + ], + "newly_poor_share": 0.017448604860594213, + "newly_nonpoor_share": 0.013590677109667515, + "overall_draw_sd": 0.000614283579310217, + "distribution": { + "mean": 5983.595486821716, + "share_zero": 0.06477683038642965, + "aggregate_billions": 855.0448175027071, + "near_line_mean": 3841.982031171648, + "median": 3678.2, + "p75": 8153.5, + "p90": 14275.5, + "p95": 19281.9, + "p99": 33146.9, + "band_means": [ + 1613.4017042778405, + 2949.148528598858, + 3446.2031241285563, + 4546.273683539886, + 5953.482764171465, + 8633.30246884922 + ], + "mean_draw_sd": 59.6027629665208 + }, + "no_allocation_subset": { + "rate": 0.124081832502952, + "diff": 0.00519425779000969, + "diff_se": 0.0009369883893189433 + }, + "n_draws": 10 + }, + "logit_threshold_hotdeck": { + "rate": [ + 0.13697340322602986, + 0.1396156131590491, + 0.1270335252368159, + 0.16473230428613586, + 0.23695623528749968, + 0.12183536540276378, + 0.1374719779418524, + 0.1778061781367272, + 0.04123334301814727, + 0.186701926157094, + 0.3348399478156611, + 0.06972871204876063, + 0.26202015908743975, + 1.0, + 0.32357982703444177, + 0.06684028188867659, + 0.01835308778771935, + 0.003692565364105016, + 0.0007961812308279958 + ], + "rate_se": [ + 0.0017948775882505672, + 0.003246392216841143, + 0.0019876196216614688, + 0.0032313270871775694, + 0.0045119754043835125, + 0.001805225792306883, + 0.0034428190748220524, + 0.0033097923148856863, + 0.0011917317467107003, + 0.007489971086042907, + 0.006838923296799728, + 0.008870998215748688, + 0.006967605417407373, + 0.0, + 0.005930325667855044, + 0.0020774144372910195, + 0.0006415166231634716, + 0.0002107046731267152, + 6.605984334044724e-05 + ], + "diff": [ + 0.00585665507967556, + 0.005702465475440843, + 0.004260631731256015, + 0.010953268309591702, + 0.01045358451308967, + 0.005211202199411569, + 0.005673658219017413, + 0.012255070018811681, + 0.001959895579250609, + 0.0036727478118027956, + 0.013699504611464797, + -0.001858163737385593, + 0.0045040977795455794, + 0.0, + 0.0745318738705267, + 0.008805084579518195, + 0.002739288087353157, + 0.00048696307598164413, + 0.00029109887791624366 + ], + "diff_se": [ + 0.0006794175225984243, + 0.001069732189337917, + 0.0007260529635477616, + 0.0013615569401235924, + 0.0018671096792058537, + 0.0007019110537027104, + 0.0011242160712046439, + 0.0014115912418256828, + 0.0007137451976704757, + 0.0032641480384659185, + 0.0018095162562039955, + 0.002788657789128847, + 0.0022989887558641235, + 0.0, + 0.007495271380460777, + 0.004971436464602843, + 0.0016725471121637968, + 0.0007252112138399993, + 0.00017413489363684505 + ], + "newly_poor_share": 0.018861674555135422, + "newly_nonpoor_share": 0.013005019475460039, + "overall_draw_sd": 0.0006676981919144081, + "distribution": { + "mean": 6154.287131518974, + "share_zero": 0.005279168755883774, + "aggregate_billions": 879.436340377347, + "near_line_mean": 4099.34977799078, + "median": 3893.4, + "p75": 8271.3, + "p90": 14399.9, + "p95": 19403.8, + "p99": 33302.0, + "band_means": [ + 1875.8936464197527, + 3218.9069321261127, + 3716.013161342137, + 4789.961982836939, + 6125.168520028066, + 8723.150028551625 + ], + "mean_draw_sd": 64.86140940233224 + }, + "no_allocation_subset": { + "rate": 0.1260954581499038, + "diff": 0.007207883436961493, + "diff_se": 0.0009350920653806055 + }, + "n_draws": 10 + }, + "logit_draw_cell_mean": { + "rate": [ + 0.13212074650669664, + 0.13439381710839018, + 0.1224075785360038, + 0.15958796018801794, + 0.23291485035790857, + 0.11698835968247441, + 0.1320999967262258, + 0.1730572933764528, + 0.036800126156612885, + 0.18152376084604177, + 0.32941327565824957, + 0.06864149295758178, + 0.25468952199155953, + 1.0, + 0.3640145736366964, + 0.010001221113481163, + 0.0, + 0.0, + 0.0 + ], + "rate_se": [ + 0.0019159051956378881, + 0.0034855844211144097, + 0.0020863389269515004, + 0.0035614726220981677, + 0.004858132470751054, + 0.0019098979400282781, + 0.0036902164971821566, + 0.003618361333627387, + 0.0013408124246821057, + 0.00775918364791198, + 0.007139178386352306, + 0.009168036744187456, + 0.0072847972851559665, + 0.0, + 0.009748022300175841, + 0.0016765935229417067, + 0.0, + 0.0, + 0.0 + ], + "diff": [ + 0.0010039983603423452, + 0.000480669424781921, + -0.00036531496955609366, + 0.005808924211473776, + 0.006412199583498562, + 0.00036419647912219744, + 0.0003016770033908156, + 0.0075061852585372735, + -0.0024733212822837733, + -0.0015054174992494285, + 0.008272832454053247, + -0.0029453828285644496, + -0.002826539316334642, + 0.0, + 0.11496662047278136, + -0.04803397619567723, + -0.015613799700366192, + -0.003205602288123372, + -0.0005050823529117521 + ], + "diff_se": [ + 0.0007772505825684602, + 0.0012620364574717916, + 0.0008111069093334946, + 0.0015901274188968618, + 0.0020757320687225375, + 0.0008115264754368732, + 0.0013281591017036637, + 0.0016327256471097775, + 0.000792296051459046, + 0.003823925435675958, + 0.0021418044997327346, + 0.002114379624585828, + 0.00252946627669586, + 0.0, + 0.009443963471884738, + 0.004972617275628578, + 0.0015372929958834547, + 0.0006807169472628365, + 0.00016539128557088874 + ], + "newly_poor_share": 0.012459071581848625, + "newly_nonpoor_share": 0.011455073221506386, + "overall_draw_sd": 0.0002070621817168782, + "distribution": { + "mean": 5988.816125525065, + "share_zero": 0.06477683038642965, + "aggregate_billions": 855.7908371955797, + "near_line_mean": 3832.4444978257015, + "median": 5857.577658358923, + "p75": 8598.156916019761, + "p90": 11435.809585492229, + "p95": 11843.143099787685, + "p99": 15585.390476190476, + "band_means": [ + 1615.8003905279763, + 2915.470013532779, + 3422.211772204054, + 4557.887923798225, + 5945.679672500067, + 8654.814102205582 + ], + "mean_draw_sd": 5.663741333820941 + }, + "no_allocation_subset": { + "rate": 0.12071737987490072, + "diff": 0.0018298051619584088, + "diff_se": 0.0010132371637605918 + }, + "n_draws": 10 + }, + "pmm": { + "rate": [ + 0.13427030256063682, + 0.13804966220979922, + 0.12463634965868274, + 0.15982791488762849, + 0.23065124752434632, + 0.11939345939049598, + 0.13590840778118202, + 0.17140487431234883, + 0.04160127888299841, + 0.1901403311770817, + 0.3244615640560235, + 0.07007638396990841, + 0.25702938396480557, + 1.0, + 0.29097602131294, + 0.06033186640202967, + 0.017308904269072382, + 0.003523517245754172, + 0.0007024324336039325 + ], + "rate_se": [ + 0.0017986573188207822, + 0.003349919517425816, + 0.0019960193611326523, + 0.003077585457601469, + 0.004660079670212162, + 0.001790901717959959, + 0.0035581776971178717, + 0.003218477754773275, + 0.0012737463640193555, + 0.007651585158299677, + 0.006932701535018352, + 0.008852882302138093, + 0.007032228668708854, + 0.0, + 0.007026196806982932, + 0.002824179476301053, + 0.0010608946897945478, + 0.00034308014192678364, + 0.00010704651364114224 + ], + "diff": [ + 0.003153554414282528, + 0.00413651452619096, + 0.0018634561531228472, + 0.006048878911084327, + 0.004148596749936312, + 0.0027692961871437632, + 0.004110088058347033, + 0.005853766194433313, + 0.002327831444101751, + 0.007111152831790507, + 0.003321120851827175, + -0.0015104918162378161, + -0.000486677343088604, + 0.0, + 0.041928068149024916, + 0.0022966690928712744, + 0.0016951045687061909, + 0.00031791495763080005, + 0.00019735008069218037 + ], + "diff_se": [ + 0.0006656177733958356, + 0.0011204494646935572, + 0.0007130394717814084, + 0.0015321059921304429, + 0.0020139739046894367, + 0.0006903773248147367, + 0.0011775163283457744, + 0.0015436465760354885, + 0.0007040806004508576, + 0.003645091395470062, + 0.0019582117081428318, + 0.002880183197418108, + 0.0022296205910069485, + 0.0, + 0.007522609147771935, + 0.00506752419306695, + 0.0017574155439823971, + 0.0007533436992944742, + 0.00020545095707093106 + ], + "newly_poor_share": 0.016296812689101266, + "newly_nonpoor_share": 0.013143258274818792, + "overall_draw_sd": 0.0007066526013457527, + "distribution": { + "mean": 5922.240125546339, + "share_zero": 0.06161216799369856, + "aggregate_billions": 846.277248939615, + "near_line_mean": 3795.798667149561, + "median": 3624.5, + "p75": 8090.9, + "p90": 14157.9, + "p95": 19157.2, + "p99": 32666.4, + "band_means": [ + 1811.4431965376484, + 2799.1295991808997, + 3451.9231436717964, + 4526.966224085245, + 5783.040535510334, + 8537.725599428359 + ], + "mean_draw_sd": 51.34901277574327 + }, + "no_allocation_subset": { + "rate": 0.12370894398263602, + "diff": 0.004821369269693704, + "diff_se": 0.0009805236692727892 + }, + "n_draws": 10 + }, + "qrf": { + "rate": [ + 0.1304418917246645, + 0.13393973543529658, + 0.12193776487928376, + 0.15282346409458544, + 0.22439268710345409, + 0.11594086535254887, + 0.1320357934178739, + 0.1647544461057615, + 0.039062485524114676, + 0.1808743102691106, + 0.32124538454334933, + 0.07150189504340257, + 0.25419635089782605, + 1.0, + 0.2545412971933657, + 0.0504193805722122, + 0.012705854518909062, + 0.0022419666268124433, + 0.0007147144606562746 + ], + "rate_se": [ + 0.0017680359660310998, + 0.0032131460408222692, + 0.0019524401343579258, + 0.0031344231764976398, + 0.004607460477414942, + 0.001785548291301098, + 0.003432182810041383, + 0.0032583328444268055, + 0.0012135916424184591, + 0.007266010679292295, + 0.0068830309228798325, + 0.008960473321918298, + 0.006858436937269146, + 0.0, + 0.006310439816679268, + 0.0022834010379751313, + 0.000742882216756691, + 0.00025402573141758997, + 0.00010302523239104148 + ], + "diff": [ + -0.0006748564216897956, + 2.6587751688317285e-05, + -0.0008351286262761326, + -0.0009555718819587167, + -0.002109963670955922, + -0.0006832978508033438, + 0.00023747369503890337, + -0.0007966620121540036, + -0.0002109619147819819, + -0.002154868076180594, + 0.00010494133915300852, + -8.498074274365663e-05, + -0.003319710410068122, + 0.0, + 0.0054933440294506475, + -0.007615816736946196, + -0.002907945181457129, + -0.0009636356613109288, + 0.00020963210774452247 + ], + "diff_se": [ + 0.0006639161495395829, + 0.0011048302288689871, + 0.000703693894937724, + 0.0013595787142926777, + 0.0018284318745864065, + 0.00068565374836129, + 0.0011947473557091912, + 0.0014616382761895189, + 0.0006960965411983834, + 0.0036166645495945703, + 0.0017074378275275107, + 0.0026009881219099897, + 0.0023046077089981837, + 0.0, + 0.007899147325533836, + 0.004780932588376095, + 0.0016651916853366157, + 0.0007323722804144055, + 0.00020184301700901504 + ], + "newly_poor_share": 0.013311703717753328, + "newly_nonpoor_share": 0.0139865601394432, + "overall_draw_sd": 0.0007480484672292954, + "distribution": { + "mean": 5797.189515730481, + "share_zero": 0.06150229455370617, + "aggregate_billions": 828.4077462160282, + "near_line_mean": 3544.1610550600444, + "median": 3623.4705882352937, + "p75": 7948.8, + "p90": 13820.329411764706, + "p95": 18425.91764705882, + "p99": 31355.53529411765, + "band_means": [ + 1426.0406603582237, + 2351.8914371191254, + 3281.5551467815335, + 4335.23025980269, + 5808.785836005858, + 8492.973953440758 + ], + "mean_draw_sd": 48.149493973807225 + }, + "no_allocation_subset": { + "rate": 0.11900340588984247, + "diff": 0.00011583117690015798, + "diff_se": 0.0008957336867377862 + }, + "n_draws": 10 + }, + "zi_qrf": { + "rate": [ + 0.13065321441816305, + 0.13456294016605738, + 0.12218243044542777, + 0.15247656505374374, + 0.22430617418538826, + 0.11618008023508287, + 0.13231570276888185, + 0.16448305112790573, + 0.03945289231662562, + 0.18260666194357905, + 0.32080684039983737, + 0.07064845364694208, + 0.2546721735327806, + 1.0, + 0.25514470392500355, + 0.052314674628251776, + 0.012975087264387288, + 0.0025954954104857857, + 0.0005457220109579402 + ], + "rate_se": [ + 0.001804718636530011, + 0.0032890607104412274, + 0.001996061056236309, + 0.003076082566497863, + 0.004638832248708925, + 0.0018025723049676116, + 0.00348411340347414, + 0.003237699638267949, + 0.0012219768458732554, + 0.007447951332274899, + 0.006879207102407378, + 0.009017159726054599, + 0.006892417035506484, + 0.0, + 0.006540448425263213, + 0.0022524504743104074, + 0.0007906954170421504, + 0.00029588639974349143, + 8.346892675942626e-05 + ], + "diff": [ + -0.0004635337281912455, + 0.0006497924824491186, + -0.0005904630601321209, + -0.0013024709228004194, + -0.002196476589021751, + -0.00044408296826933813, + 0.0005173830460468709, + -0.0010680569900097803, + 0.00017944487772895978, + -0.00042251640171214744, + -0.00033360280435895673, + -0.0009384221392041425, + -0.002843887775113574, + 0.0, + 0.006096750761088482, + -0.0057205226809066215, + -0.0026387124359789037, + -0.0006101068776375863, + 4.0639658046188044e-05 + ], + "diff_se": [ + 0.0006725128310102399, + 0.0010797774353122976, + 0.0007451719272210365, + 0.0012697236949680156, + 0.0018197819559848726, + 0.0006938649209268466, + 0.0011875625079665287, + 0.0013370642250423723, + 0.0007671266354864416, + 0.0036450531804156716, + 0.0017837674044191454, + 0.002990079981198583, + 0.0023523748274293975, + 0.0, + 0.00813524419857648, + 0.005015050593068027, + 0.0016098834904267137, + 0.000749276546613338, + 0.00018841172649408384 + ], + "newly_poor_share": 0.013685253048786874, + "newly_nonpoor_share": 0.014148786776978157, + "overall_draw_sd": 0.00044444271345847814, + "distribution": { + "mean": 5791.182914160099, + "share_zero": 0.06131936219814614, + "aggregate_billions": 827.5494138713918, + "near_line_mean": 3572.107986524211, + "median": 3610.0882352941176, + "p75": 7891.276470588235, + "p90": 13764.558823529413, + "p95": 18409.685294117633, + "p99": 31424.270588235297, + "band_means": [ + 1434.4808508689134, + 2371.6676721508657, + 3271.6450220657816, + 4388.86072395638, + 5843.438873157666, + 8441.63518059709 + ], + "mean_draw_sd": 52.29116091182236 + }, + "no_allocation_subset": { + "rate": 0.11957842269182078, + "diff": 0.0006908479788784727, + "diff_se": 0.0009744769243665113 + }, + "n_draws": 10 + }, + "no_moop": { + "rate": [ + 0.10859685755998726, + 0.11433480771362724, + 0.10420784260365511, + 0.11580074358919064, + 0.18469628260740995, + 0.09615451772207399, + 0.1127255615530577, + 0.12615228607528428, + 0.026374506666176595, + 0.1409729845709993, + 0.297760800652063, + 0.058545592740366977, + 0.23299038575738587, + 1.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ], + "rate_se": [ + 0.0017787704897606448, + 0.0032160858539244533, + 0.0019393770575135904, + 0.0029442566760378027, + 0.004634986000645334, + 0.0017425416878091306, + 0.0034265988604119876, + 0.0031026798630482146, + 0.0010597983231123106, + 0.007054781058902049, + 0.007132817383801851, + 0.008302825722358628, + 0.006990887714759059, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ], + "diff": [ + -0.02251989058636704, + -0.019578339969981026, + -0.018565050901904787, + -0.03797829238735352, + -0.04180636816700006, + -0.02046964548127822, + -0.019072758169777285, + -0.03939882204263123, + -0.012898940772720063, + -0.0420561937742919, + -0.023379642552133317, + -0.01304128304577925, + -0.024525675550508302, + 0.0, + -0.24904795316391506, + -0.0580351973091584, + -0.015613799700366192, + -0.003205602288123372, + -0.0005050823529117521 + ], + "diff_se": [ + 0.0007846074982339206, + 0.001273861322319545, + 0.0008286290599487661, + 0.0015831915498418048, + 0.0020249224383628683, + 0.000809414192678513, + 0.00133309314913825, + 0.0015710158358619069, + 0.0008288272778743621, + 0.003962644424662755, + 0.002115077900554795, + 0.0044479097021104625, + 0.002584546754914354, + 0.0, + 0.008908376186531767, + 0.004758269857449709, + 0.0015372929958834547, + 0.0006807169472628365, + 0.00016539128557088874 + ], + "newly_poor_share": 0.0, + "newly_nonpoor_share": 0.022519890586367212, + "overall_draw_sd": 0.0, + "distribution": { + "mean": 0.0, + "share_zero": 1.0, + "aggregate_billions": 0.0, + "near_line_mean": 0.0, + "median": 0.0, + "p75": 0.0, + "p90": 0.0, + "p95": 0.0, + "p99": 0.0, + "band_means": [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ] + }, + "no_allocation_subset": { + "rate": 0.1019506970734858, + "diff": -0.01693687763945652, + "diff_se": 0.0010943716814347118 + }, + "n_draws": 1 + } + }, + "data": { + "survey": "CPS ASEC 2026 public-use file", + "reference_year": 2025, + "zip_url": "https://www2.census.gov/programs-surveys/cps/datasets/2026/march/asecpub26csv.zip", + "zip_sha256": "fe819d0d2fc4470c282e76c247b8aa8b6054811619730307e9f2fe6186707d93", + "n_persons": 134729, + "n_units": 55401, + "n_households": 53344, + "population_millions": 338.55489354, + "n_units_by_fold": [ + 11076, + 11109, + 11039, + 11095, + 11082 + ], + "unit_group_counts": { + "employer": 25051, + "medicare": 19598, + "public": 5188, + "direct": 2864, + "uninsured": 2700 + }, + "share_persons_moop_allocated": 0.43122119217095056, + "share_persons_premium_allocated": 0.13230262230106363, + "share_persons_any_allocated": 0.48115105136978675, + "share_units_any_allocated": 0.5637804371762243, + "share_units_reported_zero": 0.0590242053392538, + "identity_checks": { + "resources_gap": 0.0, + "poor_mismatch": 0.0, + "moop_gap": 0.0 + } + }, + "published": { + "source": "U.S. Census Bureau, Poverty in the United States: 2025 (P60-290), September 2026", + "report_url": "https://www.census.gov/library/publications/2026/demo/p60-290.html", + "tables": { + "table_4": "https://www2.census.gov/programs-surveys/demo/tables/p60/290/table_4_spm_person.xlsx", + "table_5": "https://www2.census.gov/programs-surveys/demo/tables/p60/290/table_5_spm_program_effect_rates.xlsx", + "table_8": "https://www2.census.gov/programs-surveys/demo/tables/p60/290/table_8_spm_hist.xlsx" + }, + "accessed": "2026-11-05", + "spm_rate_percent": { + "All people": 13.11, + "Under 18": 13.4, + "18 to 64": 12.28, + "65 and older": 15.37 + }, + "spm_rate_moe90": { + "All people": 0.31, + "Under 18": 0.56, + "18 to 64": 0.34, + "65 and older": 0.53 + }, + "medical_expenses_effect_pp": { + "All people": 2.28, + "Under 18": 1.99, + "18 to 64": 1.88, + "65 and older": 3.83 + }, + "medical_expenses_effect_moe90": { + "All people": 0.13, + "Under 18": 0.21, + "18 to 64": 0.14, + "65 and older": 0.26 + }, + "number_in_poverty_thousands": 44390, + "population_thousands": 338600 + }, + "config": { + "zip_sha256": "fe819d0d2fc4470c282e76c247b8aa8b6054811619730307e9f2fe6186707d93", + "folds": 5, + "draws": 10, + "seed": 20261105, + "skip_forests": false, + "max_households": null, + "features": [ + "head_age", + "head_female", + "n_age_0_17", + "n_age_18_44", + "n_age_45_64", + "n_age_65_74", + "n_age_75_plus", + "n_female_adult", + "n_disabled", + "n_any_medicare", + "n_any_employer", + "n_any_direct", + "n_any_medicaid", + "n_any_military", + "n_uninsured", + "log_cash_income", + "cash_income_ratio", + "famtype_2", + "famtype_3", + "famtype_4", + "famtype_5", + "tenure_2", + "tenure_3", + "state_2", + "state_4", + "state_5", + "state_6", + "state_8", + "state_9", + "state_10", + "state_11", + "state_12", + "state_13", + "state_15", + "state_16", + "state_17", + "state_18", + "state_19", + "state_20", + "state_21", + "state_22", + "state_23", + "state_24", + "state_25", + "state_26", + "state_27", + "state_28", + "state_29", + "state_30", + "state_31", + "state_32", + "state_33", + "state_34", + "state_35", + "state_36", + "state_37", + "state_38", + "state_39", + "state_40", + "state_41", + "state_42", + "state_44", + "state_45", + "state_46", + "state_47", + "state_48", + "state_49", + "state_50", + "state_51", + "state_53", + "state_54", + "state_55", + "state_56" + ], + "min_cell_donors": 30, + "pmm_neighbors": 5, + "cutoff": 0.5, + "cell_levels": [ + [ + "group", + "age_band", + "income_band", + "size_band" + ], + [ + "group", + "age_band", + "income_band" + ], + [ + "group", + "age_band" + ], + [ + "group" + ], + [] + ], + "source_sha256": { + "cps_asec.py": "647857a0fa249c44", + "methods.py": "8ce97d8ddd38df23", + "evaluate.py": "09daf780b35188ad", + "run.py": "324dad8ce5a851e0" + }, + "microimpute_sha256": "e14c884a85e90486", + "environment": { + "python": "3.13.13", + "numpy": "2.5.4", + "pandas": "3.1.1", + "scipy": "1.18.1", + "scikit_learn": "1.8.0", + "quantile_forest": "1.4.2", + "microimpute": "3.1.4" + }, + "input_sha256": "5a129f5d767fff81", + "tag": "a82da82512c894e6", + "n_jobs": 4 + }, + "timings_seconds": { + "logit_family": 38.13895297050476, + "pmm": 0.5446643829345703, + "qrf": 68.79297137260437, + "zi_qrf": 185.8473162651062 + }, + "environment": { + "python": "3.13.13", + "numpy": "2.5.4", + "pandas": "3.1.1", + "scipy": "1.18.1", + "scikit_learn": "1.8.0", + "quantile_forest": "1.4.2", + "microimpute": "3.1.4" + } +} diff --git a/experiments/moop_spm/run.py b/experiments/moop_spm/run.py new file mode 100644 index 00000000..d910e96b --- /dev/null +++ b/experiments/moop_spm/run.py @@ -0,0 +1,322 @@ +"""Cross-fit MOOP imputations on the 2026 CPS ASEC and compare SPM poverty. + +Each SPM unit's reported medical out-of-pocket spending (MOOP) is +replaced by an imputation from models fit on the other folds (folds +group whole households), so every unit gets an out-of-fold imputation. +SPM poverty is then recomputed for every person under each method. + +Usage (from the repository root): + + python -m experiments.moop_spm.run --data-dir ~/data/cps_asec_2026 + +The first run downloads the Census public-use zip (139 MB) into +``--data-dir`` and checks its SHA-256. Per-fold imputations are cached in +``/cache`` so an interrupted run resumes. The cache key covers the +settings, the experiment's source files, every microimpute source file, +the versions of Python and the numerical packages, and the model inputs +with their row keys, so a change to any of them refits every fold. +""" + +from __future__ import annotations + +import argparse +import hashlib +import importlib +import json +import platform +import time +from collections.abc import Callable +from importlib.metadata import version +from pathlib import Path + +import numpy as np +import pandas as pd + +from .cps_asec import ( + REFERENCE_YEAR, + SURVEY_YEAR, + ZIP_SHA256, + ZIP_URL, + allocation_flags, + build_units, + check_spm_identity, + design_matrix, + load_persons, + sha256_of, +) +from .evaluate import evaluate +from .methods import ( + CELL_LEVELS, + DETERMINISTIC, + MIN_CELL_DONORS, + PARTICIPATION_CUTOFF, + PMM_NEIGHBORS, + FoldData, + cell_keys, + logit_family, + pmm, + qrf, + zi_qrf, +) + +HERE = Path(__file__).resolve().parent +SOURCES = ["cps_asec.py", "methods.py", "evaluate.py", "run.py"] +# Distributions whose versions can change fitted models or draws. +PACKAGES = { + "numpy": "numpy", + "pandas": "pandas", + "scipy": "scipy", + "scikit_learn": "scikit-learn", + "quantile_forest": "quantile-forest", + "microimpute": "microimpute", +} +PUBLISHED = json.loads((HERE / "published_2025.json").read_text()) + + +def assign_folds(households: pd.Series, n_folds: int, seed: int) -> np.ndarray: + """Assign whole households to folds at random, as evenly as possible.""" + unique = np.sort(households.unique()) + permuted = np.random.default_rng(seed).permutation(unique) + fold_of = pd.Series(np.arange(len(permuted)) % n_folds, index=permuted) + return households.map(fold_of).to_numpy() + + +def config_hash(config: dict) -> str: + return hashlib.sha256(json.dumps(config, sort_keys=True).encode()).hexdigest()[:16] + + +def environment() -> dict[str, str]: + """Return the Python version and the versions in ``PACKAGES``.""" + versions = {name: version(dist) for name, dist in PACKAGES.items()} + return {"python": platform.python_version(), **versions} + + +def package_sha256(package: str) -> str: + """Hash every ``.py`` file of an imported package, in path order.""" + module = importlib.import_module(package) + root = Path(module.__file__).resolve().parent + digest = hashlib.sha256() + for path in sorted(root.rglob("*.py")): + digest.update(path.relative_to(root).as_posix().encode()) + digest.update(b"\0") + digest.update(path.read_bytes()) + digest.update(b"\0") + return digest.hexdigest()[:16] + + +def frame_sha256(*frames: pd.DataFrame) -> str: + """Hash frames' column names, row keys and values.""" + digest = hashlib.sha256() + for frame in frames: + digest.update(json.dumps([str(c) for c in frame.columns]).encode()) + hashed = pd.util.hash_pandas_object(frame, index=True) + digest.update(hashed.to_numpy().tobytes()) + return digest.hexdigest()[:16] + + +def build_config( + args: argparse.Namespace, + units: pd.DataFrame, + X: pd.DataFrame, + cells: pd.DataFrame, + folds: np.ndarray, +) -> dict: + """Collect everything that can change the fold imputations.""" + fold_frame = pd.DataFrame({"fold": folds}, index=units.index) + return { + "zip_sha256": ZIP_SHA256, + "folds": args.folds, + "draws": args.draws, + "seed": args.seed, + "skip_forests": args.skip_forests, + "max_households": args.max_households, + "features": list(X.columns), + "min_cell_donors": MIN_CELL_DONORS, + "pmm_neighbors": PMM_NEIGHBORS, + "cutoff": PARTICIPATION_CUTOFF, + "cell_levels": [list(level) for level in CELL_LEVELS], + "source_sha256": {name: sha256_of(HERE / name)[:16] for name in SOURCES}, + "microimpute_sha256": package_sha256("microimpute"), + "environment": environment(), + "input_sha256": frame_sha256(units, X, cells, fold_frame), + } + + +def collect_folds( + cache: Path, + tag: str, + folds: np.ndarray, + n_folds: int, + n_draws: int, + fit_fold: Callable[[int], dict[str, np.ndarray]], +) -> tuple[dict[str, np.ndarray], dict[str, float]]: + """Load each fold's imputations from the cache, or fit and store them. + + Returns each method's draws for every unit, ``(n_draws, n_units)``, + and the fitting time of each method summed over folds. + """ + imputations: dict[str, np.ndarray] = {} + timings: dict[str, float] = {} + for fold in range(n_folds): + path = cache / f"{tag}_fold{fold}.npz" + if path.exists(): + stored = dict(np.load(path, allow_pickle=False)) + else: + start = time.time() + stored = fit_fold(fold) + np.savez_compressed(path, **stored) + print( + f"fold {fold} done in {time.time() - start:.0f}s", + flush=True, + ) + names = [str(n) for n in stored.pop("_timing_names")] + for name, seconds in zip(names, stored.pop("_timings")[0]): + timings[name] = timings.get(name, 0.0) + float(seconds) + recipient = folds == fold + for method, draws in stored.items(): + if method not in imputations: + imputations[method] = np.full((n_draws, len(folds)), np.nan) + imputations[method][:, recipient] = draws + return imputations, timings + + +def run_fold( + fold: int, + folds: np.ndarray, + units: pd.DataFrame, + X: pd.DataFrame, + cells: pd.DataFrame, + args: argparse.Namespace, +) -> dict[str, np.ndarray]: + """Fit every method on the other folds and impute this fold.""" + donor = folds != fold + recipient = folds == fold + data = FoldData( + X_donor=X.loc[donor], + moop_donor=units["moop"].to_numpy()[donor], + cells_donor=cells.loc[donor], + group_donor=units["group"].to_numpy()[donor], + X_recipient=X.loc[recipient], + cells_recipient=cells.loc[recipient], + group_recipient=units["group"].to_numpy()[recipient], + ) + rng = np.random.default_rng([args.seed, fold]) + timings: dict[str, float] = {} + out: dict[str, np.ndarray] = {} + + start = time.time() + out.update(logit_family(data, args.draws, rng)) + timings["logit_family"] = time.time() - start + + start = time.time() + out["pmm"] = pmm(data, args.draws, rng) + timings["pmm"] = time.time() - start + + if not args.skip_forests: + start = time.time() + out["qrf"] = qrf(data, args.draws, seed=args.seed + fold, n_jobs=args.n_jobs) + timings["qrf"] = time.time() - start + + start = time.time() + out["zi_qrf"] = zi_qrf(data, args.draws, seed=args.seed + fold) + timings["zi_qrf"] = time.time() - start + out["_timings"] = np.array([list(timings.values())]) + out["_timing_names"] = np.array(list(timings.keys())) + return out + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--data-dir", type=Path, required=True) + parser.add_argument("--output", type=Path, default=HERE / "results") + parser.add_argument("--folds", type=int, default=5) + parser.add_argument("--draws", type=int, default=10) + parser.add_argument("--seed", type=int, default=20261105) + parser.add_argument("--n-jobs", type=int, default=4) + parser.add_argument( + "--skip-forests", + action="store_true", + help="Skip QRF and ZI-QRF (for quick checks).", + ) + parser.add_argument( + "--max-households", + type=int, + default=None, + help="Subsample households (for quick checks).", + ) + args = parser.parse_args() + + output = args.output + cache = output / "cache" + cache.mkdir(parents=True, exist_ok=True) + + persons = load_persons(args.data_dir) + identity = check_spm_identity(persons) + if any(value != 0 for value in identity.values()): + raise ValueError(f"SPM identities fail: {identity}") + if args.max_households: + keep = ( + pd.Series(persons["PH_SEQ"].unique()) + .sample(args.max_households, random_state=args.seed) + .to_numpy() + ) + persons = persons[persons["PH_SEQ"].isin(keep)].reset_index(drop=True) + units = build_units(persons) + X = design_matrix(units) + cells = cell_keys(units) + folds = assign_folds(units["household"], args.folds, args.seed) + + config = build_config(args, units, X, cells, folds) + tag = config_hash(config) + + imputations, timings = collect_folds( + cache, + tag, + folds, + args.folds, + args.draws, + lambda fold: run_fold(fold, folds, units, X, cells, args), + ) + for method, draws in imputations.items(): + if np.isnan(draws).any(): + raise ValueError(f"{method} left units without an imputation") + if (draws < 0).any(): + raise ValueError(f"{method} produced negative MOOP") + if method in DETERMINISTIC: + # Deterministic methods: one draw carries all the information. + imputations[method] = draws[:1] + imputations["no_moop"] = np.zeros((1, len(units))) + + results = evaluate(persons, units, imputations) + allocated = allocation_flags(persons) + results["data"] = { + "survey": f"CPS ASEC {SURVEY_YEAR} public-use file", + "reference_year": REFERENCE_YEAR, + "zip_url": ZIP_URL, + "zip_sha256": ZIP_SHA256, + "n_persons": len(persons), + "n_units": len(units), + "n_households": int(units["household"].nunique()), + "population_millions": float(persons["weight"].sum() / 1e6), + "n_units_by_fold": np.bincount(folds).tolist(), + "unit_group_counts": units["group"].value_counts().to_dict(), + "share_persons_moop_allocated": float(allocated["allocated_moop"].mean()), + "share_persons_premium_allocated": float(allocated["allocated_premium"].mean()), + "share_persons_any_allocated": float(allocated.max(axis=1).mean()), + "share_units_any_allocated": float(units["any_allocated"].mean()), + "share_units_reported_zero": float((units["moop"] == 0).mean()), + "identity_checks": identity, + } + results["published"] = PUBLISHED + results["config"] = {**config, "tag": tag, "n_jobs": args.n_jobs} + results["timings_seconds"] = timings + # The cache tag covers this environment, so it describes the cached + # imputations as well as this run. + results["environment"] = config["environment"] + output.mkdir(parents=True, exist_ok=True) + (output / "results.json").write_text(json.dumps(results, indent=2) + "\n") + print(f"wrote {output / 'results.json'}") + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 6f99fa9c..cf5e77c6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -49,6 +49,7 @@ dev = [ "mypy>=1.2.3,<2.0.0", "build>=1.2.0,<2.0.0", "towncrier>=24.8.0", + "hypothesis>=6.100.0,<7.0.0", ] matching = [ diff --git a/tests/test_moop_spm_experiment.py b/tests/test_moop_spm_experiment.py new file mode 100644 index 00000000..4063c4d4 --- /dev/null +++ b/tests/test_moop_spm_experiment.py @@ -0,0 +1,477 @@ +"""Property tests for the MOOP-imputation SPM poverty experiment. + +The experiment lives in ``experiments/moop_spm``. These tests use small +synthetic data, so they run without the CPS ASEC download. They check +the invariants the published comparison relies on: + +- Fold assignment partitions households and keeps each household whole. +- Poverty is recomputed exactly from resources, MOOP and thresholds, and + it never falls when any unit's MOOP rises. +- People whose unit is poor before MOOP stay poor under any imputation. +- Weighted group rates, replicate SEs and weighted quantiles match direct + reference computations. +- Cell pools hold only same-cell donor spenders, and the nearest-donor + search agrees with a brute-force search. +- The logistic-participation variants follow their stated rules. +- Census's MOOP and Medicare premium imputations both count as allocated. +- Extracted CSVs that differ from the pinned zip are extracted again. +- The fold cache misses whenever the environment, microimpute's source + or the model inputs change. +""" + +import argparse +import os +import zipfile + +import numpy as np +import pandas as pd +import pytest +from hypothesis import HealthCheck, given, settings +from hypothesis import strategies as st + +from experiments.moop_spm import run as run_module +from experiments.moop_spm.cps_asec import allocation_flags, extract_verified +from experiments.moop_spm.evaluate import ( + RATIO_BANDS, + evaluate, + group_rates, + poverty_flags, + ratio_band_index, + replicate_se, + weighted_quantile, +) +from experiments.moop_spm.methods import ( + CELL_LEVELS, + PARTICIPATION_CUTOFF, + FoldData, + _key, + assign_cells, + fit_participation, + logit_family, + nearest_donors, + participation_probability, +) +from experiments.moop_spm.run import ( + assign_folds, + build_config, + collect_folds, + config_hash, + frame_sha256, +) + +SETTINGS = settings( + max_examples=40, + deadline=None, + suppress_health_check=[HealthCheck.too_slow], +) + + +def synthetic_spm(seed: int, n_units: int = 60): + """Persons and units frames with the columns ``evaluate`` reads.""" + rng = np.random.default_rng(seed) + unit_ids = np.arange(1, n_units + 1) * 7 + threshold = rng.uniform(20_000, 40_000, n_units) + moop = np.where(rng.random(n_units) < 0.1, 0.0, rng.gamma(1.0, 5_000, n_units)) + pre_moop = threshold * rng.uniform(0.5, 3.0, n_units) + units = pd.DataFrame( + { + "moop": moop, + "unit_weight": rng.uniform(100, 3_000, n_units), + "pre_moop_ratio": pre_moop / threshold, + "any_allocated": rng.integers(0, 2, n_units), + }, + index=pd.Index(unit_ids, name="SPM_ID"), + ) + sizes = rng.integers(1, 4, n_units) + person_unit = np.repeat(np.arange(n_units), sizes) + n_persons = len(person_unit) + persons = pd.DataFrame( + { + "SPM_ID": unit_ids[person_unit], + "A_AGE": rng.integers(0, 90, n_persons), + "PRDISFLG": rng.choice([-1, 1, 2], n_persons), + "coverage": rng.choice( + [ + "medicare", + "employer", + "direct", + "medicaid", + "military", + "uninsured", + ], + n_persons, + ), + "SPM_MEDXPNS": moop[person_unit], + "SPM_RESOURCES": (pre_moop - moop)[person_unit], + "SPM_POVTHRESHOLD": threshold[person_unit], + "weight": rng.uniform(100, 3_000, n_persons), + } + ) + persons["SPM_POOR"] = ( + persons["SPM_RESOURCES"] < persons["SPM_POVTHRESHOLD"] + ).astype(int) + replicates = persons["weight"].to_numpy()[:, None] * rng.uniform( + 0.5, 1.5, (n_persons, 160) + ) + replicate_frame = pd.DataFrame( + replicates, columns=[f"repwgt_{i}" for i in range(1, 161)] + ) + persons = pd.concat([persons, replicate_frame], axis=1) + return persons, units + + +@SETTINGS +@given( + households=st.lists(st.integers(1, 40), min_size=5, max_size=120), + n_folds=st.integers(2, 6), + seed=st.integers(0, 2**31), +) +def test_folds_partition_households(households, n_folds, seed): + series = pd.Series(households) + folds = assign_folds(series, n_folds, seed) + assert set(np.unique(folds)) <= set(range(n_folds)) + # Every household lands in exactly one fold. + assert (series.groupby(folds).nunique().sum()) == series.nunique() + per_household = pd.Series(folds).groupby(series.to_numpy()).nunique() + assert (per_household == 1).all() + # Household counts per fold differ by at most one. + household_counts = series.groupby(folds).nunique() + assert household_counts.max() - household_counts.min() <= 1 + assert (assign_folds(series, n_folds, seed) == folds).all() + + +@SETTINGS +@given(seed=st.integers(0, 2**31), bump_seed=st.integers(0, 2**31)) +def test_poverty_recomputation_and_monotonicity(seed, bump_seed): + persons, units = synthetic_spm(seed) + reported = poverty_flags(persons, units, units["moop"].to_numpy())[0] + assert (reported == (persons["SPM_POOR"] == 1)).all() + + bump = np.random.default_rng(bump_seed).uniform(0, 20_000, len(units)) + higher = poverty_flags(persons, units, units["moop"].to_numpy() + bump)[0] + # More MOOP can only push people into poverty. + assert (higher >= reported).all() + + zero = poverty_flags(persons, units, np.zeros(len(units)))[0] + assert (zero <= reported).all() + pre_moop_poor = persons["SPM_ID"].map(units["pre_moop_ratio"]).to_numpy() < 1.0 + # Units poor before MOOP stay poor under any nonnegative MOOP. + assert zero[pre_moop_poor].all() + assert higher[pre_moop_poor].all() + + +@SETTINGS +@given(seed=st.integers(0, 2**31)) +def test_group_rates_match_np_average(seed): + rng = np.random.default_rng(seed) + n = 50 + flags = rng.random(n) + weights = rng.uniform(0.1, 10, (n, 3)) + groups = rng.random((n, 4)) < 0.5 + groups[0] = True # keep every group nonempty + rates = group_rates(flags, groups, weights) + for g in range(groups.shape[1]): + for w in range(weights.shape[1]): + expected = np.average(flags[groups[:, g]], weights=weights[groups[:, g], w]) + assert rates[g, w] == pytest.approx(expected) + assert ((rates >= 0) & (rates <= 1)).all() + + +@SETTINGS +@given( + base=st.floats(-5, 5), + deviations=st.lists(st.floats(-1, 1), min_size=160, max_size=160), + scale=st.floats(-3, 3), +) +def test_replicate_se_formula(base, deviations, scale): + estimates = np.array([base] + [base + d for d in deviations]) + expected = np.sqrt(4 / 160 * np.sum(np.square(deviations))) + assert replicate_se(estimates) == pytest.approx(expected) + assert replicate_se(scale * estimates) == pytest.approx(abs(scale) * expected) + assert replicate_se(np.full(161, base)) == 0 + + +@SETTINGS +@given( + values=st.lists(st.floats(-1e6, 1e6), min_size=1, max_size=60), + q=st.floats(0.01, 0.99), +) +def test_weighted_quantile_equal_weights(values, q): + values = np.array(values) + result = weighted_quantile(values, np.ones(len(values)), q) + assert result in values + assert result == pytest.approx(np.quantile(values, q, method="inverted_cdf")) + + +def _cells(rng: np.random.Generator, n: int) -> pd.DataFrame: + return pd.DataFrame( + { + "group": rng.choice(["medicare", "employer", "public"], n), + "age_band": rng.choice(["u45", "45_64", "65p"], n), + "income_band": rng.choice(["lt1", "1_2", "2_4", "4p"], n), + "size_band": rng.choice(["1", "2", "3p"], n), + } + ) + + +@SETTINGS +@given(seed=st.integers(0, 2**31), min_donors=st.integers(1, 40)) +def test_cell_pools_are_same_cell_donor_spenders(seed, min_donors): + rng = np.random.default_rng(seed) + donors = _cells(rng, 400) + moop = np.where(rng.random(400) < 0.2, 0.0, rng.gamma(1.0, 3_000, 400)) + recipients = _cells(rng, 50) + assignment = assign_cells(donors, moop, recipients, min_donors=min_donors) + spender = moop > 0 + donor_keys = [_key(donors.loc[spender], columns) for columns in CELL_LEVELS] + recipient_keys = [_key(recipients, columns) for columns in CELL_LEVELS] + for i, (pool, level) in enumerate(zip(assignment.pools, assignment.levels)): + assert len(pool) >= min_donors + assert (pool > 0).all() + in_cell = donor_keys[level] == recipient_keys[level][i] + expected = np.sort(moop[spender][in_cell]) + assert np.array_equal(np.sort(pool), expected) + # No finer level had enough donor spenders. + for finer in range(level): + assert (donor_keys[finer] == recipient_keys[finer][i]).sum() < min_donors + assert np.allclose(assignment.means(), [p.mean() for p in assignment.pools]) + + +@SETTINGS +@given( + donor_pred=st.lists(st.floats(-50, 50, allow_nan=False), min_size=1, max_size=40), + recipient_pred=st.lists( + st.floats(-60, 60, allow_nan=False), min_size=1, max_size=20 + ), + k=st.integers(1, 8), +) +def test_nearest_donors_matches_brute_force(donor_pred, recipient_pred, k): + donor_pred = np.array(donor_pred) + recipient_pred = np.array(recipient_pred) + nearest = nearest_donors(donor_pred, recipient_pred, k) + k_eff = min(k, len(donor_pred)) + assert nearest.shape == (len(recipient_pred), k_eff) + for i, value in enumerate(recipient_pred): + assert len(set(nearest[i])) == k_eff + got = np.sort(np.abs(donor_pred[nearest[i]] - value)) + brute = np.sort(np.abs(donor_pred - value))[:k_eff] + assert np.allclose(got, brute) + + +def _fold(seed: int) -> FoldData: + rng = np.random.default_rng(seed) + n_donor, n_recipient = 600, 80 + + def features(n): + return pd.DataFrame( + { + "x1": rng.normal(size=n), + "x2": rng.normal(size=n), + "x3": rng.integers(0, 2, n).astype(float), + } + ) + + X_donor = features(n_donor) + logit = 2.0 + 1.5 * X_donor["x1"].to_numpy() + spends = rng.random(n_donor) < 1 / (1 + np.exp(-logit)) + moop = np.where(spends, rng.gamma(1.0, 4_000, n_donor), 0.0) + cells_donor = _cells(rng, n_donor) + cells_recipient = _cells(rng, n_recipient) + return FoldData( + X_donor=X_donor, + moop_donor=moop, + cells_donor=cells_donor, + group_donor=cells_donor["group"].to_numpy(), + X_recipient=features(n_recipient), + cells_recipient=cells_recipient, + group_recipient=cells_recipient["group"].to_numpy(), + ) + + +@settings(max_examples=8, deadline=None) +@given(seed=st.integers(0, 2**31)) +def test_logit_family_rules(seed): + fold = _fold(seed) + out = logit_family(fold, 4, np.random.default_rng(seed)) + prob = participation_probability( + fit_participation(fold.X_donor, fold.moop_donor), fold.X_recipient + ) + cells = assign_cells(fold.cells_donor, fold.moop_donor, fold.cells_recipient) + means = cells.means() + below = prob < PARTICIPATION_CUTOFF + for method, draws in out.items(): + assert draws.shape == (4, len(prob)) + assert (draws >= 0).all(), method + assert np.allclose(out["logit_cell_mean"][0], np.where(below, 0, means)) + assert np.allclose(out["logit_expected_value"][0], prob * means) + assert (out["logit_reg_mean"][:, below] == 0).all() + assert (out["logit_group_mean"][:, below] == 0).all() + assert (out["logit_threshold_hotdeck"][:, below] == 0).all() + for i, pool in enumerate(cells.pools): + for method in ("logit_draw_hotdeck", "logit_threshold_hotdeck"): + positive = out[method][:, i][out[method][:, i] > 0] + assert np.isin(positive, pool).all() + drawn = out["logit_draw_cell_mean"][:, i] + assert np.all((drawn == 0) | np.isclose(drawn, means[i])) + + +@settings(max_examples=10, deadline=None) +@given(seed=st.integers(0, 2**31)) +def test_evaluate_invariants(seed): + persons, units = synthetic_spm(seed, n_units=80) + bump = np.random.default_rng(seed + 1).uniform(0, 10_000, len(units)) + imputations = { + "no_moop": np.zeros((1, len(units))), + "higher": (units["moop"].to_numpy() + bump)[None, :], + "same": np.tile(units["moop"].to_numpy(), (3, 1)), + } + results = evaluate(persons, units, imputations) + groups = results["groups"] + methods = results["methods"] + nonempty = np.array(results["group_person_counts"]) > 0 + assert np.allclose(np.array(methods["same"]["diff"])[nonempty], 0) + assert methods["same"]["newly_poor_share"] == 0 + assert methods["same"]["newly_nonpoor_share"] == 0 + assert (np.array(methods["no_moop"]["diff"])[nonempty] <= 1e-12).all() + assert (np.array(methods["higher"]["diff"])[nonempty] >= -1e-12).all() + band = groups.index(["pre_moop_ratio", RATIO_BANDS[0][2]]) + if nonempty[band]: + for method in methods.values(): + assert method["rate"][band] == pytest.approx(1.0) + + +@SETTINGS +@given( + ratios=st.lists( + st.one_of( + st.floats(-5, 20, allow_nan=False), + st.sampled_from([1.0, 1.25, 1.5, 2.0, 3.0]), + ), + min_size=1, + max_size=50, + ) +) +def test_ratio_band_index_matches_band_masks(ratios): + ratios = np.array(ratios) + index = ratio_band_index(ratios) + for b, (low, high, _) in enumerate(RATIO_BANDS): + mask = (ratios >= low) & (ratios < high) + assert np.array_equal(index == b, mask) + + +def test_allocation_flags_cover_moop_and_premium(): + persons = pd.DataFrame( + { + "I_MOOP": [-1, 0, 1, 2, 3, 0, 0], + "I_MCPREM": [-1, -1, -1, 0, -1, 2, 0], + } + ) + flags = allocation_flags(persons) + assert flags["allocated_moop"].tolist() == [0, 0, 1, 1, 1, 0, 0] + assert flags["allocated_premium"].tolist() == [0, 0, 0, 0, 0, 1, 0] + + +def test_extract_verified_replaces_altered_copies(tmp_path): + zip_path = tmp_path / "release.zip" + with zipfile.ZipFile(zip_path, "w") as archive: + archive.writestr("a.csv", "x,y\n1,2\n") + archive.writestr("b.csv", "z\n3\n") + out = tmp_path / "data" + out.mkdir() + extract_verified(zip_path, out, ["a.csv", "b.csv"]) + assert (out / "a.csv").read_text() == "x,y\n1,2\n" + + # A matching copy stays in place. + os.utime(out / "b.csv", ns=(1, 1)) + extract_verified(zip_path, out, ["b.csv"]) + assert (out / "b.csv").stat().st_mtime_ns == 1 + + # Altered copies, of the same or a different size, come back. + (out / "a.csv").write_text("x,y\n9,2\n") + (out / "b.csv").write_text("z\n3\n4\n") + extract_verified(zip_path, out, ["a.csv", "b.csv"]) + assert (out / "a.csv").read_text() == "x,y\n1,2\n" + assert (out / "b.csv").read_text() == "z\n3\n" + + +def _cache_inputs(n_units: int = 12): + rng = np.random.default_rng(0) + index = pd.Index(np.arange(n_units) * 3 + 1, name="SPM_ID") + units = pd.DataFrame( + { + "moop": rng.gamma(1.0, 1_000, n_units), + "household": np.arange(n_units) // 2, + }, + index=index, + ) + X = pd.DataFrame({"x": rng.normal(size=n_units)}, index=index) + cells = pd.DataFrame({"group": ["employer"] * n_units}, index=index) + folds = assign_folds(units["household"], 3, 0) + args = argparse.Namespace( + folds=3, draws=2, seed=0, skip_forests=False, max_households=None + ) + return args, units, X, cells, folds + + +def _stub_fit(value: float, folds: np.ndarray, calls: list): + def fit(fold: int) -> dict[str, np.ndarray]: + calls.append(fold) + n = int((folds == fold).sum()) + return { + "qrf": np.full((2, n), value), + "_timings": np.array([[0.0]]), + "_timing_names": np.array(["qrf"]), + } + + return fit + + +def test_cache_misses_when_environment_or_microimpute_changes(tmp_path, monkeypatch): + args, units, X, cells, folds = _cache_inputs() + base_env = run_module.environment() + monkeypatch.setattr(run_module, "environment", lambda: base_env) + tag = config_hash(build_config(args, units, X, cells, folds)) + + calls: list = [] + first, _ = collect_folds(tmp_path, tag, folds, 3, 2, _stub_fit(17.0, folds, calls)) + assert calls == [0, 1, 2] + assert (first["qrf"] == 17.0).all() + + # Same fingerprint: every fold comes from the cache. + calls.clear() + again, _ = collect_folds(tmp_path, tag, folds, 3, 2, _stub_fit(99.0, folds, calls)) + assert calls == [] + assert (again["qrf"] == 17.0).all() + + # A new microimpute version, or edited microimpute source, refits. + newer = {**base_env, "microimpute": base_env["microimpute"] + ".post1"} + monkeypatch.setattr(run_module, "environment", lambda: newer) + new_tag = config_hash(build_config(args, units, X, cells, folds)) + assert new_tag != tag + fresh, _ = collect_folds( + tmp_path, new_tag, folds, 3, 2, _stub_fit(99.0, folds, calls) + ) + assert calls == [0, 1, 2] + assert (fresh["qrf"] == 99.0).all() + + monkeypatch.setattr(run_module, "environment", lambda: base_env) + monkeypatch.setattr(run_module, "package_sha256", lambda package: "edited") + assert config_hash(build_config(args, units, X, cells, folds)) != tag + + +def test_cache_key_tracks_model_inputs_and_row_keys(): + args, units, X, cells, folds = _cache_inputs() + tag = config_hash(build_config(args, units, X, cells, folds)) + assert config_hash(build_config(args, units, X, cells, folds)) == tag + + changed = units.copy() + changed.iloc[0, 0] += 1.0 + assert config_hash(build_config(args, changed, X, cells, folds)) != tag + + rekeyed = [frame.copy() for frame in (units, X, cells)] + for frame in rekeyed: + frame.index = frame.index + 1 + assert config_hash(build_config(args, *rekeyed, folds)) != tag + + assert config_hash(build_config(args, units, X, cells, folds[::-1])) != tag + assert frame_sha256(units) != frame_sha256(units.rename(columns={"moop": "m"}))