diff --git a/CHANGELOG.md b/CHANGELOG.md index 146140675..7f29b5c85 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,53 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- **ContinuousDiD post-fit `aggregate()` - a MIXED view/recompute adopter** + (v4 program 2(b) PR-3c; ledger row [M-025]). `ContinuousDiD.fit(aggregate=)` + is deprecated (`FutureWarning`; removed in 4.0, and the no-underscore + `"eventstudy"` spelling dies with it) in favor of post-fit + `ContinuousDiDResults.aggregate(type)` with the unified vocabulary plus + `"dose"` as this estimator's documented extra level. Unlike every prior + adopter the architecture is MIXED: the dose-response curves, the overall + binarized ATT (`ATT^{loc}` under PT; equals `ATT^{glob}` under SPT) and + `ACRT^{glob}` are ALWAYS computed by `fit()` (`aggregate="dose"` was + a fit-time no-op), so `aggregate('simple')` (2 rows, targets `att`/`acrt` + - the dual-estimand case; `n` = the disjoint treated+control unit total, + `n_kind='units'`) and `aggregate('dose')` (2N target-discriminated rows + over the dose grid) are pure VIEWS relaying stored inference verbatim and + work on ANY fit including bootstrapped ones (the stored FINITE overall + t-stat beside percentile p/CI relays through; the dose rows reproduce + `DoseResponseCurve.to_dataframe` exactly; only the df column is NaN under + bootstrap), while `aggregate('event_study')` RECOMPUTES the binarized + event study from a pruned per-cell influence payload retained on the + fit-built kit (per-(g,t) O(n_treated+n_control) IF-ingredient arrays + + unit-level arrays + the panel-level resolved survey design - on + replicate-weight designs the (n_obs x R) replicate matrix rides along; + the K-dimensional spline machinery is NOT retained, and no panel data + columns or raw unit identifiers are). Replicate-weight designs are + supported post-fit (IF-based replicate variance - no refit replay); + bootstrapped fits carry a scalars-only kit and the event-study route + fails closed (`NotImplementedError` - use the deprecated fit-time route + or re-fit with `n_bootstrap=0`; a seeded replay is a tracked follow-up). + The pre-existing fit-time value validation is retained (unknown strings + still raise `ValueError` after the warning). Post-fit-vs-fit-time + inertness is pinned at 1e-14 across plain/multi-cohort/anticipation/ + covariates/survey-TSL/zero-dose-drop/replicate(healthy,dropped, + undefined-df)/discrete/lowest-dose/not-yet-treated/universal-base + designs. This ships the FIRST heterogeneous-`target` + `AggregationResult`, with a rendering amendment: `summary()` gains a + `target` column and a neutral `estimate` heading when targets are mixed, + and `to_dataframe()` orders heterogeneous-target rows by + first-appearance target blocks (att before acrt) with labels ascending + within each block under the existing sortability guard - uniform-target + producers render byte-identically (normative rule added to + `docs/v4-design.md` section 6). `compute_honest_did` / + `compute_pretrends_power` continue to reject ContinuousDiD containers BY + DESIGN (no joint event-study covariance; the binarized bins carry no + reference-period normalization) - both terminal messages now say so. + The event-study machinery moved verbatim to the new leaf module + `diff_diff/continuous_did_aggregation.py` (shared by fit and the + post-fit kit route); the M-092 per-row-df completion row now names + ContinuousDiD alongside EfficientDiD/ImputationDiD. - **ImputationDiD + TwoStageDiD post-fit `aggregate()` via panel-backed recompute kits** (v4 program 2(b) PR-3b; ledger rows [M-021]/[M-022] flip to shimmed, new rows [M-118]/[M-119] claim the reserved balance_e diff --git a/METHODOLOGY_REVIEW.md b/METHODOLOGY_REVIEW.md index 572534929..4381399c3 100644 --- a/METHODOLOGY_REVIEW.md +++ b/METHODOLOGY_REVIEW.md @@ -705,7 +705,7 @@ and covariate-adjusted specifications.) | Field | Value | |-------|-------| -| Module | `continuous_did.py`, `continuous_did_bspline.py`, `continuous_did_results.py` | +| Module | `continuous_did.py`, `continuous_did_aggregation.py`, `continuous_did_bspline.py`, `continuous_did_results.py` | | Primary Reference | Callaway, Goodman-Bacon & Sant'Anna (2024), *Difference-in-Differences with a Continuous Treatment*, NBER WP 32117 | | R Reference | `contdid` v0.1.0 (CRAN) — two parity surfaces at relative tolerance: (a) **scalar overall ATT parity** with raw R `cont_did` / `pte_default` output at `< 0.01` (1%) on all 6 benchmarks; **scalar overall ACRT parity** with raw R `cont_did` at `< 0.01` (1%) on benchmarks 4-5; (b) **harmonized boundary-knot-normalized curve parity** with R-side ATT(d)/ACRT(d) reconstructed under `Boundary.knots = range(treated_doses)` (matching the library) at `< 0.01` max ATT(d) and `< 0.02` max ACRT(d) on benchmarks 1-3 via the benchmark harness (`_run_r_contdid` rebuilds the R-side basis under `Boundary.knots = range(treated_doses)` at `tests/test_methodology_continuous_did.py:333-367`; `_compare_with_r` orchestrates the Python-vs-R comparison at `:395-459`); benchmark 6 is event-study, scalar `overall_att` only (binarized ATT, no curve comparison and no ACRT in event-study mode). Surface (a) is direct raw-package parity; surface (b) is reconstructed-basis parity because raw `contdid` curves use `range(dvals)` instead of `range(dose)`. NOT bit-exact (`atol=1e-8`) like HAD because of the boundary-knots deviation documented below. See `tests/test_methodology_continuous_did.py::TestRBenchmark` | | Status | **Complete** | @@ -715,13 +715,13 @@ and covariate-adjusted specifications.) - [x] **PT and SPT identification** (CGBS 2024 Assumptions 1-2) — two-level parallel trends with explicit untreated-and-doses conditioning; estimands `ATT(d|d)`, `ATT(d)`, `ACRT(d)`, `ATT^{loc}`, `ATT^{glob}`, `ACRT^{glob}` defined in `docs/methodology/continuous-did.md` § 4 + REGISTRY `## ContinuousDiD` Identification block. Hand-calc coverage: `tests/test_methodology_continuous_did.py::TestLinearDoseResponse` (4 tests at `atol=1e-10` / `atol=1e-6` on no-noise linear DGP — locks the `ATT^{glob}` binarization formula `E[ΔY | D > 0] − E[ΔY | D = 0]`, the `ACRT^{glob}` plug-in average, and the `ATT(d) = 2d`, `ACRT(d) = 2` closed forms). - [x] **B-spline basis matching `splines2::bSpline`** (cubic and linear degrees, `num_knots=0` default; global boundary knots from the training-dose range, NOT per-cell) — `tests/test_methodology_continuous_did.py::TestQuadraticWithCubicBasis::test_quadratic_recovery` recovers `ATT(d) = d²` at `atol=1e-6` via a degree-3 basis (cubic spline can represent quadratic exactly). The matching basis algorithm lives in `diff_diff/continuous_did_bspline.py` (216 LoC); the boundary-knots deviation from R `contdid` is documented in the Deviations block below. - [x] **Multi-period (g,t) cell iteration with base period selection** — `TestMultiPeriodAggregation::test_multiple_groups` and `test_gt_cell_count` exercise the cohort iteration on 2-cohort staggered panels; cell counts agree with the R `ptetools`-style convention. Scalar parity with raw R `cont_did` at 1% relative further locks the staggered-aggregation surface via `TestRBenchmark::test_benchmark_4_staggered_dose` and `test_benchmark_5_not_yet_treated` (both assert overall ATT AND overall ACRT at `< 0.01`). -- [x] **Dose-response (`aggregate="dose"`) and event-study (`aggregate="eventstudy"`) aggregation** with group-proportional weights (`n_treated/n_total` per group, divided among post-treatment cells; matches R `ptetools` convention). Two R-side surfaces are exercised: (a) **scalar `overall_att`** via `TestRBenchmark::test_benchmark_1_basic_cubic` / `_2_linear` / `_3_interior_knots` / `_4_staggered_dose` / `_5_not_yet_treated` (dose mode) and `_6_event_study` (event-study mode — binarized ATT only; benchmark 6 validates the event-study code path through the scalar surface, NOT per-horizon `event_study_effects`); (b) **harmonized boundary-knot-normalized ATT(d) / ACRT(d) curves** on benchmarks 1-3 via the benchmark harness — `_run_r_contdid` at `tests/test_methodology_continuous_did.py:333-367` rebuilds the R-side basis under `Boundary.knots = range(treated_doses)` (raw `contdid` curves use `range(dvals)`, so this is reconstructed-basis parity not raw-package parity), and `_compare_with_r` orchestrates the comparison at `:395-459`. Per-benchmark tolerances: all 6 assert overall ATT at `< 0.01` (1%); benchmarks 1-3 additionally assert max ATT(d) at `< 0.01` and max ACRT(d) at `< 0.02` via the helper; benchmarks 4-5 assert overall ACRT at `< 0.01` inline. Per-horizon `event_study_effects` estimates and inference are exercised by Python-side tests at `tests/test_continuous_did.py:557-690` and `:1500-1528` (no R cross-language comparison on the per-horizon surface). Skipped if R / `contdid` not installed via `_check_r_contdid()`; benchmarks use R's `dvals` for exact evaluation-grid alignment between Python and R outputs (boundary knots are harmonized separately under surface (b) — see the `_run_r_contdid` helper's `Boundary.knots = range(treated_doses)` block at `tests/test_methodology_continuous_did.py:333-367`). +- [x] **Dose-response (always computed by `fit()`) and event-study (post-fit `results.aggregate("event_study")`; the fit-time `aggregate="eventstudy"` kwarg is deprecated, row M-025) aggregation** with group-proportional weights (`n_treated/n_total` per group, divided among post-treatment cells; matches R `ptetools` convention). Two R-side surfaces are exercised: (a) **scalar `overall_att`** via `TestRBenchmark::test_benchmark_1_basic_cubic` / `_2_linear` / `_3_interior_knots` / `_4_staggered_dose` / `_5_not_yet_treated` (dose mode) and `_6_event_study` (event-study mode — binarized ATT only; benchmark 6 validates the event-study code path through the scalar surface, NOT per-horizon `event_study_effects`); (b) **harmonized boundary-knot-normalized ATT(d) / ACRT(d) curves** on benchmarks 1-3 via the benchmark harness — `_run_r_contdid` at `tests/test_methodology_continuous_did.py:333-367` rebuilds the R-side basis under `Boundary.knots = range(treated_doses)` (raw `contdid` curves use `range(dvals)`, so this is reconstructed-basis parity not raw-package parity), and `_compare_with_r` orchestrates the comparison at `:395-459`. Per-benchmark tolerances: all 6 assert overall ATT at `< 0.01` (1%); benchmarks 1-3 additionally assert max ATT(d) at `< 0.01` and max ACRT(d) at `< 0.02` via the helper; benchmarks 4-5 assert overall ACRT at `< 0.01` inline. Per-horizon `event_study_effects` estimates and inference are exercised by Python-side tests at `tests/test_continuous_did.py:557-690` and `:1500-1528` (no R cross-language comparison on the per-horizon surface). Skipped if R / `contdid` not installed via `_check_r_contdid()`; benchmarks use R's `dvals` for exact evaluation-grid alignment between Python and R outputs (boundary knots are harmonized separately under surface (b) — see the `_run_r_contdid` helper's `Boundary.knots = range(treated_doses)` block at `tests/test_methodology_continuous_did.py:333-367`). - [x] **Multiplier bootstrap for inference** (PSU-level multiplier weights on the survey path per Phase 6) — implementation in `diff_diff/continuous_did.py`; bootstrap SE invariant on rank-deficient cells locked in `TestEdgeCasesMethodology::test_all_same_dose` (verifies `dose_response_att.se` is finite on a heterogeneous-outcome / identical-dose DGP); 80 unit tests in `tests/test_continuous_did.py` exercise the rest of the bootstrap path. - [x] **Analytical SEs via influence functions** (NOT delta method; corrected post-v3.0.0, see Corrections Made) — IF-based variance with `safe_inference()` joint-NaN consistency on all six estimand fields (`overall_att`, `overall_acrt`, dose-response, event-study). - [x] **Survey support**: weighted B-spline OLS, two-stage linearization (TSL) on influence functions, bootstrap + survey via PSU-level multiplier weights (Phase 3 + Phase 6). Boxed in REGISTRY `## ContinuousDiD` → Implementation Checklist → "Survey design support (Phase 3)" item. - [x] **`+inf` → `0` never-treated recoding** with `UserWarning` reporting the affected row count (axis-E silent-coercion fix per Phase 2 audit) — the R-style convention of `first_treat = +inf` is normalized internally but no longer absorbed silently. **Any negative `first_treat` value (including `-inf`) raises `ValueError`** with the affected row count. Locked in `tests/test_continuous_did.py`. - [x] **Zero-`first_treat` rows with nonzero `dose` force-zeroed** with `UserWarning` reporting the affected row count (axis-E silent-coercion fix per Phase 2 audit) — never-treated cells must have `D=0` for internal consistency; the previous silent zeroing is now signaled. Locked in `tests/test_continuous_did.py`. -- [x] **`bspline_derivative_design_matrix` derivative-construction failure warning** (Phase 2 axis-C #12 silent-failures audit fix) — aggregates failed basis indices into a single `UserWarning` naming them, instead of swallowing `scipy.interpolate.BSpline.ValueError` and leaving silently zeroed derivative columns. Both ACRT point estimates AND analytical/bootstrap inference read the same `dPsi` matrix (`continuous_did.py:1026-1046` and the bootstrap ACRT path at `continuous_did.py:1524-1561`), so both are biased on partial-derivative failure — the warning wording makes that explicit. The all-identical-knot degenerate case (single dose value) remains silently handled because derivatives are mathematically zero there. Locked in `tests/test_continuous_did.py::TestBSplineDerivativeDegenerateBasis` (3 tests: `test_single_dose_is_silent`, `test_valueerror_from_bspline_emits_aggregate_warning`, `test_clean_knots_emit_no_warning`); source-level aggregate-warning block at `diff_diff/continuous_did_bspline.py:150-187`. +- [x] **`bspline_derivative_design_matrix` derivative-construction failure warning** (Phase 2 axis-C #12 silent-failures audit fix) — aggregates failed basis indices into a single `UserWarning` naming them, instead of swallowing `scipy.interpolate.BSpline.ValueError` and leaving silently zeroed derivative columns. Both ACRT point estimates AND analytical/bootstrap inference read the same `dPsi` matrix (`_compute_dose_response_gt` - the `acrt_d = dPsi_eval @ beta_pred` point estimate at `continuous_did.py:1806` and the covariate-path IFs at `:1485-1526` - and the bootstrap ACRT path in `_bootstrap_gt_cell` at `continuous_did.py:2274-2340`), so both are biased on partial-derivative failure — the warning wording makes that explicit. The all-identical-knot degenerate case (single dose value) remains silently handled because derivatives are mathematically zero there. Locked in `tests/test_continuous_did.py::TestBSplineDerivativeDegenerateBasis` (3 tests: `test_single_dose_is_silent`, `test_valueerror_from_bspline_emits_aggregate_warning`, `test_clean_knots_emit_no_warning`); source-level aggregate-warning block at `diff_diff/continuous_did_bspline.py:150-187`. - [x] **Edge cases**: all-same-dose (rank-deficient design, recovers only intercept = `ATT^{glob}`, ACRT = 0 everywhere), single-treated-unit (insufficient for OLS, raises `ValueError` "No valid"), discrete-treatment (detected and warned, saturated regression deferred), rank-deficiency per cell (cell skipped under `rank_deficient_action="silent"` / `"warn"`), balanced-panel-required (matches R `contdid` v0.1.0). Locked in `TestEdgeCasesMethodology` (2 methodology tests) + rank-deficient unit tests in `test_continuous_did.py`. - [x] **Anticipation-aware not-yet-treated control mask**: when `anticipation > 0`, the not-yet-treated control mask uses `G > t + anticipation` (not just `G > t`) to exclude cohorts in the anticipation window from controls. When `anticipation=0` (default), behavior is unchanged. CHANGELOG `[3.0.x]`-era fix; locked in `test_continuous_did.py`. @@ -1464,7 +1464,7 @@ Goldens at `benchmarks/data/r_conleyreg_conley_golden.json`; generator at `bench **Outstanding Concerns — cross-estimator survey coverage boundary** (intentional, fail-closed `NotImplementedError` deferrals, not bugs; line refs current as of this review — re-grep before relying on them): - **Conley + survey_design** (open methodological question — no canonical weighted spatial-HAC under probability sampling): `conley.py:298`, `linalg.py:1367` / `linalg.py:3501`, `spillover.py:3246`. -- **Replicate-weight designs** (use TSL strata/PSU/FPC, or `n_bootstrap=0`): `synthetic_did.py:437`, `continuous_did.py:1408`, `staggered_triple_diff.py:701`, `spillover.py:2400`, `bacon.py:526`, `had.py:1751` (+ HAD pretests), `wooldridge.py:89`, `trop.py:444`, `staggered.py:2228` (CS bootstrap), `efficient_did.py:1175`, `chaisemartin_dhaultfoeuille.py:2810`. +- **Replicate-weight designs** (use TSL strata/PSU/FPC, or `n_bootstrap=0`): `synthetic_did.py:437`, `continuous_did.py:2181`, `staggered_triple_diff.py:701`, `spillover.py:2400`, `bacon.py:526`, `had.py:1751` (+ HAD pretests), `wooldridge.py:89`, `trop.py:444`, `staggered.py:2228` (CS bootstrap), `efficient_did.py:1175`, `chaisemartin_dhaultfoeuille.py:2810`. - **Survey + non-HC1 vcov** — HC2 / HC2-BM / classical **explicitly raise `NotImplementedError`** under `survey_design=` (a fail-closed guard: the survey TSL / replicate-refit variance would otherwise silently discard the requested sandwich family): `stacked_did.py:428`, `sun_abraham.py:751`, `wooldridge.py:702`, `twfe.py:252`. - **Survey + user `cluster=`** — **explicitly raises** (a fail-closed guard, not a silent drop): the survey TSL / replicate-refit variance would otherwise ignore `cluster=`, so the combination is rejected at construction: `efficient_did.py:518`, `staggered.py:1719`, `imputation.py:314`, `two_stage.py:1439`, `triple_diff.py:674`. - **SyntheticControl** — no survey support yet (`synthetic_control.py:335`). diff --git a/README.md b/README.md index 2e8728dab..a020d47d3 100644 --- a/README.md +++ b/README.md @@ -91,7 +91,7 @@ For rigorous DiD analysis, follow these 8 steps. Skipping diagnostic steps produ 4. **Choose estimator** - staggered adoption -> CS/SA/BJS (NOT plain TWFE); few treated units -> SDiD; factor confounding -> TROP; simple 2x2 -> DiD. Run `BaconDecomposition` to diagnose TWFE bias. 5. **Estimate** - `estimator.fit(data, ...)`. Always print the cluster count first and choose inference method based on the result (cluster-robust if >= 50 clusters, wild bootstrap if fewer). 6. **Sensitivity analysis** - `compute_honest_did(results)` for bounds under PT violations (MultiPeriodDiD, CS, or dCDH natively; a StackedDiD `results.aggregate('event_study')` container also admits - needs `kappa_pre >= 2`), `run_all_placebo_tests()` for 2x2 falsification, specification comparisons for staggered designs. -7. **Heterogeneity** - CS: `results.aggregate('group')`/`'event_study'` (post-fit, no refit); SA: `results.event_study_effects` / `to_dataframe(level='cohort')`; Stacked: `results.aggregate('event_study')`/`'simple'` post-fit views (surface always computed since 3.9); EDiD: `results.aggregate(...)` post-fit from retained EIFs (3.9); ImputationDiD/TwoStageDiD: `results.aggregate(...)` post-fit from panel-backed kits (3.9); subgroup re-estimation. +7. **Heterogeneity** - CS: `results.aggregate('group')`/`'event_study'` (post-fit, no refit); SA: `results.event_study_effects` / `to_dataframe(level='cohort')`; Stacked: `results.aggregate('event_study')`/`'simple'` post-fit views (surface always computed since 3.9); EDiD: `results.aggregate(...)` post-fit from retained EIFs (3.9); ImputationDiD/TwoStageDiD: `results.aggregate(...)` post-fit from panel-backed kits (3.9); ContinuousDiD: `results.aggregate('dose'/'simple'/'event_study')` post-fit (3.9; dose/simple are views, event_study recomputes); subgroup re-estimation. 8. **Robustness** - compare 2-3 estimators (CS vs SA vs BJS), report with and without covariates (shows whether conditioning drives identification), present pre-trends and sensitivity bounds. Full guide: `diff_diff.get_llm_guide("practitioner")`. diff --git a/TODO.md b/TODO.md index fc72170d8..9cf78393b 100644 --- a/TODO.md +++ b/TODO.md @@ -23,8 +23,9 @@ Related tracking surfaces: |-------|----------|--------|--------|----------| | EfficientDiD `aggregate()` on bootstrapped fits fails closed (M-023); wiring `BootstrapReplaySpec` (or retaining the n_bootstrap x n_gt draw matrix materialized at fit) would enable exact post-fit replay of percentile inference | `diff_diff/efficient_did_results.py`, `diff_diff/aggregation.py` | 2(b) PR-3a | Mid | Low | | ImputationDiD/TwoStageDiD `aggregate()` on bootstrapped fits fails closed (M-021/M-022); ImputationDiD's per-target psi machinery makes seeded replay tractable (the panel-backed kit retains everything the psi precompute reads), TwoStageDiD's per-level GMM scores are function-locals and would need retention | `diff_diff/imputation_results.py`, `diff_diff/two_stage_results.py`, `diff_diff/aggregation.py` | 2(b) PR-3b | Mid | Low | +| ContinuousDiD `aggregate('event_study')` on bootstrapped fits fails closed (M-025); a seeded post-fit bootstrap-ES replay is tractable - the multiplier draws are seeded (`np.random.default_rng(self.seed)`) - but needs the FULL per-cell `_bootstrap_info` (bread/ee_treated/Psi_eval/dPsi_*/beta_pred) the pruned kit deliberately drops, so shipping it means a kit-payload change with its own memory contract | `diff_diff/continuous_did_aggregation.py`, `diff_diff/continuous_did_results.py` | 2(b) PR-3c | Mid | Low | | diagnostic_report's ES-gated checks read the raw `event_study_effects` field, which post-fit `results.aggregate()` never populates - their remediation strings steer users to the deprecated fit-time kwarg (qualified "deprecated but functional until 4.0" since 2(b) PR-3b); teach the checks to consume a post-fit container (or recompute via the kit) before 4.0 removes the kwarg | `diff_diff/diagnostic_report.py` | 2(b) PR-3b | Mid | Medium | -| EfficientDiD and ImputationDiD are the outstanding M-092 event-study df-provenance holes: the container's per-row df is all-NaN even on survey fits where a finite `_survey_df` governed the p-values (the container-level scalar `df_survey` IS exposed - the hole is the PER-ROW column only; no event_study_df/df_inference field; pre-existing, NOT a regression of the M-023 PR - today's builder output is identical). The kits now retain the scalar (ImputationDiD's since 2(b) PR-3b - same shape: scalar `df_survey` exposed, per-row column all-NaN, identical to its fit-time surface); threading it into the per-row channel is a contained follow-up | `diff_diff/efficient_did_results.py`, `diff_diff/imputation_results.py`, `diff_diff/results_base.py` | 2(b) PR-3a | Quick | Low | +| EfficientDiD, ImputationDiD and ContinuousDiD are the outstanding M-092 event-study df-provenance holes: the container's per-row df is all-NaN even on survey fits where a finite `_survey_df` governed the p-values (the container-level scalar `df_survey` IS exposed - the hole is the PER-ROW column only; no event_study_df/df_inference field; pre-existing, NOT a regression of the M-023 PR - today's builder output is identical). The kits now retain the scalar (ImputationDiD's since 2(b) PR-3b, ContinuousDiD's since 2(b) PR-3c - same shape: scalar `df_survey` exposed, per-row column all-NaN, identical to each fit-time surface); threading it into the per-row channel is a contained follow-up | `diff_diff/efficient_did_results.py`, `diff_diff/imputation_results.py`, `diff_diff/continuous_did_results.py`, `diff_diff/results_base.py` | 2(b) PR-3a | Quick | Low | | practitioner `step_name="heterogeneity"` producer-side collisions: three OTHER estimators' advice steps reuse the key with non-heterogeneity labels (`:975` ContinuousDiD dose-response, `:1022` Triple placebo-group, `:1413` LPDiD WAS arrays), so DiagnosticReport's heterogeneity completion silently drops that unrelated advice from `next_steps` via `_filter_steps` - the same latent collision fixed for StackedDiD in M-024 (renamed to `sub_experiment_balance`). Renaming these changes those estimators' report output; audit + rename with per-estimator pins. | `diff_diff/practitioner.py` | 2(b) PR-2 review R9 | Quick | Low | | PreTrendsPower `violation='linear'` on CS `base_period='varying'` input targets the wrong alternative: `δ_pre = M · \|t\|` assumes level coefficients against a common reference, but varying-base pre-treatment effects are consecutive-period comparisons (constant increments under a linear trend). Both CS-sourced routes now WARN (REGISTRY PreTrendsPower Note), and universal-base GAPPED grids fail closed via the `reference_event_times` common-reference guard; what remains is the varying-base resolution - either transforming the violation vector through each coefficient's actual base mapping (needs per-horizon base provenance) or requiring `base_period='universal'` for the linear benchmark - a per-estimator methodology decision with a hand-calculated linear-violation gate | `diff_diff/pretrends.py` | 2(b) PR-1 R5 | Mid | Medium | | `run_placebo_test`/`run_all_placebo_tests` expose ONE overloaded `time` param (calendar for timing/group tests, 0/1 post dummy forwarded into the M-137/M-138 callees - works only because a two-period 0/1 calendar column is both); redesign the wrapper signatures (dual `time=`/`post=` params or a split) so each callee gets an honestly-named argument | `diff_diff/diagnostics.py` | 2(c)-ii PR-B | Mid | Medium | diff --git a/benchmarks/speed_review/bench_dose_response.py b/benchmarks/speed_review/bench_dose_response.py index 9dd387651..bbfbecb1d 100644 --- a/benchmarks/speed_review/bench_dose_response.py +++ b/benchmarks/speed_review/bench_dose_response.py @@ -1,9 +1,10 @@ """ Scenario 6: Pricing dose-response with ContinuousDiD cubic spline. -Chains: CDiD fit with aggregate='dose' (overall ATT + ACRT + dose-response -curves + bootstrap 199) -> dataframe extraction -> event-study pre-trend -> -binarized-DiD comparison -> spline sensitivity (degree=1, num_knots=2). +Chains: CDiD fit (overall ATT + ACRT + dose-response curves are always +computed; bootstrap 199) -> dataframe extraction -> post-fit event-study +pre-trend (results.aggregate('event_study'), row M-025) -> binarized-DiD +comparison -> spline sensitivity (degree=1, num_knots=2). Data shape: 500 stores x 6 quarterly periods, 1 cohort at period 3, log-normal dose. Matches Tutorial 14 scaled from 200 to 500 units. @@ -50,15 +51,14 @@ def cdid_cubic_fit(): cdid = ContinuousDiD( degree=3, num_knots=1, n_bootstrap=199, seed=123, ) - results["cubic"] = cdid.fit(**fit_kwargs, aggregate="dose") + results["cubic"] = cdid.fit(**fit_kwargs) def extract_curves(): - # The cubic fit used aggregate="dose", so only dose-response and - # group-time levels are available on the result. Event-study is - # extracted separately in the dedicated pretrend phase below. - # NB: ContinuousDiD uses 'eventstudy' for fit(aggregate=...) but - # 'event_study' for to_dataframe(level=...). Two different - # spellings within one estimator - flagged in performance-plan.md. + # Dose-response and group-time levels are always available on the + # result (fit computes them unconditionally). The event study is + # produced post-fit in the dedicated pretrend phase below via + # results.aggregate('event_study') - the unified spelling; the + # fit-time 'eventstudy' kwarg is deprecated (row M-025). r = results["cubic"] out = {} for level in ("dose_response", "group_time"): @@ -66,12 +66,13 @@ def extract_curves(): results["curves"] = out def cdid_event_study(): + # Post-fit route (M-025): the phase still times the event-study + # computation - the analytical fit plus the kit recompute. cdid = ContinuousDiD( degree=3, num_knots=1, n_bootstrap=0, seed=123, ) - results["event_study"] = cdid.fit( - **fit_kwargs, aggregate="eventstudy", - ) + res = cdid.fit(**fit_kwargs) + results["event_study"] = res.aggregate("event_study") def binarized_comparison(): # Derive post from the actual first_treat cohort in the data so @@ -92,13 +93,13 @@ def spline_sensitivity_linear(): cdid = ContinuousDiD( degree=1, num_knots=0, n_bootstrap=199, seed=123, ) - results["linear"] = cdid.fit(**fit_kwargs, aggregate="dose") + results["linear"] = cdid.fit(**fit_kwargs) def spline_sensitivity_more_knots(): cdid = ContinuousDiD( degree=3, num_knots=2, n_bootstrap=199, seed=123, ) - results["many_knots"] = cdid.fit(**fit_kwargs, aggregate="dose") + results["many_knots"] = cdid.fit(**fit_kwargs) phases = [ ("1_cdid_cubic_spline_bootstrap199", cdid_cubic_fit), diff --git a/diff_diff/aggregation.py b/diff_diff/aggregation.py index e9f55599e..828cb8427 100644 --- a/diff_diff/aggregation.py +++ b/diff_diff/aggregation.py @@ -274,7 +274,12 @@ def to_dataframe(self) -> pd.DataFrame: Rows are ordered by ``label`` when the labels are homogeneously sortable, and in producer order otherwise (mixed-type cohort labels - cannot be ordered without raising). + cannot be ordered without raising). Heterogeneous-``target`` + containers (ContinuousDiD's att/acrt, row M-025) order by + FIRST-APPEARANCE target blocks instead - producer order, NOT + lexicographic, which would put ``"acrt"`` before ``"att"`` - with + labels ascending within each block under the same sortability + guard. """ data: Dict[str, Any] = { "level": self.level, @@ -291,7 +296,18 @@ def to_dataframe(self) -> pd.DataFrame: "df": self.df, } frame = pd.DataFrame(data, columns=list(AGGREGATION_SCHEMA)) - if len(frame) > 1 and _sortable(self.label): + distinct_targets = list(dict.fromkeys(self.target)) + if len(distinct_targets) > 1 and len(frame) > 1: + rank = {t: i for i, t in enumerate(distinct_targets)} + target_rank = np.array([rank[t] for t in self.target]) + if _sortable(self.label): + # lexsort: LAST key is primary - target blocks first, + # labels ascending within each block. + order = np.lexsort((self.label, target_rank)) + else: + order = np.argsort(target_rank, kind="stable") + frame = frame.iloc[order].reset_index(drop=True) + elif len(frame) > 1 and _sortable(self.label): order = np.argsort(self.label, kind="stable") frame = frame.iloc[order].reset_index(drop=True) return frame @@ -348,16 +364,36 @@ def summary(self, alpha: Optional[float] = None) -> str: return "\n".join(lines) n_label = "n" if self.n_kind is None else f"n[{self.n_kind}]" - lines.append(f"{'label':>14} {'ATT':>11} {'SE':>10} {'t':>8} {'p':>8} {n_label:>10}") - lines.append("-" * 64) frame = self.to_dataframe() - for _, row in frame.iterrows(): - n_disp = "" if not np.isfinite(row["n"]) else f"{row['n']:.0f}" + if len(dict.fromkeys(self.target)) > 1: + # Heterogeneous targets (ContinuousDiD's att/acrt, row M-025): + # a target column disambiguates the duplicate labels and the + # estimate heading goes neutral - the hard-coded 'ATT' would + # mislabel every acrt row. Uniform-target containers render + # exactly as before (byte-stable). lines.append( - f"{str(row['label']):>14} {row['att']:>11.4f} {row['se']:>10.4f} " - f"{row['t_stat']:>8.3f} {row['p_value']:>8.4f} {n_disp:>10}" + f"{'label':>14} {'target':>8} {'estimate':>11} {'SE':>10} " + f"{'t':>8} {'p':>8} {n_label:>10}" ) - lines.append("-" * 64) + lines.append("-" * 73) + for _, row in frame.iterrows(): + n_disp = "" if not np.isfinite(row["n"]) else f"{row['n']:.0f}" + lines.append( + f"{str(row['label']):>14} {str(row['target']):>8} " + f"{row['att']:>11.4f} {row['se']:>10.4f} " + f"{row['t_stat']:>8.3f} {row['p_value']:>8.4f} {n_disp:>10}" + ) + lines.append("-" * 73) + else: + lines.append(f"{'label':>14} {'ATT':>11} {'SE':>10} {'t':>8} {'p':>8} {n_label:>10}") + lines.append("-" * 64) + for _, row in frame.iterrows(): + n_disp = "" if not np.isfinite(row["n"]) else f"{row['n']:.0f}" + lines.append( + f"{str(row['label']):>14} {row['att']:>11.4f} {row['se']:>10.4f} " + f"{row['t_stat']:>8.3f} {row['p_value']:>8.4f} {n_disp:>10}" + ) + lines.append("-" * 64) lines.append(f"Confidence intervals at alpha={self.alpha}.") if self.weight is None: lines.append("Per-row aggregation weights are not defined for this level.") @@ -389,6 +425,21 @@ class AggregationKit: by design and the exclusion above applies to everything OUTSIDE the enumerated bookkeeping payload. + PRUNED-PAYLOAD VARIANT (ContinuousDiD [M-025]): its event-study + recompute needs per-(g, t) IF INGREDIENTS rather than a per-unit EIF + dict or the panel - ``bookkeeping`` retains a pruned + ``_bootstrap_info`` subset (treated/control indices, ``delta_y_treated``, + ``ee_control``, masses, the covariate-path ``if_att_glob``), O(n_treated + + n_control) per cell, plus unit-level arrays and - on survey fits - + the PANEL-LEVEL ``ResolvedSurveyDesign`` (the recompute performs the + unit collapse itself; on replicate designs this carries the + (n_obs x R) replicate matrix). The K-dimensional spline machinery + (bread, ``ee_treated``, ``Psi_eval``, ``dPsi_*``) is NOT retained; + bootstrap fits retain scalars only (their event-study route fails + closed). ``influence`` is empty by design here too. The ``simple`` / + ``dose`` levels are pure views over stored public results fields and + read no kit at all. + Attributes ---------- bookkeeping : dict diff --git a/diff_diff/continuous_did.py b/diff_diff/continuous_did.py index d361f80d0..38284cc3b 100644 --- a/diff_diff/continuous_did.py +++ b/diff_diff/continuous_did.py @@ -9,6 +9,7 @@ inference. """ +import dataclasses import warnings from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple @@ -16,11 +17,13 @@ import pandas as pd from diff_diff._base import BaseEstimator -from diff_diff._deprecation import warn_deprecated_kwarg +from diff_diff._deprecation import NOT_SUPPLIED, warn_deprecated_kwarg +from diff_diff.aggregation import AggregationKit from diff_diff.bootstrap_utils import ( compute_effect_bootstrap_stats, generate_bootstrap_weights_batch, ) +from diff_diff.continuous_did_aggregation import _ContinuousDiDAggregationMixin from diff_diff.continuous_did_bspline import ( SATURATED_TOL, bspline_derivative_design_matrix, @@ -51,7 +54,111 @@ __all__ = ["ContinuousDiD", "ContinuousDiDResults", "DoseResponseCurve"] -class ContinuousDiD(BaseEstimator): +#: Pruned per-cell payload keys the post-fit event-study recompute reads +#: (row M-025). Exactly the ``_bootstrap_info`` subset +#: ``_compute_event_study_inference`` consumes - the K-dimensional spline +#: machinery (bread, ee_treated, Psi_eval, dPsi_*, beta_pred) is +#: deliberately NOT retained and dies with fit(). ``w_treated``/ +#: ``w_control``/``w_treated_arr`` are copied only when present because +#: the consumer's survey-mass branch keys on ``"w_treated" in b_info``. +_ES_PAYLOAD_KEYS = ( + "treated_indices", + "control_indices", + "n_treated", + "n_control", + "att_glob", + "mu_0", + "delta_y_treated", + "ee_control", + "w_treated", + "w_control", + "w_treated_arr", +) + + +def _build_continuous_aggregation_kit( + estimator: "ContinuousDiD", + gt_results: Dict[Tuple, Dict], + gt_bootstrap_info: Dict[Tuple, Dict], + precomp: Dict[str, Any], + resolved_survey: Optional["ResolvedSurveyDesign"], + has_post_cells: bool, + survey_df: Optional[int], + survey_metadata: Optional[Any], +) -> AggregationKit: + """Build the post-fit aggregation kit for ContinuousDiD (row M-025). + + ``influence`` is empty BY DESIGN: the event-study recompute reads the + pruned per-cell payload + unit-level arrays in ``bookkeeping``, not a + per-unit EIF dict on the kit's influence contract; the ``simple`` / + ``dose`` levels are pure views over stored public results fields and + never read the kit at all. + + On bootstrap fits (``n_bootstrap > 0``) the kit is SCALARS-ONLY: + ``aggregate('event_study')`` fails closed before reading any payload + and the views never read the kit, so a populated payload there would + be pure dead retention. The scalars still distinguish the + bootstrap-NotImplementedError gate from the legacy-pickle no-kit + ValueError. + """ + is_bootstrap = estimator.n_bootstrap > 0 + gt_summary: Dict[Tuple, Dict[str, Any]] = {} + gt_es_payload: Dict[Tuple, Dict[str, Any]] = {} + if not is_bootstrap: + for gt, r in gt_results.items(): + gt_summary[gt] = { + "att_glob": float(r["att_glob"]), + "n_treated": int(r["n_treated"]), + } + b_info = gt_bootstrap_info.get(gt, {}) + if not b_info: + gt_es_payload[gt] = {} + continue + pruned = {k: b_info[k] for k in _ES_PAYLOAD_KEYS if k in b_info} + cov_if = b_info.get("cov_if") + pruned["cov_if"] = ( + { + "cell_indices": cov_if["cell_indices"], + "if_att_glob": cov_if["if_att_glob"], + } + if cov_if is not None + else None + ) + gt_es_payload[gt] = pruned + bookkeeping: Dict[str, Any] = { + "gt_summary": gt_summary, + "gt_es_payload": gt_es_payload, + "n_units": None if is_bootstrap else precomp["n_units"], + "unit_cohorts": None if is_bootstrap else precomp["unit_cohorts"], + "unit_survey_weights": (None if is_bootstrap else precomp.get("unit_survey_weights")), + "unit_first_panel_row": (None if is_bootstrap else precomp["unit_first_panel_row"]), + # PANEL-LEVEL design ref (locked decision: the unit-level collapse + # stays inside the verbatim recompute body; on replicate designs + # this carries the (n_obs x R) replicate matrix - documented in + # the REGISTRY memory contract). + "resolved_survey": None if is_bootstrap else resolved_survey, + "has_post_cells": has_post_cells, + "survey_df": survey_df, + "n_bootstrap": estimator.n_bootstrap, + "base_period": estimator.base_period, + # Fit-final COPY - the ES carrier metadata source; never the + # mutable public field (post-fit mutation of replicate_method / + # df_survey must not reach post-fit provenance). + "survey_metadata": ( + dataclasses.replace(survey_metadata) if survey_metadata is not None else None + ), + } + return AggregationKit( + bookkeeping=bookkeeping, + influence={}, + alpha=estimator.alpha, + anticipation=estimator.anticipation, + cband=False, + bootstrap=None, + ) + + +class ContinuousDiD(_ContinuousDiDAggregationMixin, BaseEstimator): """ Continuous Difference-in-Differences estimator. @@ -142,8 +249,9 @@ class ContinuousDiD(BaseEstimator): >>> est = ContinuousDiD(n_bootstrap=199, seed=42) >>> results = est.fit(data, outcome="outcome", unit="unit", ... time="period", first_treat="first_treat", - ... dose="dose", aggregate="dose") + ... dose="dose") >>> results.overall_att # doctest: +SKIP + >>> results.aggregate("dose") # doctest: +SKIP """ _VALID_CONTROL_GROUPS = {"never_treated", "not_yet_treated", "lowest_dose"} @@ -261,7 +369,7 @@ def fit( time: str, first_treat: str, dose: str, - aggregate: Optional[str] = None, + aggregate: Any = NOT_SUPPLIED, survey_design: Optional["SurveyDesign"] = None, covariates: Optional[List[str]] = None, ) -> ContinuousDiDResults: @@ -283,8 +391,16 @@ def fit( dose : str Continuous dose column. aggregate : str, optional - ``"dose"`` for dose-response aggregation, ``"eventstudy"`` for - binarized event study. + DEPRECATED (row M-025, removed in 4.0) - aggregate as a + post-fit step instead: ``results.aggregate('event_study')`` + for the binarized event study (underscored - the + ``"eventstudy"`` spelling dies with this parameter), or + ``results.aggregate('dose')`` / ``results.aggregate('simple')`` + views. The dose-response curves and overall ATT/ACRT are + always computed by ``fit()``, so ``aggregate="dose"`` was + already a no-op. Supplying ANY value (including ``None``) + warns ``FutureWarning``; supplied values still run the legacy + routing unchanged until 4.0. survey_design : SurveyDesign, optional Survey design specification for design-based inference. Supports weighted estimation and Taylor series linearization @@ -300,6 +416,31 @@ def fit( ------- ContinuousDiDResults """ + # M-025 deprecation shim: a plain fit() never warns; supplying + # aggregate= with ANY value (None included) warns once, then the + # legacy routing below runs unchanged - "eventstudy" still + # computes the fit-time surface and invalid strings still reach + # the pre-existing ValueError. Only the SENTINEL normalizes to + # None (it would otherwise fail the _VALID_AGGREGATES check on + # every plain fit). The post-fit successor validates its own + # (unified) vocabulary. + if aggregate is not NOT_SUPPLIED: + warnings.warn( + "ContinuousDiD.fit(aggregate=) is deprecated and will be " + "removed in 4.0. Fit once, then aggregate as a post-fit " + "step: results = ContinuousDiD().fit(...); " + "results.aggregate('event_study') (note the underscore - " + "the 'eventstudy' spelling dies with this parameter) / " + ".aggregate('dose') / .aggregate('simple'). The " + "dose-response curves and overall ATT/ACRT are always " + "computed by fit(), so aggregate='dose' was already " + "redundant.", + FutureWarning, + stacklevel=2, + ) + else: + aggregate = None + # 1. Validate & prepare _VALID_AGGREGATES = (None, "dose", "eventstudy") if aggregate not in _VALID_AGGREGATES: @@ -969,116 +1110,24 @@ def fit( acrt_d_ci_lower[idx] = ci[0] acrt_d_ci_upper[idx] = ci[1] - # Event study analytical SEs + # Event study analytical SEs - the body lives in + # continuous_did_aggregation._compute_event_study_inference, + # shared verbatim with the post-fit aggregate('event_study') + # recompute (row M-025). fit passes its full gt_results / + # gt_bootstrap_info locals; the method reads only the + # kit-compatible key subset. if event_study_effects is not None: - n_units = precomp["n_units"] - unit_sw = precomp.get("unit_survey_weights") - - # Build unit-level ResolvedSurveyDesign once (reused per bin) - unit_resolved_es = None - if resolved_survey is not None: - row_idx = precomp["unit_first_panel_row"] - unit_resolved_es = resolved_survey.subset_to_units_by_row_idx( - row_idx, unit_weights=precomp.get("unit_survey_weights") - ) - - for e_val, info_e in event_study_effects.items(): - # Collect (g,t) cells for this event-time bin - e_gts = [gt for gt in gt_results if gt[1] - gt[0] == e_val] - if not e_gts: - continue - # Weights within this bin: survey-weighted mass or n_treated - if unit_sw is not None: - unit_cohorts = precomp["unit_cohorts"] - ns = np.array( - [float(np.sum(unit_sw[unit_cohorts == gt[0]])) for gt in e_gts], - dtype=float, - ) - else: - ns = np.array( - [gt_results[gt]["n_treated"] for gt in e_gts], - dtype=float, - ) - total_n = ns.sum() - if total_n == 0: - continue - ws = ns / total_n - - # Build per-unit IF for this event-time bin - if_es = np.zeros(n_units) - for idx_cell, gt in enumerate(e_gts): - b_info = gt_bootstrap_info.get(gt, {}) - if not b_info: - continue - w = ws[idx_cell] - # Covariate path: the binarized event-study effect is - # att_glob, whose per-unit cell IF is precomputed. - cov_if = b_info.get("cov_if") - if cov_if is not None: - np.add.at( - if_es, - cov_if["cell_indices"], - w * cov_if["if_att_glob"], - ) - continue - treated_idx = b_info["treated_indices"] - control_idx = b_info["control_indices"] - n_t = b_info["n_treated"] - n_c = b_info["n_control"] - # Use survey-weighted masses when available - if "w_treated" in b_info: - n_t = b_info["w_treated"] - n_c = b_info["w_control"] - n_total_gt = n_t + n_c - p_1 = n_t / n_total_gt - p_0 = n_c / n_total_gt - att_glob_gt = b_info["att_glob"] - mu_0 = b_info["mu_0"] - delta_y_treated = b_info["delta_y_treated"] - ee_control = b_info["ee_control"] - sw_treated = b_info.get("w_treated_arr") - - for k, uid in enumerate(treated_idx): - score_k = delta_y_treated[k] - att_glob_gt - mu_0 - if sw_treated is not None: - score_k = sw_treated[k] * score_k - if_es[uid] += w * score_k / p_1 / n_total_gt - for k, uid in enumerate(control_idx): - if_es[uid] -= w * ee_control[k] / p_0 / n_total_gt - - # Compute SE: survey-aware TSL or standard sqrt(sum(IF^2)) - if unit_resolved_es is not None: - if unit_resolved_es.uses_replicate_variance: - from diff_diff.survey import compute_replicate_if_variance - - # Score-scale: psi = w * if_es (matches TSL bread) - psi_es = unit_resolved_es.weights * if_es - variance, _nv = compute_replicate_if_variance( - psi_es, unit_resolved_es - ) - es_se = ( - float(np.sqrt(max(variance, 0.0))) - if np.isfinite(variance) - else np.nan - ) - else: - X_ones_es = np.ones((n_units, 1)) - tsl_scale_es = float(unit_resolved_es.weights.sum()) - if_es_tsl = if_es * tsl_scale_es - vcov_es = compute_survey_vcov( - X_ones_es, if_es_tsl, unit_resolved_es - ) - es_se = float(np.sqrt(np.abs(vcov_es[0, 0]))) - else: - es_se = float(np.sqrt(np.sum(if_es**2))) - - t_stat, p_val, ci_es = safe_inference( - info_e["effect"], es_se, self.alpha, df=_survey_df - ) - info_e["se"] = es_se - info_e["t_stat"] = t_stat - info_e["p_value"] = p_val - info_e["conf_int"] = ci_es + self._compute_event_study_inference( + event_study_effects, + gt_summary=gt_results, + gt_es_payload=gt_bootstrap_info, + n_units=precomp["n_units"], + unit_cohorts=precomp["unit_cohorts"], + unit_survey_weights=precomp.get("unit_survey_weights"), + unit_first_panel_row=precomp["unit_first_panel_row"], + resolved_survey=resolved_survey, + survey_df=_survey_df, + ) # 6. Assemble results dose_response_att = DoseResponseCurve( @@ -1123,7 +1172,7 @@ def fit( n_control_units_out = n_control reference_dose_out = None - return ContinuousDiDResults( + results = ContinuousDiDResults( dose_response_att=dose_response_att, dose_response_acrt=dose_response_acrt, overall_att=overall_att, @@ -1163,6 +1212,19 @@ def fit( event_study_effects=event_study_effects, survey_metadata=survey_metadata, ) + # Post-fit aggregation kit (row M-025): attached on EVERY fit; + # scalars-only on bootstrap fits (the ES route fails closed there). + results._aggregation_kit = _build_continuous_aggregation_kit( + self, + gt_results=gt_results, + gt_bootstrap_info=gt_bootstrap_info, + precomp=precomp, + resolved_survey=resolved_survey, + has_post_cells=len(post_gt) > 0, + survey_df=_survey_df, + survey_metadata=survey_metadata, + ) + return results # ------------------------------------------------------------------ # Internal helpers @@ -1880,49 +1942,9 @@ def _deriv(z: np.ndarray) -> np.ndarray: "_bootstrap_info": bootstrap_info, } - def _aggregate_event_study( - self, - gt_results: Dict[Tuple, Dict], - gt_bootstrap_info: Dict[Tuple, Dict] = None, - unit_survey_weights: Optional[np.ndarray] = None, - unit_cohorts: Optional[np.ndarray] = None, - anticipation: int = 0, - ) -> Dict[int, Dict[str, Any]]: - """Aggregate binarized ATT_glob by relative period.""" - effects_by_e: Dict[int, List[Tuple[float, float, Tuple]]] = {} - - for (g, t), r in gt_results.items(): - e = t - g - if anticipation > 0 and e < -anticipation: - continue - if e not in effects_by_e: - effects_by_e[e] = [] - # Compute weight for this (g,t) cell - if unit_survey_weights is not None and unit_cohorts is not None: - # Survey-weighted: sum of survey weights for treated units in group g - g_mask = unit_cohorts == g - cell_weight = float(np.sum(unit_survey_weights[g_mask])) - else: - cell_weight = float(r["n_treated"]) - effects_by_e[e].append((r["att_glob"], cell_weight, (g, t))) - - result = {} - for e, entries in sorted(effects_by_e.items()): - effects = np.array([x[0] for x in entries]) - weights = np.array([x[1] for x in entries]) - if np.sum(weights) > 0: - w = weights / np.sum(weights) - agg = float(np.sum(w * effects)) - else: - agg = np.nan - result[e] = { - "effect": agg, - "se": np.nan, - "t_stat": np.nan, - "p_value": np.nan, - "conf_int": (np.nan, np.nan), - } - return result + # _aggregate_event_study moved verbatim to + # continuous_did_aggregation._ContinuousDiDAggregationMixin (row M-025); + # fit-time call sites resolve via the mixin base. def _compute_analytical_se( self, diff --git a/diff_diff/continuous_did_aggregation.py b/diff_diff/continuous_did_aggregation.py new file mode 100644 index 000000000..90860fcba --- /dev/null +++ b/diff_diff/continuous_did_aggregation.py @@ -0,0 +1,213 @@ +""" +Shared event-study aggregation machinery for ContinuousDiD (row M-025). + +Leaf module: hosts ``_ContinuousDiDAggregationMixin``, the binarized +event-study aggregation shared between fit-time ``fit(aggregate= +"eventstudy")`` (deprecated) and the post-fit +``ContinuousDiDResults.aggregate("event_study")`` recompute. The import +DAG forces the split: ``continuous_did.py`` imports +``continuous_did_results.py``, so the results module cannot reach the +estimator module - both instead import this leaf (numpy + diff_diff.utils +/ diff_diff.survey only). + +The mixin's two methods are the verbatim bodies previously private to +``continuous_did.py``: + +- ``_aggregate_event_study`` - the per-relative-period reweighting of + per-(g, t) binarized ``att_glob`` values (cohort survey mass or + ``n_treated`` weights). +- ``_compute_event_study_inference`` - fit's formerly-inline analytical + ES-SE block: per-bin per-unit influence functions from the pruned + per-cell payload, with plain / TSL-survey / replicate-weight variance + and ``safe_inference`` at the fit's ``alpha`` and survey df. + +The post-fit route calls both on a throwaway +``_ContinuousKitAggregator`` host (``continuous_did_results.py``) whose +inputs come exclusively from the fit-built ``AggregationKit``; fit() +calls them with its locals. Neither method writes ``self`` state. +""" + +from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple + +import numpy as np + +from diff_diff.survey import compute_survey_vcov +from diff_diff.utils import safe_inference + +if TYPE_CHECKING: + from diff_diff.survey import ResolvedSurveyDesign + + +class _ContinuousDiDAggregationMixin: + """Binarized event-study aggregation shared by fit-time and post-fit. + + Host attribute contract (the throwaway kit aggregator sets exactly + this; ``ContinuousDiD`` carries it as an ``__init__`` param): + + - ``alpha`` - significance level for ``safe_inference``. + """ + + alpha: float + + def _aggregate_event_study( + self, + gt_results: Dict[Tuple, Dict], + gt_bootstrap_info: Optional[Dict[Tuple, Dict]] = None, + unit_survey_weights: Optional[np.ndarray] = None, + unit_cohorts: Optional[np.ndarray] = None, + anticipation: int = 0, + ) -> Dict[int, Dict[str, Any]]: + """Aggregate binarized ATT_glob by relative period.""" + effects_by_e: Dict[int, List[Tuple[float, float, Tuple]]] = {} + + for (g, t), r in gt_results.items(): + e = t - g + if anticipation > 0 and e < -anticipation: + continue + if e not in effects_by_e: + effects_by_e[e] = [] + # Compute weight for this (g,t) cell + if unit_survey_weights is not None and unit_cohorts is not None: + # Survey-weighted: sum of survey weights for treated units in group g + g_mask = unit_cohorts == g + cell_weight = float(np.sum(unit_survey_weights[g_mask])) + else: + cell_weight = float(r["n_treated"]) + effects_by_e[e].append((r["att_glob"], cell_weight, (g, t))) + + result = {} + for e, entries in sorted(effects_by_e.items()): + effects = np.array([x[0] for x in entries]) + weights = np.array([x[1] for x in entries]) + if np.sum(weights) > 0: + w = weights / np.sum(weights) + agg = float(np.sum(w * effects)) + else: + agg = np.nan + result[e] = { + "effect": agg, + "se": np.nan, + "t_stat": np.nan, + "p_value": np.nan, + "conf_int": (np.nan, np.nan), + } + return result + + def _compute_event_study_inference( + self, + event_study_effects: Dict[int, Dict[str, Any]], + gt_summary: Dict[Tuple, Dict], + gt_es_payload: Dict[Tuple, Dict], + n_units: int, + unit_cohorts: np.ndarray, + unit_survey_weights: Optional[np.ndarray], + unit_first_panel_row: np.ndarray, + resolved_survey: Optional["ResolvedSurveyDesign"], + survey_df: Optional[int], + ) -> None: + """Fill analytical se/t/p/CI on the binarized event-study rows. + + The verbatim body of fit's former inline "Event study analytical + SEs" block: mutates ``event_study_effects`` rows in place. The + fit-time caller passes its ``gt_results`` / ``gt_bootstrap_info`` + locals; the post-fit caller passes the kit's ``gt_summary`` / + pruned ``gt_es_payload`` (key-compatible by construction). Only + runs when the fit was analytical with post-treatment cells - the + callers own that gating. + """ + unit_sw = unit_survey_weights + + # Build unit-level ResolvedSurveyDesign once (reused per bin) + unit_resolved_es = None + if resolved_survey is not None: + row_idx = unit_first_panel_row + unit_resolved_es = resolved_survey.subset_to_units_by_row_idx( + row_idx, unit_weights=unit_survey_weights + ) + + for e_val, info_e in event_study_effects.items(): + # Collect (g,t) cells for this event-time bin + e_gts = [gt for gt in gt_summary if gt[1] - gt[0] == e_val] + if not e_gts: + continue + # Weights within this bin: survey-weighted mass or n_treated + if unit_sw is not None: + ns = np.array( + [float(np.sum(unit_sw[unit_cohorts == gt[0]])) for gt in e_gts], + dtype=float, + ) + else: + ns = np.array( + [gt_summary[gt]["n_treated"] for gt in e_gts], + dtype=float, + ) + total_n = ns.sum() + if total_n == 0: + continue + ws = ns / total_n + + # Build per-unit IF for this event-time bin + if_es = np.zeros(n_units) + for idx_cell, gt in enumerate(e_gts): + b_info = gt_es_payload.get(gt, {}) + if not b_info: + continue + w = ws[idx_cell] + # Covariate path: the binarized event-study effect is + # att_glob, whose per-unit cell IF is precomputed. + cov_if = b_info.get("cov_if") + if cov_if is not None: + np.add.at( + if_es, + cov_if["cell_indices"], + w * cov_if["if_att_glob"], + ) + continue + treated_idx = b_info["treated_indices"] + control_idx = b_info["control_indices"] + n_t = b_info["n_treated"] + n_c = b_info["n_control"] + # Use survey-weighted masses when available + if "w_treated" in b_info: + n_t = b_info["w_treated"] + n_c = b_info["w_control"] + n_total_gt = n_t + n_c + p_1 = n_t / n_total_gt + p_0 = n_c / n_total_gt + att_glob_gt = b_info["att_glob"] + mu_0 = b_info["mu_0"] + delta_y_treated = b_info["delta_y_treated"] + ee_control = b_info["ee_control"] + sw_treated = b_info.get("w_treated_arr") + + for k, uid in enumerate(treated_idx): + score_k = delta_y_treated[k] - att_glob_gt - mu_0 + if sw_treated is not None: + score_k = sw_treated[k] * score_k + if_es[uid] += w * score_k / p_1 / n_total_gt + for k, uid in enumerate(control_idx): + if_es[uid] -= w * ee_control[k] / p_0 / n_total_gt + + # Compute SE: survey-aware TSL or standard sqrt(sum(IF^2)) + if unit_resolved_es is not None: + if unit_resolved_es.uses_replicate_variance: + from diff_diff.survey import compute_replicate_if_variance + + # Score-scale: psi = w * if_es (matches TSL bread) + psi_es = unit_resolved_es.weights * if_es + variance, _nv = compute_replicate_if_variance(psi_es, unit_resolved_es) + es_se = float(np.sqrt(max(variance, 0.0))) if np.isfinite(variance) else np.nan + else: + X_ones_es = np.ones((n_units, 1)) + tsl_scale_es = float(unit_resolved_es.weights.sum()) + if_es_tsl = if_es * tsl_scale_es + vcov_es = compute_survey_vcov(X_ones_es, if_es_tsl, unit_resolved_es) + es_se = float(np.sqrt(np.abs(vcov_es[0, 0]))) + else: + es_se = float(np.sqrt(np.sum(if_es**2))) + + t_stat, p_val, ci_es = safe_inference(info_e["effect"], es_se, self.alpha, df=survey_df) + info_e["se"] = es_se + info_e["t_stat"] = t_stat + info_e["p_value"] = p_val + info_e["conf_int"] = ci_es diff --git a/diff_diff/continuous_did_results.py b/diff_diff/continuous_did_results.py index e609f32cf..6d85090f9 100644 --- a/diff_diff/continuous_did_results.py +++ b/diff_diff/continuous_did_results.py @@ -5,18 +5,34 @@ and aggregated estimation results. """ +import dataclasses from dataclasses import dataclass, field -from typing import Any, Dict, List, Optional, Tuple +from typing import Any, ClassVar, Dict, List, Optional, Tuple import numpy as np import pandas as pd +from diff_diff.aggregation import AggregationMixin, AggregationResult +from diff_diff.continuous_did_aggregation import _ContinuousDiDAggregationMixin from diff_diff.results import _format_survey_block, _get_significance_stars -from diff_diff.results_base import BaseResults +from diff_diff.results_base import BaseResults, build_event_study_surface +from diff_diff.utils import safe_inference __all__ = ["ContinuousDiDResults", "DoseResponseCurve"] +class _ContinuousKitAggregator(_ContinuousDiDAggregationMixin): + """Throwaway host for the post-fit event-study recompute (row M-025). + + A fresh instance runs each ``aggregate('event_study')`` call so the + recompute can never read or write estimator/results state. Sets + exactly the mixin's host-attribute contract: ``alpha``. + """ + + def __init__(self, alpha: float) -> None: + self.alpha = alpha + + @dataclass class DoseResponseCurve: """ @@ -79,7 +95,7 @@ def to_dataframe(self) -> pd.DataFrame: @dataclass -class ContinuousDiDResults(BaseResults): +class ContinuousDiDResults(BaseResults, AggregationMixin): """ Results from Continuous Difference-in-Differences estimation. @@ -159,6 +175,16 @@ class ContinuousDiDResults(BaseResults): event_study_effects: Optional[Dict[int, Dict[str, Any]]] = field(default=None) # Survey design metadata (SurveyMetadata instance from diff_diff.survey) survey_metadata: Optional[Any] = field(default=None) + # Post-fit aggregation kit (row M-025), attached by ContinuousDiD.fit(). + # Declared LAST for positional-__init__ compatibility. Only the + # 'event_study' recompute reads it; 'simple'/'dose' are views. + _aggregation_kit: Optional[Any] = field(default=None, repr=False, compare=False) + + # Post-fit aggregation routing (M-122 contract). ContinuousDiD's extra + # 'dose' level is documented in the ledger row and v4-design section 6; + # no level takes balance_e (the estimator has no balance_e machinery). + _AGGREGATE_SUPPORTED: ClassVar[Tuple[str, ...]] = ("simple", "event_study", "dose") + _AGGREGATE_BALANCE_E_TYPES: ClassVar[Tuple[str, ...]] = () # --- Inference-field aliases (balance/external-adapter compatibility) --- # ATT-side is the headline contract; ACRT remains accessible via overall_acrt_*. @@ -449,7 +475,14 @@ def to_dataframe(self, level: str = "dose_response") -> pd.DataFrame: return pd.DataFrame(rows) elif level == "event_study": if self.event_study_effects is None: - raise ValueError("Event study effects not computed. Use aggregate='eventstudy'.") + raise ValueError( + "Event study effects not computed. Call " + "results.aggregate('event_study') for the unified " + "post-fit container (on a bootstrapped fit, re-fit " + "with n_bootstrap=0 or use the deprecated fit-time " + "aggregate='eventstudy'); a result unpickled from an " + "older release must be re-fit with diff-diff >= 3.9." + ) rows = [] for rel_t, data in sorted(self.event_study_effects.items()): rows.append( @@ -478,3 +511,188 @@ def is_significant(self) -> bool: def significance_stars(self) -> str: """Significance stars for overall ATT.""" return _get_significance_stars(self.overall_att_p_value) + + # ------------------------------------------------------------------ + # Post-fit aggregation (row M-025, on the M-122 contract). + # MIXED architecture: 'simple' and 'dose' are pure VIEWS over stored + # public fields (the dCDH precedent - nothing recomputed, so they + # work on ANY fit including bootstrap fits and legacy pickles, + # relaying the stored inference verbatim); 'event_study' is a KIT + # RECOMPUTE (the EfficientDiD class) from the pruned per-cell IF + # payload, failing closed on bootstrap fits. + # ------------------------------------------------------------------ + + def _stored_inference_df(self) -> float: + """The df the STORED overall inference actually used (NaN = none). + + Bootstrap fits carry percentile p/CI - no df governs them. + Otherwise ``dose_response_att.df_survey`` is the stored provenance + channel for fit's ``_survey_df`` (the same value every + ``safe_inference`` call received); the replicate-undefined + 0-sentinel and ``None`` both report NaN in the df COLUMN, while + the view relays still pass the RAW stored value into their own + ``safe_inference`` derivations (``DoseResponseCurve.to_dataframe`` + parity). + """ + if self.n_bootstrap > 0: + return float("nan") + df_survey = self.dose_response_att.df_survey + if df_survey is not None and np.isfinite(df_survey) and df_survey > 0: + return float(df_survey) + return float("nan") + + def _aggregate_compute( + self, level: str, *, weights: Optional[str], balance_e: Optional[int] + ) -> Any: + if level == "simple": + # 2-row VIEW of the stored overall estimands: ContinuousDiD's + # headline parameters are the binarized overall ATT (ATT^{loc} + # under PT; equals ATT^{glob} under SPT) AND ACRT^{glob}, so the + # target column discriminates two "overall" rows (the + # container spec's dual-estimand case). Relays are strictly + # bit-exact - on bootstrap fits the stored quintet includes a + # FINITE safe_inference t beside the percentile p/CI and it + # relays through unchanged; only the df column is NaN there. + att_ci = self.overall_att_conf_int + acrt_ci = self.overall_acrt_conf_int + n_total = float(self.n_treated_units + self.n_control_units) + df_val = self._stored_inference_df() + return AggregationResult( + level="simple", + label=np.array(["overall", "overall"], dtype=object), + target=np.array(["att", "acrt"], dtype=object), + att=np.array([self.overall_att, self.overall_acrt], dtype=float), + se=np.array([self.overall_att_se, self.overall_acrt_se], dtype=float), + t_stat=np.array([self.overall_att_t_stat, self.overall_acrt_t_stat], dtype=float), + p_value=np.array( + [self.overall_att_p_value, self.overall_acrt_p_value], + dtype=float, + ), + conf_int_lower=np.array([att_ci[0], acrt_ci[0]], dtype=float), + conf_int_upper=np.array([att_ci[1], acrt_ci[1]], dtype=float), + # Treated and control unit sets are DISJOINT for this + # estimator (unlike Imputation/TwoStage), so the CS + # disjoint-total convention applies. + n=np.array([n_total, n_total], dtype=float), + df=np.array([df_val, df_val], dtype=float), + alpha=self.alpha, + n_kind="units", + weight=np.array([1.0, 1.0], dtype=float), + estimator="ContinuousDiD", + ) + + if level == "dose": + # 2N-row VIEW of the stored dose-response curves: att block + # then acrt block (first-appearance target order). t/p + # reproduce each DoseResponseCurve.to_dataframe exactly - + # including the bootstrap branch (stored p, NaN t) and the + # raw stored df_survey (0-sentinel included) fed to + # safe_inference on the analytical branch. + blocks = [] + for curve, target in ( + (self.dose_response_att, "att"), + (self.dose_response_acrt, "acrt"), + ): + n_grid = len(curve.effects) + if curve.n_bootstrap > 0 and curve.p_value is not None: + t_stat = np.full(n_grid, np.nan) + p_value = np.asarray(curve.p_value, dtype=float) + else: + t_stat = np.full(n_grid, np.nan) + p_value = np.full(n_grid, np.nan) + for i in range(n_grid): + t_i, p_i, _ = safe_inference( + curve.effects[i], curve.se[i], df=curve.df_survey + ) + t_stat[i] = t_i + p_value[i] = p_i + blocks.append((curve, target, t_stat, p_value)) + df_val = self._stored_inference_df() + return AggregationResult( + level="dose", + label=np.concatenate( + [np.asarray(c.dose_grid, dtype=object) for c, _, _, _ in blocks] + ), + target=np.array( + ["att"] * len(blocks[0][0].effects) + ["acrt"] * len(blocks[1][0].effects), + dtype=object, + ), + att=np.concatenate([c.effects for c, _, _, _ in blocks]).astype(float), + se=np.concatenate([c.se for c, _, _, _ in blocks]).astype(float), + t_stat=np.concatenate([t for _, _, t, _ in blocks]), + p_value=np.concatenate([p for _, _, _, p in blocks]), + conf_int_lower=np.concatenate([c.conf_int_lower for c, _, _, _ in blocks]).astype( + float + ), + conf_int_upper=np.concatenate([c.conf_int_upper for c, _, _, _ in blocks]).astype( + float + ), + # Grid evaluation points carry no count and no aggregation + # mass - inventing either would be a fabricated number. + n=np.full(2 * len(blocks[0][0].effects), np.nan), + df=np.full(2 * len(blocks[0][0].effects), df_val), + alpha=self.alpha, + n_kind=None, + weight=None, + estimator="ContinuousDiD", + ) + + # level == "event_study": the kit recompute. + kit = self._aggregation_kit + if kit is None: + raise ValueError( + "This ContinuousDiDResults has no aggregation kit - it is " + "attached by ContinuousDiD.fit(); a result unpickled from " + "an older release will not have one. Re-fit with " + "diff-diff >= 3.9 to enable post-fit aggregate()." + ) + bk = kit.bookkeeping + if bk["n_bootstrap"] > 0: + raise NotImplementedError( + "aggregate('event_study') on a bootstrapped ContinuousDiD " + "fit is not implemented - the fit-time event study used " + "multiplier-bootstrap inference whose per-cell draws are " + "not retained, and an analytical recompute would silently " + "differ. Until 4.0 the deprecated fit-time " + "aggregate='eventstudy' still computes the bootstrap " + "surface, or re-fit with n_bootstrap=0 for the analytical " + "post-fit route." + ) + host = _ContinuousKitAggregator(alpha=kit.alpha) + es = host._aggregate_event_study( + bk["gt_summary"], + gt_bootstrap_info=None, + unit_survey_weights=bk["unit_survey_weights"], + unit_cohorts=bk["unit_cohorts"], + anticipation=kit.anticipation, + ) + if bk["has_post_cells"]: + host._compute_event_study_inference( + es, + gt_summary=bk["gt_summary"], + gt_es_payload=bk["gt_es_payload"], + n_units=bk["n_units"], + unit_cohorts=bk["unit_cohorts"], + unit_survey_weights=bk["unit_survey_weights"], + unit_first_panel_row=bk["unit_first_panel_row"], + resolved_survey=bk["resolved_survey"], + survey_df=bk["survey_df"], + ) + # else: fit-faithful empty-post_gt quirk - the fit-time surface + # also leaves ES rows at NaN inference when no post-treatment + # cells exist. + meta = bk["survey_metadata"] + meta = dataclasses.replace(meta) if meta is not None else None + carrier = dataclasses.replace( + self, + event_study_effects=es, + survey_metadata=meta, + alpha=kit.alpha, + anticipation=kit.anticipation, + # _provenance_kwargs reads base_period off the carrier - it + # rides the kit like its siblings alpha/anticipation so + # post-fit mutation of the public field cannot reach + # recomputed provenance. + base_period=bk["base_period"], + ) + return build_event_study_surface(carrier) diff --git a/diff_diff/guides/llms-autonomous.txt b/diff_diff/guides/llms-autonomous.txt index 9100a132c..086c98e69 100644 --- a/diff_diff/guides/llms-autonomous.txt +++ b/diff_diff/guides/llms-autonomous.txt @@ -219,8 +219,8 @@ view. Every field below appears as a top-level key in that dict. `first_treat == 0`, force-zeroes nonzero `dose` on those rows with a `UserWarning`, drops units where `first_treat > 0` AND `dose == 0`, and rejects negative dose only among treated units - where `first_treat > 0`; see `continuous_did.py:276-327` and - `:348-360`). + where `first_treat > 0`; see `continuous_did.py:528-593` and + `:632-651`). In the canonical `ContinuousDiD` setup (Callaway, Goodman-Bacon, Sant'Anna 2024), the dose `D_i` is **time-invariant per unit** @@ -297,7 +297,7 @@ view. Every field below appears as a top-level key in that dict. treated-dose requirement. A continuous panel with negative non-zero doses (e.g. `dose_min == -1.5`) labeled as `first_treat > 0` would be rejected at fit time - (``continuous_did.py:287-294``); the same negative-dose units + (``continuous_did.py:586-593``); the same negative-dose units labeled as `first_treat == 0` would be coerced to dose=0 with a `UserWarning` instead. See §5.2 for the standard-workflow walkthrough. diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index b3ff97543..07baa2375 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -767,7 +767,10 @@ cdid.fit( time: str, first_treat: str, dose: str, # Column with continuous treatment dose - aggregate: str = None, # None, "dose", "eventstudy" + aggregate: Any = NOT_SUPPLIED, # DEPRECATED (M-025) - removed in 4.0; use post-fit results.aggregate(); + # legacy values None/"dose"/"eventstudy" still route until 4.0 (any + # supplied value warns; "dose" was always a no-op - the curves are + # always computed) covariates: list[str] = None, # Conditional parallel trends (X). None = unconditional ) -> ContinuousDiDResults ``` @@ -779,8 +782,16 @@ from diff_diff import ContinuousDiD est = ContinuousDiD(n_bootstrap=199, seed=42) results = est.fit(data, outcome='outcome', unit='unit', time='period', - first_treat='first_treat', dose='dose', aggregate='dose') + first_treat='first_treat', dose='dose') results.print_summary() +dose_table = results.aggregate('dose') # ATT(d) + ACRT(d) rows (view; any fit) +overall = results.aggregate('simple') # att + acrt rows (view; any fit) +# Binarized event study: analytical fits only (bootstrapped fits raise - +# re-fit with n_bootstrap=0, or the deprecated fit-time aggregate='eventstudy' +# until 4.0): +es = ContinuousDiD(seed=42).fit(data, outcome='outcome', unit='unit', + time='period', first_treat='first_treat', + dose='dose').aggregate('event_study') ``` ### HeterogeneousAdoptionDiD @@ -1785,7 +1796,7 @@ Each event study effect dict contains: `effect`, `se`, `t_stat`, `p_value`, `con | `n_clusters` | `int | None` | Number of effective clusters; `None` under survey designs and `None` under EfficientDiD's default unclustered fit | | `df_convention` | `str | None` | On the knob-carrying containers (`StackedDiDResults`, `ImputationDiDResults`, `WooldridgeDiDResults`, `LPDiDResults`, `SunAbrahamResults`): the configured df convention (3.9 / M-127); `StackedDiDResults` additionally carries `inference_df` (the overall-ATT df actually used) | -**Methods:** `summary()`, `print_summary()`, `to_dataframe()`, `to_dict()` (flat dict of headline aliases + `vcov_type` + conditional `cluster_name`/`n_clusters`/`n_bootstrap`/`inference_method`); `aggregate(type, weights=None, *, balance_e=None)` on the shipped post-fit adopters (`StackedDiDResults` views since 3.9/M-024; `EfficientDiDResults` recomputes from retained EIFs since 3.9/M-023; `ImputationDiDResults` and `TwoStageDiDResults` recompute from their PANEL-BACKED kits since 3.9/M-021/M-022 - bootstrapped fits raise on all three recomputing adopters) +**Methods:** `summary()`, `print_summary()`, `to_dataframe()`, `to_dict()` (flat dict of headline aliases + `vcov_type` + conditional `cluster_name`/`n_clusters`/`n_bootstrap`/`inference_method`); `aggregate(type, weights=None, *, balance_e=None)` on the shipped post-fit adopters (`StackedDiDResults` views since 3.9/M-024; `EfficientDiDResults` recomputes from retained EIFs since 3.9/M-023; `ImputationDiDResults` and `TwoStageDiDResults` recompute from their PANEL-BACKED kits since 3.9/M-021/M-022 - bootstrapped fits raise on all three recomputing adopters; `ContinuousDiDResults` is MIXED since 3.9/M-025 - `'simple'`/`'dose'` are views over stored fields that work on any fit incl. bootstrapped, `'event_study'` recomputes from a pruned per-cell IF kit and raises on bootstrapped fits) ### ContinuousDiDResults @@ -1810,7 +1821,7 @@ Each event study effect dict contains: `effect`, `se`, `t_stat`, `p_value`, `con | `n_obs` | `int` | Number of observations | | `n_treated_units` | `int` | Treated units | | `n_control_units` | `int` | Control units | -| `event_study_effects` | `dict[int, dict] or None` | Event study effects (if `aggregate="eventstudy"`) | +| `event_study_effects` | `dict[int, dict] or None` | Event study effects (populated only by the DEPRECATED fit-time `aggregate="eventstudy"`; the post-fit successor `results.aggregate('event_study')` returns a unified `EventStudyResults` container instead of populating this field) | **DoseResponseCurve** sub-dataclass: @@ -1823,7 +1834,7 @@ Each event study effect dict contains: `effect`, `se`, `t_stat`, `p_value`, `con | `conf_int_upper` | `np.ndarray` | Upper CI bound | | `target` | `str` | `"att"` or `"acrt"` | -**Methods:** `summary()`, `print_summary()`, `to_dataframe()` +**Methods:** `summary()`, `print_summary()`, `to_dataframe()`, `aggregate(type, ...)` (M-025 MIXED adopter: `'simple'` = 2-row att+acrt view, `'dose'` = 2N-row target-discriminated curve view - both work on any fit; `'event_study'` = kit recompute, analytical fits only) ### HeterogeneousAdoptionDiDResults diff --git a/diff_diff/guides/llms-practitioner.txt b/diff_diff/guides/llms-practitioner.txt index 5434bf73d..2198e4875 100644 --- a/diff_diff/guides/llms-practitioner.txt +++ b/diff_diff/guides/llms-practitioner.txt @@ -433,6 +433,19 @@ print(results.aggregate('event_study', balance_e=2).to_dataframe()) results = cs.fit(data, ..., aggregate='all') ``` +### For ContinuousDiD (MIXED post-fit `aggregate()`, row M-025) +```python +# Dose-response curves + overall ATT/ACRT are ALWAYS computed by fit(); +# 'dose' and 'simple' are views that work on ANY fit (bootstrapped too): +results = cdid.fit(data, ...) # no aggregate= kwarg (deprecated) +print(results.aggregate('dose').to_dataframe()) # ATT(d) + ACRT(d) rows +print(results.aggregate('simple').to_dataframe()) # att + acrt overall rows +# The binarized event study RECOMPUTES from the retained kit - analytical +# fits only (a bootstrapped fit raises; re-fit with n_bootstrap=0 or use +# the deprecated fit-time aggregate='eventstudy' until 4.0): +print(results.aggregate('event_study').to_dataframe()) +``` + ### For SunAbraham (no `aggregate` parameter) SA computes event-study and cohort effects automatically during `fit()`: ```python diff --git a/diff_diff/guides/llms.txt b/diff_diff/guides/llms.txt index 4c1bb284b..cf92172b8 100644 --- a/diff_diff/guides/llms.txt +++ b/diff_diff/guides/llms.txt @@ -21,7 +21,7 @@ diagnostic steps produces unreliable results. 4. **Choose estimator** — staggered adoption → CS/SA/BJS (NOT plain TWFE); few treated units → SDiD; factor confounding → TROP; simple 2x2 → DiD. Run `BaconDecomposition` to diagnose TWFE bias. 5. **Estimate** — `estimator.fit(data, ...)`. Always print the cluster count first and choose inference method based on the result (cluster-robust if >= 50 clusters, wild bootstrap if fewer). 6. **Sensitivity analysis** — `compute_honest_did(results)` for bounds under PT violations (MultiPeriodDiD, CS, or dCDH natively; a StackedDiD `results.aggregate('event_study')` container also admits - needs `kappa_pre >= 2` so estimated pre-periods exist), `run_all_placebo_tests()` for 2x2 falsification, specification comparisons for staggered designs. -7. **Heterogeneity** — CS: `results.aggregate('group')`/`.aggregate('event_study')` post-fit, no refit (fit-time `aggregate=`/`balance_e=` are deprecated since 3.9, removed in 4.0; `compute_honest_did` / `compute_pretrends_power` / `plot_event_study` all accept the post-fit `results.aggregate('event_study')` container directly; EXCEPTION: a BOOTSTRAPPED CS fit's `aggregate()` raises — use the fit-time aggregation there); dCDH: `results.aggregate('event_study')`/`.aggregate('simple')` post-fit views (bootstrap fits included — pure views); SA: `results.event_study_effects`/`to_dataframe(level='cohort')`; Stacked: `results.aggregate('event_study')`/`.aggregate('simple')` post-fit views (the surface is ALWAYS computed at fit since 3.9 - row M-024 - and the container admits into `compute_honest_did`/`compute_pretrends_power` with `kappa_pre >= 2`); EDiD: `results.aggregate('event_study')`/`.aggregate('group')`/`.aggregate('simple')` post-fit, RECOMPUTED from retained EIFs (3.9, row M-023; fit-time `aggregate=`/`balance_e=` deprecated; bootstrapped EDiD fits raise - use the fit-time aggregation there; EDiD containers are NOT admitted into honest/pretrends - no joint ES covariance); BJS/Gardner: `results.aggregate('event_study')`/`.aggregate('group')`/`.aggregate('simple')` post-fit on ImputationDiD and TwoStageDiD too (3.9, rows M-021/M-022; recomputed from panel-backed kits, `balance_e=` on `aggregate('event_study')`; bootstrapped fits raise - use the deprecated fit-time aggregation there; their containers are not admitted into honest/pretrends - Imputation by design, TwoStage deferred pending a normalization derivation); subgroup re-estimation. +7. **Heterogeneity** — CS: `results.aggregate('group')`/`.aggregate('event_study')` post-fit, no refit (fit-time `aggregate=`/`balance_e=` are deprecated since 3.9, removed in 4.0; `compute_honest_did` / `compute_pretrends_power` / `plot_event_study` all accept the post-fit `results.aggregate('event_study')` container directly; EXCEPTION: a BOOTSTRAPPED CS fit's `aggregate()` raises — use the fit-time aggregation there); dCDH: `results.aggregate('event_study')`/`.aggregate('simple')` post-fit views (bootstrap fits included — pure views); SA: `results.event_study_effects`/`to_dataframe(level='cohort')`; Stacked: `results.aggregate('event_study')`/`.aggregate('simple')` post-fit views (the surface is ALWAYS computed at fit since 3.9 - row M-024 - and the container admits into `compute_honest_did`/`compute_pretrends_power` with `kappa_pre >= 2`); EDiD: `results.aggregate('event_study')`/`.aggregate('group')`/`.aggregate('simple')` post-fit, RECOMPUTED from retained EIFs (3.9, row M-023; fit-time `aggregate=`/`balance_e=` deprecated; bootstrapped EDiD fits raise - use the fit-time aggregation there; EDiD containers are NOT admitted into honest/pretrends - no joint ES covariance); BJS/Gardner: `results.aggregate('event_study')`/`.aggregate('group')`/`.aggregate('simple')` post-fit on ImputationDiD and TwoStageDiD too (3.9, rows M-021/M-022; recomputed from panel-backed kits, `balance_e=` on `aggregate('event_study')`; bootstrapped fits raise - use the deprecated fit-time aggregation there; their containers are not admitted into honest/pretrends - Imputation by design, TwoStage deferred pending a normalization derivation); CGBS continuous: ContinuousDiD is a MIXED adopter (3.9, row M-025) - `results.aggregate('dose')` (ATT(d)+ACRT(d) rows) and `.aggregate('simple')` (att+acrt rows) are views over the always-computed curves and work on ANY fit incl. bootstrapped, while `.aggregate('event_study')` recomputes the binarized event study from a pruned per-cell IF kit and raises on bootstrapped fits (use the deprecated fit-time `aggregate='eventstudy'` there until 4.0; its container is not admitted into honest/pretrends - no joint ES covariance and no reference normalization); subgroup re-estimation. 8. **Robustness** — compare 2-3 estimators (CS vs SA vs BJS), MUST report with and without covariates (shows whether conditioning drives identification), present pre-trends and sensitivity bounds. After estimation, call `practitioner_next_steps(results)` for context-aware diff --git a/diff_diff/honest_did.py b/diff_diff/honest_did.py index 7cc4d5fe5..45bea41fc 100644 --- a/diff_diff/honest_did.py +++ b/diff_diff/honest_did.py @@ -609,10 +609,13 @@ def _extract_container_params( "the native results object where supported " "(MultiPeriodDiDResults, CallawaySantAnnaResults, or " "ChaisemartinDHaultfoeuilleResults); " - "EfficientDiDResults and ImputationDiDResults containers are " - "rejected BY DESIGN (their surfaces carry no joint event-study " - "covariance - per-horizon SEs only; see the REGISTRY " - "EfficientDiD and ImputationDiD Notes); TwoStageDiDResults " + "EfficientDiDResults, ImputationDiDResults and " + "ContinuousDiDResults containers are rejected BY DESIGN " + "(their surfaces carry no joint event-study covariance - " + "per-horizon SEs only - and ContinuousDiD's binarized bins " + "additionally carry no reference-period normalization at " + "all; see the REGISTRY EfficientDiD, ImputationDiD and " + "ContinuousDiD Notes); TwoStageDiDResults " "container admission is DEFERRED pending a normalization " "derivation - analytical fits carry the joint Gardner-GMM " "covariance, but the pre-period coefficients are stage-1 " diff --git a/diff_diff/practitioner.py b/diff_diff/practitioner.py index 1cbd45855..0d6897a24 100644 --- a/diff_diff/practitioner.py +++ b/diff_diff/practitioner.py @@ -1079,6 +1079,27 @@ def _handle_continuous(results: Any): priority="medium", step_name="sensitivity", ), + _step( + baker_step=7, + label="Aggregate post-fit (event study / dose / simple)", + why=( + "The fit-time aggregate= kwarg is deprecated (row M-025; " + "removed in 4.0). Aggregate as a post-fit step instead: " + "aggregate('event_study') recomputes the binarized event " + "study from the retained kit (analytical fits only - on a " + "bootstrapped fit it raises; re-fit with n_bootstrap=0 or " + "use the deprecated fit-time aggregate='eventstudy' until " + "4.0), while aggregate('dose') and aggregate('simple') are " + "views over the always-computed curves and overall " + "ATT/ACRT, available on any fit." + ), + code=( + "es = results.aggregate('event_study') # analytical fits\n" + "dose_table = results.aggregate('dose')\n" + "overall = results.aggregate('simple') # att + acrt rows" + ), + step_name="aggregation", + ), ] warnings = _check_nan_att(results) return steps, warnings @@ -1259,8 +1280,9 @@ def _handle_had(results: Any): "cdid = ContinuousDiD()\n" "cdid_results = cdid.fit(\n" " data_cdid, outcome='y', unit='unit', time='t',\n" - " first_treat='first_treat', dose='d',\n" - " aggregate='dose')" + " first_treat='first_treat', dose='d')\n" + "# Dose-response curves are always computed:\n" + "# cdid_results.dose_response_att / .aggregate('dose')" ), step_name="estimator_selection", ), @@ -1412,12 +1434,13 @@ def _handle_had_event_study(results: Any): "uses never-treated by default (or " "control_group='lowest_dose' for Remark 3.1 when " "P(D=0)=0). Two ContinuousDiD aggregation " - "surfaces are relevant and distinct: `aggregate='dose'` " - "(the default) produces the per-dose ATT(d) / ACRT(d) " - "curves on `results.dose_response_att` / " - "`results.dose_response_acrt`; `aggregate='eventstudy'` " - "produces a binarized event-study of `att_glob` on " - "`results.event_study_effects` (NOT per-dose by horizon). " + "surfaces are relevant and distinct: the per-dose ATT(d) / " + "ACRT(d) curves are ALWAYS computed by fit() (on " + "`results.dose_response_att` / `results.dose_response_acrt`, " + "or as a table via `results.aggregate('dose')`); the " + "binarized event-study of `att_glob` comes from post-fit " + "`results.aggregate('event_study')` (NOT per-dose by " + "horizon). " "Pick the aggregation that matches the estimand you " "actually want. HAD itself remains valid even with a " "small share of never-treated units (paper compatibility); " @@ -1433,9 +1456,10 @@ def _handle_had_event_study(results: Any): ), code=( "# HAD reports per-event-time WAS at the dose boundary.\n" - "# For per-dose ATT(d)/ACRT(d) curves, use ContinuousDiD with\n" - "# the DEFAULT aggregate='dose' (NOT 'eventstudy' - that gives\n" - "# binarized event-study of att_glob, not per-dose curves).\n" + "# For per-dose ATT(d)/ACRT(d) curves, use ContinuousDiD -\n" + "# fit() always computes them (post-fit aggregate('event_study')\n" + "# gives the binarized event-study of att_glob instead, which\n" + "# is NOT per-dose).\n" "from diff_diff import ContinuousDiD\n" "# ContinuousDiD requires a TIME-INVARIANT per-unit dose.\n" "# Re-prepare the panel (e.g. collapse each unit's positive\n" @@ -1443,10 +1467,10 @@ def _handle_had_event_study(results: Any): "cdid = ContinuousDiD()\n" "cdid_res = cdid.fit(\n" " data_cdid, outcome='y', unit='unit', time='t',\n" - " first_treat='first_treat', dose='d',\n" - " aggregate='dose')\n" + " first_treat='first_treat', dose='d')\n" "# Per-dose curves live here:\n" - "# cdid_res.dose_response_att / .dose_response_acrt" + "# cdid_res.dose_response_att / .dose_response_acrt\n" + "# (tabular: cdid_res.aggregate('dose'))" ), step_name="estimator_selection", ), diff --git a/diff_diff/pretrends.py b/diff_diff/pretrends.py index d24618c07..e517f77b0 100644 --- a/diff_diff/pretrends.py +++ b/diff_diff/pretrends.py @@ -1445,11 +1445,15 @@ def _extract_container_pre_period_params( f"(got source={surface.source!r}). For other estimators " "pass the native results object where supported " "(MultiPeriodDiDResults, CallawaySantAnnaResults, or " - "SunAbrahamResults); EfficientDiDResults and " - "ImputationDiDResults containers are rejected BY DESIGN " + "SunAbrahamResults); EfficientDiDResults, " + "ImputationDiDResults and ContinuousDiDResults containers " + "are rejected BY DESIGN " "(their surfaces carry no joint event-study covariance - " - "per-horizon SEs only; see the REGISTRY EfficientDiD and " - "ImputationDiD Notes); TwoStageDiDResults container " + "per-horizon SEs only - and ContinuousDiD's binarized " + "bins additionally carry no reference-period " + "normalization at all; see the REGISTRY EfficientDiD, " + "ImputationDiD and ContinuousDiD Notes); " + "TwoStageDiDResults container " "admission is DEFERRED pending a normalization derivation " "- analytical fits carry the joint Gardner-GMM covariance, " "but the pre-period coefficients are stage-1 residual " diff --git a/diff_diff/profile.py b/diff_diff/profile.py index 5b0548faf..33f5db300 100644 --- a/diff_diff/profile.py +++ b/diff_diff/profile.py @@ -98,17 +98,17 @@ class TreatmentDoseShape: ``ContinuousDiD`` out — for graded-adoption panels with dose changes use ``HeterogeneousAdoptionDiD``. 3. ``PanelProfile.is_balanced == True``. Actual fit-time gate - (``continuous_did.py:329-338``); not ``first_treat``-dependent. + (``continuous_did.py:653-661``); not ``first_treat``-dependent. 4. Absence of the ``duplicate_unit_time_rows`` alert. The precompute path silently resolves duplicate ``(unit, time)`` - cells via last-row-wins (``continuous_did.py:818-823``); + cells via last-row-wins (``continuous_did.py:1256-1261``); **not** a fit-time raise. The agent must deduplicate before fit because ``ContinuousDiD`` will otherwise overwrite silently. 5. ``treatment_dose.dose_min > 0`` (over non-zero doses). Predicts ``ContinuousDiD.fit()``'s strictly-positive-treated- dose requirement (raises ``ValueError`` on negative dose for - ``first_treat > 0`` units, ``continuous_did.py:287-294``). + ``first_treat > 0`` units, ``continuous_did.py:586-593``). Failure means some treated units have negative dose; see routing notes below. diff --git a/diff_diff/results_base.py b/diff_diff/results_base.py index a07130a93..af1c3d0cc 100644 --- a/diff_diff/results_base.py +++ b/diff_diff/results_base.py @@ -626,7 +626,10 @@ def summary(self, alpha: Optional[float] = None) -> str: "StackedDiDResults": "re-fit with diff-diff >= 3.9, which always computes the surface", "StaggeredTripleDiffResults": "refit with aggregate='event_study' (or 'all')", "EfficientDiDResults": "call results.aggregate('event_study') (on a bootstrapped fit, re-fit with n_bootstrap=0 or the deprecated fit-time aggregate=)", - "ContinuousDiDResults": "refit with aggregate='eventstudy' (or 'all')", + "ContinuousDiDResults": ( + "call results.aggregate('event_study') (on a bootstrapped fit, " + "re-fit with n_bootstrap=0 or the deprecated fit-time aggregate=)" + ), "WooldridgeDiDResults": "call results.aggregate('event_study') first", "SpilloverDiDResults": "refit with event_study=True", "ChaisemartinDHaultfoeuilleResults": "refit with L_max >= 1", diff --git a/docs/api/_autosummary/diff_diff.ContinuousDiDResults.rst b/docs/api/_autosummary/diff_diff.ContinuousDiDResults.rst index c2a10b072..8d496a7d0 100644 --- a/docs/api/_autosummary/diff_diff.ContinuousDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.ContinuousDiDResults.rst @@ -12,6 +12,7 @@ .. autosummary:: ~ContinuousDiDResults.__init__ + ~ContinuousDiDResults.aggregate ~ContinuousDiDResults.print_summary ~ContinuousDiDResults.summary ~ContinuousDiDResults.to_dataframe diff --git a/docs/api/continuous_did.rst b/docs/api/continuous_did.rst index ad442ff01..5b762132a 100644 --- a/docs/api/continuous_did.rst +++ b/docs/api/continuous_did.rst @@ -78,6 +78,7 @@ Results container for Continuous DiD estimation. ~ContinuousDiDResults.summary ~ContinuousDiDResults.print_summary ~ContinuousDiDResults.to_dataframe + ~ContinuousDiDResults.aggregate DoseResponseCurve ----------------- @@ -108,10 +109,10 @@ Basic usage:: est = ContinuousDiD(n_bootstrap=199, seed=42) results = est.fit(data, outcome='outcome', unit='unit', time='period', first_treat='first_treat', - dose='dose', aggregate='dose') + dose='dose') results.print_summary() -Accessing dose-response curves:: +Accessing dose-response curves (always computed by ``fit()``):: # ATT(d) dose-response curve as DataFrame att_df = results.dose_response_att.to_dataframe() @@ -124,13 +125,22 @@ Accessing dose-response curves:: print(f"Overall ATT: {results.overall_att:.3f} (SE: {results.overall_att_se:.3f})") print(f"Overall ACRT: {results.overall_acrt:.3f} (SE: {results.overall_acrt_se:.3f})") -Event study aggregation:: - - # Dynamic effects (binarized ATT by relative period) - results_es = est.fit(data, outcome='outcome', unit='unit', - time='period', first_treat='first_treat', - dose='dose', aggregate='eventstudy') - es_df = results_es.to_dataframe(level='event_study') + # Or as post-fit aggregation tables (views - work on any fit): + dose_table = results.aggregate('dose') # ATT(d) + ACRT(d) rows + overall = results.aggregate('simple') # att + acrt overall rows + +Event study aggregation (post-fit, row M-025):: + + # Dynamic effects (binarized ATT by relative period). The event-study + # recompute needs an ANALYTICAL fit - on a bootstrapped fit + # aggregate('event_study') raises; the deprecated fit-time + # aggregate='eventstudy' remains available for that case until 4.0. + est_es = ContinuousDiD(seed=42) + results_es = est_es.fit(data, outcome='outcome', unit='unit', + time='period', first_treat='first_treat', + dose='dose') + es = results_es.aggregate('event_study') + es_df = es.to_dataframe() Comparison with CallawaySantAnna -------------------------------- diff --git a/docs/choosing_estimator.rst b/docs/choosing_estimator.rst index cdd1a6d6f..2623bf7e6 100644 --- a/docs/choosing_estimator.rst +++ b/docs/choosing_estimator.rst @@ -373,9 +373,9 @@ Use :class:`~diff_diff.ContinuousDiD` when: est = ContinuousDiD(n_bootstrap=199, seed=42) results = est.fit(data, outcome='outcome', unit='unit', time='period', first_treat='first_treat', - dose='dose', aggregate='dose') + dose='dose') - # Overall effect and dose-response curve + # Overall effect and dose-response curve (always computed by fit) print(f"Overall ATT: {results.overall_att:.3f}") att_curve = results.dose_response_att.to_dataframe() diff --git a/docs/dev-status.md b/docs/dev-status.md index 7ec61d460..5ea31b45b 100644 --- a/docs/dev-status.md +++ b/docs/dev-status.md @@ -33,7 +33,7 @@ Target: ideally < 1000 lines per module; modules ≥3000 lines are candidates fo | `synthetic_control.py` | 2526 | Monitor | | `prep_dgp.py` | 2524 | Monitor | | `estimators.py` | 2441 | Monitor | -| `continuous_did.py` | 2431 | Monitor | +| `continuous_did.py` | 2459 | Monitor | | `sun_abraham.py` | 2314 | Monitor | | `triple_diff.py` | 2231 | Monitor | | `wooldridge.py` | 2192 | Monitor | diff --git a/docs/doc-deps.yaml b/docs/doc-deps.yaml index e38a6342d..67ce44bac 100644 --- a/docs/doc-deps.yaml +++ b/docs/doc-deps.yaml @@ -57,6 +57,7 @@ groups: - diff_diff/chaisemartin_dhaultfoeuille_results.py continuous_did: - diff_diff/continuous_did.py + - diff_diff/continuous_did_aggregation.py - diff_diff/continuous_did_bspline.py - diff_diff/continuous_did_results.py stacked_did: diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 5f5244520..d1bd3451b 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -1291,7 +1291,7 @@ there (analytical `ACRT(d_1)` SE matches the bootstrap). Otherwise the not-yet-treated controls. When `anticipation=0` (default), behavior is unchanged. - **Boundary knots**: Knots are built once from all treated doses (global, not per-cell) to ensure a common basis across (g,t) cells for aggregation. Evaluation grid is clamped to training-dose boundary knots (`range(dose)`). R's `contdid` v0.1.0 has an inconsistency where `splines2::bSpline(dvals)` uses `range(dvals)` instead of `range(dose)`, which can produce extrapolation artifacts at dose grid extremes. Our approach avoids extrapolation and is methodologically sound. -- **Note:** `bspline_derivative_design_matrix` previously swallowed `ValueError` from `scipy.interpolate.BSpline` in the per-basis derivative loop, leaving affected columns of the derivative design matrix as zero with no user-facing signal. It now aggregates the failed basis indices and emits ONE `UserWarning` naming them. Both ACRT point estimates and analytical/bootstrap inference read the same `dPsi` matrix (see `continuous_did.py:1026-1046` and the bootstrap ACRT path at `continuous_did.py:1524-1561`), so both are biased on a partial derivative-construction failure — the warning wording makes that explicit. The all-identical-knot degenerate case (single dose value) remains silently handled — derivatives there are mathematically zero. Axis-C finding #12 in the Phase 2 silent-failures audit. +- **Note:** `bspline_derivative_design_matrix` previously swallowed `ValueError` from `scipy.interpolate.BSpline` in the per-basis derivative loop, leaving affected columns of the derivative design matrix as zero with no user-facing signal. It now aggregates the failed basis indices and emits ONE `UserWarning` naming them. Both ACRT point estimates and analytical/bootstrap inference read the same `dPsi` matrix (see `_compute_dose_response_gt` - the `acrt_d = dPsi_eval @ beta_pred` point estimate at `continuous_did.py:1806` and the covariate-path IFs at `:1485-1526` - and the bootstrap ACRT path in `_bootstrap_gt_cell` at `continuous_did.py:2274-2340`), so both are biased on a partial derivative-construction failure — the warning wording makes that explicit. The all-identical-knot degenerate case (single dose value) remains silently handled — derivatives there are mathematically zero. Axis-C finding #12 in the Phase 2 silent-failures audit. ### Deviations from the paper / from R / library extensions @@ -1315,6 +1315,8 @@ labels.* 6. **Note (discrete-treatment saturated regression — library extension beyond `contdid` v0.1.0):** `treatment_type="discrete"` estimates the dose-response by a **saturated regression** (CGBS 2024 Eq. 4.1) — one indicator per distinct dose level, so `beta_j = mean_{D=d_j}(ΔY − control) = ATT(d_j)` (a per-level 2×2 DiD) — instead of the B-spline sieve. `ACRT(d_j)` is the paper's **backward difference** on the grid `{d_0 = 0, d_1, …, d_J}` (Eq. 4.1 makes `d_0 = 0` the omitted category with `ATT(0) = 0`): `ACRT(d_j) = [ATT(d_j) − ATT(d_{j-1})]/(d_j − d_{j-1})` for `j ≥ 2`, and at the lowest positive level it references the zero-dose baseline, `ACRT(d_1) = [ATT(d_1) − 0]/(d_1 − 0) = ATT(d_1)/d_1`. So a single positive dose (`J = 1`, e.g. binary `D ∈ {0,1}`) yields `ACRT(d_1) = ATT(d_1)/d_1`, and for `d_1 = 1` the documented binary identity `ACRT = ATT` holds exactly. This is a **library extension**: `contdid` v0.1.0 accepts `treatment_type` in its signature but **does not implement the discrete path** (documented "Discrete treatment not yet implemented"), so there is **no external R anchor**. It is instead an *exact* basis swap of the B-spline design/evaluation/derivative trio for an indicator/identity/finite-difference trio; every downstream quantity is linear in `beta`, so the analytical-SE / multiplier-bootstrap / covariate (reg,dr) / survey machinery is reused unchanged and reduces *analytically* to the per-level 2×2 DiD (`bread @ psi_bar = ones(J)`; the common control mean cancels in the `j ≥ 2` adjacent differences whose `L`-rows sum to 0). **reg vs dr:** the constant DR augmentation `η̄_cont` cancels in the `j ≥ 2` differences, so `ACRT(d_j)` point AND SE are identical for `reg`/`dr` there; but `ACRT(d_1) = ATT(d_1)/d_1` references the fixed baseline `ATT(0) = 0` (not shifted by `η̄_cont`), so `reg` and `dr` genuinely **differ at `ACRT(d_1)` by `η̄_cont/d_1`** (and correspondingly in `ACRT^glob` via the `d_1` mass) — the dr influence function carries the augmentation variance at `d_1` (validated: analytical `ACRT(d_1)` SE matches the multiplier bootstrap). Validation (R-free, in CI): exact hand-calc of `ATT(d_j)`/`ACRT`/`overall_att` and the analytical SE against a direct per-level 2×2 reconstruction (`~1e-12`/`~1e-10`), DGP recovery, and MC coverage for analytical + bootstrap (`tests/test_methodology_continuous_did.py::TestDiscreteSaturated`, `tests/test_continuous_did.py::TestDiscreteSaturatedAPI`). **Fail-closed policies (no-silent-failures):** (i) multi-cohort fits with **heterogeneous dose support** across cohorts raise `NotImplementedError` — an absent global level yields a dropped zero column (`att_d[level]=0`) that the plain-sum dose aggregation would bias toward zero (support-aware aggregation is deferred; single-cohort, 2-period, and shared-support multi-cohort are supported); (ii) a requested `dvals` value that is not an observed dose level raises `ValueError` (a saturated model cannot be evaluated off-support); (iii) an over-parameterized fit (`< 2` treated units per level, or `J > n_treated/2`) warns (degenerate per-level SE); (iv) with `survey_design=`, any dose level with **zero effective treated mass in a `(g,t)` cell** raises `ValueError` — a per-cell check (not just the global positive-weight check), so a level that survey/subpopulation weights zero out for one cohort while another cohort keeps it cannot silently drop to a zero-coefficient saturated column. Cross-references `docs/methodology/continuous-did.md` § 5.1. 7. **Note (lowest-dose-as-control, Remark 3.1 — library extension beyond `contdid` v0.1.0):** `control_group="lowest_dose"` implements CGBS 2024 Remark 3.1 for settings with no untreated group (`P(D=0) = 0`): the lowest-dose group `d_L` becomes the comparison and the estimand is `ATT(d) − ATT(d_L)` (SPT), with `ATT(d_L) = 0` the omitted reference. Mechanically it is a **control-group swap** — the D=0 control pool is replaced by the `d_L` group; the entire linear influence-function / bootstrap / event-study / survey machinery is control-group-generic and reused unchanged (`ee_control` already carries the reference-group variance, so **no new SE plumbing**). On the discrete saturated basis the backward-difference operator's reference shifts from `0` to `d_L` (`ACRT(d_1) = ATT(d_1)/(d_1 − d_L)`); on the continuous B-spline path the reference shifts only `μ_0` (the level), leaving `ACRT = spline'` unchanged. `contdid` v0.1.0 does **not** implement Remark 3.1, so there is **no external R anchor**; validation (R-free, in CI): an **exact `d_L → 0` equivalence** anchor (relabelling a `never_treated` panel's D=0 group as a tiny common dose `d_L = ε` reproduces the `never_treated` ATT and SE exactly, for any ε), a discrete hand-calc of `ATT(d)−ATT(d_L)`/`ACRT`/`overall_att`/`overall_acrt` and the per-level 2×2 SE (`~1e-10`), continuous mass-point DGP recovery, analytical-vs-bootstrap SE agreement, a pre-period placebo, and MC coverage (`tests/test_methodology_continuous_did.py::TestLowestDose`, `tests/test_continuous_did.py::TestLowestDoseAPI`). The continuous path requires a genuine **mass point** at the minimum dose (`>= 2` units at `d_L`, i.e. `P(D=d_L) > 0`) — the Remark 3.1 identification condition; a singleton minimum fails closed. **Fail-closed policies (no-silent-failures):** (i) never-treated units present with `lowest_dose` → `ValueError` (they would be silently dropped); (ii) singleton `d_L` (no mass point) → `ValueError`; (iii) no treated dose above `d_L` → `ValueError`; (iv) user `dvals ≤ d_L` → `ValueError` (`d_L` is the omitted reference); (v) survey/subpopulation weighting that leaves the `d_L` group with `< 2` positive-weight units → `ValueError` (a single positive-weight reference unit gives `ee_control = 0`, i.e. zero control-side variance — the effective-`>= 2` analogue of the raw mass-point guard, applied after weighting); (vi) a boundary gap `d_1 − d_L` that is a tiny fraction of the dose range warns (huge boundary ACRT/SE). **Deferred (fail-closed `NotImplementedError` + TODO):** multi-cohort `lowest_dose` (needs a within-cohort reference + support-aware cross-cohort aggregation) and `covariates=` × `lowest_dose` (conditional-PT-relative-to-`d_L` estimand). Cross-references `docs/methodology/continuous-did.md` § 5.6. +8. **Note (post-fit aggregate() - rows M-025/M-122):** `fit(aggregate=)` is deprecated in 3.9 (removed in 4.0; the no-underscore `"eventstudy"` spelling dies with it) in favor of post-fit `results.aggregate(type)` on the unified vocabulary + `'dose'` as this estimator's documented extra level; the PRE-EXISTING fit-time value validation is retained (unknown strings still raise `ValueError` after the deprecation warning - unlike the EfficientDiD/Imputation shims, which never validated). (a) **MIXED view/recompute architecture** (unique among the aggregate-postfit adopters): the dose-response curves and the overall binarized ATT (ATT^{loc} under PT; equals ATT^{glob} under SPT) plus ACRT^{glob} are ALWAYS computed by `fit()` (`aggregate="dose"` was a fit-time no-op), so `aggregate('simple')` (2 rows, targets att/acrt - the dual-estimand case the `AggregationResult.target` column exists for; `n` = the DISJOINT treated+control unit total, `n_kind='units'`) and `aggregate('dose')` (2N target-discriminated rows; labels = the dose grid twice; no count/mass per row) are pure VIEWS over stored public fields, PERMITTED on bootstrap fits (the dCDH deviation from CallawaySantAnna's fail-closed rule, applied per level) - they relay stored inference verbatim, including the FINITE `safe_inference` t-stat fit stores beside the percentile p/CI on bootstrapped overall rows and the `DoseResponseCurve.to_dataframe`-exact NaN-t derivation on dose rows; only the df column is uniformly NaN under bootstrap. (b) **`aggregate('event_study')` recomputes** the binarized event study from a pruned per-cell IF payload retained on the fit-built kit: per-(g,t) treated/control positional indices, `delta_y_treated`, `ee_control`, masses and the covariate-path `if_att_glob` (O(n_treated+n_control) per cell), unit-level arrays, and - on survey fits - the PANEL-LEVEL `ResolvedSurveyDesign` ref (the recompute performs the unit collapse itself, keeping the moved body verbatim; on replicate designs the (n_obs x R) replicate matrix rides along - the documented memory cost). The K-dimensional spline machinery (bread, `ee_treated`, `Psi_eval`, `dPsi_*`) is NOT retained; no panel data columns and no raw unit identifiers are retained. Replicate-weight designs ARE supported post-fit (IF-based `compute_replicate_if_variance` - no refit replay). Bootstrap fits carry a SCALARS-ONLY kit and the event-study route fails closed (`NotImplementedError`; the deprecated fit-time `aggregate='eventstudy'` computes the bootstrap surface until 4.0, or re-fit with `n_bootstrap=0`) - a seeded post-fit bootstrap replay is the TODO.md row. (c) **Fit-faithful quirk:** when no post-treatment (g,t) cells exist, event-study rows keep NaN inference on BOTH routes (the fit-time surface never fills them). (d) **df provenance:** the stored `dose_response_att.df_survey` channel (the value every fit-time `safe_inference` received) drives the views' df column (finite-and-positive else NaN - the replicate-undefined 0 sentinel reports NaN in the column but feeds the t/p derivation raw); the post-fit event-study container exposes the scalar `df_survey` channel only (all-NaN per-row df - the M-092 completion hole, tracked in TODO.md). (e) **Rendering:** this is the FIRST heterogeneous-`target` `AggregationResult`; `summary()`/`to_dataframe()` gained the target column / first-appearance target-block ordering amendment (normative rule in `docs/v4-design.md` section 6; uniform-target producers byte-stable). (f) **Consumer admission:** `compute_honest_did`/`compute_pretrends_power` reject ContinuousDiD containers BY DESIGN - no joint event-study covariance exists (per-bin IF SEs only) and the binarized bins carry no reference-period normalization at all (see M-093). Warning stacklevels in the moved bodies remain tuned for the fit-time frame depth, so post-fit-route warnings attribute to a library frame (the shipped EfficientDiD convention). + ### Implementation Checklist - [x] B-spline basis construction matching R's `splines2::bSpline` (global knots from all treated doses; boundary knots use training-dose range; see deviation note above) diff --git a/docs/performance-plan.md b/docs/performance-plan.md index dba55237d..c09e7e655 100644 --- a/docs/performance-plan.md +++ b/docs/performance-plan.md @@ -835,12 +835,15 @@ scenario development. None are silent-failures and none belong in this PR or in the silent-failures audit; logging here for awareness. 1. **`aggregate` / `level` parameter naming is inconsistent.** CS accepts - `aggregate="event_study"`; ContinuousDiD requires + `aggregate="event_study"`; ContinuousDiD historically required `aggregate="eventstudy"` on `fit()` **but** `level="event_study"` on `to_dataframe()`. Two different spellings within one estimator plus a third cross-estimator spelling. Surfaced when the P1 exit-propagation fix stopped silently swallowing the resulting `ValueError` in the - dose-response benchmark. Route: API-consistency cleanup, minor. + dose-response benchmark. RESOLVED by row M-025: the fit-time kwarg is + deprecated (removed in 4.0, its `"eventstudy"` spelling with it) and + the post-fit successor `results.aggregate("event_study")` uses the + unified underscored vocabulary. 2. **`generate_survey_did_data(panel=True)` `treated` column.** Row-level active-treatment indicator that is zero in pre-periods, which makes it quietly incompatible with `check_parallel_trends` (expects unit-level diff --git a/docs/performance-scenarios.md b/docs/performance-scenarios.md index 938b564a0..69f7e444a 100644 --- a/docs/performance-scenarios.md +++ b/docs/performance-scenarios.md @@ -300,18 +300,19 @@ serves a different purpose: R-parity accuracy). They complement it. ```python ContinuousDiD(degree=3, num_knots=1, n_bootstrap=199).fit( data, outcome="y", unit="unit", time="period", first_treat="first_treat", - dose="dose", aggregate="dose", + dose="dose", ) ``` -- **Operation chain.** (1) CDiD fit with `aggregate="dose"` - produces - overall ATT, overall ACRT, and the dose-response curves; (2) extract +- **Operation chain.** (1) CDiD fit - produces + overall ATT, overall ACRT, and the dose-response curves + unconditionally (the fit-time `aggregate=` kwarg is deprecated, row + M-025); (2) extract `results.to_dataframe(level="dose_response")` and - `level="group_time"` (event-study is not populated by a dose-only - fit, so it is extracted in a separate step); (3) a second CDiD fit - with `aggregate="eventstudy"` for pre-trend diagnostics (note the - spelling: `fit(aggregate="eventstudy")` with no underscore, but - `to_dataframe(level="event_study")` with underscore - see the - correctness-adjacent observations in `performance-plan.md`); + `level="group_time"`; (3) a second analytical CDiD fit followed by + post-fit `results.aggregate("event_study")` for pre-trend + diagnostics (the unified underscored spelling - the fit-time + no-underscore `"eventstudy"` value dies with the deprecated kwarg in + 4.0); (4) compare to a binarized DiD fit on the same data to quantify information loss from binarizing; (5) alternate `degree=1` (linear) and (6) `num_knots=2` refits for spline-sensitivity. The dose-curve diff --git a/docs/tutorials/14_continuous_did.ipynb b/docs/tutorials/14_continuous_did.ipynb index 73008923c..8f1a5219c 100644 --- a/docs/tutorials/14_continuous_did.ipynb +++ b/docs/tutorials/14_continuous_did.ipynb @@ -177,7 +177,7 @@ "source": [ "## Basic Estimation\n", "\n", - "The `ContinuousDiD` estimator follows an sklearn-like API. Call `fit()` with column names for outcome, unit, time, first treatment period, and dose. Setting `aggregate=\"dose\"` computes the dose-response curve and global summary parameters." + "The `ContinuousDiD` estimator follows an sklearn-like API. Call `fit()` with column names for outcome, unit, time, first treatment period, and dose. The dose-response curve and global summary parameters are always computed (the fit-time `aggregate=` keyword is deprecated; aggregation is a post-fit step via `results.aggregate(...)`)." ] }, { @@ -195,7 +195,6 @@ " time=\"period\",\n", " first_treat=\"first_treat\",\n", " dose=\"dose\",\n", - " aggregate=\"dose\",\n", ")\n", "\n", "print(results.summary())" @@ -316,7 +315,9 @@ "\n", "An event study aggregation binarizes treatment (treated vs. untreated) and shows effects by relative time period. This is useful for **pre-trends diagnostics**: pre-treatment coefficients should be near zero if parallel trends holds.\n", "\n", - "Event study works best with multiple cohorts (to have richer variation in event-time). A single cohort can still produce event-time estimates, but diagnostics are more informative with multiple cohorts. We generate a new dataset with two treatment cohorts." + "Event study works best with multiple cohorts (to have richer variation in event-time). A single cohort can still produce event-time estimates, but diagnostics are more informative with multiple cohorts. We generate a new dataset with two treatment cohorts.\n", + "\n", + "The event study is a **post-fit** aggregation: call `results.aggregate(\"event_study\")` on an analytical fit (the fit-time `aggregate=\"eventstudy\"` keyword is deprecated). On a bootstrapped fit the post-fit event study raises - re-fit with `n_bootstrap=0`, or keep the deprecated fit-time keyword until 4.0 if you need bootstrap event-study inference." ] }, { @@ -338,7 +339,7 @@ " seed=42,\n", ")\n", "\n", - "est_es = ContinuousDiD(n_bootstrap=199, seed=42)\n", + "est_es = ContinuousDiD(seed=42)\n", "results_es = est_es.fit(\n", " data_es,\n", " outcome=\"outcome\",\n", @@ -346,12 +347,12 @@ " time=\"period\",\n", " first_treat=\"first_treat\",\n", " dose=\"dose\",\n", - " aggregate=\"eventstudy\",\n", ")\n", "\n", - "# Event study table\n", - "es_df = results_es.to_dataframe(level=\"event_study\")\n", - "es_df['pre_period'] = es_df['relative_period'] < 0\n", + "# Post-fit event study (unified EventStudyResults container)\n", + "es = results_es.aggregate(\"event_study\")\n", + "es_df = es.to_dataframe()\n", + "es_df['pre_period'] = es_df['event_time'] < 0\n", "print(\"Event Study Effects (binarized ATT by relative period):\")\n", "print(es_df.to_string(index=False))" ] @@ -366,8 +367,8 @@ "if HAS_MATPLOTLIB:\n", " fig, ax = plt.subplots(figsize=(8, 5))\n", "\n", - " rel = es_df['relative_period'].values\n", - " att = es_df['att_glob'].values\n", + " rel = es_df['event_time'].values\n", + " att = es_df['att'].values\n", " ci_lo = es_df['conf_int_lower'].values\n", " ci_hi = es_df['conf_int_upper'].values\n", "\n", @@ -432,8 +433,7 @@ "print(\"-\" * 58)\n", "for cfg in configs:\n", " est_cfg = ContinuousDiD(degree=cfg[\"degree\"], num_knots=cfg[\"num_knots\"], seed=42)\n", - " res_cfg = est_cfg.fit(data_log, \"outcome\", \"unit\", \"period\", \"first_treat\", \"dose\",\n", - " aggregate=\"dose\")\n", + " res_cfg = est_cfg.fit(data_log, \"outcome\", \"unit\", \"period\", \"first_treat\", \"dose\")\n", " spline_results[cfg[\"label\"]] = res_cfg\n", " print(f\"{cfg['label']:<35} {res_cfg.overall_att:>10.4f} {res_cfg.overall_acrt:>11.4f}\")" ] @@ -503,8 +503,7 @@ "print(\"-\" * 63)\n", "for cg in [\"never_treated\", \"not_yet_treated\"]:\n", " est_cg = ContinuousDiD(control_group=cg, seed=42)\n", - " res_cg = est_cg.fit(data_mc, \"outcome\", \"unit\", \"period\", \"first_treat\", \"dose\",\n", - " aggregate=\"dose\")\n", + " res_cg = est_cg.fit(data_mc, \"outcome\", \"unit\", \"period\", \"first_treat\", \"dose\")\n", " print(f\"{cg:<20} {res_cg.overall_att:>10.4f} {res_cg.overall_att_se:>10.4f} \"\n", " f\"{res_cg.overall_acrt:>11.4f} {res_cg.overall_acrt_se:>10.4f}\")" ] @@ -535,13 +534,11 @@ "source": [ "# Analytical SEs (default)\n", "est_ana = ContinuousDiD(seed=42)\n", - "res_ana = est_ana.fit(data, \"outcome\", \"unit\", \"period\", \"first_treat\", \"dose\",\n", - " aggregate=\"dose\")\n", + "res_ana = est_ana.fit(data, \"outcome\", \"unit\", \"period\", \"first_treat\", \"dose\")\n", "\n", "# Bootstrap SEs\n", "est_boot = ContinuousDiD(n_bootstrap=199, seed=42)\n", - "res_boot = est_boot.fit(data, \"outcome\", \"unit\", \"period\", \"first_treat\", \"dose\",\n", - " aggregate=\"dose\")\n", + "res_boot = est_boot.fit(data, \"outcome\", \"unit\", \"period\", \"first_treat\", \"dose\")\n", "\n", "print(f\"{'Method':<15} {'ATT_glob':>10} {'SE':>10} {'95% CI':>25}\")\n", "print(\"-\" * 62)\n", @@ -649,7 +646,7 @@ "\n", "- **`level=\"dose_response\"`** (default): The dose-response curve with ATT(d) and ACRT(d)\n", "- **`level=\"group_time\"`**: Underlying group-time cell estimates\n", - "- **`level=\"event_study\"`**: Event study effects (only available when fitted with `aggregate=\"eventstudy\"`)" + "- **`level=\"event_study\"`**: Event study effects (legacy surface - populated only by the deprecated fit-time `aggregate=\"eventstudy\"`; prefer the post-fit `results.aggregate(\"event_study\")` container and its `to_dataframe()`)" ] }, { @@ -700,7 +697,7 @@ "| `n_bootstrap` | 0 | Bootstrap iterations (0 = analytical SEs only) |\n", "| `bootstrap_weights` | `\"rademacher\"` | `\"rademacher\"`, `\"mammen\"`, or `\"webb\"` |\n", "| `seed` | `None` | Random seed for reproducibility |\n", - "| `aggregate` | `None` | `\"dose\"` for dose-response, `\"eventstudy\"` for event study |\n", + "| `aggregate` (fit) | deprecated | Fit-time aggregation is deprecated (row M-025); use post-fit `results.aggregate(\"event_study\"/\"dose\"/\"simple\")` - the dose curves are always computed |\n", "\n", "### Reference\n", "\n", diff --git a/docs/v4-deprecations.yaml b/docs/v4-deprecations.yaml index e6921203c..fec662f08 100644 --- a/docs/v4-deprecations.yaml +++ b/docs/v4-deprecations.yaml @@ -291,11 +291,12 @@ rows: introduced_in: "3.9" deprecated_in: "3.9" removed_in: "4.0" - status: planned - phase: 2 + status: shimmed + phase: 5 warning: FutureWarning - code_refs: [diff_diff/continuous_did.py] - notes: "The no-underscore 'eventstudy' spelling dies with this param; aggregate() accepts only the unified vocabulary (+ 'dose' as this estimator's extra level)." + test_ref: tests/test_aggregate_contract.py + code_refs: [diff_diff/continuous_did.py, diff_diff/continuous_did_aggregation.py, diff_diff/continuous_did_results.py, diff_diff/aggregation.py, diff_diff/results_base.py, diff_diff/honest_did.py, diff_diff/pretrends.py, diff_diff/practitioner.py, diff_diff/guides/llms-practitioner.txt] + notes: "Shimmed in 3.9: fit(aggregate=) warns via the shared NOT_SUPPLIED sentinel (a plain fit() never warns; supplying ANY value incl. None warns once, then the legacy routing runs unchanged). Unlike the EfficientDiD/Imputation shims, fit-time VALUE VALIDATION pre-existed and STAYS: unknown strings still raise ValueError after the warning ((None, 'dose', 'eventstudy') only). The no-underscore 'eventstudy' spelling dies with the param in 4.0; aggregate() accepts only the unified vocabulary + 'dose' as this estimator's documented extra level. MIXED VIEW/RECOMPUTE architecture (unique in the register): 'simple' and 'dose' are pure VIEWS over stored public fields - the dose curves and the overall binarized ATT (ATT^{loc} under PT; equals ATT^{glob} under SPT) plus ACRT^{glob} are ALWAYS computed by fit (aggregate='dose' was a fit-time no-op) - so both levels are PERMITTED on bootstrap fits (the dCDH deviation from CallawaySantAnna's fail-closed rule, applied PER-LEVEL), relaying stored inference verbatim: the overall rows carry the FINITE safe_inference t fit stores beside percentile p/CI, the dose rows reproduce DoseResponseCurve.to_dataframe (NaN t under bootstrap), and only the df column is uniformly NaN there. 'event_study' is a PRUNED-IF-PAYLOAD kit recompute (see the AggregationKit docstring variant): per-(g,t) O(n_treated+n_control) IF-ingredient arrays + unit-level arrays + the PANEL-LEVEL resolved survey design (on replicate designs the (n_obs x R) replicate matrix rides along - a unit-level collapse was reviewed and declined for verbatim-move safety); no panel data columns and no raw unit identifiers are retained; bootstrap fits get a SCALARS-ONLY kit and the ES route fails closed (NotImplementedError naming the fit-time route / n_bootstrap=0 refit; ContinuousDiDResults has no bootstrap_results field, so the config gate is the honest one). Replicate-weight designs ARE supported post-fit (IF-based compute_replicate_if_variance - no refit replay, the contrast with M-021/M-022). Fit-faithful quirk: empty-post_gt fits leave ES rows at NaN inference on both routes (has_post_cells flag). CONTAINER SHAPES: simple = 2 rows (targets att/acrt - the dual-estimand case the AggregationResult target column exists for; n = disjoint treated+control units total, n_kind='units'); dose = 2N target-discriminated rows (labels = the dose grid twice; n NaN / n_kind None / weight None - grid evaluation points carry no count or mass); df from the stored dose_response_att.df_survey channel (finite-and->0 else NaN; the raw stored value incl. the replicate 0-sentinel feeds the to_dataframe-exact t/p derivation). Ships the FIRST heterogeneous-target container, with the AggregationResult summary()/to_dataframe rendering amendment (target column + neutral estimate heading when targets mixed; FIRST-APPEARANCE target-block ordering, labels ascending within block under the _sortable guard; uniform-target producers byte-stable) - normative rule in v4-design section 6. Supported subset simple|event_study|dose - group/calendar fail closed via the mixin; balance_e applies to no level (empty _AGGREGATE_BALANCE_E_TYPES); weights= rejected. Admission: honest/pretrends containers rejected BY DESIGN (no joint ES covariance; bins not reference-normalized - see M-093). Bootstrap-ES post-fit replay is the TODO.md row." - id: M-026 kind: param group: aggregate-postfit @@ -988,7 +989,7 @@ rows: phase: 2 test_ref: tests/test_event_study_surface.py code_refs: [diff_diff/results_base.py, diff_diff/__init__.py] - notes: "Phase 2 unified event-study representation (spec section 5): EventStudyResults container + builders for the 14 producers (CallawaySantAnna, SunAbraham, ImputationDiD, TwoStageDiD, StackedDiD, SpilloverDiD, ContinuousDiD, EfficientDiD, WooldridgeDiD, StaggeredTripleDifference, MultiPeriodDiD, LPDiD, ChaisemartinDHaultfoeuille, HeterogeneousAdoptionDiD). Canonical quintet columns, explicit is_reference marking (successor to the retiring sentinels [M-093]), vcov+vcov_index ordering, cband columns, event_time_convention metadata. Public exposure rides aggregate(type='event_study') in Phase 2 PR (b); merged TWFE returns it in Phase 3 [M-010]. introduced_in gates the 3.9 cut, mirroring [M-091]. Born done in this introducing diff (builder is package-internal; the class is exported). Amended pre-cut (same test_ref, introduced_in 3.9 not yet released): df became PER-ROW (one entry per event time, the df each stored p/CI actually used; joins the pinned to_dataframe schema) and StackedDiD/TwoStageDiD persist their internal full ES VCVs (event_study_vcov/_index/_df container fields; mode-gated for TwoStageDiD bootstrap/replicate). Completed pre-cut by the remaining producer channels: SunAbraham (per-event dict) and de Chaisemartin-D'Haultfoeuille (scalar) event_study_df, plus LPDiD pooled_df for the headline pre/post windows - every producer whose inference records a df now exposes it. Amended pre-cut a second time (2(b) PR-1, with M-026): three optional PROVENANCE fields appended last - base_period, anticipation, and df_survey (the fit's resolved SCALAR inference df beside the per-row channel: survey_metadata.df_survey with replicate-undefined mapping to the 0.0 fail-closed sentinel, else df_inference, else None - the per-row df column cannot encode that sentinel because __post_init__ NaNs it wherever p is non-finite) - threaded by the builders (the _empty_surface early return included) so the container consumers (HonestDiD's universal-base check, PreTrendsPower's anticipation cutoff, the honest df extraction) read fit-faithful values instead of dropping them. Amended pre-cut a third time (same PR): a fourth provenance field reference_event_times (also new on CallawaySantAnnaResults, computed at fit under base_period=universal) - the DISTINCT per-cohort positional-base event times, the common-reference signal is_reference cannot carry on gapped grids where a cohort's base overlaps another cohort's estimated horizon; HonestDiD and PreTrendsPower fail closed on more than one entry, on BOTH input routes (REGISTRY HonestDiD common-reference-guard Note). SCOPE QUALIFIER (2(b) PR-3a, with M-023): the 'every producer whose inference records a df now exposes it' sentence is satisfied on EfficientDiD through the SCALAR df_survey channel only - its newly public aggregate('event_study') container has no per-row df source (no event_study_df/df_inference field; the per-row column is all-NaN, contract-permitted for a producer that records none); threading the retained kit scalar into the per-row channel is the TODO.md M-092-completion row. The qualifier EXTENDS to ImputationDiD (2(b) PR-3b, with M-021): its newly public container likewise exposes only the scalar df_survey channel (no per-row df source; all-NaN per-row column, identical to its fit-time surface) - the same TODO.md completion row names it." + notes: "Phase 2 unified event-study representation (spec section 5): EventStudyResults container + builders for the 14 producers (CallawaySantAnna, SunAbraham, ImputationDiD, TwoStageDiD, StackedDiD, SpilloverDiD, ContinuousDiD, EfficientDiD, WooldridgeDiD, StaggeredTripleDifference, MultiPeriodDiD, LPDiD, ChaisemartinDHaultfoeuille, HeterogeneousAdoptionDiD). Canonical quintet columns, explicit is_reference marking (successor to the retiring sentinels [M-093]), vcov+vcov_index ordering, cband columns, event_time_convention metadata. Public exposure rides aggregate(type='event_study') in Phase 2 PR (b); merged TWFE returns it in Phase 3 [M-010]. introduced_in gates the 3.9 cut, mirroring [M-091]. Born done in this introducing diff (builder is package-internal; the class is exported). Amended pre-cut (same test_ref, introduced_in 3.9 not yet released): df became PER-ROW (one entry per event time, the df each stored p/CI actually used; joins the pinned to_dataframe schema) and StackedDiD/TwoStageDiD persist their internal full ES VCVs (event_study_vcov/_index/_df container fields; mode-gated for TwoStageDiD bootstrap/replicate). Completed pre-cut by the remaining producer channels: SunAbraham (per-event dict) and de Chaisemartin-D'Haultfoeuille (scalar) event_study_df, plus LPDiD pooled_df for the headline pre/post windows - every producer whose inference records a df now exposes it. Amended pre-cut a second time (2(b) PR-1, with M-026): three optional PROVENANCE fields appended last - base_period, anticipation, and df_survey (the fit's resolved SCALAR inference df beside the per-row channel: survey_metadata.df_survey with replicate-undefined mapping to the 0.0 fail-closed sentinel, else df_inference, else None - the per-row df column cannot encode that sentinel because __post_init__ NaNs it wherever p is non-finite) - threaded by the builders (the _empty_surface early return included) so the container consumers (HonestDiD's universal-base check, PreTrendsPower's anticipation cutoff, the honest df extraction) read fit-faithful values instead of dropping them. Amended pre-cut a third time (same PR): a fourth provenance field reference_event_times (also new on CallawaySantAnnaResults, computed at fit under base_period=universal) - the DISTINCT per-cohort positional-base event times, the common-reference signal is_reference cannot carry on gapped grids where a cohort's base overlaps another cohort's estimated horizon; HonestDiD and PreTrendsPower fail closed on more than one entry, on BOTH input routes (REGISTRY HonestDiD common-reference-guard Note). SCOPE QUALIFIER (2(b) PR-3a, with M-023): the 'every producer whose inference records a df now exposes it' sentence is satisfied on EfficientDiD through the SCALAR df_survey channel only - its newly public aggregate('event_study') container has no per-row df source (no event_study_df/df_inference field; the per-row column is all-NaN, contract-permitted for a producer that records none); threading the retained kit scalar into the per-row channel is the TODO.md M-092-completion row. The qualifier EXTENDS to ImputationDiD (2(b) PR-3b, with M-021): its newly public container likewise exposes only the scalar df_survey channel (no per-row df source; all-NaN per-row column, identical to its fit-time surface) - the same TODO.md completion row names it. The qualifier EXTENDS to ContinuousDiD (2(b) PR-3c, with M-025): its newly public container exposes only the scalar df_survey channel via the carrier's survey_metadata (no event_study_df/df_inference field, so _from_relative_dict publishes an all-NaN per-row df column on survey fits whose ES rows received a finite _survey_df - identical on the fit-time and post-fit routes); the same TODO.md completion row names it." - id: M-093 kind: behavior group: results-contract @@ -999,7 +1000,7 @@ rows: status: planned phase: 5 code_refs: [diff_diff/staggered_results.py, diff_diff/sun_abraham.py, diff_diff/imputation_results.py, diff_diff/two_stage_results.py, diff_diff/stacked_did_results.py, diff_diff/efficient_did_results.py, diff_diff/continuous_did_results.py, diff_diff/wooldridge_results.py, diff_diff/chaisemartin_dhaultfoeuille_results.py, diff_diff/lpdid_results.py, diff_diff/staggered_triple_diff_results.py, diff_diff/results.py, diff_diff/had.py, diff_diff/visualization/_event_study.py, diff_diff/honest_did.py, diff_diff/pretrends.py] - notes: "4.0 sentinel retirement + schema enforcement (spec section 5): the n_groups==0 / n_obs==0 reference-row sentinels retire; every estimator's to_dataframe(level='event_study') emits the [M-092] column schema; the plotter / HonestDiD / PreTrendsPower consume the unified surface. Thirteen ES-carrying source modules enumerated file-by-file (results.py covers MultiPeriodDiD + SpilloverDiD) plus the plotter. behavior-at-done requires test_ref. PARTIAL PRE-CUT DELIVERY (2(b) PR-1, with M-026): the consumer half shipped for CS-SOURCED containers - compute_honest_did, compute_pretrends_power (both with SOURCE-SCOPED admission: source == CallawaySantAnnaResults only, rejecting dCDH l1 containers BY DESIGN since their placebo semantics need honest_did's native branch, and every other producer pending its own aggregate() migration) and plot_event_study / plot_honest_event_study (no source scoping - plotting is label-faithful). Admission widening is each later shim PR's methodology decision, not automatic. SECOND PRE-CUT AMENDMENT (2(b) PR-2, with M-024): admission widened to STACKED-SOURCED containers in compute_honest_did + compute_pretrends_power (source in {CallawaySantAnnaResults, StackedDiDResults}; kappa_pre >= 2 required for estimated pre-periods; honest additionally needs a non-singular full retained event-study covariance (pre+post sub-block); withheld-inference rows admitted with a source-scoped warning; the seven producer-derived guard messages now derive the producer from surface.source). Source scoping continues for all other producers. THIRD PRE-CUT AMENDMENT (2(b) PR-3a, with M-023): admission NOT widened to EfficientDiD-sourced containers BY DESIGN (the dCDH precedent) - the PRIMARY ground is the absent joint event-study covariance (container vcov=None, all-NaN per-row df); reference semantics are regime-dependent (PT-All: no reference row, universal first-period baseline; PT-Post: a materialized mechanical zero anchor at e=-1-anticipation when estimated, marked via the membership-gated reference_period property - never fabricated when absent); both terminal TypeErrors now state the EfficientDiD rejection explicitly; source scoping continues to reject. FOURTH PRE-CUT AMENDMENT (2(b) PR-3b, with M-021/M-022): admission NOT widened by the Imputation/TwoStage migrations - ImputationDiD rejected BY DESIGN (no joint ES covariance; per-horizon Theorem-3 conservative SEs only); TwoStageDiD DEFERRED despite carrying the joint Gardner-GMM covariance on analytical fits (M-092), because its pre-period coefficients are stage-1 residual means (the reference horizon is dropped from the no-intercept Stage-2 design; the zero anchor row is appended mechanically), not contrasts against the advertised reference, while HonestDiD's Delta^RM/Delta^SD arithmetic hard-codes the delta_0=0 normalization into its boundary/bridge constraints - admission awaits a normalization derivation (DEFERRED.md paper-gated row); both terminal TypeErrors state the per-producer grounds. Source scoping continues to reject all other producers. This row stays planned: its transition remains the 4.0 sentinel retirement + full-producer schema enforcement." + notes: "4.0 sentinel retirement + schema enforcement (spec section 5): the n_groups==0 / n_obs==0 reference-row sentinels retire; every estimator's to_dataframe(level='event_study') emits the [M-092] column schema; the plotter / HonestDiD / PreTrendsPower consume the unified surface. Thirteen ES-carrying source modules enumerated file-by-file (results.py covers MultiPeriodDiD + SpilloverDiD) plus the plotter. behavior-at-done requires test_ref. PARTIAL PRE-CUT DELIVERY (2(b) PR-1, with M-026): the consumer half shipped for CS-SOURCED containers - compute_honest_did, compute_pretrends_power (both with SOURCE-SCOPED admission: source == CallawaySantAnnaResults only, rejecting dCDH l1 containers BY DESIGN since their placebo semantics need honest_did's native branch, and every other producer pending its own aggregate() migration) and plot_event_study / plot_honest_event_study (no source scoping - plotting is label-faithful). Admission widening is each later shim PR's methodology decision, not automatic. SECOND PRE-CUT AMENDMENT (2(b) PR-2, with M-024): admission widened to STACKED-SOURCED containers in compute_honest_did + compute_pretrends_power (source in {CallawaySantAnnaResults, StackedDiDResults}; kappa_pre >= 2 required for estimated pre-periods; honest additionally needs a non-singular full retained event-study covariance (pre+post sub-block); withheld-inference rows admitted with a source-scoped warning; the seven producer-derived guard messages now derive the producer from surface.source). Source scoping continues for all other producers. THIRD PRE-CUT AMENDMENT (2(b) PR-3a, with M-023): admission NOT widened to EfficientDiD-sourced containers BY DESIGN (the dCDH precedent) - the PRIMARY ground is the absent joint event-study covariance (container vcov=None, all-NaN per-row df); reference semantics are regime-dependent (PT-All: no reference row, universal first-period baseline; PT-Post: a materialized mechanical zero anchor at e=-1-anticipation when estimated, marked via the membership-gated reference_period property - never fabricated when absent); both terminal TypeErrors now state the EfficientDiD rejection explicitly; source scoping continues to reject. FOURTH PRE-CUT AMENDMENT (2(b) PR-3b, with M-021/M-022): admission NOT widened by the Imputation/TwoStage migrations - ImputationDiD rejected BY DESIGN (no joint ES covariance; per-horizon Theorem-3 conservative SEs only); TwoStageDiD DEFERRED despite carrying the joint Gardner-GMM covariance on analytical fits (M-092), because its pre-period coefficients are stage-1 residual means (the reference horizon is dropped from the no-intercept Stage-2 design; the zero anchor row is appended mechanically), not contrasts against the advertised reference, while HonestDiD's Delta^RM/Delta^SD arithmetic hard-codes the delta_0=0 normalization into its boundary/bridge constraints - admission awaits a normalization derivation (DEFERRED.md paper-gated row); both terminal TypeErrors state the per-producer grounds. FIFTH PRE-CUT AMENDMENT (2(b) PR-3c, with M-025): admission NOT widened by the ContinuousDiD migration - rejected BY DESIGN on two independent grounds: no joint event-study covariance exists (per-bin IF SEs only), and the binarized bins carry NO reference-period normalization at all (no reference row exists in the surface; each bin is a raw reweighted binarized ATT level), so HonestDiD's delta_0=0 arithmetic has no anchor even in principle without new methodology; both terminal TypeErrors name the ContinuousDiD rejection alongside EfficientDiD/Imputation. Source scoping continues to reject all other producers. This row stays planned: its transition remains the 4.0 sentinel retirement + full-producer schema enforcement." # ---- Behavior policies (schema-tracked, spec-governed; no reality probe) - - id: M-080 diff --git a/docs/v4-design.md b/docs/v4-design.md index a67c874bc..e48b523a6 100644 --- a/docs/v4-design.md +++ b/docs/v4-design.md @@ -443,6 +443,19 @@ EXISTING post-fit `aggregate(type=)` - the emfx-style prior art for this section's pattern - plus `summary(aggregation=)` and `to_dataframe(aggregation=)` [M-044] [M-086] [M-087]. +**Heterogeneous-`target` rendering** (added with [M-025], the first +producer of a mixed-target `AggregationResult` - ContinuousDiD's att/acrt +dual estimand): when a container carries more than one distinct `target`, +`summary()` renders a `target` column and a neutral `estimate` heading +(the uniform-target `ATT` heading would mislabel the other estimand's +rows), and `to_dataframe()` orders rows by FIRST-APPEARANCE target blocks +(producer order - NOT lexicographic, which would invert att/acrt) with +labels ascending within each block under the same `_sortable` guard as +the uniform path (mixed-type labels keep producer order, never raise). +Uniform-target containers render byte-identically to the pre-[M-025] +output; the machine surfaces (`to_dataframe` columns, `to_dict`) carried +`target` from birth and are unchanged. + **Semantics.** `aggregate()` re-aggregates WITHOUT refitting, from influence functions retained on the results object. **Correction (Phase 2b PR 1):** this section previously said "CallawaySantAnna already stores them". It does not - diff --git a/tests/test_aggregate_contract.py b/tests/test_aggregate_contract.py index c48937173..c76cbc60d 100644 --- a/tests/test_aggregate_contract.py +++ b/tests/test_aggregate_contract.py @@ -3,7 +3,9 @@ The ``test_ref`` for every aggregate-postfit ledger row this file pins: M-020/M-023/M-021/M-022 (the CS / EfficientDiD / Imputation / TwoStage ``fit(aggregate=)`` shims), M-024/M-026 (the Stacked / dCDH shims + view -relays), M-117/M-120/M-118/M-119 (``balance_e`` moves onto ``aggregate()``) +relays), M-025 (the ContinuousDiD shim + MIXED view/recompute aggregate() - +'simple'/'dose' views, 'event_study' pruned-IF-payload kit), +M-117/M-120/M-118/M-119 (``balance_e`` moves onto ``aggregate()``) and M-122 (``AggregationResult``). The headline gate is NUMERICAL INERTNESS: for every supported type, @@ -1116,6 +1118,7 @@ def test_mixin_hooks_are_not_dataclass_fields(self): assert "EfficientDiDResults" in checked assert "ImputationDiDResults" in checked assert "TwoStageDiDResults" in checked + assert "ContinuousDiDResults" in checked # --------------------------------------------------------------------------- # @@ -2819,3 +2822,565 @@ def test_empty_balance_window(self, twostage_fitted): es = twostage_fitted.aggregate("event_study", balance_e=100) df_ = es.to_dataframe() assert (df_["is_reference"] | ~np.isfinite(df_["att"])).all() + + +# --------------------------------------------------------------------------- # +# ContinuousDiD (row M-025): fit(aggregate=) shim + the MIXED view/recompute +# aggregate() - 'simple'/'dose' are pure VIEWS over stored fields (permitted +# on bootstrap fits), 'event_study' is a pruned-IF-payload kit recompute +# (bootstrap fails closed; replicate designs supported, no refit replay) +# --------------------------------------------------------------------------- # + +CONT_KW = dict( + outcome="outcome", unit="unit", time="period", first_treat="first_treat", dose="dose" +) + + +def _cont_panel(seed=19, n_units=110, n_periods=6, cohort_periods=None): + from diff_diff import generate_continuous_did_data + + return generate_continuous_did_data( + n_units=n_units, + n_periods=n_periods, + cohort_periods=cohort_periods or [3, 4], + seed=seed, + ) + + +def _cont_covariate_panel(seed=23): + d = _cont_panel(seed=seed) + rng = np.random.default_rng(seed) + xmap = {u: rng.normal() for u in d["unit"].unique()} + d["x1"] = d["unit"].map(xmap) + return d + + +def _cont_survey_panel(seed=29, replicate=False, degenerate=None, zero_dose_treated=False): + d = _cont_panel(seed=seed) + rng = np.random.default_rng(seed) + wmap = {u: rng.uniform(0.5, 2.0) for u in d["unit"].unique()} + d["w"] = d["unit"].map(wmap) + d["strata"] = d["unit"] % 4 + d["psu"] = d["unit"] + rep_cols = [] + if replicate: + n_rep = 8 + for r in range(n_rep): + col = f"rw{r}" + rep_cols.append(col) + if degenerate == "dropped" and r >= n_rep - 2: + d[col] = 0.0 + elif degenerate == "undefined" and r >= 1: + d[col] = 0.0 + else: + jitter = {u: rng.uniform(0.1, 2.0) for u in d["unit"].unique()} + d[col] = d["unit"].map(jitter) * d["w"] + if zero_dose_treated: + # Treated units with dose == 0: fires the fit-time drop (with + # UserWarning) AND the post-drop survey re-resolution - the exact + # resolved_survey object the kit snapshots. + treated_units = sorted(d.loc[d["first_treat"] > 0, "unit"].unique())[:4] + d.loc[d["unit"].isin(treated_units), "dose"] = 0.0 + return d, rep_cols + + +def _cont_survey_design(rep_cols=None, tsl=True): + from diff_diff import SurveyDesign + + if rep_cols: + return SurveyDesign(weights="w", replicate_weights=rep_cols, replicate_method="JK1") + if tsl: + return SurveyDesign(weights="w", strata="strata", psu="psu") + return SurveyDesign(weights="w") + + +def _cont_discrete_panel(seed=31): + d = _cont_panel(seed=seed, cohort_periods=[3]) + rng = np.random.default_rng(seed) + treated = d["first_treat"] > 0 + level_map = {u: float(rng.choice([1.0, 2.0, 3.0])) for u in d.loc[treated, "unit"].unique()} + d.loc[treated, "dose"] = d.loc[treated, "unit"].map(level_map) + return d + + +def _cont_lowest_dose_panel(seed=37, n_units=60): + """Single cohort, no never-treated units, mass point at d_L.""" + rng = np.random.default_rng(seed) + rows = [] + for u in range(n_units): + d_u = 0.5 if u < 10 else float(rng.uniform(1.0, 3.0)) + g = 3 + for t in range(1, 6): + treated = t >= g + rows.append( + { + "unit": u, + "period": t, + "first_treat": g, + "dose": d_u, + "outcome": u / 25 + + 0.3 * t + + (1.0 + 0.8 * d_u if treated else 0.0) + + rng.normal(0, 0.3), + } + ) + return pd.DataFrame(rows) + + +def _cont_empty_post_panel(seed=41, n_units=50): + """All treatment starts AFTER the last observed period: pre-period + (g,t) cells exist but post_gt is empty (the fit-time warning path; + overall/dose fields are all-NaN and ES rows keep NaN inference).""" + rng = np.random.default_rng(seed) + rows = [] + for u in range(n_units): + g = 0 if u % 3 == 0 else 5 + d_u = 0.0 if g == 0 else float(rng.uniform(0.5, 2.0)) + for t in range(1, 5): + rows.append( + { + "unit": u, + "period": t, + "first_treat": g, + "dose": d_u, + "outcome": u / 20 + 0.4 * t + rng.normal(0, 0.3), + } + ) + return pd.DataFrame(rows) + + +def _fit_cont(data, *, est_kw=None, **fit_kw): + from diff_diff import ContinuousDiD + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return ContinuousDiD(**(est_kw or {})).fit(data, **CONT_KW, **fit_kw) + + +@pytest.fixture(scope="module") +def cont_panel(): + return _cont_panel() + + +@pytest.fixture(scope="module") +def cont_fitted(cont_panel): + """Plain fit - the pruned-payload kit powers aggregate('event_study').""" + return _fit_cont(cont_panel) + + +@pytest.fixture(scope="module") +def cont_fit_time(cont_panel): + """Deprecated fit-time aggregate='eventstudy' - the inertness reference.""" + return _fit_cont(cont_panel, aggregate="eventstudy") + + +@pytest.fixture(scope="module") +def cont_bootstrap(cont_panel): + return _fit_cont(cont_panel, est_kw=dict(n_bootstrap=30, seed=3)) + + +class TestContinuousShim: + def test_plain_fit_does_not_warn(self, cont_panel): + from diff_diff import ContinuousDiD + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + ContinuousDiD().fit(cont_panel, **CONT_KW) + assert [w for w in caught if issubclass(w.category, FutureWarning)] == [] + + def test_aggregate_kwarg_warns_even_at_none(self, cont_panel): + from diff_diff import ContinuousDiD + + with pytest.warns(FutureWarning, match=r"ContinuousDiD\.fit\(aggregate=\)"): + ContinuousDiD().fit(cont_panel, **CONT_KW, aggregate=None) + + def test_dose_value_warns_and_is_inert(self, cont_panel, cont_fitted): + """aggregate='dose' was always a fit-time no-op - the deprecated + call warns (the message says so) and produces the identical fit.""" + from diff_diff import ContinuousDiD + + with pytest.warns(FutureWarning, match="already[\\s\\n ]*redundant"): + res = ContinuousDiD().fit(cont_panel, **CONT_KW, aggregate="dose") + assert res.overall_att == cont_fitted.overall_att + assert res.overall_acrt_se == cont_fitted.overall_acrt_se + np.testing.assert_array_equal( + res.dose_response_att.effects, cont_fitted.dose_response_att.effects + ) + assert res.event_study_effects is None + + def test_warn_and_still_work(self, cont_fit_time): + """The deprecated 'eventstudy' value still computes the fit-time + surface (legacy routing unchanged until 4.0).""" + assert cont_fit_time.event_study_effects is not None + assert any(np.isfinite(v["se"]) for v in cont_fit_time.event_study_effects.values()) + + def test_invalid_value_warns_then_raises(self, cont_panel): + """The PRE-EXISTING fit-time value validation survives the shim + (a delta vs the EfficientDiD/Imputation no-validation shims): + unknown strings still raise ValueError - after the warning.""" + from diff_diff import ContinuousDiD + + with pytest.warns(FutureWarning): + with pytest.raises(ValueError, match="Invalid aggregate"): + ContinuousDiD().fit(cont_panel, **CONT_KW, aggregate="event_study") + + +class TestContinuousAggregate: + # ---------------- event_study: kit-recompute inertness ---------------- + + def test_event_study_inert(self, cont_fitted, cont_fit_time): + es = cont_fitted.aggregate("event_study") + _assert_es_container_matches_fit_time(es, cont_fit_time, "cont/es") + # n_kind is a container ATTRIBUTE (not an EVENT_STUDY_SCHEMA + # column): ContinuousDiD's ES dict carries no count key. + assert es.n_kind is None + + @pytest.mark.parametrize( + "arm", + [ + "multi_cohort", + "anticipation", + "covariates", + "survey_tsl", + "survey_zero_dose_drop", + "replicate_healthy", + "replicate_dropped", + "replicate_undefined", + "discrete", + "lowest_dose", + "not_yet_treated", + "nondefault_config", + ], + ) + def test_event_study_inert_across_designs(self, arm): + est_kw, fit_kw = {}, {} + if arm == "multi_cohort": + data = _cont_panel(seed=43, cohort_periods=[3, 4, 5]) + elif arm == "anticipation": + data = _cont_panel(seed=47) + est_kw = dict(anticipation=1) + elif arm == "covariates": + data = _cont_covariate_panel() + fit_kw = dict(covariates=["x1"]) + elif arm == "survey_tsl": + data, _ = _cont_survey_panel() + fit_kw = dict(survey_design=_cont_survey_design()) + elif arm == "survey_zero_dose_drop": + data, _ = _cont_survey_panel(zero_dose_treated=True) + fit_kw = dict(survey_design=_cont_survey_design()) + elif arm == "replicate_healthy": + data, rep = _cont_survey_panel(replicate=True) + fit_kw = dict(survey_design=_cont_survey_design(rep_cols=rep)) + elif arm == "replicate_dropped": + data, rep = _cont_survey_panel(replicate=True, degenerate="dropped") + fit_kw = dict(survey_design=_cont_survey_design(rep_cols=rep)) + elif arm == "replicate_undefined": + data, rep = _cont_survey_panel(replicate=True, degenerate="undefined") + fit_kw = dict(survey_design=_cont_survey_design(rep_cols=rep)) + elif arm == "discrete": + data = _cont_discrete_panel() + est_kw = dict(treatment_type="discrete") + elif arm == "lowest_dose": + data = _cont_lowest_dose_panel() + est_kw = dict(control_group="lowest_dose") + elif arm == "not_yet_treated": + data = _cont_panel(seed=53, cohort_periods=[3, 5]) + est_kw = dict(control_group="not_yet_treated") + else: # nondefault_config + data = _cont_panel(seed=59) + est_kw = dict(alpha=0.10, base_period="universal") + plain = _fit_cont(data, est_kw=est_kw, **fit_kw) + ref = _fit_cont(data, est_kw=est_kw, aggregate="eventstudy", **fit_kw) + es = plain.aggregate("event_study") + _assert_es_container_matches_fit_time(es, ref, f"cont/{arm}") + + def test_nondefault_provenance_discriminates_defaults(self): + """alpha=0.10 / base_period='universal' must SURFACE - a kit that + hard-codes the defaults would pass mutation-isolation alone.""" + data = _cont_panel(seed=59) + res = _fit_cont(data, est_kw=dict(alpha=0.10, base_period="universal")) + es = res.aggregate("event_study") + assert es.alpha == 0.10 + assert es.base_period == "universal" + s = res.aggregate("simple") + assert s.alpha == 0.10 + + def test_empty_post_gt_all_levels(self): + """Fit-faithful empty-post_gt quirk: ES rows keep NaN inference on + BOTH routes; the views relay the stored all-NaN fields.""" + data = _cont_empty_post_panel() + plain = _fit_cont(data) + ref = _fit_cont(data, aggregate="eventstudy") + es = plain.aggregate("event_study") + _assert_es_container_matches_fit_time(es, ref, "cont/empty_post") + frame = es.to_dataframe() + assert not np.isfinite(frame["se"]).any() + s = plain.aggregate("simple") + assert np.isnan(s.att).all() and np.isnan(s.se).all() + d = plain.aggregate("dose") + assert np.isnan(d.att).all() and np.isnan(d.se).all() + + # ---------------- dose view: per-curve to_dataframe parity ---------------- + + def test_dose_view_per_curve_parity(self, cont_fitted): + agg = cont_fitted.aggregate("dose") + n = len(cont_fitted.dose_grid) + assert list(agg.target) == ["att"] * n + ["acrt"] * n + for block, curve in ( + (slice(0, n), cont_fitted.dose_response_att), + (slice(n, 2 * n), cont_fitted.dose_response_acrt), + ): + frame = curve.to_dataframe() + np.testing.assert_array_equal(agg.att[block], frame["effect"].to_numpy()) + np.testing.assert_array_equal(agg.se[block], frame["se"].to_numpy()) + np.testing.assert_array_equal( + agg.conf_int_lower[block], frame["conf_int_lower"].to_numpy() + ) + np.testing.assert_allclose( + agg.t_stat[block], frame["t_stat"].to_numpy(), rtol=0, atol=0, equal_nan=True + ) + np.testing.assert_allclose( + agg.p_value[block], frame["p_value"].to_numpy(), rtol=0, atol=0, equal_nan=True + ) + assert agg.n_kind is None and agg.weight is None + assert np.isnan(agg.n).all() + + def test_dose_view_survey_df(self): + """Survey fits: finite df threads into the derived t/p AND the df + column; the per-curve to_dataframe oracle still holds.""" + data, _ = _cont_survey_panel() + res = _fit_cont(data, survey_design=_cont_survey_design()) + assert res.dose_response_att.df_survey is not None + agg = res.aggregate("dose") + n = len(res.dose_grid) + frame = res.dose_response_att.to_dataframe() + np.testing.assert_allclose( + agg.p_value[:n], frame["p_value"].to_numpy(), rtol=0, atol=0, equal_nan=True + ) + finite = np.isfinite(agg.p_value) + assert np.isfinite(agg.df[finite]).all() + + def test_dose_view_replicate_undefined_sentinel(self): + """The replicate-undefined 0-sentinel: NaN t/p via safe_inference + (the raw stored value feeds the derivation) and NaN df column.""" + data, rep = _cont_survey_panel(replicate=True, degenerate="undefined") + res = _fit_cont(data, survey_design=_cont_survey_design(rep_cols=rep)) + assert res.dose_response_att.df_survey == 0 + agg = res.aggregate("dose") + assert np.isnan(agg.t_stat).all() and np.isnan(agg.p_value).all() + assert np.isnan(agg.df).all() + + def test_dose_view_bootstrap_relays_stored_p(self, cont_bootstrap): + agg = cont_bootstrap.aggregate("dose") + n = len(cont_bootstrap.dose_grid) + assert np.isnan(agg.t_stat).all() + np.testing.assert_array_equal( + agg.p_value[:n], np.asarray(cont_bootstrap.dose_response_att.p_value, dtype=float) + ) + assert np.isnan(agg.df).all() + + # ---------------- simple view ---------------- + + def test_simple_view_bit_exact(self, cont_fitted): + s = cont_fitted.aggregate("simple") + assert list(s.label) == ["overall", "overall"] + assert list(s.target) == ["att", "acrt"] + assert s.att[0] == cont_fitted.overall_att + assert s.att[1] == cont_fitted.overall_acrt + assert s.se[0] == cont_fitted.overall_att_se + assert s.se[1] == cont_fitted.overall_acrt_se + assert s.t_stat[0] == cont_fitted.overall_att_t_stat + assert s.p_value[1] == cont_fitted.overall_acrt_p_value + assert s.conf_int_lower[0] == cont_fitted.overall_att_conf_int[0] + assert s.conf_int_upper[1] == cont_fitted.overall_acrt_conf_int[1] + # Disjoint treated/control unit sets -> the CS total convention. + assert s.n_kind == "units" + expected_n = float(cont_fitted.n_treated_units + cont_fitted.n_control_units) + assert (s.n == expected_n).all() + np.testing.assert_array_equal(s.weight, [1.0, 1.0]) + + def test_simple_view_bootstrap_finite_t_relays(self, cont_bootstrap): + """fit() stores a FINITE safe_inference t beside the percentile + p/CI on bootstrap fits - the relay carries it through verbatim + (bit-exact relay, NOT NaN); only the df column is NaN.""" + s = cont_bootstrap.aggregate("simple") + assert s.t_stat[0] == cont_bootstrap.overall_att_t_stat + assert np.isfinite(s.t_stat[0]) + assert s.p_value[0] == cont_bootstrap.overall_att_p_value + assert np.isnan(s.df).all() + + def test_simple_view_survey_df(self): + data, _ = _cont_survey_panel() + res = _fit_cont(data, survey_design=_cont_survey_design()) + s = res.aggregate("simple") + assert np.isfinite(s.df[0]) + assert s.df[0] == float(res.dose_response_att.df_survey) + + # ---------------- heterogeneous-target rendering ---------------- + + def test_summary_renders_target_column(self, cont_fitted): + s_text = cont_fitted.aggregate("simple").summary() + assert "target" in s_text and "estimate" in s_text + assert "acrt" in s_text + # The uniform-target 'ATT' heading must NOT appear as a column head. + header_line = [ln for ln in s_text.splitlines() if "estimate" in ln][0] + assert "ATT" not in header_line + + def test_uniform_target_summary_byte_stable(self, fitted): + """A uniform-target producer's summary() renders EXACTLY as before + the heterogeneous-target amendment (no target column, ATT head).""" + s_text = fitted.aggregate("simple").summary() + assert "target" not in s_text + assert "ATT" in s_text + + def test_dose_ordering_att_block_first_with_unsorted_dvals(self): + """FIRST-APPEARANCE target blocks (att then acrt - NOT lexicographic, + which would invert) with labels ascending WITHIN each block; the + custom dvals grid is unsorted so within-block sorting is actually + exercised (the default grid is ascending by construction).""" + data = _cont_panel(seed=61) + dvals = np.array([2.0, 1.0, 1.5]) + res = _fit_cont(data, est_kw=dict(dvals=dvals)) + frame = res.aggregate("dose").to_dataframe() + n = len(dvals) + assert list(frame["target"]) == ["att"] * n + ["acrt"] * n + att_labels = frame["label"][:n].astype(float).to_numpy() + acrt_labels = frame["label"][n:].astype(float).to_numpy() + np.testing.assert_array_equal(att_labels, np.sort(dvals)) + np.testing.assert_array_equal(acrt_labels, np.sort(dvals)) + + def test_mixed_type_labels_preserve_producer_order(self): + """The _sortable fallback survives the heterogeneous-target branch: + mixed-type labels keep producer order per block, never raise.""" + agg = AggregationResult( + level="dose", + label=np.array(["b", 2, "a", 1], dtype=object), + target=np.array(["att", "att", "acrt", "acrt"], dtype=object), + att=np.zeros(4), + se=np.ones(4), + t_stat=np.zeros(4), + p_value=np.ones(4), + conf_int_lower=np.zeros(4), + conf_int_upper=np.zeros(4), + n=np.full(4, np.nan), + df=np.full(4, np.nan), + ) + frame = agg.to_dataframe() + assert list(frame["label"]) == ["b", 2, "a", 1] + assert list(frame["target"]) == ["att", "att", "acrt", "acrt"] + assert "target" in agg.summary() + + # ---------------- bootstrap gating + kit shape ---------------- + + def test_bootstrap_event_study_fails_closed(self, cont_bootstrap): + with pytest.raises(NotImplementedError, match="bootstrap"): + cont_bootstrap.aggregate("event_study") + + def test_bootstrap_views_still_work(self, cont_bootstrap): + assert cont_bootstrap.aggregate("simple") is not None + assert cont_bootstrap.aggregate("dose") is not None + + def test_bootstrap_kit_is_scalars_only(self, cont_bootstrap): + """Dead-retention guarantee: the ES payload can never be consumed + on a bootstrap fit, so it is not retained.""" + bk = cont_bootstrap._aggregation_kit.bookkeeping + assert bk["gt_es_payload"] == {} + assert bk["gt_summary"] == {} + assert bk["n_units"] is None + assert bk["unit_cohorts"] is None + assert bk["resolved_survey"] is None + assert bk["n_bootstrap"] == 30 + + # ---------------- isolation (ES route only - views are views) ---------------- + + def test_es_isolation_from_public_field_mutation(self, cont_panel): + res = _fit_cont(cont_panel) + baseline = res.aggregate("event_study").to_dataframe() + res.event_study_effects = {99: {"effect": 1.0}} + res.group_time_effects = {} + res.groups = [] + res.alpha = 0.5 + res.anticipation = 7 + res.base_period = "mutated" + if res.survey_metadata is not None: + res.survey_metadata.df_survey = -1 + again = res.aggregate("event_study").to_dataframe() + pd.testing.assert_frame_equal(baseline, again) + + def test_bootstrap_gate_reads_kit_not_field(self, cont_panel): + """The sharpest isolation arm: mutating res.n_bootstrap = 0 must + NOT bypass the fail-closed gate (it reads the kit).""" + res = _fit_cont(cont_panel, est_kw=dict(n_bootstrap=30, seed=3)) + res.n_bootstrap = 0 + with pytest.raises(NotImplementedError, match="bootstrap"): + res.aggregate("event_study") + + def test_repeated_calls_idempotent(self, cont_fitted): + a = cont_fitted.aggregate("event_study").to_dataframe() + b = cont_fitted.aggregate("event_study").to_dataframe() + pd.testing.assert_frame_equal(a, b) + + def test_aggregate_does_not_mutate_survey_metadata(self): + data, rep = _cont_survey_panel(replicate=True) + res = _fit_cont(data, survey_design=_cont_survey_design(rep_cols=rep)) + before = copy.deepcopy(res.survey_metadata.__dict__) + res.aggregate("event_study") + res.aggregate("simple") + res.aggregate("dose") + assert res.survey_metadata.__dict__ == before + + # ---------------- no-kit legacy + vocabulary ---------------- + + def test_no_kit_es_raises_views_work(self, cont_panel): + res = _fit_cont(cont_panel) + object.__setattr__(res, "_aggregation_kit", None) + with pytest.raises(ValueError, match="aggregation kit"): + res.aggregate("event_study") + # Views need no kit - deliberately still work on legacy pickles. + assert res.aggregate("simple") is not None + assert res.aggregate("dose") is not None + + @pytest.mark.parametrize("bad", ["group", "calendar", "all", "nonsense"]) + def test_unsupported_types_fail_closed(self, cont_fitted, bad): + with pytest.raises(ValueError, match="Unsupported aggregation type"): + cont_fitted.aggregate(bad) + + @pytest.mark.parametrize("level", ["simple", "event_study", "dose"]) + def test_balance_e_rejected_empty_vocabulary(self, cont_fitted, level): + with pytest.raises(ValueError, match="no aggregation type on this estimator"): + cont_fitted.aggregate(level, balance_e=1) + + def test_weights_rejected(self, cont_fitted): + with pytest.raises(ValueError, match="does not accept a weights selector"): + cont_fitted.aggregate("simple", weights="cell") + + # ---------------- pickle + retention ---------------- + + def test_pickle_round_trip(self, cont_fitted): + clone = pickle.loads(pickle.dumps(cont_fitted)) + a = cont_fitted.aggregate("event_study").to_dataframe() + b = clone.aggregate("event_study").to_dataframe() + pd.testing.assert_frame_equal(a, b) + pd.testing.assert_frame_equal( + cont_fitted.aggregate("dose").to_dataframe(), + clone.aggregate("dose").to_dataframe(), + ) + + def test_no_raw_unit_identifiers_are_retained(self, cont_fitted): + """The kit stores positional indices and first_treat cohort values + only - never the raw unit identifier column.""" + bk = cont_fitted._aggregation_kit.bookkeeping + + def _walk(obj): + if isinstance(obj, dict): + for v in obj.values(): + yield from _walk(v) + elif isinstance(obj, (list, tuple)): + for v in obj: + yield from _walk(v) + else: + yield obj + + for leaf in _walk({k: v for k, v in bk.items() if k != "survey_metadata"}): + assert not isinstance(leaf, pd.DataFrame), "kit retains a DataFrame" + assert not isinstance(leaf, pd.Series), "kit retains a Series" diff --git a/tests/test_continuous_did.py b/tests/test_continuous_did.py index 0afb1dc56..509f09c89 100644 --- a/tests/test_continuous_did.py +++ b/tests/test_continuous_did.py @@ -18,6 +18,13 @@ from diff_diff.continuous_did_results import ContinuousDiDResults from diff_diff.prep_dgp import generate_continuous_did_data +# M-025: sites deliberately kept on the deprecated fit-time aggregate= +# (bootstrap event studies - the aggregated bootstrap surface is +# fit-time-only until replay ships - and legacy-surface pins) run under +# this module-scoped suppression; the shim behavior itself is pinned in +# tests/test_aggregate_contract.py::TestContinuousShim. +pytestmark = pytest.mark.filterwarnings(r"ignore:ContinuousDiD\.fit\(aggregate=\):FutureWarning") + # ============================================================================= # B-Spline Basis Tests # ============================================================================= @@ -375,10 +382,19 @@ def test_invalid_aggregate_raises(self): } ) est = ContinuousDiD() - with pytest.raises(ValueError, match="Invalid aggregate"): - est.fit( - data, "outcome", "unit", "period", "first_treat", "dose", aggregate="event_study" - ) + # M-025: the shim warns on ANY supplied value BEFORE the surviving + # value validation raises - both must fire, in that order. + with pytest.warns(FutureWarning, match=r"ContinuousDiD\.fit\(aggregate=\)"): + with pytest.raises(ValueError, match="Invalid aggregate"): + est.fit( + data, + "outcome", + "unit", + "period", + "first_treat", + "dose", + aggregate="event_study", + ) def test_no_never_treated_error(self): data = pd.DataFrame( @@ -529,7 +545,6 @@ def test_multi_period_aggregation(self): "period", "first_treat", "dose", - aggregate="dose", ) # With linear DGP (ATT(d) = 1 + 2d) and degree=1, should recover well # ACRT should be close to 2.0 @@ -550,7 +565,6 @@ def test_single_cohort_aggregation(self): "period", "first_treat", "dose", - aggregate="dose", ) assert len(results.groups) == 1 assert np.isfinite(results.overall_att) @@ -615,22 +629,13 @@ def test_event_study_not_yet_treated(self): noise_sd=0.5, ) est = ContinuousDiD(control_group="not_yet_treated", n_bootstrap=0) - results = est.fit( - data, - "outcome", - "unit", - "period", - "first_treat", - "dose", - aggregate="eventstudy", - ) - assert results.event_study_effects is not None - rel_periods = sorted(results.event_study_effects.keys()) + results = est.fit(data, "outcome", "unit", "period", "first_treat", "dose") + es = results.aggregate("event_study").to_dataframe() + rel_periods = sorted(es["event_time"]) assert min(rel_periods) < 0 # Pre-treatment assert max(rel_periods) >= 0 # Post-treatment - for e, info in results.event_study_effects.items(): - assert np.isfinite(info["effect"]), f"effect is NaN for e={e}" - assert np.isfinite(info["se"]), f"SE is NaN for e={e}" + assert np.isfinite(es["att"]).all(), "effect is NaN for some bin" + assert np.isfinite(es["se"]).all(), "SE is NaN for some bin" def test_event_study_universal_base_period(self): """Event study with base_period='universal' and analytic SE.""" @@ -642,22 +647,13 @@ def test_event_study_universal_base_period(self): noise_sd=0.5, ) est = ContinuousDiD(base_period="universal", n_bootstrap=0) - results = est.fit( - data, - "outcome", - "unit", - "period", - "first_treat", - "dose", - aggregate="eventstudy", - ) - assert results.event_study_effects is not None - rel_periods = sorted(results.event_study_effects.keys()) + results = est.fit(data, "outcome", "unit", "period", "first_treat", "dose") + es = results.aggregate("event_study").to_dataframe() + rel_periods = sorted(es["event_time"]) assert min(rel_periods) < 0 # Pre-treatment assert max(rel_periods) >= 0 # Post-treatment - for e, info in results.event_study_effects.items(): - assert np.isfinite(info["effect"]), f"effect is NaN for e={e}" - assert np.isfinite(info["se"]), f"SE is NaN for e={e}" + assert np.isfinite(es["att"]).all(), "effect is NaN for some bin" + assert np.isfinite(es["se"]).all(), "SE is NaN for some bin" def test_event_study_not_yet_treated_bootstrap(self, ci_params): """Event study with not_yet_treated control group and bootstrap SE.""" @@ -1260,23 +1256,15 @@ def test_anticipation_event_study(self): seed=42, ) est = ContinuousDiD(anticipation=1, n_bootstrap=0) - results = est.fit( - data, - "outcome", - "unit", - "period", - "first_treat", - "dose", - aggregate="eventstudy", - ) - assert results.event_study_effects is not None + results = est.fit(data, "outcome", "unit", "period", "first_treat", "dose") + es = results.aggregate("event_study").to_dataframe() # With anticipation=1 and g=3, post-treatment starts at t=2 (g - anticipation). # Relative times e = t - g, so t=2 → e=-1 (the anticipation period). - rel_times = sorted(results.event_study_effects.keys()) + rel_times = sorted(es["event_time"]) assert ( -1 in rel_times ), f"Anticipation period e=-1 missing from event study; got {rel_times}" - assert np.isfinite(results.event_study_effects[-1]["effect"]) + assert np.isfinite(es.loc[es["event_time"] == -1, "att"]).all() def test_anticipation_event_study_excludes_contaminated_periods(self): """With anticipation=2, event study should not contain e < -2.""" @@ -1315,17 +1303,9 @@ def test_anticipation_event_study_excludes_contaminated_periods(self): data = pd.DataFrame(rows) est = ContinuousDiD(anticipation=2, n_bootstrap=0) - results = est.fit( - data, - "outcome", - "unit", - "period", - "first_treat", - "dose", - aggregate="eventstudy", - ) - assert results.event_study_effects is not None - for e in results.event_study_effects.keys(): + results = est.fit(data, "outcome", "unit", "period", "first_treat", "dose") + es = results.aggregate("event_study").to_dataframe() + for e in es["event_time"]: assert e >= -2, f"Found relative period e={e} with anticipation=2; " f"expected e >= -2" def test_anticipation_not_yet_treated_excludes_anticipation_window(self): @@ -1512,23 +1492,17 @@ def test_event_study_analytical_se_finite(self): noise_sd=0.5, ) est = ContinuousDiD(n_bootstrap=0) - results = est.fit( - data, - "outcome", - "unit", - "period", - "first_treat", - "dose", - aggregate="eventstudy", - ) - assert results.event_study_effects is not None - for e, info in results.event_study_effects.items(): - assert np.isfinite(info["se"]), f"SE is NaN for e={e}" - assert info["se"] > 0, f"SE is non-positive for e={e}" - assert np.isfinite(info["t_stat"]), f"t_stat is NaN for e={e}" - assert np.isfinite(info["p_value"]), f"p_value is NaN for e={e}" - assert 0 <= info["p_value"] <= 1, f"p_value out of range for e={e}" - lo, hi = info["conf_int"] + results = est.fit(data, "outcome", "unit", "period", "first_treat", "dose") + es = results.aggregate("event_study").to_dataframe() + assert len(es) > 0 + for _, row in es.iterrows(): + e = row["event_time"] + assert np.isfinite(row["se"]), f"SE is NaN for e={e}" + assert row["se"] > 0, f"SE is non-positive for e={e}" + assert np.isfinite(row["t_stat"]), f"t_stat is NaN for e={e}" + assert np.isfinite(row["p_value"]), f"p_value is NaN for e={e}" + assert 0 <= row["p_value"] <= 1, f"p_value out of range for e={e}" + lo, hi = row["conf_int_lower"], row["conf_int_upper"] assert np.isfinite(lo) and np.isfinite(hi), f"conf_int contains NaN for e={e}" @@ -1619,7 +1593,6 @@ def test_ipw_with_covariates_raises(self): "period", "first_treat", "dose", - aggregate="dose", covariates=["x1"], ) @@ -1628,7 +1601,7 @@ def test_ipw_without_covariates_ok(self): # break the unconditional path. data = _cov_data() est = ContinuousDiD(estimation_method="ipw") - res = est.fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose") + res = est.fit(data, "outcome", "unit", "period", "first_treat", "dose") assert np.isfinite(res.overall_att) def test_survey_with_covariates_raises(self): @@ -1645,7 +1618,6 @@ def test_survey_with_covariates_raises(self): "period", "first_treat", "dose", - aggregate="dose", covariates=["x1"], survey_design=SurveyDesign(weights="w"), ) @@ -1661,7 +1633,6 @@ def test_missing_covariate_column_raises(self): "period", "first_treat", "dose", - aggregate="dose", covariates=["not_a_col"], ) @@ -1679,7 +1650,6 @@ def test_missing_covariate_values_raise(self): "period", "first_treat", "dose", - aggregate="dose", covariates=["x1"], ) @@ -1711,7 +1681,6 @@ def test_covariate_metadata_on_results(self): "period", "first_treat", "dose", - aggregate="dose", covariates=["x1"], ) assert res.covariates == ["x1"] @@ -1744,7 +1713,6 @@ def test_covariate_eventstudy_and_bootstrap(self): "period", "first_treat", "dose", - aggregate="dose", covariates=["x1"], ) boot = ContinuousDiD(estimation_method=method, n_bootstrap=199, seed=3).fit( @@ -1754,7 +1722,6 @@ def test_covariate_eventstudy_and_bootstrap(self): "period", "first_treat", "dose", - aggregate="dose", covariates=["x1"], ) assert np.isfinite(boot.overall_att_se) @@ -1765,24 +1732,22 @@ def test_clone_refit_idempotent(self): data = _cov_data() with pytest.warns(FutureWarning, match=r"\(covariates=\) is deprecated"): est = ContinuousDiD(covariates=["x1"], estimation_method="dr", seed=1) - r1 = est.fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose") + r1 = est.fit(data, "outcome", "unit", "period", "first_treat", "dose") # raw-keep: the clone re-warns because the config still carries the # deprecated ctor covariates (M-084, documented). with pytest.warns(FutureWarning, match=r"\(covariates=\) is deprecated"): clone = ContinuousDiD(**est.get_params()) - r2 = clone.fit(data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose") + r2 = clone.fit(data, "outcome", "unit", "period", "first_treat", "dose") assert abs(float(r1.overall_att) - float(r2.overall_att)) < 1e-12 assert abs(float(r1.overall_att_se) - float(r2.overall_att_se)) < 1e-12 def test_no_covariate_path_unchanged(self): """Passing covariates=None runs the unchanged unconditional path.""" data = _cov_data() - base = ContinuousDiD(seed=1).fit( - data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose" - ) + base = ContinuousDiD(seed=1).fit(data, "outcome", "unit", "period", "first_treat", "dose") # estimation_method has no effect without covariates alt = ContinuousDiD(estimation_method="reg", seed=1).fit( - data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose" + data, "outcome", "unit", "period", "first_treat", "dose" ) assert float(base.overall_att) == float(alt.overall_att) assert float(base.overall_att_se) == float(alt.overall_att_se) @@ -1830,7 +1795,6 @@ def add(ft, d): time="period", first_treat="first_treat", dose="dose", - aggregate="dose", ) @@ -2008,10 +1972,10 @@ def test_continuous_default_matches_explicit(self): """Default treatment_type is 'continuous'; explicit value gives identical output.""" data = generate_continuous_did_data(n_units=120, n_periods=3, seed=13) r_default = ContinuousDiD(n_bootstrap=0).fit( - data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose" + data, "outcome", "unit", "period", "first_treat", "dose" ) r_explicit = ContinuousDiD(treatment_type="continuous", n_bootstrap=0).fit( - data, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose" + data, "outcome", "unit", "period", "first_treat", "dose" ) np.testing.assert_allclose( r_default.dose_response_att.effects, r_explicit.dose_response_att.effects @@ -2098,7 +2062,6 @@ def test_continuous_no_mass_point_raises(self): time="period", first_treat="first_treat", dose="dose", - aggregate="dose", ) def test_continuous_one_dose_above_warns(self): @@ -2121,7 +2084,6 @@ def test_continuous_one_dose_above_warns(self): time="period", first_treat="first_treat", dose="dose", - aggregate="dose", ) assert np.all(np.isfinite(res.dose_response_att.se)) assert np.allclose(res.dose_response_acrt.effects, 0.0) @@ -2168,12 +2130,12 @@ def test_event_study_lowest_dose(self): {1.0: 0.5, 2.0: 1.5}, n_control=0, cohorts=(2,), n_periods=4, noise=0.0, seed=6 ) res = ContinuousDiD(control_group="lowest_dose", treatment_type="discrete").fit( - df, "outcome", "unit", "period", "first_treat", "dose", aggregate="eventstudy" + df, "outcome", "unit", "period", "first_treat", "dose" ) - assert res.event_study_effects is not None + es = res.aggregate("event_study").to_dataframe() # Pre-period event bins (e < 0) difference out to ~0 (both groups untreated). - pre = {e: v for e, v in res.event_study_effects.items() if e < 0} - assert pre and all(abs(v["effect"]) < 1e-9 for v in pre.values()) + pre = es[es["event_time"] < 0] + assert len(pre) > 0 and (pre["att"].abs() < 1e-9).all() def test_survey_zeroed_dL_group_raises(self): """A survey design zeroing the entire d_L reference group fails closed.""" diff --git a/tests/test_event_study_consumers.py b/tests/test_event_study_consumers.py index 1c16e68cb..7897a1706 100644 --- a/tests/test_event_study_consumers.py +++ b/tests/test_event_study_consumers.py @@ -2015,3 +2015,35 @@ def test_real_twostage_container_rejected_deferred(self): compute_honest_did(surface, M=1.0) with pytest.raises(TypeError, match="DEFERRED pending a normalization"): compute_pretrends_power(surface, M=1.0) + + +class TestContinuousContainerRejection: + """A REAL ContinuousDiD post-fit container is rejected BY DESIGN. + + Two independent grounds, both named in the clause: no joint event-study + covariance (per-bin IF SEs only), and the binarized bins carry NO + reference-period normalization at all (no reference row exists). The + match strings pin the NEW clause text, not the class name (the + got source={...!r} interpolation would match "ContinuousDiDResults" + even without the message edit). + """ + + def test_real_continuous_container_rejected_by_design(self): + from diff_diff import ContinuousDiD, generate_continuous_did_data + + d = generate_continuous_did_data(n_units=90, n_periods=6, cohort_periods=[3, 4], seed=17) + res = ContinuousDiD().fit( + d, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + dose="dose", + ) + surface = res.aggregate("event_study") + assert surface.source == "ContinuousDiDResults" + assert surface.vcov is None # no joint ES covariance exists + with pytest.raises(TypeError, match="no reference-period[\\s\\n ]*normalization"): + compute_honest_did(surface, M=1.0) + with pytest.raises(TypeError, match="no reference-period[\\s\\n ]*normalization"): + compute_pretrends_power(surface, M=1.0) diff --git a/tests/test_guides.py b/tests/test_guides.py index 4fb13bbd1..1522cc35e 100644 --- a/tests/test_guides.py +++ b/tests/test_guides.py @@ -67,7 +67,7 @@ def test_content_stability_autonomous_fingerprints(): "must be removed before fitting." ) # ContinuousDiD also requires strictly positive treated doses - # (`continuous_did.py:287-294` raises on negative dose support). + # (`continuous_did.py:586-593` raises on negative dose support). # The autonomous guide must list `dose_min > 0` so an agent reading # `treatment_dose.dose_min == -1.5` knows to route the panel away # from ContinuousDiD before paying for the failed fit. @@ -897,7 +897,10 @@ def _section(self, header): nxt = text.index("\n### ", start + 1) return text[start:nxt] - def _assert_shim_lines(self, section, fit_call, agg_row, bal_row): + def _assert_shim_lines(self, section, fit_call, agg_row, bal_row=None): + """``bal_row=None`` skips the balance_e-line lookup - required for + estimators without a balance_e param (ContinuousDiD, M-025), where + the unconditional ``next(...)`` would raise StopIteration.""" fit_start = section.index(fit_call) fit_block = section[fit_start : section.index("\n)", fit_start)] agg_line = next( @@ -906,11 +909,12 @@ def _assert_shim_lines(self, section, fit_call, agg_row, bal_row): assert "NOT_SUPPLIED" in agg_line assert f"DEPRECATED ({agg_row})" in agg_line assert "results.aggregate()" in agg_line - bal_line = next( - line for line in fit_block.splitlines() if line.strip().startswith("balance_e") - ) - assert "NOT_SUPPLIED" in bal_line - assert f"DEPRECATED ({bal_row})" in bal_line + if bal_row is not None: + bal_line = next( + line for line in fit_block.splitlines() if line.strip().startswith("balance_e") + ) + assert "NOT_SUPPLIED" in bal_line + assert f"DEPRECATED ({bal_row})" in bal_line def test_llms_full_imputation_fit_aggregate_line_documents_shim(self): # The M-021/M-118 twin of the EfficientDiD pin above - no other @@ -919,3 +923,7 @@ def test_llms_full_imputation_fit_aggregate_line_documents_shim(self): def test_llms_full_two_stage_fit_aggregate_line_documents_shim(self): self._assert_shim_lines(self._section("### TwoStageDiD"), ".fit(", "M-022", "M-119") + + def test_llms_full_continuous_fit_aggregate_line_documents_shim(self): + # M-025: no balance_e twin exists - bal_row stays None. + self._assert_shim_lines(self._section("### ContinuousDiD"), "cdid.fit(", "M-025") diff --git a/tests/test_methodology_continuous_did.py b/tests/test_methodology_continuous_did.py index 9263f06dc..c0f373b3c 100644 --- a/tests/test_methodology_continuous_did.py +++ b/tests/test_methodology_continuous_did.py @@ -433,14 +433,10 @@ def _compare_with_r( staggered=staggered, ) - # Map R aggregation names to Python aggregate parameter - py_aggregate = None - if aggregation == "dose": - py_aggregate = "dose" - elif aggregation == "eventstudy": - py_aggregate = "eventstudy" - - # Python estimation using R's dvals for exact grid match + # Python estimation using R's dvals for exact grid match. M-025: + # the dose curves are always computed by fit(); the event-study + # route (R aggregation "eventstudy") is exercised POST-FIT so the + # benchmark keeps validating that code path. dvals = np.array(r_out["dvals"]) est = ContinuousDiD( degree=degree, @@ -455,8 +451,9 @@ def _compare_with_r( "period", "first_treat", "dose", - aggregate=py_aggregate, ) + if aggregation == "eventstudy": + results.aggregate("event_study") # Compare overall ATT r_overall_att = r_out["overall_att"] @@ -607,7 +604,6 @@ def test_benchmark_4_staggered_dose(self): "period", "first_treat", "dose", - aggregate="dose", ) # Overall ATT @@ -705,7 +701,6 @@ def test_benchmark_5_not_yet_treated(self): "period", "first_treat", "dose", - aggregate="dose", ) att_diff = abs(results.overall_att - r_out["overall_att"]) / ( @@ -795,8 +790,11 @@ def test_benchmark_6_event_study(self): "period", "first_treat", "dose", - aggregate="eventstudy", ) + # M-025: keep validating the event-study code path (this + # benchmark's documented purpose) via the post-fit route - + # overall_att itself is aggregate-independent. + results.aggregate("event_study") # Compare overall ATT (binarized) att_diff = abs(results.overall_att - r_out["overall_att"]) / ( @@ -947,7 +945,6 @@ def _fit_cov(df, method, **kw): "period", "first_treat", "dose", - aggregate="dose", covariates=["x1", "x2"], ) @@ -1077,9 +1074,7 @@ def test_covariate_dgp_recovery(self): with warnings.catch_warnings(): warnings.simplefilter("ignore") - uncond = ContinuousDiD().fit( - df, "outcome", "unit", "period", "first_treat", "dose", aggregate="dose" - ) + uncond = ContinuousDiD().fit(df, "outcome", "unit", "period", "first_treat", "dose") assert abs(float(uncond.overall_att) - tau) > 0.5 # unconditional biased for method in ("reg", "dr"): res = _fit_cov(df, method) @@ -1268,7 +1263,6 @@ class TestDiscreteSaturated: time="period", first_treat="first_treat", dose="dose", - aggregate="dose", ) def test_hand_calc_att_acrt_overall(self): @@ -1388,7 +1382,6 @@ class TestLowestDose: time="period", first_treat="first_treat", dose="dose", - aggregate="dose", ) def test_dL_to_zero_exact_equivalence(self): diff --git a/tests/test_practitioner.py b/tests/test_practitioner.py index 7f72f6c56..6541850da 100644 --- a/tests/test_practitioner.py +++ b/tests/test_practitioner.py @@ -642,6 +642,49 @@ def _results(self): return r +class TestContinuousAggregationStep: + """M-025 post-fit aggregation guidance for ContinuousDiD. + + Deliberately NOT on ``_AggregationStepMixin``: that mixin pins + ``aggregate('group')`` (not in ContinuousDiD's supported set), + simulates bootstrap via a ``bootstrap_results`` field + (ContinuousDiDResults keys on ``n_bootstrap`` and has no such field), + and expects an ``aggregate='all'`` fallback (never a valid + ContinuousDiD value). ContinuousDiD's step is a single unconditional + step whose wording carries the bootstrap carve-out. + """ + + @staticmethod + def _agg_step(output): + return [s for s in output["next_steps"] if "Aggregate post-fit" in s["label"]] + + def _results(self): + r = ContinuousDiDResults.__new__(ContinuousDiDResults) + r.overall_att = 0.6 + r.overall_att_se = 0.15 + return r + + def test_aggregation_step_present_with_all_routes(self): + output = practitioner_next_steps(self._results(), verbose=False) + steps = self._agg_step(output) + assert len(steps) == 1 + assert "M-025" in steps[0]["why"] + assert "results.aggregate('event_study')" in steps[0]["code"] + assert "results.aggregate('dose')" in steps[0]["code"] + assert "results.aggregate('simple')" in steps[0]["code"] + # The bootstrap carve-out is carried in the wording. + assert "n_bootstrap=0" in steps[0]["why"] + + def test_aggregation_step_name_is_non_steps_key(self): + from diff_diff.practitioner import STEPS + + assert "aggregation" not in STEPS + output = practitioner_next_steps( + self._results(), completed_steps=["heterogeneity"], verbose=False + ) + assert len(self._agg_step(output)) == 1 + + # --------------------------------------------------------------------------- # Tests: unknown result type fallback # --------------------------------------------------------------------------- diff --git a/tests/test_profile_panel.py b/tests/test_profile_panel.py index 601f7b9b7..95ca202d5 100644 --- a/tests/test_profile_panel.py +++ b/tests/test_profile_panel.py @@ -1235,7 +1235,7 @@ def test_treatment_dose_min_flags_negative_dose_continuous_panels(): consistent with the dose column on per-unit treated/untreated status (negative-dose units labeled `first_treat > 0`), `ContinuousDiD.fit()` would raise `ValueError` at - `continuous_did.py:287-294` ("Dose must be strictly positive for + `continuous_did.py:586-593` ("Dose must be strictly positive for treated units (D > 0)"). `ContinuousDiD` as currently implemented does not apply on this panel. `HeterogeneousAdoptionDiD` is also NOT a routing alternative here: HAD requires non-negative diff --git a/tests/test_survey_phase3.py b/tests/test_survey_phase3.py index acc83ec34..eaa7f3703 100644 --- a/tests/test_survey_phase3.py +++ b/tests/test_survey_phase3.py @@ -1421,8 +1421,11 @@ def test_efficient_did_weights_only_se(self, staggered_survey_data): assert r.overall_se > 0 assert r.overall_se > 0.01 # Not artificially tiny + @pytest.mark.filterwarnings(r"ignore:ContinuousDiD\.fit\(aggregate=\):FutureWarning") def test_continuous_did_eventstudy_survey(self, continuous_survey_data): - """ContinuousDiD aggregate=eventstudy should work with survey design.""" + """The deprecated fit-time aggregate=eventstudy still works with a + survey design (M-025 legacy routing; the post-fit successor's + survey arms are pinned in tests/test_aggregate_contract.py).""" from diff_diff import ContinuousDiD sd = SurveyDesign(weights="weight", strata="stratum")