From 653c25795e8927b6aebfcac4ed5b59bed58a16fe Mon Sep 17 00:00:00 2001 From: igerber Date: Wed, 26 Aug 2026 12:10:24 -0400 Subject: [PATCH 1/2] feat(dml): DMLDiD panel=False repeated cross sections - Chang (2020) Case 2 (DML PR-B2) Ships Chang (2020) Case 2 (repeated cross sections) as an RCS lane inside DMLDiD, selected by a new `panel: bool = True` constructor parameter (mirrors CallawaySantAnna). Declared RCS only: one row per unit, unique row IDs required; no allow_unbalanced_panel RC-routing. Score family (diff_diff/_dr_scores.py): - chang_rcs_score: Eq 3.2 uncentered summand (D - g)/(p lam (1-lam)(1-g)) * ((T - lam) y - l2_hat), with the single control-only (T - lam) Y regression (Chang's I_kz^c) as the outcome nuisance and global within-cell p_hat = mean(D), lam_hat = mean(T). - chang_rcs_lambda_slope: sample-analogue G2_lambda from the closed-form d/d-lambda of psi_2 (the paper prints no estimator - documented implementation decision, REGISTRY Note). - chang_rcs_score_augmented: Theorem 2 psi_bar = summand - D theta/p + G2_lambda (T - lam); SE = sqrt(mean(psi_bar^2)/n_cell). The lambda-correction is mandatory (review warning quoted in docstring) and pinned by a regression test that recomputes the no-lambda SE. - Shared fail-closed validator (strict-binary D AND T, strictly interior p_hat/lam_hat, ps in [0,1)). Estimator lane (diff_diff/dml_did.py): - panel=False validation: unique unit IDs, covariates required, stationary-sampling UserWarning (Assumption 2.3) emitted only after the declared-RCS structure validates. - _precompute_rcs: per-row cohorts, is_panel False, canonical_size = n_obs; agg_cohort_masses deliberately unset (float-keyed lookups collide on >2^53 int64 labels; the bincount fallback is numerically identical under unique row IDs). - _compute_dml_rcs_gt: pooled two-period cells on level outcomes, FOUR-group guard (any empty treated/control x period group -> zero_treated_control), D x T 4-class stratified folds (singleton stratum -> cross_fit_degenerate), propensity clip-never-drop, lam_hat extremeness warning, per-observation IF payload psi_bar/n_cell, diagnostics gain lam_hat and g2_lambda. - RCS aggregation weights are the FIXED cohort row masses via per-cell agg_weight (CS-RCS convention, WIF-consistent SEs); aggregate('total') fails closed on RCS fits. - Results/reports design-aware: summary Design line + obs labels, BusinessReport cites Assumption 2.3 + 3.2(h) on RCS, practitioner snippet carries panel=False, target-parameter text names the cohort-mass weighting, DiagnosticReport skips the Goodman-Bacon check on RCS fits (auto-refit path only; precomputed passthrough honored). Validation (no Case 2 parity oracle exists - DoubleMLDIDCSBinary implements the Sant'Anna-Zhao 4-regression score and omits the lambda-correction): - benchmarks/doubleml/chang_rcs_characterization.py: CHARACTERIZATION spike (not parity) with pasted transcript; Part 1 documents nonzero Chang-vs-DoubleML gaps on shared folds, Part 2 golden literals pin public DMLDiD(panel=False) == hand Eq 3.2/Thm 2 pipeline at 0.0. - Equation-level fixtures at 1e-15/1e-12, finite-difference checks of G2_lambda and d/dp psi_2, oracle-nuisance recovery (true propensity + true l20), hand-pipeline replay at rtol 1e-14, DR in BOTH nuisance directions, ATT recovery on a fixed RCS DGP (theta_0 = 3), lambda-correction regression guard, native golden reproduction with live-fit gap pins, slow-lane Monte Carlo coverage, payload/idx contracts, degenerate handling with assert_nan_inference, >2^53 int64 cohort-label aggregation + bootstrap replay, variance- conventions dml_did_rcs row (table regenerated). Docs: REGISTRY DMLDiD section covers both cases (Case 2 equations, 8 new Notes, Case 2 assumption citations 3.2(a)/3.2(h) + Assumption 2.3 bullet, characterization reference block); dml_did.rst full rework (Case 2 methodology, declared-RCS restrictions, loud no-survey-weights note); chang-2020-review.md checklist flips; llms*.txt guides; README/index one-liners; choosing_estimator, migration-4.0, survey-roadmap, practitioner_decision_tree, REPORTING.md, variance-conventions prose, references.rst, doc-deps.yaml, CHANGELOG. ROADMAP DML section removed (Case 3 stays in DEFERRED.md); new TODO row tracks Chang section-4 RCS DGP replication. --- CHANGELOG.md | 14 + DEFERRED.md | 2 +- README.md | 2 +- ROADMAP.md | 10 - TODO.md | 3 +- .../doubleml/chang_rcs_characterization.py | 265 +++++++++ diff_diff/_dr_scores.py | 185 ++++++- diff_diff/_reporting_helpers.py | 24 +- diff_diff/business_report.py | 26 +- diff_diff/diagnostic_report.py | 18 + diff_diff/dml_did.py | 514 ++++++++++++++++-- diff_diff/dml_did_results.py | 40 +- diff_diff/guides/llms-autonomous.txt | 10 +- diff_diff/guides/llms-full.txt | 7 +- diff_diff/guides/llms-practitioner.txt | 5 +- diff_diff/guides/llms.txt | 2 +- diff_diff/mmm.py | 9 +- diff_diff/practitioner.py | 4 +- docs/api/dml_did.rst | 76 ++- docs/choosing_estimator.rst | 27 +- docs/doc-deps.yaml | 5 +- docs/index.rst | 2 +- docs/methodology/REGISTRY.md | 102 +++- docs/methodology/REPORTING.md | 12 + docs/methodology/papers/chang-2020-review.md | 20 +- docs/methodology/variance-conventions.md | 7 +- docs/migration-4.0.md | 2 +- docs/practitioner_decision_tree.rst | 7 +- docs/references.rst | 2 +- docs/survey-roadmap.md | 2 +- tests/test_dml_did.py | 513 +++++++++++++++++ tests/test_dr_scores.py | 98 ++++ tests/test_methodology_dml_did.py | 362 +++++++++++- tests/test_methodology_dr_scores.py | 169 +++++- tests/test_profile_panel.py | 3 +- tests/test_variance_conventions.py | 25 + 36 files changed, 2423 insertions(+), 151 deletions(-) create mode 100644 benchmarks/doubleml/chang_rcs_characterization.py diff --git a/CHANGELOG.md b/CHANGELOG.md index fae464658..0033bb570 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 HonestDiD / PreTrendsPower / DiagnosticReport / BusinessReport / MMM totals integration. See the REGISTRY "DMLDiD" section for equations and documented implementation choices. +- **`DMLDiD(panel=False)` — Chang (2020) Case 2 repeated cross sections**: + the same staggered per-(g,t)-cell architecture on DECLARED cross-sectional + data (one observation per row, row-unique unit IDs): pooled two-period + cells on level outcomes, the single control-only `(T - λ̂)·Y` outcome + nuisance (`chang_rcs_score`), D×T-stratified folds, and the λ-corrected + Theorem 2 variance (`chang_rcs_score_augmented` with the explicit + `Ĝ₂λ(T − λ̂)` term). Per-observation influence functions; aggregation + weights are fixed cohort row masses (WIF-consistent SEs); + `aggregate('total')` fails closed on RCS fits; no survey weights (weighted + RCS belongs to `CallawaySantAnna(panel=False, survey_design=...)`). + Validated by equation-level fixtures, oracle closed forms, + derivative-identity checks, double robustness in both directions, and a + committed `DoubleMLDIDCSBinary` characterization spike (no parity oracle + exists — DoubleML's RCS score differs and omits the λ term). ## [3.10.0] - 2026-08-22 diff --git a/DEFERRED.md b/DEFERRED.md index b6e35d19f..f0840b8ee 100644 --- a/DEFERRED.md +++ b/DEFERRED.md @@ -156,4 +156,4 @@ decisions (refactor waivers, perf trade-offs, test-infrastructure calls) are rec | **`ImputationDiD` SE vcov is already rank-guarded upstream.** Excluded from the structural rank-guard sweep: the lead/effect vcov comes from `solve_ols(..., return_vcov=True, rank_deficient_action=...)` at the OLS fit (`imputation.py:~2316`), which already drops rank-deficient columns. The only raw inverse (`solve(V_gamma, gamma)`, `imputation.py:~2530`) is the pretrends **Wald F-test statistic** with a safe `NaN` fallback — a test statistic, not a sandwich bread — so there is no garbage-SE exposure. No structural rank-guard needed. | `imputation.py` | structural-rank-guard / 2026-06-28 | | **TWFE HC2/HC2-BM full-dummy dedup: drift-prone duplication already resolved by the shared builder; full delegation waived.** The former Actionable row (origin: follow-up review, citing pre-#655 line numbers) asked to extract a shared dummy-construction helper or delegate TWFE's HC2/HC2-BM path to DiD's `fixed_effects=` branch. The shared-helper half SHIPPED in #655: both sites now delegate dummy construction, drop-first convention, FE column naming, and the duplicate-term backstop to the single `build_fe_dummy_blocks` (`utils.py`) + `validate_design_term_names` implementation (`twfe.py::fit` full-dummy branch; `estimators.py::DifferenceInDifferences.fit` `fixed_effects=` branch) — the FE-naming / survey-behavior drift risk the row targeted is gone. What remains per site is ~4 lines of genuinely estimator-specific design-matrix assembly (TWFE stacks `const`/`ATT`/covariates; DiD stacks its formula terms), which is not drift-prone duplication. The remaining full-delegation option — routing `TWFE.fit` through DiD machinery with TWFE-specific cluster-default threading — would touch TWFE's user-visible result surface (coefficient-dict keys, cluster-label conventions, warning text) for near-zero residual benefit; waived on cost/benefit. | `twfe.py::fit`, `estimators.py::DifferenceInDifferences.fit`, `utils.py::build_fe_dummy_blocks` | #655 / 2026-07-10 | | **Survey TSL SE intentionally counts genuine-subpopulation zero-weight PSUs (matches R, NOT a bug).** Recorded as the REGISTRY § "Subpopulation Analysis" TSL-meat Note (Lumley 2004 §3.4 full-design domain convention; R `survey::svyrecvar(subset())` parity); regression-locked by `tests/test_survey.py::TestZeroWeightPsuConventionWaiver`. | `survey.py` (`_compute_stratified_psu_meat`) | PR-B / 2026-06-30 | -| DMLDiD survey/cluster support: Chang (2020) assumes i.i.d. sampling, so `DMLDiD.fit()` accepts no `survey_design=`/`cluster=` (bare TypeError). The per-unit augmented-score influence function IS unit-level clustering (REGISTRY DMLDiD M-080 Note), so the 4.0 auto-cluster default flip is inert; what is missing is the COARSER-than-unit CR1 surface CS exposes and any design-based (survey) variance. The CS mixins were verified safe with survey keys absent (every read is `.get` with panel fallback), so the extension is additive: survey keys in `_precompute` + kit label threading (`_BOOTSTRAP_LABEL` reaches the survey-bootstrap <2-PSU warning) + a clustered-score variance derivation consistent with the cross-fitting. Needs a methodology decision on clustered cross-fitting (cluster-level folds vs unit folds with clustered scores). | `diff_diff/dml_did.py` | DML PR-B1 | Low | +| DMLDiD survey/cluster support: Chang (2020) assumes i.i.d. sampling, so `DMLDiD.fit()` accepts no `survey_design=`/`cluster=` (bare TypeError). The per-unit augmented-score influence function IS unit-level clustering (REGISTRY DMLDiD M-080 Note), so the 4.0 auto-cluster default flip is inert; what is missing is the COARSER-than-unit CR1 surface CS exposes and any design-based (survey) variance. The CS mixins were verified safe with survey keys absent (every read is `.get` with panel fallback), so the extension is additive: survey keys in `_precompute` AND `_precompute_rcs` (both design lanes) + kit label threading (`_BOOTSTRAP_LABEL` reaches the survey-bootstrap <2-PSU warning) + a clustered-score variance derivation consistent with the cross-fitting — per-UNIT influence functions on the panel lane, per-OBSERVATION on the declared-RCS lane (where clustering would group rows into design clusters). Needs a methodology decision on clustered cross-fitting (cluster-level folds vs unit folds with clustered scores). Note RCS data is typically survey data (BRFSS/ACS/CPS), so the RCS lane raises the priority of the survey half. | `diff_diff/dml_did.py` | DML PR-B1 | Low | diff --git a/README.md b/README.md index fc4ad2e5d..db347ed3b 100644 --- a/README.md +++ b/README.md @@ -121,7 +121,7 @@ Full guide: `diff_diff.get_llm_guide("practitioner")`. - [LPDiD](https://diff-diff.readthedocs.io/en/stable/api/lpdid.html) - Dube, Girardi, Jorda & Taylor (2025) Local Projections DiD: per-horizon long-difference event study on clean controls (no negative weighting), variance- or equally-weighted ATT, for absorbing or non-absorbing (reversible) treatment - [ChangesInChanges](https://diff-diff.readthedocs.io/en/stable/api/changes_in_changes.html) - Athey & Imbens (2006) nonlinear/distributional DiD for the 2x2 design: full counterfactual distribution and quantile treatment effects via CDF transformation, plus the QDiD comparison estimator via `method="qdid"`; bootstrap inference; R qte parity. Alias `CiC` - [LWDiD](https://diff-diff.readthedocs.io/en/stable/api/lwdid.html) - Lee & Wooldridge (2025, 2026) rolling-transformation DiD: unit-specific demean/detrend converts panel to cross-section, staggered adoption, `estimation_method` in `reg`/`ipw`/`dr`/`psm` (the papers' RA/IPW/IPWRA plus propensity-score matching), exact small-N inference on the classical collapsed regression -- [DMLDiD](https://diff-diff.readthedocs.io/en/stable/api/dml_did.html) - Chang (2020) double/debiased machine learning DiD: staggered ATT(g,t) with cross-fitted ML nuisance learners (DML2) and Neyman-orthogonal scores, for flexible/high-dimensional covariate adjustment under conditional parallel trends +- [DMLDiD](https://diff-diff.readthedocs.io/en/stable/api/dml_did.html) - Chang (2020) double/debiased machine learning DiD: staggered ATT(g,t) with cross-fitted ML nuisance learners (DML2) and Neyman-orthogonal scores, for flexible/high-dimensional covariate adjustment under conditional parallel trends; panel or declared repeated cross sections (`panel=False`) - [BaconDecomposition](https://diff-diff.readthedocs.io/en/stable/api/bacon.html) - Goodman-Bacon (2021) decomposition for diagnosing TWFE bias in staggered settings ## Diagnostics & Sensitivity diff --git a/ROADMAP.md b/ROADMAP.md index e4856dae6..efd956d3b 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -105,16 +105,6 @@ Unified framework encompassing synthetic control and regression approaches via l **Reference**: Athey et al. (2021), *Journal of the American Statistical Association*. -### Double / Debiased ML for DiD — repeated cross sections - -The panel lane shipped as `DMLDiD` (Chang 2020 Case 1, staggered ATT(g,t)). -Remaining scope: Chang's Case 2 repeated-cross-section score (Equation 3.2) -with the lambda-corrected variance (a proposed `chang_rcs_score` helper -alongside the existing panel scores in `diff_diff/_dr_scores.py`), and its -staggered RCS lane. - -**Reference**: Chang, N.-C. (2020), *The Econometrics Journal* 23(2), 177-191. - ### Alternative Inference Methods - **Randomization inference**: exact p-values for small samples. diff --git a/TODO.md b/TODO.md index 3b295879a..fde3148c0 100644 --- a/TODO.md +++ b/TODO.md @@ -82,7 +82,8 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m |-------|----------|--------|--------|----------| | `CallawaySantAnnaResults.summary(alpha=...)` (and sibling staggered summaries relaying the parent renderer) relabels the confidence-interval header at the requested alpha while always printing the FIT-TIME `overall_conf_int` — silent mislabeling on any non-fit alpha (bootstrap percentile intervals cannot be reconstructed from the SE at all). DMLDiDResults now rejects non-fit alpha (the EventStudyResults.summary precedent, results_base.py:643-660); apply the same guard (or a real recomputation on analytical fits) to the CS-family summaries | `diff_diff/staggered_results.py` | DML PR-B1 review | Quick | Medium | | `plot_event_study`'s pointwise CI reconstruction (`_event_study.py:302-304`) gates on `np.isfinite(std_err)` but not `> 0` — a zero-SE period draws a finite zero-width interval while its stored inference is all-NaN (the CS-surface twin of the `plot_group_effects` alternate-alpha gate fixed in DML PR-B1; pre-existing main behavior, cross-surface-twins audit) — apply the same `se > 0` NaN-gate | `diff_diff/visualization/_event_study.py` | DML PR-B1 review | Quick | Medium | -| DMLDiD tutorial notebook (CONTRIBUTING's new-estimator checklist requires a tutorial; deferred to its own PR per the numbers-locked notebook protocol — prototype in scripts, lock numbers, assemble + execute once, register in `docs/tutorials/index.rst` with a toctree short label + group card): staggered DGP with nonlinear covariate confounding, learner comparison (linear/ridge/sieve/sklearn object), post-fit aggregation + HonestDiD via the event-study container, seed/reproducibility note | `docs/tutorials/`, `docs/tutorials/index.rst` | DML PR-B1 | Mid | Medium | +| DMLDiD tutorial notebook (CONTRIBUTING's new-estimator checklist requires a tutorial; deferred to its own PR per the numbers-locked notebook protocol — prototype in scripts, lock numbers, assemble + execute once, register in `docs/tutorials/index.rst` with a toctree short label + group card): staggered DGP with nonlinear covariate confounding, learner comparison (linear/ridge/sieve/sklearn object), post-fit aggregation + HonestDiD via the event-study container, seed/reproducibility note, and a `panel=False` repeated-cross-section example | `docs/tutorials/`, `docs/tutorials/index.rst` | DML PR-B1 | Mid | Medium | +| Replicate Chang (2020) §4's own RCS simulation DGPs (pp. 17-21, "fully specified" per the paper review) as recovery/coverage fixtures for the `DMLDiD(panel=False)` lane — the shipped tests use a library-authored RCS design (documented in the REGISTRY checklist caveat); needs the paper PDF to extract the parameterization | `tests/test_methodology_dml_did.py`, `docs/methodology/papers/chang-2020-review.md` | DML PR-B2 | Mid | Low | | Optional scheduled end-to-end execution gate for the MMM tutorials (29/30): a cron-only workflow (or extension of `mmm-interop.yml`) that executes both notebooks in isolated exact-pin environments, so a stale/invalid committed posterior cannot stay green indefinitely - today the hybrid posture (deliberate: notebooks execute locally with committed outputs; CI smoke-tests the exporters without sampling; drift tests pin source + committed-output needles) leaves the MCMC claims un-re-executed in CI | `.github/workflows/mmm-interop.yml`, `docs/tutorials/29_mmm_calibration_pymc.ipynb`, `docs/tutorials/30_mmm_calibration_meridian.ipynb` | mmm-interop | Mid | Low | | Committed `fixest::feols` event-study golden for TWFE `event_study=True` (within + pooled specs, unbalanced + covariate panels, matched CR1 cluster convention, per-period effects + vcov block) - the in-suite gates are shared-core cross-checks (TWFE-within == MPD-absorb, pooled == MPD bit-exact), so a defect common to the shared core would pass; the live-R harness (`benchmarks/R/benchmark_multiperiod.R`, `feols(y ~ treated * time_f \| unit)`) validated the within design in `docs/benchmarks.rst` but is not a committed regression test - follow the `fixest_did_twfe_golden.json` committed-golden pattern (pytest.skip when absent) | `tests/test_fixest_did_twfe_parity.py`, `benchmarks/R/` | 3(a) R2 | Mid | Medium | | Type-blind `n_bootstrap` acceptance in already-validated estimators - HAD bool (`isinstance(..., int)` passes `True`, runs as 1 replicate), dCDH bool+float (its bare `< 0` check passes both `True` and `2.5`), TROP float (`2.5` passes the `>= 2` floor), SyntheticDiD float under all three variance methods + bool/negative under jackknife (its floor check is skipped there) - align these local checks with the `utils.validate_n_bootstrap` type guard (M-081 kept them out of the sweep: it scoped to previously-UNvalidated estimators only) | `diff_diff/had.py`, `diff_diff/chaisemartin_dhaultfoeuille.py`, `diff_diff/trop.py`, `diff_diff/synthetic_did.py` | 2(d) PR-B | Quick | Low | diff --git a/benchmarks/doubleml/chang_rcs_characterization.py b/benchmarks/doubleml/chang_rcs_characterization.py new file mode 100644 index 000000000..6c8a70e93 --- /dev/null +++ b/benchmarks/doubleml/chang_rcs_characterization.py @@ -0,0 +1,265 @@ +"""DoubleML RCS CHARACTERIZATION spike for the DMLDiD panel=False lane. + +CHARACTERIZATION, NOT PARITY. ``DoubleMLDIDCSBinary`` (doubleml 0.11.4; +``score="observational"``, ``in_sample_normalization=False``, trimming 0.01) +implements the Sant'Anna-Zhao repeated-cross-section score — FOUR +treatment-by-period outcome regressions — while Chang (2020) Equation 3.2 +uses ONE control-only regression of ``(T - lam) * y`` on X, and +DoubleML's variance omits Chang's Theorem 2 ``G_2lambda * (T - lam)`` +correction entirely. The two estimators are consistent for the same +estimand but carry different first-order influence functions and +finite-sample values: ATT gaps are O_p(N^{-1/2}) by construction and SE +gaps additionally carry the lambda term. This spike DOCUMENTS both gaps per cell +(and reports the lambda-omitted Chang SE beside them); the machine-precision +assertions here are SELF-parity — the public ``DMLDiD(panel=False)`` +estimator vs the hand-rolled Chang pipeline under DMLDiD's own +reconstructed folds (Part 2), whose 12-decimal literals are the goldens +consumed by tests/test_methodology_dml_did.py::TestRCSGoldenCharacterization +(native reproduction swaps sklearn's lbfgs logit for the library IRLS +solver; tolerance atol 2e-4 ATT / 1e-5 SE, the B0/B1 precedent). + +Environment (side venv; doubleml/sklearn are NEVER diff-diff dependencies): + + python -m venv .venv-doubleml + .venv-doubleml/bin/pip install "doubleml==0.11.4" scikit-learn + .venv-doubleml/bin/python benchmarks/doubleml/chang_rcs_characterization.py + +Observed transcript (2026-08-26, doubleml 0.11.4, sklearn 1.9.0, macOS arm64): + + Part 1 — Chang Eq 3.2 vs DoubleMLDIDCSBinary (shared folds; gaps EXPECTED): + cell (g=3, t=3, base=2) CHANG ATT = 1.875325609339 SE(full) = 0.243971253081 SE(no-lam) = 0.249877072301 | DML ATT = 2.121300323064 SE = 0.109099376634 | ATT gap -2.460e-01 SE gap +1.349e-01 SE(no-lam) gap +1.408e-01 + cell (g=3, t=4, base=2) CHANG ATT = 2.360914148301 SE(full) = 0.252067309584 SE(no-lam) = 0.258532386882 | DML ATT = 2.293055793734 SE = 0.106740749269 | ATT gap +6.786e-02 SE gap +1.453e-01 SE(no-lam) gap +1.518e-01 + cell (g=4, t=3, base=2) CHANG ATT = -0.075262310611 SE(full) = 0.187206946073 SE(no-lam) = 0.187207274463 | DML ATT = -0.057963388953 SE = 0.106181054127 | ATT gap -1.730e-02 SE gap +8.103e-02 SE(no-lam) gap +8.103e-02 + cell (g=4, t=4, base=3) CHANG ATT = 1.799124069195 SE(full) = 0.260544978066 SE(no-lam) = 0.265671872127 | DML ATT = 1.995045506760 SE = 0.105732614154 | ATT gap -1.959e-01 SE gap +1.548e-01 SE(no-lam) gap +1.599e-01 + CHARACTERIZATION OK (nonzero, bounded gaps — different scores, same estimand) + + Part 2 — public DMLDiD(panel=False) vs hand Chang (DMLDiD's own folds): + cell (g=3, t=3, base=2) DMLDiD ATT = 1.867260790882 SE = 0.244040011648 | diff vs hand ATT +0.000e+00 SE +0.000e+00 + cell (g=3, t=4, base=2) DMLDiD ATT = 2.347579854514 SE = 0.251882814022 | diff vs hand ATT +0.000e+00 SE +0.000e+00 + cell (g=4, t=3, base=2) DMLDiD ATT = -0.089752999600 SE = 0.187258740736 | diff vs hand ATT +0.000e+00 SE +0.000e+00 + cell (g=4, t=4, base=3) DMLDiD ATT = 1.822066292980 SE = 0.260472487961 | diff vs hand ATT +0.000e+00 SE +0.000e+00 + SELF-PARITY OK (public DMLDiD cells within 1e-10 of the hand Chang pipeline) + +Interpretation. The ATT gaps (1.7e-2 to 2.5e-1) are finite-sample +differences between two DIFFERENT orthogonal scores estimating the same +estimand — the pre-registered outcome for a characterization, and exactly +why the REGISTRY records DoubleMLDIDCSBinary as "not an oracle" for +Case 2. The SE columns are NOT directly comparable beyond direction: on +top of the missing lambda term, DoubleMLDIDCSBinary scatters its +per-cell score to the FULL frame with zero fill (did_cs_binary.py's +_set_id_positions), so its variance scaling factor is the full N rather +than the cell's pooled row count — the transcript records the observed +values rather than claiming a decomposition of the gap. The lambda +term's own effect is isolated in the SE(full)-vs-SE(no-lam) columns +(same pipeline, one term removed): up to ~2.6% relative on these cells, +and near-zero exactly where lam_hat ~ 0.5 (cell (4,3): the (1-2*lam) +factor vanishes) — matching the closed form. Part 2's zero diffs are the +real precision anchor: the shipped estimator IS the hand-rolled +Equation 3.2 / Theorem 2 pipeline under identical folds and learners. +""" + +import numpy as np +import pandas as pd +from doubleml.data import DoubleMLPanelData +from doubleml.did import DoubleMLDIDCSBinary +from sklearn.linear_model import LinearRegression, LogisticRegression + +SEED = 7 +N_ROWS = 4000 +PERIODS = [1, 2, 3, 4] +COHORTS = [0, 3, 4] # never-treated + two staggered cohorts +K = 5 +TRIM = 1e-2 + +rng = np.random.default_rng(SEED) +cohort = rng.choice(COHORTS, size=N_ROWS, p=[0.5, 0.25, 0.25]) +tt = rng.choice(PERIODS, size=N_ROWS) +X_row = rng.standard_normal((N_ROWS, 2)) +y = ( + 0.8 * X_row[:, 0] + - 0.4 * X_row[:, 1] + + 0.3 * tt + + 0.5 * X_row[:, 0] * tt / 4 # covariate-dependent trend (X matters) + + rng.standard_normal(N_ROWS) +) +post = (cohort > 0) & (tt >= cohort) +y = y + post * (2.0 + 0.2 * (tt - cohort)) # heterogeneous dynamic effect + +df = pd.DataFrame( + { + "id": np.arange(N_ROWS), + "t": tt, + "g": cohort, + "y": y, + "x1": X_row[:, 0], + "x2": X_row[:, 1], + } +) +panel = DoubleMLPanelData(df, y_col="y", d_cols="g", t_col="t", id_col="id", x_cols=["x1", "x2"]) + +CELLS = [(3, 3, 2), (3, 4, 2), (4, 3, 2), (4, 4, 3)] # (g, t_eval, base) + + +def chang_cell(idx, t_eval, g_val, smpls): + """Hand-rolled Chang Eq 3.2 cell under the given fold splits. + + Returns (theta, se_full, se_no_lambda) — se_full carries the Theorem 2 + G_2lambda correction, se_no_lambda deliberately omits it (the + 'plausible implementation bug' the review warns against), so the + transcript records how much the lambda term moves the SE. + """ + D = (cohort[idx] == g_val).astype(float) + T = (tt[idx] == t_eval).astype(float) + yv = y[idx] + X = X_row[idx] + n = idx.shape[0] + p_hat = D.mean() + lam = T.mean() + oof_g = np.empty(n) + oof_l = np.empty(n) + r = (T - lam) * yv + for tr, te in smpls: + lg = LogisticRegression(penalty=None, solver="lbfgs", max_iter=1000).fit(X[tr], D[tr]) + oof_g[te] = np.clip(lg.predict_proba(X[te])[:, 1], TRIM, 1 - TRIM) + ctrl = tr[D[tr] == 0] + oof_l[te] = LinearRegression().fit(X[ctrl], r[ctrl]).predict(X[te]) + w = (D - oof_g) / (p_hat * lam * (1 - lam) * (1 - oof_g)) + summand = w * ((T - lam) * yv - oof_l) + theta = summand.mean() + odds = (D - oof_g) / (1 - oof_g) + g2 = np.mean( + -((1 - 2 * lam) / (lam**2 * (1 - lam) ** 2)) * (odds / p_hat) * ((T - lam) * yv - oof_l) + - (yv / (p_hat * lam * (1 - lam))) * odds + ) + psi_full = summand - D * theta / p_hat + g2 * (T - lam) + psi_no_lam = summand - D * theta / p_hat + se_full = np.sqrt(np.mean(psi_full**2) / n) + se_no_lam = np.sqrt(np.mean(psi_no_lam**2) / n) + return float(theta), float(se_full), float(se_no_lam) + + +# --------------------------------------------------------------------------- +# Part 1 — characterization vs DoubleMLDIDCSBinary under SHARED folds. +# --------------------------------------------------------------------------- +print("Part 1 — Chang Eq 3.2 vs DoubleMLDIDCSBinary (shared folds; gaps EXPECTED):") +gaps_nonzero = True +for g_val, t_eval, base in CELLS: + in_cell = ((cohort == g_val) | (cohort == 0)) & ((tt == t_eval) | (tt == base)) + idx = np.flatnonzero(in_cell) + n = idx.shape[0] + D = (cohort[idx] == g_val).astype(float) + T = (tt[idx] == t_eval).astype(float) + + cell_rng = np.random.default_rng(1000 + g_val * 10 + t_eval) + strata = (D + 2 * T).astype(int) + test_folds = [[] for _ in range(K)] + cursor = 0 + for s_val in np.unique(strata): + members = np.flatnonzero(strata == s_val) + members = cell_rng.permutation(members) + for m_i in members: + test_folds[cursor % K].append(m_i) + cursor += 1 + test_folds = [np.sort(np.array(f, dtype=int)) for f in test_folds] + smpls = [(np.setdiff1d(np.arange(n), te), te) for te in test_folds] + + m = DoubleMLDIDCSBinary( + panel, + g_value=g_val, + t_value_pre=base, + t_value_eval=t_eval, + ml_g=LinearRegression(), + ml_m=LogisticRegression(penalty=None, solver="lbfgs", max_iter=1000), + control_group="never_treated", + n_folds=K, + n_rep=1, + score="observational", + in_sample_normalization=False, + trimming_threshold=TRIM, + draw_sample_splitting=False, + ) + m.set_sample_splitting([smpls]) + m.fit() + att_dml, se_dml = float(m.coef[0]), float(m.se[0]) + + theta, se_full, se_no_lam = chang_cell(idx, t_eval, g_val, smpls) + d_att = theta - att_dml + print( + f"cell (g={g_val}, t={t_eval}, base={base}) CHANG ATT = {theta:.12f} " + f"SE(full) = {se_full:.12f} SE(no-lam) = {se_no_lam:.12f} | " + f"DML ATT = {att_dml:.12f} SE = {se_dml:.12f} | " + f"ATT gap {d_att:+.3e} SE gap {se_full - se_dml:+.3e} " + f"SE(no-lam) gap {se_no_lam - se_dml:+.3e}" + ) + # Characterization honesty: the scores are DIFFERENT — a zero gap would + # mean the two implementations coincide, which they must not. + gaps_nonzero &= abs(d_att) > 1e-12 + # Sanity band: both estimate the same estimand. + gaps_nonzero &= abs(d_att) < 0.5 + +print( + "CHARACTERIZATION OK (nonzero, bounded gaps — different scores, same estimand)" + if gaps_nonzero + else "CHARACTERIZATION FAILED" +) + +# --------------------------------------------------------------------------- +# Part 2 — SELF-parity: public DMLDiD(panel=False) vs the hand-rolled Chang +# pipeline under DMLDiD's OWN reconstructed folds (sklearn learners on both +# sides — machine precision expected). The printed 12-decimal CHANG-side +# numbers are the goldens for TestRCSGoldenCharacterization. +# --------------------------------------------------------------------------- +import sys # noqa: E402 + +sys.path.insert(0, ".") +from diff_diff import DMLDiD # noqa: E402 +from diff_diff._crossfit import assign_folds # noqa: E402 + +SEED_DML = 11 +est = DMLDiD( + propensity_learner=LogisticRegression(penalty=None, solver="lbfgs", max_iter=1000), + outcome_learner=LinearRegression(), + seed=SEED_DML, + n_folds=K, + pscore_trim=TRIM, + panel=False, +) +import warnings # noqa: E402 + +with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = est.fit(df, outcome="y", unit="id", time="t", first_treat="g", covariates=["x1", "x2"]) + +sorted_cohorts = [3, 4] +sorted_periods = PERIODS +ok2 = True +print("\nPart 2 — public DMLDiD(panel=False) vs hand Chang (DMLDiD's own folds):") +for g_val, t_eval, base in CELLS: + entry = res.group_time_effects[(g_val, t_eval)] + in_cell = ((cohort == g_val) | (cohort == 0)) & ((tt == t_eval) | (tt == base)) + idx = np.flatnonzero(in_cell) + n = idx.shape[0] + D = (cohort[idx] == g_val).astype(float) + T = (tt[idx] == t_eval).astype(float) + g_idx = sorted_cohorts.index(g_val) + t_idx = sorted_periods.index(t_eval) + rng2 = np.random.default_rng(np.random.SeedSequence(entropy=SEED_DML, spawn_key=(g_idx, t_idx))) + folds = assign_folds(n, K, rng=rng2, stratify=D + 2.0 * T) + smpls = [ + (np.flatnonzero(folds.fold_ids != k), np.flatnonzero(folds.fold_ids == k)) for k in range(K) + ] + theta, se_full, _ = chang_cell(idx, t_eval, g_val, smpls) + d_att = float(entry["effect"]) - theta + d_se = float(entry["se"]) - se_full + print( + f"cell (g={g_val}, t={t_eval}, base={base}) DMLDiD ATT = {entry['effect']:.12f} " + f"SE = {entry['se']:.12f} | diff vs hand ATT {d_att:+.3e} SE {d_se:+.3e}" + ) + ok2 &= abs(d_att) < 1e-10 and abs(d_se) < 1e-10 + +print( + "SELF-PARITY OK (public DMLDiD cells within 1e-10 of the hand Chang pipeline)" + if ok2 + else "SELF-PARITY FAILED" +) +raise SystemExit(0 if (gaps_nonzero and ok2) else 1) diff --git a/diff_diff/_dr_scores.py b/diff_diff/_dr_scores.py index b3a9a2f7e..043ce1320 100644 --- a/diff_diff/_dr_scores.py +++ b/diff_diff/_dr_scores.py @@ -1,6 +1,6 @@ -"""Shared doubly-robust panel DiD scores (private DML infrastructure). +"""Shared doubly-robust DiD scores (private DML infrastructure). -Two distinct score families live here — the repo's methodology review +Three distinct score families live here — the repo's methodology review (``docs/methodology/papers/chang-2020-review.md``) pins that they are NOT interchangeable: @@ -15,6 +15,16 @@ cross-fitted (DML) nuisances, normalized by the unconditional treated share ``p_hat``, with the augmented variant carrying the finite-dimensional treated-share variance correction (``G_1p = -theta/p_hat``). +- :func:`chang_rcs_score` / :func:`chang_rcs_lambda_slope` / + :func:`chang_rcs_score_augmented` — the Chang (2020) Case 2 (repeated + cross sections) score on LEVEL outcomes with the post-period sampling share + ``lam_hat``, whose Theorem 2 variance carries BOTH finite-dimensional + corrections (the treated-share fold-in plus an EXPLICIT + ``G_2lambda * (T - lam_hat)`` term). The Case 2 outcome nuisance is a + SINGLE control-only regression of ``(T - lam_hat) * Y`` on X — deliberately + different from Sant'Anna-Zhao/DoubleML's four treatment-by-period outcome + regressions, so ``doubleml.DoubleMLDIDCSBinary`` is a characterization + anchor, not a parity oracle. References ---------- @@ -35,6 +45,9 @@ "drdid_panel_inf_func", "chang_panel_score", "chang_panel_score_augmented", + "chang_rcs_score", + "chang_rcs_lambda_slope", + "chang_rcs_score_augmented", ] @@ -243,3 +256,171 @@ def chang_panel_score_augmented( "(p_hat = 1 means no comparison population; the estimand is unidentified)" ) return summand - D * theta / p_hat + + +def _validate_chang_rcs_inputs( + y: np.ndarray, + D: np.ndarray, + T: np.ndarray, + m2_hat: np.ndarray, + ps: np.ndarray, + p_hat: float, + lam_hat: float, + context: str, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Shared validation for the Chang Case 2 (repeated cross sections) scores.""" + y = np.asarray(y, dtype=np.float64) + D = np.asarray(D, dtype=np.float64) + T = np.asarray(T, dtype=np.float64) + m2_hat = np.asarray(m2_hat, dtype=np.float64) + ps = np.asarray(ps, dtype=np.float64) + for name, arr in (("y", y), ("D", D), ("T", T), ("m2_hat", m2_hat), ("ps", ps)): + if arr.ndim != 1: + raise ValueError(f"{context}: {name} must be 1-dimensional, got ndim={arr.ndim}") + n = y.shape[0] + if n == 0: + raise ValueError(f"{context}: inputs are empty (n=0); the score is undefined") + for name, arr in (("y", y), ("D", D), ("T", T), ("m2_hat", m2_hat), ("ps", ps)): + if arr.shape[0] != n: + raise ValueError( + f"{context}: {name} has length {arr.shape[0]}, expected {n} (length of y)" + ) + if not np.all(np.isfinite(arr)): + raise ValueError(f"{context}: {name} contains non-finite values") + if not np.all((D == 0.0) | (D == 1.0)): + raise ValueError(f"{context}: D must be strictly binary 0/1") + if not np.all((T == 0.0) | (T == 1.0)): + raise ValueError(f"{context}: T must be strictly binary 0/1") + if not np.isfinite(p_hat) or not (0.0 < p_hat < 1.0): + raise ValueError( + f"{context}: p_hat must satisfy 0 < p_hat < 1 strictly, got {p_hat!r} " + "(p_hat = 1 means no comparison population; the estimand is unidentified)" + ) + if not np.isfinite(lam_hat) or not (0.0 < lam_hat < 1.0): + raise ValueError( + f"{context}: lam_hat must satisfy 0 < lam_hat < 1 strictly, got " + f"{lam_hat!r} (lam_hat = 0 or 1 means all observations lie in one " + "period; 1/(lam*(1-lam)) is undefined and the estimand is " + "unidentified)" + ) + if np.any(ps < 0.0) or np.any(ps >= 1.0): + raise ValueError( + f"{context}: ps must lie in [0, 1) strictly; clip/trim propensity " + "scores before calling (values >= 1 would divide by zero)" + ) + return y, D, T, m2_hat, ps + + +def chang_rcs_score( + y: np.ndarray, + D: np.ndarray, + T: np.ndarray, + m2_hat: np.ndarray, + ps: np.ndarray, + p_hat: float, + lam_hat: float, +) -> np.ndarray: + """Chang (2020) Case 2 per-observation UNCENTERED score summand. + + Returns ``summand_i = (D_i - ps_i) / (p_hat * lam_hat * (1 - lam_hat) * + (1 - ps_i)) * ((T_i - lam_hat) * y_i - m2_hat_i)``, whose sample mean is + the ATT (Equation 3.2's ``psi_2`` equals ``summand - theta``). ``m2_hat`` + is the SINGLE Case 2 outcome nuisance ``l_20(X) = E[(T - lam) * Y | X, + D=0]`` — one cross-fitted regression of ``(T - lam_hat) * Y`` on X trained + on control observations only (Chang's ``I_kz^c``), NOT the four + treatment-by-period regressions of the Sant'Anna-Zhao/DoubleML RCS score. + ``ps`` is cross-fitted and trimmed by the CALLER's policy. + + ``p_hat`` and ``lam_hat`` are caller-supplied; the library convention + (REGISTRY "DMLDiD" Notes) is the FULL-SAMPLE-within-cell treated share and + post-period sampling share (``mean(D)`` / ``mean(T)``), mirroring the + Case 1 global-``p_hat`` convention and DoubleML's ``t.mean()``. + """ + y, D, T, m2_hat, ps = _validate_chang_rcs_inputs( + y, D, T, m2_hat, ps, p_hat, lam_hat, "chang_rcs_score" + ) + weight = (D - ps) / (p_hat * lam_hat * (1.0 - lam_hat) * (1.0 - ps)) + return weight * ((T - lam_hat) * y - m2_hat) + + +def chang_rcs_lambda_slope( + y: np.ndarray, + D: np.ndarray, + T: np.ndarray, + m2_hat: np.ndarray, + ps: np.ndarray, + p_hat: float, + lam_hat: float, +) -> float: + """Chang (2020) Case 2 lambda-slope estimator ``G_2lambda``. + + Sample mean of the closed-form derivative ``d/d(lambda) psi_2`` evaluated + at the plug-in nuisances (recovered from the proof of Theorem 2, p. 55 + display; the paper prints NO explicit ``G_2lambda`` estimator — this + natural sample analogue is a documented implementation decision, REGISTRY + "DMLDiD" Note):: + + d_lam psi_2_i = -((1 - 2*lam) / (lam**2 * (1-lam)**2)) + * ((D_i - ps_i) / (p_hat * (1 - ps_i))) + * ((T_i - lam) * y_i - m2_hat_i) + - (y_i / (p_hat * lam * (1-lam))) + * ((D_i - ps_i) / (1 - ps_i)) + + Only consistency is required of this estimator (Theorem 2 imposes no + rate). The first term equals ``-((1 - 2*lam) / (lam*(1-lam))) * + summand_i`` — an algebraic identity the test suite cross-checks. + """ + y, D, T, m2_hat, ps = _validate_chang_rcs_inputs( + y, D, T, m2_hat, ps, p_hat, lam_hat, "chang_rcs_lambda_slope" + ) + odds = (D - ps) / (1.0 - ps) + term1 = ( + -((1.0 - 2.0 * lam_hat) / (lam_hat**2 * (1.0 - lam_hat) ** 2)) + * (odds / p_hat) + * ((T - lam_hat) * y - m2_hat) + ) + term2 = -(y / (p_hat * lam_hat * (1.0 - lam_hat))) * odds + return float(np.mean(term1 + term2)) + + +def chang_rcs_score_augmented( + summand: np.ndarray, + D: np.ndarray, + T: np.ndarray, + y: np.ndarray, + m2_hat: np.ndarray, + ps: np.ndarray, + theta: float, + p_hat: float, + lam_hat: float, +) -> np.ndarray: + """Chang (2020) Case 2 augmented score with BOTH finite-dim corrections. + + Returns ``psi_bar_i = summand_i - D_i * theta / p_hat + G_2lambda * + (T_i - lam_hat)`` — Theorem 2's combined score: the treated-share + correction ``G_2p = -theta/p_hat`` folds into the score exactly as in + Case 1, while the lambda-correction stays an EXPLICIT extra term + (``G_2lambda`` computed internally via :func:`chang_rcs_lambda_slope`). + The variance estimator is ``SE = sqrt(mean(psi_bar**2) / N)``. + + Per the methodology review: "Omitting the λ-correction term is a + plausible implementation bug the proof structure warns against" — the + bare ``psi_2`` squared is NOT the Theorem 2 estimator, and no DoubleML + parity anchor exists for this object (``DoubleMLDIDCSBinary``'s variance + omits the lambda term; see the committed characterization spike). + """ + context = "chang_rcs_score_augmented" + summand = np.asarray(summand, dtype=np.float64) + if summand.ndim != 1: + raise ValueError(f"{context}: summand must be 1-dimensional") + if summand.shape[0] == 0: + raise ValueError(f"{context}: inputs are empty (n=0); the score is undefined") + if not np.all(np.isfinite(summand)): + raise ValueError(f"{context}: summand contains non-finite values") + if not np.isfinite(theta): + raise ValueError(f"{context}: theta must be finite, got {theta!r}") + y, D, T, m2_hat, ps = _validate_chang_rcs_inputs(y, D, T, m2_hat, ps, p_hat, lam_hat, context) + if summand.shape[0] != y.shape[0]: + raise ValueError(f"{context}: summand has length {summand.shape[0]}, expected {y.shape[0]}") + g2_lambda = chang_rcs_lambda_slope(y, D, T, m2_hat, ps, p_hat, lam_hat) + return summand - D * theta / p_hat + g2_lambda * (T - lam_hat) diff --git a/diff_diff/_reporting_helpers.py b/diff_diff/_reporting_helpers.py index 6473b50fe..456d46c49 100644 --- a/diff_diff/_reporting_helpers.py +++ b/diff_diff/_reporting_helpers.py @@ -150,16 +150,28 @@ def describe_target_parameter(results: Any) -> Dict[str, Any]: } if name == "DMLDiDResults": + if getattr(results, "panel", True) is False: + weight_name = "cohort-mass-weighted" + weight_clause = ( + "cell weights are FIXED cohort row masses (the CS-RCS " + "convention, WIF-consistent; the per-cell complete-case " + "``n_treated`` is display-only; REGISTRY DMLDiD RCS Note)" + ) + else: + weight_name = "valid-treated-count-weighted" + weight_clause = ( + "cell weights are the per-cell " + "complete-case ``n_treated`` (full cohort masses enter only " + "the aggregation weight influence function; REGISTRY DMLDiD " + "complete-case Note)" + ) return { - "name": "overall ATT (valid-treated-count-weighted average of ATT(g,t))", + "name": f"overall ATT ({weight_name} average of ATT(g,t))", "definition": ( - "A per-cell valid-treated-count-weighted average of " + f"A per-cell {weight_name} average of " "group-time ATTs ``ATT(g, t)`` across post-anticipation " "cells (``t >= g - anticipation``; ``t >= g`` when " - "anticipation=0) — cell weights are the per-cell " - "complete-case ``n_treated`` (full cohort masses enter only " - "the aggregation weight influence function; REGISTRY DMLDiD " - "complete-case Note) — " + "anticipation=0) — " + weight_clause + " — " "where each cell is a conditional-on-covariates ATT " "estimated by Chang (2020)'s cross-fitted Neyman-orthogonal " "score (DML2). ``overall_att`` is the simple-aggregation " diff --git a/diff_diff/business_report.py b/diff_diff/business_report.py index e1e2e3328..7d740da1b 100644 --- a/diff_diff/business_report.py +++ b/diff_diff/business_report.py @@ -1575,6 +1575,26 @@ def _describe_assumption(estimator_name: str, results: Any = None) -> Dict[str, ), } if estimator_name == "DMLDiDResults": + is_rcs = getattr(results, "panel", True) is False + if is_rcs: + tail = ( + "cross-fitting (DML2) removing own-observation overfitting. " + "The repeated-cross-section design (Chang Case 2) " + "additionally assumes STATIONARY cross-sectional sampling " + "(Assumption 2.3: each wave samples the same target " + "population, so the composition of (D, X) is stable across " + "waves while outcomes are the period-specific potential " + "outcomes — not data-checkable), and valid " + "normal inference requires the nuisance learners to satisfy " + "Chang (2020)'s Case 2 rate conditions (Assumption 3.2(h))." + ) + else: + tail = ( + "cross-fitting (DML2) removing own-observation overfitting; " + "valid normal inference additionally requires the nuisance " + "learners to satisfy Chang (2020)'s rate conditions " + "(Assumption 3.1(f))." + ) return { "parallel_trends_variant": "conditional_on_covariates", "no_anticipation": True, @@ -1585,11 +1605,7 @@ def _describe_assumption(estimator_name: str, results: Any = None) -> Dict[str, "supplied covariates, per treatment cohort and period " "(group-time ATT), plus no anticipation. The " "Neyman-orthogonal score makes the estimate first-order " - "insensitive to machine-learning regularization bias, with " - "cross-fitting (DML2) removing own-observation overfitting; " - "valid normal inference additionally requires the nuisance " - "learners to satisfy Chang (2020)'s rate conditions " - "(Assumption 3.1(f))." + "insensitive to machine-learning regularization bias, with " + tail ), } if estimator_name in { diff --git a/diff_diff/diagnostic_report.py b/diff_diff/diagnostic_report.py index 22ad9a20b..367d3f89f 100644 --- a/diff_diff/diagnostic_report.py +++ b/diff_diff/diagnostic_report.py @@ -1265,6 +1265,24 @@ def _instance_skip_reason(self, check: str) -> Optional[str]: "``precomputed={'bacon': ...}`` with a survey-aware " "decomposition." ) + # Declared repeated-cross-section fits (panel=False on any + # panel-attribute-bearing producer — DMLDiD and CS alike): + # BaconDecomposition refits a TWFE decomposition that needs a + # real panel. On one-observation-per-unit data it does not + # raise — treatment is collinear with the absorbed unit fixed + # effects and it returns a meaningless twfe_estimate of 0.0 — + # so skip with an explicit reason. Placed AFTER the precomputed + # passthrough above: a user-supplied precomputed Bacon result + # (e.g. computed against a real panel elsewhere) stays honored. + if getattr(r, "panel", True) is False: + return ( + "Goodman-Bacon decomposition requires panel data; this " + "fit used a declared repeated cross section " + "(panel=False), where the TWFE replay is degenerate " + "(treatment collinear with unit fixed effects). Supply " + "``precomputed={'bacon': ...}`` computed on a real " + "panel if one exists." + ) return None if check == "heterogeneity": # Needs multiple group or event-study effects. Use len() rather than diff --git a/diff_diff/dml_did.py b/diff_diff/dml_did.py index c3e98ddca..a0de7a09d 100644 --- a/diff_diff/dml_did.py +++ b/diff_diff/dml_did.py @@ -1,21 +1,27 @@ -"""DMLDiD: Double/Debiased Machine Learning DiD (Chang 2020), staggered panel. - -Implements Chang (2020, The Econometrics Journal 23(2), 177-191) Case 1 -(panel) as a staggered ATT(g,t) estimator: each Callaway-Sant'Anna style -(g, t) cell is a 2-period Chang problem — cross-fitted nuisances -(propensity ``g_hat`` and control outcome-change regression ``m_hat``), -the Neyman-orthogonal score ``psi_1`` (``chang_panel_score``), and the -augmented-score plug-in variance (``chang_panel_score_augmented``; -``SE = sqrt(mean(psi_bar**2)/n)``). The classic 2-period design is the -degenerate single-cell case. Covariates are REQUIRED (Chang's estimator -exists for the high-dimensional-X setting; use CallawaySantAnna without -covariates). +"""DMLDiD: Double/Debiased Machine Learning DiD (Chang 2020), staggered. + +Implements Chang (2020, The Econometrics Journal 23(2), 177-191) as a +staggered ATT(g,t) estimator: each Callaway-Sant'Anna style (g, t) cell is +a 2-period Chang problem. ``panel=True`` (default) runs Case 1 (repeated +outcomes) — cross-fitted nuisances (propensity ``g_hat`` and control +outcome-change regression ``m_hat``), the Neyman-orthogonal score +``psi_1`` (``chang_panel_score``), and the augmented-score plug-in +variance (``chang_panel_score_augmented``; +``SE = sqrt(mean(psi_bar**2)/n)``). ``panel=False`` runs Case 2 (declared +repeated cross sections) — level outcomes, the single control-only +``(T - lam_hat) * Y`` outcome nuisance (``chang_rcs_score``), and the +λ-corrected Theorem 2 variance (``chang_rcs_score_augmented``). The +classic 2-period design is the degenerate single-cell case. Covariates +are REQUIRED (Chang's estimator exists for the high-dimensional-X +setting; use CallawaySantAnna without covariates). ``DMLDiD`` writes the CallawaySantAnna per-(g,t) ``influence_func_info`` -payload (per-unit entries ``psi_bar_i / n_cell``, so ``sqrt(sum(if**2))`` -IS the cell SE) and inherits the CS aggregation + multiplier-bootstrap -mixins: event study with sup-t bands, group/simple/total aggregation, and -post-fit ``results.aggregate()`` with bootstrap replay. +payload (per-sampling-unit entries ``psi_bar_i / n_cell`` — per unit on +panel fits, per observation on RCS fits — so ``sqrt(sum(if**2))`` IS the +cell SE) and inherits the CS aggregation + multiplier-bootstrap mixins: +event study with sup-t bands, group/simple aggregation (plus total on +panel fits only; RCS fits fail ``total`` closed), and post-fit +``results.aggregate()`` with bootstrap replay. See docs/methodology/REGISTRY.md "DMLDiD" for equations, implementation Notes (global p-hat, D-stratified folds, pooled fold weighting, trimming, @@ -32,7 +38,13 @@ from diff_diff._base import BaseEstimator from diff_diff._crossfit import DegenerateFoldError, assign_folds, cross_fit_predict -from diff_diff._dr_scores import chang_panel_score, chang_panel_score_augmented +from diff_diff._dr_scores import ( + chang_panel_score, + chang_panel_score_augmented, + chang_rcs_lambda_slope, + chang_rcs_score, + chang_rcs_score_augmented, +) from diff_diff._learners import ( _CLASSIFIER_NAMES, _REGRESSOR_NAMES, @@ -221,14 +233,23 @@ def _validate_seed(seed: Any) -> Optional[int]: class DMLDiD(CallawaySantAnnaBootstrapMixin, CallawaySantAnnaAggregationMixin, BaseEstimator): - """Chang (2020) DML DiD: staggered panel ATT(g,t) with cross-fitted ML nuisances. + """Chang (2020) DML DiD: staggered ATT(g,t) with cross-fitted ML nuisances. + + ``panel=True`` (default) estimates Case 1 on panel data; ``panel=False`` + estimates Case 2 on declared repeated cross sections (level outcomes, + λ-corrected variance). - Per (g, t) cell: D-stratified K-fold cross-fitting of the propensity + Per (g, t) cell: K-fold cross-fitting of the propensity (``propensity_learner``, out-of-fold ``predict_proba``) and the - control-only outcome-change regression (``outcome_learner``, trained on - the cell's controls, Chang's ``I_kz^c``), then the pooled orthogonal - score mean and the augmented-score plug-in SE. Aggregation (event study, - group, simple, total) is POST-FIT via ``results.aggregate()``. + control-only outcome regression (``outcome_learner``, trained on the + cell's controls, Chang's ``I_kz^c``), then the pooled orthogonal score + mean and the augmented-score plug-in SE. On panel fits the folds are + D-stratified and the outcome nuisance is the outcome-change regression + ``E[dY | X, D=0]``; on RCS fits the folds are D x T stratified and the + nuisance is the level regression ``E[(T - lam)Y | X, D=0]`` with the + lambda-corrected variance. Aggregation (event study, group, simple; + plus total on panel fits — RCS fits fail ``total`` closed) is POST-FIT + via ``results.aggregate()``. Parameters ---------- @@ -271,6 +292,20 @@ class DMLDiD(CallawaySantAnnaBootstrapMixin, CallawaySantAnnaAggregationMixin, B Propensity clip bound: fitted propensities are clipped to ``[pscore_trim, 1 - pscore_trim]`` after the extremeness warning (clip, never drop; Chang's paper gives no trimming rule). + panel : bool, default True + ``True`` estimates Chang's Case 1 (repeated outcomes) on panel data + (one row per unit-period, cell score on outcome CHANGES). ``False`` + estimates Case 2 (repeated cross sections) on DECLARED + cross-sectional data — one observation per row with a row-unique + ``unit`` ID, cell score on outcome LEVELS with the post-period + sampling share λ̂ and the λ-corrected Theorem 2 variance. Case 2 + additionally assumes stationary cross-sectional sampling (Chang + Assumption 2.3: each wave samples the SAME target population — + the composition of ``(D, X)`` is stable across waves, while + outcomes are the period-specific potential outcomes, so trends + and treatment effects are expected, not violations), which is not + data-checkable — ``fit()`` warns. ``aggregate('total')`` is unavailable on RCS fits (fails + closed, the library-wide RC convention). """ _BOOTSTRAP_LABEL: ClassVar[str] = "DMLDiD" @@ -289,6 +324,7 @@ def __init__( base_period: str = "varying", cband: bool = True, pscore_trim: float = 0.01, + panel: bool = True, ) -> None: # Raw assignment, then ONE shared validator (also re-run at the top # of fit() as the direct-mutation defense) validates and normalizes @@ -305,6 +341,7 @@ def __init__( self.base_period = base_period self.cband = cband self.pscore_trim = pscore_trim + self.panel = panel self._revalidate_config() # Fitted-state lifecycle (house convention; not inherited under this MRO). @@ -363,6 +400,13 @@ def _revalidate_config(self) -> None: ) self.cband = bool(self.cband) self.pscore_trim = _validate_pscore_trim(self.pscore_trim) + if not isinstance(self.panel, (bool, np.bool_)): + raise ValueError( + f"panel must be a bool, got {self.panel!r} (type " + f"{type(self.panel).__name__}) — truthy strings like 'False' " + "would silently select the panel lane" + ) + self.panel = bool(self.panel) # ------------------------------------------------------------------ # Input validation @@ -534,6 +578,38 @@ def _validate_and_prepare( ) df[col] = converted.astype(np.float64) + # Structural design guards. Panel: one row per (unit, time) and a + # time-invariant cohort label. Declared RCS: one row per unit — with + # row-unique IDs the two panel guards below are subsumed (duplicate + # (unit, time) is impossible and per-row first_treat is trivially + # unit-constant). + if not self.panel: + if df[unit].duplicated().any(): + raise ValueError( + "panel=False requires unique unit IDs (one observation per " + "unit). Found duplicate unit IDs. If your data is a panel, " + "use panel=True." + ) + self._check_control_group_availability(df, unit, first_treat) + # Emit only AFTER the declared-RCS structure has fully + # validated — a failing call raises without the misleading + # suggestion that estimation began. + warnings.warn( + "panel=False uses Chang (2020) Case 2 repeated-cross-section " + "scores, which assume stationary cross-sectional sampling " + "(Assumption 2.3): each wave samples the SAME target " + "population — conditional on the period, rows are i.i.d. " + "draws from the distribution of (Y(0), D, X) (pre) or " + "(Y(1), D, X) (post), so the composition of (D, X) is stable " + "across waves while outcomes are the period-specific " + "potential outcomes (trends and treatment effects are " + "expected, not violations). This assumption is not " + "data-checkable.", + UserWarning, + stacklevel=3, + ) + return df, covariates + # Duplicate (unit, time) rows. dup_mask = df.duplicated(subset=[unit, time], keep=False) if dup_mask.any(): @@ -552,7 +628,18 @@ def _validate_and_prepare( "must be constant per unit" ) - # Control-group availability (CS parity, incl. the NESTING). + self._check_control_group_availability(df, unit, first_treat) + + return df, covariates + + def _check_control_group_availability( + self, df: pd.DataFrame, unit: str, first_treat: str + ) -> None: + """Control-group availability (CS parity, incl. the NESTING). + + Shared by both design lanes: on declared RCS the groupby-first + degenerates correctly to per-row under the unique-row-ID guard. + """ unit_cohorts_series = df.groupby(unit)[first_treat].first() n_never = int((unit_cohorts_series == 0).sum()) n_cohorts = int(unit_cohorts_series[unit_cohorts_series > 0].nunique()) @@ -568,8 +655,6 @@ def _validate_and_prepare( "cohorts when there are no never-treated units." ) - return df, covariates - def _normalize_label_columns( self, time_vals: np.ndarray, @@ -839,6 +924,63 @@ def _precompute( "n_units": n_units, } + def _precompute_rcs( + self, + df: pd.DataFrame, + outcome: str, + unit: str, + time: str, + first_treat: str, + covariates: List[str], + ) -> Dict[str, Any]: + """Declared-RCS bookkeeping: rows ARE the sampling units. + + Mirrors the minimal CS RCS precompute contract: ``all_units`` = + observation positions, ``unit_to_idx`` = None, per-ROW + ``unit_cohorts``, ``canonical_size`` = n_obs — the aggregation and + bootstrap mixins branch off exactly these (per-row multiplier + weights; per-obs cohort bincount pg basis; ``aggregate('total')`` + fails closed on ``is_panel: False``). + + ``agg_cohort_masses`` is deliberately NOT set: under unique row IDs + the aggregation cache's bincount fallback yields exactly the fixed + per-cohort row masses (never-treated included, denominator = n_obs), + so the WIF pg basis is already the fixed-cohort-mass one and the + per-cell ``agg_weight`` (from ``rcs_cohort_masses``, int-keyed by the + canonical labels) aligns the point-estimate weights with it. + Supplying the key would be a numeric no-op that routes lookups + through float()-keyed dict access, where distinct int64 cohorts + above 2**53 (admissible through 2**62 by the label pipeline) + collide. ``obs_per_unit`` and survey keys likewise omitted (true RCS + is one row per unit — a non-None obs_per_unit would divide the WIF). + """ + n_obs = len(df) + unit_cohorts = df[first_treat].to_numpy() + treatment_groups = sorted(g for g in df[first_treat].unique() if g > 0) + time_periods = sorted(df[time].unique()) + period_to_col = {t: i for i, t in enumerate(time_periods)} + return { + "all_units": np.arange(n_obs), + "unit_to_idx": None, + "unit_cohorts": unit_cohorts, + "obs_time": df[time].to_numpy(), + "obs_outcome": df[outcome].to_numpy(), + "obs_covariates": df[covariates].to_numpy(dtype=np.float64), + "cohort_masks": {g: unit_cohorts == g for g in treatment_groups}, + # +inf cohorts were recoded to 0 upstream, so == 0 is complete. + "never_treated_mask": unit_cohorts == 0, + "period_to_col": period_to_col, + "observed_sorted": sorted(period_to_col), + "time_periods": time_periods, + "treatment_groups": treatment_groups, + "is_panel": False, + "canonical_size": n_obs, + "n_units": n_obs, + "rcs_cohort_masses": { + g: int(np.count_nonzero(unit_cohorts == g)) for g in treatment_groups + }, + } + def _sanitize_learner_error(self, exc: BaseException) -> str: """Persisted error text for cross_fit_diagnostics / to_dict exports. @@ -1082,6 +1224,273 @@ def _compute_dml_gt( } return gt_entry, if_entry, diagnostics + def _compute_dml_rcs_gt( + self, + precomputed: Dict[str, Any], + g: Any, + t: Any, + g_idx: int, + t_idx: int, + root_entropy: int, + dropped_units_out: Optional[set] = None, + ) -> Tuple[Dict[str, Any], Optional[Dict[str, Any]], Optional[Dict[str, Any]]]: + """One RCS (g, t) cell — Chang Case 2. Same return contract as the + panel cell: (gt_entry, if_entry|None, diagnostics|None). + + Pooled two-period cell on LEVEL outcomes: four disjoint row groups + (treated/control x {t, base}), one cross-fit propensity on the pooled + rows, and the SINGLE Case 2 outcome nuisance — a control-only (both + periods, Chang's I_kz^c) regression of (T - lam_hat) * y on X. + """ + observed_sorted = precomputed["observed_sorted"] + period_to_col = precomputed["period_to_col"] + base = _select_base_period_impl(self.base_period, self.anticipation, g, t, observed_sorted) + if base is None or base not in period_to_col or t not in period_to_col: + return _nan_gt_entry(skip_reason="missing_period"), None, None + + unit_cohorts = precomputed["unit_cohorts"] + treated_mask = precomputed["cohort_masks"][g] + if self.control_group == "never_treated": + control_mask = precomputed["never_treated_mask"] + else: # not_yet_treated (CS semantics: untreated at max(t, base) + k) + nyt_threshold = max(t, base) + self.anticipation + control_mask = precomputed["never_treated_mask"] | ( + (unit_cohorts > nyt_threshold) & (unit_cohorts != g) + ) + + obs_time = precomputed["obs_time"] + y_obs = precomputed["obs_outcome"] + X_obs = precomputed["obs_covariates"] + at_t = obs_time == t + at_base = obs_time == base + # Per-ROW complete cases (RCS covariates live on the row itself; no + # base-period X exists). + valid = np.isfinite(y_obs) & np.all(np.isfinite(X_obs), axis=1) + + treated_t = treated_mask & at_t & valid + treated_b = treated_mask & at_base & valid + control_t = control_mask & at_t & valid + control_b = control_mask & at_base & valid + if dropped_units_out is not None: + dropped_units_out.update( + np.flatnonzero((treated_mask | control_mask) & (at_t | at_base) & ~valid).tolist() + ) + n_treated = int(np.sum(treated_t | treated_b)) + n_control = int(np.sum(control_t | control_b)) + # FOUR-group guard (CS RCS precedent): the Case 2 estimand needs + # treated AND control rows in BOTH periods; skip vocabulary reuses + # zero_treated_control for any empty group. + if ( + min( + int(treated_t.sum()), + int(treated_b.sum()), + int(control_t.sum()), + int(control_b.sum()), + ) + == 0 + ): + return ( + _nan_gt_entry( + n_treated=n_treated, + n_control=n_control, + skip_reason="zero_treated_control", + ), + None, + None, + ) + + cell_mask = treated_t | treated_b | control_t | control_b + cell_idx = np.flatnonzero(cell_mask) + n_cell = cell_idx.shape[0] + D_cell = treated_mask[cell_idx].astype(np.float64) + T_cell = at_t[cell_idx].astype(np.float64) + y_cell = y_obs[cell_idx] + X_cell = X_obs[cell_idx] + # Global-within-cell shares (REGISTRY convention: mirrors the Case 1 + # global p_hat and DoubleML's d.mean()/t.mean()); strictly interior + # BY the four-group guard. + p_hat = float(D_cell.mean()) + lam_hat = float(T_cell.mean()) + if min(p_hat, 1.0 - p_hat) < self.pscore_trim: + warnings.warn( + f"DMLDiD cell (g={g}, t={t}): empirical treated share " + f"p_hat={p_hat:.4f} over the pooled two-period rows " + f"({n_treated} treated / {n_control} control) is extreme " + f"(min(p, 1-p) < pscore_trim={self.pscore_trim}); the Chang " + "score and variance scale with powers of 1/p_hat, so this " + "cell's estimate may be unstable.", + UserWarning, + stacklevel=3, + ) + if min(lam_hat, 1.0 - lam_hat) < self.pscore_trim: + warnings.warn( + f"DMLDiD cell (g={g}, t={t}): post-period sampling share " + f"lam_hat={lam_hat:.4f} is extreme (min(lam, 1-lam) < " + f"pscore_trim={self.pscore_trim}); Chang's Case 2 bounds " + "carry powers of 1/(lam*(1-lam)) up to cubes, so this cell's " + "estimate may be unstable.", + UserWarning, + stacklevel=3, + ) + + # Per-cell fold draw seeded exactly like the panel lane; strata are + # the FOUR D x T classes (DoubleML's d + 2t encoding — REGISTRY + # deviation Note): every training complement then carries control + # rows in both periods by construction, Chang's fold-composition + # requirement for fitting l_2 on I_kz^c. + seed_seq = np.random.SeedSequence(entropy=root_entropy, spawn_key=(g_idx, t_idx)) + rng = np.random.default_rng(seed_seq) + diagnostics: Dict[str, Any] = { + "propensity": None, + "outcome": None, + "p_hat": float(p_hat), + "lam_hat": float(lam_hat), + "n_clipped_ps": None, + "fold_seed": {"entropy": int(root_entropy), "spawn_key": [int(g_idx), int(t_idx)]}, + } + + try: + folds = assign_folds(n_cell, self.n_folds, rng=rng, stratify=D_cell + 2.0 * T_cell) + except ValueError as exc: + # Cell smaller than n_folds, or a singleton D x T stratum (e.g. + # ONE treated row in the base period cannot be cross-fitted). + diagnostics["skip_reason"] = "cross_fit_degenerate" + diagnostics["error"] = str(exc) + return ( + _nan_gt_entry( + n_treated=n_treated, + n_control=n_control, + skip_reason="cross_fit_degenerate", + ), + None, + diagnostics, + ) + + context = f"DMLDiD (g={g}, t={t})" + try: + with np.errstate(over="ignore", invalid="ignore"): + ps_res = cross_fit_predict( + make_learner(self.propensity_learner, kind="classifier"), + X_cell, + D_cell, + folds, + predict_method="predict_proba", + context_label=f"{context} propensity", + ) + r_cell = (T_cell - lam_hat) * y_cell + or_res = cross_fit_predict( + make_learner(self.outcome_learner, kind="regressor"), + X_cell, + r_cell, + folds, + predict_method="predict", + fit_mask=(D_cell == 0.0), + context_label=f"{context} outcome", + ) + except DegenerateFoldError as exc: + diagnostics["skip_reason"] = "cross_fit_degenerate" + diagnostics["error"] = self._sanitize_learner_error(exc) + return ( + _nan_gt_entry( + n_treated=n_treated, + n_control=n_control, + skip_reason="cross_fit_degenerate", + ), + None, + diagnostics, + ) + + ps_raw = ps_res.oof_predictions + _check_propensity_diagnostics(ps_raw, self.pscore_trim) + ps = np.clip(ps_raw, self.pscore_trim, 1.0 - self.pscore_trim) + n_clipped = int(np.sum((ps_raw < self.pscore_trim) | (ps_raw > 1.0 - self.pscore_trim))) + m2_hat = or_res.oof_predictions + + diagnostics["propensity"] = { + "fold_losses": [float(v) for v in ps_res.fold_losses], + "n_fit_per_fold": [int(v) for v in ps_res.n_fit_per_fold], + } + diagnostics["outcome"] = { + "fold_losses": [float(v) for v in or_res.fold_losses], + "n_fit_per_fold": [int(v) for v in or_res.n_fit_per_fold], + } + diagnostics["n_clipped_ps"] = n_clipped + + try: + with np.errstate(over="ignore", invalid="ignore"): + summand = chang_rcs_score(y_cell, D_cell, T_cell, m2_hat, ps, p_hat, lam_hat) + theta = float(np.mean(summand)) + psi_bar = chang_rcs_score_augmented( + summand, D_cell, T_cell, y_cell, m2_hat, ps, theta, p_hat, lam_hat + ) + g2_lambda = chang_rcs_lambda_slope( + y_cell, D_cell, T_cell, m2_hat, ps, p_hat, lam_hat + ) + se = float(np.sqrt(np.mean(psi_bar**2) / n_cell)) + except ValueError as exc: + diagnostics["skip_reason"] = "non_finite_score" + diagnostics["error"] = self._sanitize_learner_error(exc) + return ( + _nan_gt_entry( + n_treated=n_treated, + n_control=n_control, + skip_reason="non_finite_score", + ), + None, + diagnostics, + ) + if not ( + np.isfinite(theta) + and np.isfinite(se) + and np.isfinite(g2_lambda) + and np.all(np.isfinite(psi_bar)) + ): + diagnostics["skip_reason"] = "non_finite_score" + return ( + _nan_gt_entry( + n_treated=n_treated, + n_control=n_control, + skip_reason="non_finite_score", + ), + None, + diagnostics, + ) + diagnostics["g2_lambda"] = float(g2_lambda) + + t_stat, p_value, conf_int = safe_inference(theta, se, alpha=self.alpha) + gt_entry = { + "effect": theta, + "se": se, + "t_stat": t_stat, + "p_value": p_value, + "conf_int": conf_int, + # DISPLAY counts: pooled two-period valid rows. Aggregation + # weights come from agg_weight below (fixed cohort row mass, the + # CS-RCS convention) — never from these counts. + "n_treated": n_treated, + "n_control": n_control, + "skip_reason": None, + "agg_weight": precomputed["rcs_cohort_masses"][g], + } + + # Payload: per-OBS entries psi_bar_i / n_cell over BOTH periods' + # rows, so sqrt(sum(if^2)) IS the cell SE. cell_idx is a flatnonzero + # of one mask => strictly increasing and duplicate-free, and the + # D-partition keeps treated_idx/control_idx disjoint (the fancy-+= + # scatter contract in the aggregation layer). + n_units = precomputed["n_units"] + inf_full = np.zeros(n_units) + inf_full[cell_idx] = psi_bar / n_cell + treated_idx = cell_idx[D_cell == 1.0].astype(np.int64) + control_idx = cell_idx[D_cell == 0.0].astype(np.int64) + if_entry = { + "treated_idx": treated_idx, + "control_idx": control_idx, + "treated_inf": inf_full[treated_idx], + "control_inf": inf_full[control_idx], + } + return gt_entry, if_entry, diagnostics + # ------------------------------------------------------------------ # Fit # ------------------------------------------------------------------ @@ -1099,7 +1508,12 @@ def fit( df, covariates = self._validate_and_prepare( data, outcome, unit, time, first_treat, covariates ) - precomputed = self._precompute(df, outcome, unit, time, first_treat, covariates) + if self.panel: + precomputed = self._precompute(df, outcome, unit, time, first_treat, covariates) + cell_fn = self._compute_dml_gt + else: + precomputed = self._precompute_rcs(df, outcome, unit, time, first_treat, covariates) + cell_fn = self._compute_dml_rcs_gt treatment_groups = precomputed["treatment_groups"] time_periods = precomputed["time_periods"] observed_sorted = precomputed["observed_sorted"] @@ -1119,7 +1533,7 @@ def fit( self.base_period, self.anticipation, g, time_periods, observed_sorted ): t_idx = time_periods.index(t) - gt_entry, if_entry, diagnostics = self._compute_dml_gt( + gt_entry, if_entry, diagnostics = cell_fn( precomputed, g, t, g_idx, t_idx, root_entropy, dropped_units ) group_time_effects[(g, t)] = gt_entry @@ -1168,17 +1582,29 @@ def fit( # the sole user-visible trace). Accumulated per cell during the # estimation loop — no second O(n_units x n_cells) sweep. if dropped_units: - warnings.warn( - f"{len(dropped_units)} unit(s) were excluded from at least one " - "(group, time) cell they would otherwise join, due to a " - "missing or NON-FINITE outcome, a non-finite covariate at the " - "cell's base period, or an outcome difference overflowing to " - "non-finite. DMLDiD estimates each cell on its complete cases " - "(point weights use per-cell valid counts; aggregation cohort " - "masses use full cohorts — see REGISTRY.md).", - UserWarning, - stacklevel=2, - ) + if self.panel: + warnings.warn( + f"{len(dropped_units)} unit(s) were excluded from at least one " + "(group, time) cell they would otherwise join, due to a " + "missing or NON-FINITE outcome, a non-finite covariate at the " + "cell's base period, or an outcome difference overflowing to " + "non-finite. DMLDiD estimates each cell on its complete cases " + "(point weights use per-cell valid counts; aggregation cohort " + "masses use full cohorts — see REGISTRY.md).", + UserWarning, + stacklevel=2, + ) + else: + warnings.warn( + f"{len(dropped_units)} observation(s) were excluded from a " + "(group, time) cell they would otherwise join, due to a " + "missing or NON-FINITE outcome or a non-finite covariate " + "on the row. DMLDiD estimates each cell on its complete " + "cases (aggregation weights use fixed cohort row masses — " + "see REGISTRY.md).", + UserWarning, + stacklevel=2, + ) # Universal base period: per-cohort zero reference cells (full # nine-key CS dict; the kit hard-reads effect AND n_treated, the @@ -1194,7 +1620,7 @@ def fit( cohort_mass = float(np.count_nonzero(unit_cohorts == g)) if cohort_mass <= 0: continue - group_time_effects[(g, base)] = { + ref_entry: Dict[str, Any] = { "effect": 0.0, "se": np.nan, "t_stat": np.nan, @@ -1205,6 +1631,13 @@ def fit( "skip_reason": None, "is_reference": True, } + if not self.panel: + # Keep the RCS pg basis uniform: reference cells carry the + # same fixed cohort row mass as estimated cells. + ref_entry["agg_weight"] = precomputed["rcs_cohort_masses"].get( + g, int(round(cohort_mass)) + ) + group_time_effects[(g, base)] = ref_entry influence_func_info[(g, base)] = { "treated_idx": np.array([], dtype=np.int64), "control_idx": np.array([], dtype=np.int64), @@ -1313,6 +1746,7 @@ def fit( n_bootstrap=self.n_bootstrap, bootstrap_weights=self.bootstrap_weights, cband=self.cband, + panel=self.panel, ) results._aggregation_kit = _build_aggregation_kit( cast(Any, self), # duck-typed host contract (alpha/anticipation/cband) diff --git a/diff_diff/dml_did_results.py b/diff_diff/dml_did_results.py index 01f23f01e..3707f8d79 100644 --- a/diff_diff/dml_did_results.py +++ b/diff_diff/dml_did_results.py @@ -1,18 +1,20 @@ -"""Results container for the DMLDiD estimator (Chang 2020, staggered panel). +"""Results container for the DMLDiD estimator (Chang 2020, staggered). ``DMLDiDResults`` subclasses :class:`~diff_diff.staggered_results.CallawaySantAnnaResults` and inherits its kit-based post-fit ``aggregate()`` machinery (event study -with sup-t bands, group/simple/total aggregation, bootstrap replay). The -subclass adds the cross-fitting provenance the parent has no concept of: -learner specs, fold count, per-cell cross-fit diagnostics, and the -inference-provenance fields (``seed``/``n_bootstrap``/``bootstrap_weights``/ -``cband``) that move estimates or inference. +with sup-t bands, group/simple aggregation, bootstrap replay; ``total`` on +panel fits — repeated-cross-section fits fail it closed). The subclass adds +the cross-fitting provenance the parent has no concept of: learner specs, +fold count, per-cell cross-fit diagnostics, and the inference-provenance +fields (``seed``/``n_bootstrap``/``bootstrap_weights``/``cband``) that move +estimates or inference. ``vcov_type`` stays at the inherited ``"hc1"``: in this library ``hc1`` with -``cluster=None`` IS the per-unit influence-function variance by definition -(REGISTRY.md "IF-based variance estimators..." — the default), and DMLDiD's -augmented-score SE ``sqrt(mean(psi_bar**2)/n)`` is exactly a per-unit IF -variance, so ``"hc1"`` is the family-consistent, truthful value. +``cluster=None`` IS the per-sampling-unit influence-function variance by +definition (REGISTRY.md "IF-based variance estimators..." — the default), +and DMLDiD's augmented-score SE ``sqrt(mean(psi_bar**2)/n)`` is exactly +that — per UNIT on panel fits, per OBSERVATION on repeated-cross-section +fits (rows are the sampling units there). """ from dataclasses import dataclass, field @@ -53,11 +55,12 @@ def _convert(value: Any) -> Any: @dataclass class DMLDiDResults(CallawaySantAnnaResults): - """Results from DMLDiD (Chang 2020 DML DiD, staggered panel ATT(g,t)). + """Results from DMLDiD (Chang 2020 DML DiD, staggered ATT(g,t)). Inherits the full Callaway-Sant'Anna results surface — ``att``/``se`` aliases, ``to_dataframe``, post-fit ``aggregate()`` (simple / event_study - / group / total) with bootstrap replay — and adds the DML provenance + / group, plus total on panel fits; repeated-cross-section fits fail + ``total`` closed) with bootstrap replay — and adds the DML provenance fields below. Every inherited CS-only field that DMLDiD never populates (``epv_*``, ``pscore_fallback``, ``allow_unbalanced_panel``, ``used_rc_on_unbalanced_panel``, ``cluster_name``, ``n_clusters``, @@ -66,8 +69,10 @@ class DMLDiDResults(CallawaySantAnnaResults): ``event_study_df``) stays at its inherited default and is inert. ``group_effects`` stays ``None`` permanently — ``aggregate("group")`` returns a separate container and never mutates this object. ``panel`` - is NOT inert (report rendering reads it for the sample-count label); - its inherited ``True`` is correct for the panel-only design. + is NOT inert: it carries the fit's DECLARED design — ``True`` for the + Case 1 panel lane, ``False`` for the Case 2 repeated-cross-section lane + (``DMLDiD(panel=False)``) — and report rendering plus the aggregation + ``n_kind`` read it for "units" vs "observations" semantics. Attributes ---------- @@ -169,7 +174,12 @@ def summary(self, alpha: Optional[float] = None) -> str: for entry in self.cross_fit_diagnostics.values() if entry.get("skip_reason") is not None ) - header = [ + header = [] + if self.panel is False: + # Conditional line only (panel output stays byte-stable): the + # declared design changes score, variance, and count semantics. + header.append(f"{'Design:':<30} {'repeated cross sections':>10}") + header += [ f"{'Propensity learner:':<30} {self.propensity_learner!r:>10}", f"{'Outcome learner:':<30} {self.outcome_learner!r:>10}", f"{'Cross-fitting folds (K):':<30} {self.n_folds:>10}", diff --git a/diff_diff/guides/llms-autonomous.txt b/diff_diff/guides/llms-autonomous.txt index aab49a748..e7cc3c58f 100644 --- a/diff_diff/guides/llms-autonomous.txt +++ b/diff_diff/guides/llms-autonomous.txt @@ -349,7 +349,7 @@ supported / out of scope; `warn` supported but with documented caveats; | `StackedDiD` | ✓ | ✓ | ✗ | ✗ | ✗ | ✓ | ✗ | ✗ | ✓ | | `WooldridgeDiD` (ETWFE) | ✓ | ✓ | ✗ | ✗ | ✗ | ✓ | ✗ | ✗ | ✓ | | `LWDiD` | ✓ | ✓ | ✗ | ✗ | ✓ | warn | partial | ✗ | ✓ | -| `DMLDiD` | ✓ | ✓ | ✗ | ✗ | partial | ✓ (REQUIRED) | ✗ | ✗ | ✗ (per-unit IF variance IS unit clustering; no coarser cluster=) | +| `DMLDiD` | ✓ | ✓ | ✗ | ✗ | partial | ✓ (REQUIRED) | ✗ | ✗ | ✗ (per-sampling-unit IF variance IS unit clustering; no coarser cluster=; panel=False = declared RCS) | | `EfficientDiD` | ✓ | ✓ | ✗ | ✗ | partial | ✓ | ✗ | ✗ | ✓ | | `SyntheticDiD` | ✓ | ✗ | ✗ | ✗ | ✓ | ✓ | ✓ | ✗ | partial | | `TROP` | ✓ | ✓ | ✗ | ✗ | ✗ | ✗ | ✓ | ✗ | partial | @@ -763,6 +763,14 @@ Explicit RCS support in this library: - `CallawaySantAnna(panel=False)` - repeated-cross-section mode per REGISTRY.md §CallawaySantAnna; use this variant on RCS data. +- `DMLDiD(panel=False)` - Chang (2020) Case 2 declared-RCS mode + (row-unique unit IDs; level-outcome orthogonal scores with the + lambda-corrected variance; REGISTRY.md §DMLDiD). CAVEATS: it carries + NO survey weights even though RCS data is typically survey data + (weighted RCS -> CallawaySantAnna(panel=False, survey_design=...)), + and NO cluster= - the generic cluster-on-the-unit-proxy instruction + below CANNOT be followed on DMLDiD (its per-observation IF variance + is the only inference surface); aggregate('total') fails closed. - `TripleDifference` - DDD cross-sectional use cases are documented in `docs/choosing_estimator.rst`; the two-period DDD estimator does not require within-unit tracking when the third comparison axis diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index cffe01b30..5acaa51d7 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -1397,7 +1397,7 @@ Caveats and contracts: ### DMLDiD -Chang (2020, The Econometrics Journal 23(2)) double/debiased machine learning DiD, extended to staggered adoption: each Callaway-Sant'Anna style (g, t) cell is a 2-period Chang Case 1 problem with cross-fitted ML nuisances (DML2, per-cell D-stratified K-fold) and the Neyman-orthogonal score, so the ATT is first-order insensitive to nuisance regularization bias. Covariates are REQUIRED (conditional parallel trends per Abadie 2005; for unconditional staggered DiD use CallawaySantAnna). Per-cell SE = sqrt(mean(psi_bar^2)/n) from Chang's Theorem-2 augmented score; the analytical anchor is DoubleML (committed machine-precision parity spikes, doubleml==0.11.4). +Chang (2020, The Econometrics Journal 23(2)) double/debiased machine learning DiD, extended to staggered adoption: each Callaway-Sant'Anna style (g, t) cell is a 2-period Chang problem with cross-fitted ML nuisances (DML2) and a Neyman-orthogonal score, so the ATT is first-order insensitive to nuisance regularization bias. panel=True (default) runs Case 1 on panel data (outcome-change score, D-stratified folds); panel=False runs Case 2 on DECLARED repeated cross sections (level outcomes, one row per unit, D-x-period-stratified folds, the single control-only (T-lambda)*Y outcome nuisance, and the lambda-corrected Theorem 2 variance). Covariates are REQUIRED (conditional parallel trends per Abadie 2005; for unconditional staggered DiD use CallawaySantAnna). Per-cell SE = sqrt(mean(psi_bar^2)/n) from Chang's Theorem-2 augmented score; the panel anchor is DoubleML (committed machine-precision parity spikes, doubleml==0.11.4), while the RCS lane has NO parity oracle (DoubleMLDIDCSBinary uses a different score and omits the lambda term — a committed characterization spike documents the divergence). ```python DMLDiD( @@ -1413,15 +1413,16 @@ DMLDiD( base_period="varying", # or "universal" (CS semantics) cband=True, # sup-t bands on the event-study replay pscore_trim=0.01, # clip fitted propensities (never drop) + panel=True, # False = declared repeated cross sections (Chang Case 2) ).fit(data, outcome, unit, time, first_treat, covariates) -> DMLDiDResults ``` Key contracts: -- Aggregation is POST-FIT ONLY: `results.aggregate('event_study'/'group'/'simple'/'total')` (the fit-time `event_study_effects` surface is never populated by design; `plot_event_study(results.aggregate('event_study'))`). Bootstrapped fits replay the fit-time bootstrap on the recompute levels. +- Aggregation is POST-FIT ONLY: `results.aggregate('event_study'/'group'/'simple')`, plus `'total'` on panel fits — RCS fits fail `aggregate('total')` closed (the fit-time `event_study_effects` surface is never populated by design; `plot_event_study(results.aggregate('event_study'))`). Bootstrapped fits replay the fit-time bootstrap on the recompute levels. - HonestDiD / PreTrendsPower consume the aggregate('event_study') container (admitted source); the native `compute_honest_did(results)` route raises with the container instruction. Varying-base containers warn; use base_period='universal' for clean Rambachan-Roth interpretation. - `SieveLearner(k_max=None, criterion="bic")` is the exported configurable learner (adaptive polynomial degree by IC); any sklearn-style estimator object also plugs in (seed stochastic learners yourself — the library seed pins folds, not learner internals). - Per-cell complete cases (one consolidated unbalanced-input warning); degenerate cells (fewer members than folds, singleton treated stratum, fail-closed learner) become NaN cells with machine-readable skip_reason and a consolidated warning; surviving cells still aggregate. -- Panel only; no survey_design=/cluster= (Chang assumes i.i.d.; per-unit IF variance IS unit-level clustering). Reproducibility: set seed (fold draws move point estimates). +- No survey_design=/cluster= on EITHER design (Chang assumes i.i.d.; the per-sampling-unit IF variance IS unit-level clustering — per unit on panel, per observation on RCS). panel=False requires row-unique unit IDs and assumes stationary cross-sectional sampling (Assumption 2.3, warned, not data-checkable) — note RCS data is typically survey data (BRFSS/ACS/CPS) and DMLDiD carries NO survey weights; use CallawaySantAnna(panel=False, survey_design=...) for weighted RCS. Reproducibility: set seed (fold draws move point estimates). ### TROP diff --git a/diff_diff/guides/llms-practitioner.txt b/diff_diff/guides/llms-practitioner.txt index 7e1224062..34048fca6 100644 --- a/diff_diff/guides/llms-practitioner.txt +++ b/diff_diff/guides/llms-practitioner.txt @@ -227,7 +227,10 @@ Is treatment adoption staggered (multiple cohorts, different timing)? | | learners + orthogonal scores per (g,t) | | cell; covariates REQUIRED, aggregation | | post-fit, set seed= for reproducible -| | fold draws +| | fold draws; panel=False = declared +| | repeated cross sections (Case 2, +| | lambda-corrected variance; no survey +| | weights) | \-- WooldridgeDiD (ETWFE) -- nonlinear outcomes (logit/Poisson) or saturated OLS | |-- NO, simple 2x2 design: diff --git a/diff_diff/guides/llms.txt b/diff_diff/guides/llms.txt index c892964de..9443ddea6 100644 --- a/diff_diff/guides/llms.txt +++ b/diff_diff/guides/llms.txt @@ -81,7 +81,7 @@ The site is organized into 5 sections, each with a landing page: - [ChangesInChanges](https://diff-diff.readthedocs.io/en/stable/api/changes_in_changes.html): Athey & Imbens (2006) nonlinear/distributional DiD for the 2x2 design: recovers the treated group's full counterfactual outcome distribution and quantile treatment effects (ATT + QTE grid) via the CDF transformation `F_10(F_00^{-1}(F_01(y)))`; invariant to monotone outcome transformations (unconditional fits; the covariate QR branch is not); bootstrap inference (panel or repeated cross-section resampling); point parity with R `qte::CiC()`, including its covariate branch (`covariates=` -> per-cell linear quantile regression, Melly-Santangelo-style conditional CiC). Continuous outcomes, numeric covariates. Alias `CiC`. - [QDiD](https://diff-diff.readthedocs.io/en/stable/api/changes_in_changes.html): **Deprecated 3.9, removed 4.0 - use `ChangesInChanges(method="qdid")`.** Athey & Imbens (2006) quantile DiD comparison estimator (additive quantile-by-quantile DiD, matching R `qte::QDiD()` including its covariate branch via `covariates=`); same bootstrap machinery as ChangesInChanges. The paper recommends CiC over QDiD (scale-dependent model with testable restrictions; a non-monotonicity warning fires when violated - unconditional fits only, the covariate-path counterfactual quantile curve is monotone by construction). - [LWDiD](https://diff-diff.readthedocs.io/en/stable/api/lwdid.html): Lee & Wooldridge (2025, 2026) rolling-transformation DiD — unit-specific demean/detrend converts panel to cross-section, supports staggered adoption with flexible control groups. Signature: `LWDiD(rolling='demean', estimation_method='reg', vcov_type='hc1', cluster=None, control_group='not_yet_treated', alpha=0.05, n_bootstrap=0, seed=None, pscore_trim=0.01, n_neighbors=1, caliper=None, with_replacement=True, n_jobs=1).fit(data, outcome, unit, time, treatment, first_treat=None, covariates=None)`. `estimation_method` values: `reg` (papers' RA), `ipw`, `dr` (papers' IPWRA, doubly robust), `psm`; `vcov_type` values: `classical`/`hc1`/`hc2`/`hc3` for `reg`; `ipw`/`dr` accept `hc1` only (influence-function variance); `psm` accepts `hc1` as configuration only - PSM inference is unavailable (NaN) pending an Abadie-Imbens matching variance; cluster-robust inference via the constructor's `cluster=` column (hc1/CR1 only, not a `vcov_type` value; rejected for `psm`). Per-period effects: post-fit `results.aggregate('event_study')`. -- [DMLDiD](https://diff-diff.readthedocs.io/en/stable/api/dml_did.html): Chang (2020) double/debiased machine learning DiD — staggered ATT(g,t) with cross-fitted ML nuisances (DML2) and Neyman-orthogonal scores; covariates REQUIRED (conditional parallel trends). Signature: `DMLDiD(propensity_learner='logit', outcome_learner='linear', n_folds=5, control_group='never_treated', anticipation=0, alpha=0.05, n_bootstrap=0, bootstrap_weights=None, seed=None, base_period='varying', cband=True, pscore_trim=0.01).fit(data, outcome, unit, time, first_treat, covariates)`. Learners: string names (`linear`/`ridge`/`sieve` regressors, `logit` classifier) or any object with fit/predict(_proba) (sklearn-compatible); `SieveLearner(k_max, criterion)` is exported for adaptive polynomial nuisances. Aggregation is POST-FIT: `results.aggregate('event_study'/'group'/'simple'/'total')` (sup-t bands via bootstrap replay). With seed=None point estimates vary across fits (random folds); set seed for reproducibility. +- [DMLDiD](https://diff-diff.readthedocs.io/en/stable/api/dml_did.html): Chang (2020) double/debiased machine learning DiD — staggered ATT(g,t) with cross-fitted ML nuisances (DML2) and Neyman-orthogonal scores; covariates REQUIRED (conditional parallel trends). Signature: `DMLDiD(propensity_learner='logit', outcome_learner='linear', n_folds=5, control_group='never_treated', anticipation=0, alpha=0.05, n_bootstrap=0, bootstrap_weights=None, seed=None, base_period='varying', cband=True, pscore_trim=0.01, panel=True).fit(data, outcome, unit, time, first_treat, covariates)`. `panel=False` = declared repeated cross sections (Chang Case 2: level outcomes, row-unique unit IDs, lambda-corrected variance; no survey weights — weighted RCS belongs to CallawaySantAnna(panel=False, survey_design=...)). Learners: string names (`linear`/`ridge`/`sieve` regressors, `logit` classifier) or any object with fit/predict(_proba) (sklearn-compatible); `SieveLearner(k_max, criterion)` is exported for adaptive polynomial nuisances. Aggregation is POST-FIT: `results.aggregate('event_study'/'group'/'simple')`, plus `'total'` on panel fits (RCS fits fail 'total' closed); sup-t bands via bootstrap replay. With seed=None point estimates vary across fits (random folds); set seed for reproducibility. - [BaconDecomposition](https://diff-diff.readthedocs.io/en/stable/api/bacon.html): Goodman-Bacon (2021) decomposition for diagnosing TWFE bias in staggered settings ## Diagnostics and Sensitivity Analysis diff --git a/diff_diff/mmm.py b/diff_diff/mmm.py index 98e4609f3..38c39a59a 100644 --- a/diff_diff/mmm.py +++ b/diff_diff/mmm.py @@ -123,9 +123,12 @@ "contributing (g, t) cells (n_kind='cells'). Use " "results.aggregate('total') instead - its single row is the " "estimator-owned total incremental outcome over the per-cell " - "complete-case treated units and needs no scale; DMLDiD fits are " - "panel-only with no survey support, so the total route is available " - "on every fit (bootstrapped fits replay). If you need a different " + "complete-case treated units and needs no scale. " + "PANEL fits (the default) always support the total route (bootstrapped fits replay); repeated-cross-section fits " + "(panel=False) fail aggregate('total') closed on EVERY adopter " + "(CallawaySantAnna's RCS fits included) - estimator-owned totals need " + "per-unit tracking; on RCS fits pass a caller-defined numeric scale, " + "or use a panel fit when genuine longitudinal data exists. If you need a different " "estimand (e.g. full-cohort exposure on a panel with missing cells), " "pass a numeric scale." ), diff --git a/diff_diff/practitioner.py b/diff_diff/practitioner.py index ebeb65fac..c0ae286fa 100644 --- a/diff_diff/practitioner.py +++ b/diff_diff/practitioner.py @@ -540,7 +540,9 @@ def _handle_dml_did(results: Any): ), code=( "# Refit with alternative nuisance learners:\n" - "alt = DMLDiD(outcome_learner='sieve', seed=0).fit(\n" + "alt = DMLDiD(outcome_learner='sieve', seed=0" + + (", panel=False" if getattr(results, "panel", True) is False else "") + + ").fit(\n" " df, outcome=..., unit=..., time=..., first_treat=...,\n" " covariates=[...])\n" "print(alt.att, results.att) # should be close" diff --git a/docs/api/dml_did.rst b/docs/api/dml_did.rst index 8ae3e8873..ace59d885 100644 --- a/docs/api/dml_did.rst +++ b/docs/api/dml_did.rst @@ -4,19 +4,26 @@ DMLDiD — Double/Debiased Machine Learning DiD Chang (2020)'s double/debiased machine learning (DML) estimator for Difference-in-Differences with covariates, extended to staggered adoption: each Callaway-Sant'Anna style group-time cell :math:`(g, t)` is estimated -as a 2-period Chang Case 1 problem with cross-fitted machine-learning -nuisance functions and a Neyman-orthogonal score, so the ATT estimate is +as a 2-period Chang problem with cross-fitted machine-learning nuisance +functions and a Neyman-orthogonal score, so the ATT estimate is first-order insensitive to the regularization bias of the nuisance -learners. The classic 2-period design is the degenerate single-cell case. +learners. ``panel=True`` (the default) runs Case 1 (repeated outcomes) on +panel data; ``panel=False`` runs Case 2 (repeated cross sections) on +declared cross-sectional data — one observation per row. The classic +2-period design is the degenerate single-cell case of either lane. Identification follows Abadie (2005): CONDITIONAL parallel trends — the untreated potential-outcome trend is parallel across treated and control units only after conditioning on covariates :math:`X`. Chang's contribution is valid :math:`\sqrt{n}` inference when the two nuisance -functions (the propensity score :math:`g_0(X) = P(D=1|X)` and the control -outcome-change regression :math:`\ell_0(X) = E[\Delta Y | X, D=0]`) are +functions — the propensity score :math:`g_0(X) = P(D=1|X)` plus, on panel +fits, the control outcome-change regression +:math:`\ell_0(X) = E[\Delta Y | X, D=0]` (Case 1) or, on RCS fits, the +control level regression +:math:`\ell_{20}(X) = E[(T - \lambda) Y | X, D=0]` (Case 2) — are estimated by machine learning under DML2 cross-fitting — PROVIDED the -learners satisfy Chang's rate conditions (Assumption 3.1(f): each +learners satisfy Chang's rate conditions (Assumption 3.1(f) for Case 1 / +3.2(h) for Case 2: each nuisance at :math:`o_p(n^{-1/4})` in the :math:`L_2` norm, plus a product remainder bound). Cross-fitting removes overfitting bias but cannot substitute for the rate conditions, and the reported fold losses do not @@ -55,10 +62,34 @@ error is the plug-in :math:`\sqrt{\overline{\bar\psi^2} / n}`. This exact object was matched to DoubleML at machine precision in the committed parity spikes (``benchmarks/doubleml/``). +**Case 2 — repeated cross sections (Chang 2020, Equation 3.2;** +``panel=False``\ **).** The cell pools the two periods' rows (post indicator +:math:`T_i = 1\{\text{time}_i = t\}`) and scores LEVEL outcomes: + +.. math:: + + \text{summand}_i = \frac{D_i - \hat g(X_i)}{\hat p\,\hat\lambda(1-\hat\lambda)\,(1 - \hat g(X_i))}\,\bigl((T_i - \hat\lambda)Y_i - \hat\ell_2(X_i)\bigr) + +with :math:`\hat\lambda = \text{mean}(T)` the post-period sampling share and +:math:`\hat\ell_2` the SINGLE cross-fitted control-only regression of +:math:`(T - \hat\lambda)Y` on :math:`X` (Chang's :math:`I_{kz}^c`) — one +regression, deliberately different from the Sant'Anna-Zhao/DoubleML +four-regression RCS score. The Theorem 2 variance carries BOTH +finite-dimensional corrections: the treated-share fold-in plus an explicit +:math:`\hat G_{2\lambda}(T_i - \hat\lambda)` term (the λ-correction the +paper's proof structure warns is easy to omit), with +:math:`\hat G_{2\lambda}` the sample mean of the closed-form +:math:`\partial_\lambda \psi_2`. Folds are stratified on the four +:math:`D \times T` classes so every training fold carries control rows in +both periods. RCS aggregation weights are FIXED cohort row masses (the +CS-RCS convention, keeping the variance the influence function of the +reported aggregate); ``aggregate('total')`` is unavailable on RCS fits. + **Aggregation.** ``DMLDiD`` writes the CallawaySantAnna per-cell influence-function payload and inherits the CS aggregation and -multiplier-bootstrap machinery: event-study / group / simple / total -aggregations are produced **post-fit** via ``results.aggregate(...)``, +multiplier-bootstrap machinery: event-study / group / simple aggregations +(plus total on panel fits — RCS fits fail ``total`` closed) are produced +**post-fit** via ``results.aggregate(...)``, with sup-t uniform bands and bootstrap replay on bootstrapped fits. See ``docs/methodology/REGISTRY.md`` "DMLDiD" for the full equations, @@ -135,18 +166,33 @@ Restrictions - **Covariates required** — ``fit(covariates=None)`` or an empty list raises, directing to :class:`~diff_diff.CallawaySantAnna`. -- **Panel only** — repeated cross sections are not supported (Chang's - Case 2 score is a planned follow-up). +- **Declared designs only** — ``panel=True`` needs one row per + (unit, period); ``panel=False`` needs ROW-UNIQUE unit IDs (one + observation per row) and additionally assumes STATIONARY cross-sectional + sampling (Chang Assumption 2.3: each wave samples the same target + population — the composition of :math:`(D, X)` is stable across waves + while outcomes are the period-specific potential outcomes; warned at + fit, not data-checkable). + ``aggregate('total')`` is unavailable on RCS fits (fails closed). + NOTE: repeated-cross-section data is typically SURVEY data (BRFSS / ACS / + CPS), and ``DMLDiD`` carries NO survey weights — use + :class:`~diff_diff.CallawaySantAnna` ``(panel=False, survey_design=...)`` + for weighted RCS designs. - **No survey/cluster support** — Chang (2020) assumes i.i.d. sampling; - ``fit()`` accepts no ``survey_design=`` or ``cluster=``. The per-unit - influence-function variance IS unit-level clustering; coarser + ``fit()`` accepts no ``survey_design=`` or ``cluster=``. The + per-sampling-unit influence-function variance — per unit on panel fits, + per observation on repeated-cross-section fits — IS unit-level + clustering; coarser clustering is a tracked follow-up (``DEFERRED.md``). - **Propensity clipping, never dropping** — fitted propensities are clipped to ``[pscore_trim, 1 - pscore_trim]`` after an extremeness warning (the paper gives no trimming rule). -- **Per-cell complete cases** — a unit with a missing/non-finite outcome - or a non-finite base-period covariate at a cell is excluded from that - cell only (one consolidated ``UserWarning`` reports the drops). +- **Per-cell complete cases** — on panel fits, a unit with a + missing/non-finite outcome or a non-finite base-period covariate at a + cell is excluded from that cell only; on RCS fits the exclusion is + row-level (a non-finite outcome or covariate on the row — there is no + base-period covariate). One consolidated ``UserWarning`` reports the + drops. - **Degenerate cells skip loudly** — a cell that cannot be cross-fitted (fewer members than folds, a singleton treated/control stratum, a fail-closed learner error) is recorded as a NaN cell with a diff --git a/docs/choosing_estimator.rst b/docs/choosing_estimator.rst index b7a08808f..4bedb1391 100644 --- a/docs/choosing_estimator.rst +++ b/docs/choosing_estimator.rst @@ -656,32 +656,41 @@ identification strategy with analytical (non-bootstrap) inference. DMLDiD (Chang 2020) ~~~~~~~~~~~~~~~~~~~ -**When to use**: Staggered (or 2-period) panel designs where parallel trends -is plausible only CONDITIONAL on covariates and the covariate adjustment +**When to use**: Staggered (or 2-period) designs — panel data or declared +repeated cross sections (``panel=False``) — where parallel trends is +plausible only CONDITIONAL on covariates and the covariate adjustment must be flexible or high-dimensional — double/debiased machine learning (DML2 cross-fitting + Neyman-orthogonal scores) makes the ATT first-order insensitive to the nuisance learners' regularization bias. **Key features**: -- Per-(g, t) cell Chang (2020) Case 1 estimation on the Callaway-Sant'Anna - cell architecture; the 2-period design is the degenerate single-cell case +- Per-(g, t) cell Chang (2020) estimation on the Callaway-Sant'Anna cell + architecture — Case 1 (outcome changes) on panel data, Case 2 (level + outcomes with the λ-corrected variance) on declared repeated cross + sections; the 2-period design is the degenerate single-cell case - Configurable nuisance learners: string names (``'logit'``; ``'linear'``, ``'ridge'``, ``'sieve'``) or any object with ``fit``/``predict`` (``predict_proba``) — sklearn estimators plug in directly; :class:`~diff_diff.SieveLearner` is the exported adaptive-degree option - Covariates are REQUIRED (the no-covariates case routes to :class:`~diff_diff.CallawaySantAnna`) -- Post-fit aggregation (``results.aggregate('event_study'/'group'/'simple'/ - 'total')``) with sup-t bands and bootstrap replay; HonestDiD / +- Post-fit aggregation (``results.aggregate('event_study'/'group'/ + 'simple')``, plus ``'total'`` on panel fits — RCS fits fail ``total`` + closed) with sup-t bands and bootstrap replay; HonestDiD / PreTrendsPower consume the event-study container -- Analytical augmented-score inference anchored to DoubleML at machine - precision (committed parity spikes) +- Analytical augmented-score inference: the panel lane is anchored to + DoubleML at machine precision (committed parity spikes); the RCS lane + has no DoubleML oracle (different score) and is validated by equation + fixtures + a committed characterization spike **vs Callaway-Sant'Anna**: same cell architecture and aggregation surface; DMLDiD replaces CS's parametric nuisances with cross-fitted ML learners — prefer it when the covariate relationship is nonlinear/high-dimensional, -prefer CS otherwise (fewer moving parts, survey/cluster support). +prefer CS otherwise (fewer moving parts, survey/cluster support). Both +handle declared repeated cross sections via ``panel=False``, but only CS +carries survey weights there — RCS data is typically survey data, so +weighted RCS belongs to CS. **Example**:: diff --git a/docs/doc-deps.yaml b/docs/doc-deps.yaml index 882e091f2..ba5eb0179 100644 --- a/docs/doc-deps.yaml +++ b/docs/doc-deps.yaml @@ -979,6 +979,9 @@ sources: - path: docs/methodology/REGISTRY.md section: "DMLDiD" type: methodology + - path: docs/index.rst + section: "Supported Estimators (one-line catalog row)" + type: user_guide diff_diff/dml_did_results.py: drift_risk: low @@ -1452,7 +1455,7 @@ sources: - path: docs/methodology/REGISTRY.md section: "Cross-fitting, DR-score, and ridge infrastructure (DML)" type: methodology - note: "Shared DR panel scores: drdid_panel_inf_func (Sant'Anna-Zhao locally efficient IF, relocated verbatim from ContinuousDiD._dr_cell_inf_func, oracle-pinned in tests/test_dr_scores.py) and the Chang (2020) Case 1 orthogonal score pair (chang_panel_score / chang_panel_score_augmented) consumed by the DMLDiD estimator. The two families are NOT interchangeable (Chang normalizes by the unconditional treated share; SZ is self-normalized)." + note: "Shared DR panel scores: drdid_panel_inf_func (Sant'Anna-Zhao locally efficient IF, relocated verbatim from ContinuousDiD._dr_cell_inf_func, oracle-pinned in tests/test_dr_scores.py) and the Chang (2020) Case 1 orthogonal score pair (chang_panel_score / chang_panel_score_augmented) and Case 2 RCS family (chang_rcs_score / chang_rcs_lambda_slope / chang_rcs_score_augmented, lambda-corrected Theorem 2 variance) consumed by the DMLDiD estimator (panel=True/False). The two families are NOT interchangeable (Chang normalizes by the unconditional treated share; SZ is self-normalized)." - path: docs/methodology/papers/chang-2020-review.md type: methodology note: "The paper review's 'Relation to Existing diff-diff Estimators' section names drdid_panel_inf_func as the DRDID-parity score location; keep the cross-reference current if functions move." diff --git a/docs/index.rst b/docs/index.rst index 4e1c4b76c..dc6744f54 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -186,7 +186,7 @@ Supported Estimators * - :class:`~diff_diff.LWDiD` - Lee & Wooldridge (2025, 2026) rolling-transformation DiD; ``rolling='detrend'`` handles heterogeneous linear trends * - :class:`~diff_diff.DMLDiD` - - Chang (2020) double/debiased ML DiD; staggered ATT(g,t) with cross-fitted nuisance learners + - Chang (2020) double/debiased ML DiD; staggered ATT(g,t) with cross-fitted nuisance learners (panel or declared repeated cross sections) * - :class:`~diff_diff.QDiD` - Quantile DiD comparison estimator applying DiD quantile-by-quantile (deprecated 3.9 - use :class:`~diff_diff.ChangesInChanges` with ``method="qdid"``) * - :class:`~diff_diff.RegressionDiscontinuity` diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 1541ba57e..bb9483fca 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -859,8 +859,11 @@ double-weight). `doubleml.DoubleMLDIDBinary` via `benchmarks/doubleml/chang_staggered_parity.py` (same pin); `DoubleMLDIDMulti` orchestrates the same per-cell estimator over staggered - timing. Not oracles: `DoubleMLDIDCS` (different RCS score), - `DoubleMLDIDMulti` for Chang Cases 2-3. + timing. Not oracles: `DoubleMLDIDCS`/`DoubleMLDIDCSBinary` (different RCS + score, no λ-correction — the committed characterization spike + `benchmarks/doubleml/chang_rcs_characterization.py` documents the + divergence for the shipped Case 2 lane), `DoubleMLDIDMulti` for Chang + Case 3. - R: none for Chang's estimator; `DRDID::drdid_panel` for the SZ score. --- @@ -2791,17 +2794,34 @@ control-variable dimension `d` may exceed the sample size `N`, by constructing **Neyman-orthogonal scores** (each = Abadie's score + a mean-zero adjustment term) and estimating with the **Chernozhukov et al. (2018) DML cross-fitting algorithm (DML2 variant)**. The shipped `DMLDiD` estimator implements **Case 1 -(repeated outcomes / panel)** as a STAGGERED ATT(g,t) estimator: each -Callaway-Sant'Anna style `(g, t)` cell is a 2-period Chang problem on the -cell's `ΔY = Y_t − Y_base`, treated indicator `D = 1{cohort = g}`, and -base-period covariates; the classic 2-period design is the degenerate -single-cell case. Case 2 (repeated cross sections; Equation 3.2 score + -λ-corrected variance) is a planned follow-up; Case 3 (multilevel treatment -intensity) is deferred (`DEFERRED.md`). +(repeated outcomes / panel; `panel=True`, the default)** and **Case 2 +(repeated cross sections; `panel=False`)** as a STAGGERED ATT(g,t) estimator: +each Callaway-Sant'Anna style `(g, t)` cell is a 2-period Chang problem. On +the panel lane the cell score runs on `ΔY = Y_t − Y_base`, treated indicator +`D = 1{cohort = g}`, and base-period covariates. On the declared-RCS lane +(one observation per row, row-unique unit IDs) the cell pools the two +periods' rows and runs Equation 3.2 on LEVEL outcomes with the post-period +sampling share `λ̂` and the λ-corrected Theorem 2 variance; the classic +2-period design is the degenerate single-cell case of either lane. Case 3 +(multilevel treatment intensity) is deferred (`DEFERRED.md`). + +**Case 2 equations (as implemented).** Per cell with pooled two-period rows, +post indicator `T = 1{time = t}`, `p̂ = mean(D)` and `λ̂ = mean(T)` (global +within cell):: + + summand_i = (D_i − ĝ(X_i)) / (p̂ λ̂ (1−λ̂)(1−ĝ(X_i))) · ((T_i − λ̂) Y_i − ℓ̂₂(X_i)) + θ̂ = mean(summand) (Equation 3.2: ψ₂ = summand − θ) + ψ̄_i = summand_i − D_i θ̂ / p̂ + Ĝ₂λ (T_i − λ̂) (Theorem 2, BOTH corrections) + SE = sqrt(mean(ψ̄²) / n_cell) + +where `ℓ̂₂` is the SINGLE cross-fitted control-only regression of +`(T − λ̂)·Y` on X (Chang's `I_kz^c`) and `Ĝ₂λ` is the sample mean of the +closed-form `∂λψ₂` (`chang_rcs_lambda_slope`). - **Note:** Staggered extension framing. Chang (2020) is a 2-period, common-timing paper (treatment occurs only at the second period). - `DMLDiD` applies the Case 1 score PER (g, t) CELL under the + `DMLDiD` applies the corresponding Case 1 (panel) or Case 2 (RCS) score + PER (g, t) CELL under the Callaway-Sant'Anna cell architecture (positional base periods, R `did` parity; never-treated / not-yet-treated control semantics; per-cell complete cases), which is a library extension of the paper's design, not a @@ -2812,8 +2832,9 @@ intensity) is deferred (`DEFERRED.md`). *Assumption checks / warnings:* - **Conditional parallel trends (Assumption 2.1, Abadie 2005):** `E[Y^0(1) - Y^0(0) | X, D = 1] = E[Y^0(1) - Y^0(0) | X, D = 0]` per cell. Untestable; the estimator adds **no identification assumptions beyond Abadie (2005)** (Section 2, p. 8). -- **Overlap (Assumption 2.2):** `P(D = 1) > 0` and `P(D = 1 | X) < 1` a.s. Regularity Assumptions 3.1(a) strengthen this to **strict overlap**: `Pr(κ ≤ g_0(X) ≤ 1 - κ) = 1` for some fixed `κ > 0`, imposed on the **estimated** propensity too (pp. 28, 37). The theory does not cover fitted propensities approaching 0 or 1 — see the trimming Note below. -- **First-stage rate condition (Assumptions 3.1(f) + Theorem 1):** BOTH conditions hold jointly — the bundle envelope `‖η̂_k - η_0‖_{P,2} ≤ ε_N` with `ε_N = o(N^{-1/4})` (each nuisance component must meet the rate; a fast learner cannot compensate a slow one), AND the product bound `‖ĝ - g_0‖²_{P,2} + ‖ĝ - g_0‖_{P,2}·‖ℓ̂ - ℓ_0‖_{P,2} ≤ ε_N²`. Not detectable at runtime; documented. +- **Overlap (Assumption 2.2):** `P(D = 1) > 0` and `P(D = 1 | X) < 1` a.s. Regularity Assumptions 3.1(a) (Case 1) / 3.2(a) (Case 2) strengthen this to **strict overlap**: `Pr(κ ≤ g_0(X) ≤ 1 - κ) = 1` for some fixed `κ > 0`, imposed on the **estimated** propensity too (pp. 28, 37). The theory does not cover fitted propensities approaching 0 or 1 — see the trimming Note below. +- **Stationary RCS sampling (Assumption 2.3, `panel=False` only):** conditional on `T = 0` (resp. `T = 1`), rows are i.i.d. draws from the distribution of `(Y(0), D, X)` (resp. `(Y(1), D, X)`) — each wave samples the same target population, so the composition of `(D, X)` is stable across waves while outcomes are the period-specific potential outcomes (trends and treatment effects are expected, not violations). Not data-checkable; surfaced as a fit-time `UserWarning` after the declared-RCS structure validates. +- **First-stage rate condition (Assumptions 3.1(f) (Case 1) / 3.2(h) (Case 2) + Theorem 1):** BOTH conditions hold jointly — the bundle envelope `‖η̂_k - η_0‖_{P,2} ≤ ε_N` with `ε_N = o(N^{-1/4})` (each nuisance component must meet the rate; a fast learner cannot compensate a slow one), AND the product bound `‖ĝ - g_0‖²_{P,2} + ‖ĝ - g_0‖_{P,2}·‖ℓ̂ - ℓ_0‖_{P,2} ≤ ε_N²`. Not detectable at runtime; documented. *Estimator equation (Equation 3.1, per cell, as implemented — `chang_panel_score`):* @@ -2941,6 +2962,51 @@ the finite-dimensional `p_0` is handled by the variance correction below. COARSER-than-unit clustering, a surface DMLDiD ships without: the paper assumes i.i.d. sampling, and the survey/cluster follow-up is the tracked `DEFERRED.md` row. +- **Note:** Global λ̂ convention (Case 2) — `λ̂` is the FULL-SAMPLE-within-cell + post-period sampling share `mean(T)` over the pooled two-period rows, + mirroring the global-p̂ convention above and DoubleML's `t.mean()` (the + paper's per-fold `λ̂_k` printing carries the same contradiction as `p̂_k` — + see Gaps in the paper review). +- **Note:** Ĝ₂λ estimator (Case 2) — the paper prints NO explicit estimator + for the λ-slope `G_2λ0`; `chang_rcs_lambda_slope` uses the natural sample + analogue — the mean of the closed-form `∂λψ₂` (recovered from the Theorem 2 + proof, p. 55 display) at the plug-in nuisances. Theorem 2 requires only + consistency of this estimator, no rate. +- **Note:** λ-correction is MANDATORY in the Case 2 variance — the review + warns verbatim that "omitting the λ-correction term is a plausible + implementation bug the proof structure warns against"; + `TestLambdaCorrectionRegression` pins at the estimator level that the + reported SE differs from the λ-term-omitted recomputation. +- **Note:** D×T fold stratification (Case 2) — RCS cells stratify the fold + draw on the FOUR D×period classes (`D + 2T`, DoubleML's `d + 2t` encoding), + a deviation from Chang's plain random partition that structurally + guarantees control rows in BOTH periods in every training complement + (Chang's `I_kz^c` fold-composition requirement). Consequence: any singleton + D×T stratum (e.g. ONE treated row in the base period) dies as + `cross_fit_degenerate`. +- **Note:** Case 2 four-group guard — a cell needs treated AND control rows + in BOTH periods; any empty group skips as `zero_treated_control` (the + vocabulary is reused, not widened). A new λ̂-extremeness warning mirrors + the empirical-p̂ one (Case 2 bounds carry `1/(λ(1−λ))` up to cubes). +- **Note:** RCS aggregation weights (Case 2) — event-study/group/simple + aggregation weights each estimated cell by its FIXED cohort row mass + (per-cell `agg_weight`, the CS-RCS convention), so the aggregation WIF + variance is the influence function of the reported aggregate (the pg basis + is the per-row cohort bincount, identical under unique row IDs); + `n_treated`/`n_control` on RCS cells are pooled two-period valid-row + DISPLAY counts only, never weights. `agg_cohort_masses` is deliberately + NOT set in the precompute: it would be a numeric no-op that routes lookups + through `float()`-keyed dict access, where distinct int64 cohorts above + 2^53 (admissible through 2^62 by the label pipeline) collide. +- **Note:** `aggregate('total')` fails closed on RCS fits (the library-wide + repeated-cross-section convention; `is_panel: False` kit bookkeeping); + panel fits keep the total route. Per-observation IF entries make + `vcov_type="hc1"` the per-SAMPLING-UNIT variance on RCS (rows are the + sampling units). +- **Note:** Case-2-only moment conditions — Assumption 3.2's level-outcome + bounds (`E[Y²|X] ≤ C`, `|E[YU]| ≤ C`) belong to the repeated-cross-section + case ONLY and are NOT imposed on the panel path; stationary sampling + (Assumption 2.3) is warned at fit (not data-checkable). *Edge cases:* - Propensity near 0/1: clipped per the trimming Note (error bounds blow up as @@ -2966,9 +3032,17 @@ the finite-dimensional `p_0` is handled by the variance correction below. every cell's ATT and SE at <1e-15 under shared folds and pinned config). `DoubleMLDIDMulti` orchestrates the same per-cell estimator over staggered timing (and remains NOT an oracle for Chang Cases 2-3). +- Case 2 (RCS): NO parity oracle exists — `DoubleMLDIDCSBinary` implements + the Sant'Anna-Zhao four-regression RCS score and its variance omits + Chang's λ-correction. The committed CHARACTERIZATION spike + `benchmarks/doubleml/chang_rcs_characterization.py` (same doubleml==0.11.4 + pin) documents the per-cell divergence under shared folds, isolates the + λ-term's SE effect, and anchors the shipped estimator by SELF-parity + (public `DMLDiD(panel=False)` == the hand-rolled Eq 3.2 / Thm 2 pipeline + at 0.0 observed diff under identical folds and sklearn learners). - R / Stata: none — the paper ships no companion package. -**Requirements checklist (shipped Case 1 / staggered lane):** +**Requirements checklist (shipped Case 1 panel + Case 2 RCS staggered lanes):** - [x] Neyman-orthogonal Case 1 score (3.1) implemented exactly (Abadie score + mean-zero adjustment; `chang_panel_score`) - [x] DML2 cross-fitting: per-cell K-fold partition, nuisances fit on fold complements, never on the evaluation fold - [x] Outcome nuisance `ℓ̂` fit on the UNTREATED subsample of the training complement only (`I_kz^c` → `fit_mask=(D==0)`) @@ -2978,7 +3052,7 @@ the finite-dimensional `p_0` is handled by the variance correction below. - [x] Per-cell degenerate guards (zero treated/control, cell < K, singleton stratum, empty untreated complement) — closed skip vocabulary, consolidated warning - [x] Normal-approximation inference via `safe_inference()` - [x] Validation: 2-period DoubleMLDID + staggered per-cell DoubleMLDIDBinary parity spikes (version-pinned, committed, golden literals consumed in-tests) + oracle-nuisance closed-form equivalence + degenerate-cell hand-pipeline equivalence (rtol 1e-14 — BLAS reduces differently-laid-out inputs in platform-dependent order, so bit identity does not hold cross-platform) + Monte Carlo coverage sanity -- [ ] Case 2 (repeated cross sections): Equation 3.2 score + λ-corrected variance — planned follow-up (ROADMAP) +- [x] Case 2 (repeated cross sections): Equation 3.2 score + λ-corrected Theorem 2 variance — SHIPPED as `DMLDiD(panel=False)` (`chang_rcs_score` / `chang_rcs_lambda_slope` / `chang_rcs_score_augmented`; equation-level fixtures, oracle closed forms, derivative-identity checks, DR both directions, characterization spike, MC coverage). The paper's own §4 RCS simulation DGPs are NOT replicated (tracked TODO row — needs the paper PDF pp. 17-21); the shipped recovery tests use a library-authored RCS design. - [ ] Case 3 (multilevel treatment): deferred (`DEFERRED.md`; implementation-required overlap conditions per the paper review's Case 3 caution) --- diff --git a/docs/methodology/REPORTING.md b/docs/methodology/REPORTING.md index 0b4a86577..223ef9084 100644 --- a/docs/methodology/REPORTING.md +++ b/docs/methodology/REPORTING.md @@ -278,6 +278,18 @@ a library setting. weighted estimate; users can alternatively pass `precomputed={'bacon': ...}` with a survey-aware result. +- **Note:** Declared repeated-cross-section fits skip Bacon. Any + panel-attribute-bearing producer fitted with `panel=False` (DMLDiD + and CallawaySantAnna alike) skips the Goodman-Bacon check with an + explicit reason instead of replaying: `BaconDecomposition` on + one-observation-per-unit data is degenerate (treatment collinear + with the absorbed unit fixed effects — it returns a meaningless + `twfe_estimate` of 0.0 rather than raising). The gate sits AFTER the + precomputed-Bacon passthrough, so `precomputed={'bacon': ...}` + computed on a real panel stays honored. BusinessReport's DMLDiD + identification narrative is likewise design-aware: RCS fits cite the + stationary-sampling Assumption 2.3 and the Case 2 rate package. + The simple 2x2 parallel-trends helper (`utils.check_parallel_trends`) has no survey-aware variant. On a survey-backed `DiDResults` the check is skipped **unconditionally**, regardless of whether diff --git a/docs/methodology/papers/chang-2020-review.md b/docs/methodology/papers/chang-2020-review.md index d25d14254..b53b81976 100644 --- a/docs/methodology/papers/chang-2020-review.md +++ b/docs/methodology/papers/chang-2020-review.md @@ -54,7 +54,7 @@ The conventional plug-in of (2.1) fails with ML first stages (p. 8): the score h Gateaux derivative in `g_0`, and ML nuisances converge slower than `N^{-1/2}` due to regularization bias. -*Estimator equation — orthogonal scores (Equations 3.1-3.3 in paper; Case 1 SHIPPED as the staggered `DMLDiD` estimator — see REGISTRY.md "DMLDiD"; Cases 2-3 remain unimplemented):* +*Estimator equation — orthogonal scores (Equations 3.1-3.3 in paper; Cases 1-2 SHIPPED as the staggered `DMLDiD` estimator, `panel=True`/`False` — see REGISTRY.md "DMLDiD"; Case 3 remains unimplemented):* Each score = Abadie's score + a mean-zero adjustment term (`c_1`, `c_2`, `c_w`), so the same ATT is identified. @@ -177,22 +177,22 @@ the same variance estimators remain consistent under kernel first stages. - Python (validation oracle — **panel lane only**): `doubleml.DoubleMLDID` implements a Chang/Zimmert-style orthogonal panel score and is the closest parity anchor for Case 1 — but only under a specific configuration (its `in_sample_normalization` option changes the score's normalization; pin the config and verify score equivalence at the equation level before treating any run as an oracle). **Version caveat:** `DoubleMLDID` is deprecated upstream ("will be removed with version 0.12.0. Please use DoubleMLDIDBinary instead", verified 2026-08-22) — pin the exact DoubleML version used for golden-fixture generation, archive the fixtures in-repo, and verify `DoubleMLDIDBinary`'s score/normalization equivalence separately before adopting it as the replacement anchor. **Scope caveat:** treat DoubleML as an **equation-level score oracle**, not a full-estimator finite-sample oracle, unless the fixture supplies identical fold assignments and scalar-nuisance normalization: upstream uses the global treated share and treatment-stratified sample splitting, while Chang specifies random folds and leaves the fold-level `p̂_k` convention ambiguous (see Gaps). Full-estimator comparisons under differing conventions are approximate/asymptotic — a finite-sample mismatch there does not falsify a paper-faithful implementation, and exact-parity fixtures must not silently import DoubleML's normalization as if it were Chang's. Add independent fixtures for whichever `p̂_k` convention is selected. - **Not oracles:** `doubleml.DoubleMLDIDCS` uses a Sant'Anna-Zhao-style repeated-cross-section score with four treatment-by-period outcome regressions and different normalizations — a *related* estimator, not an implementation of Chang's single-`ℓ_20` Equation 3.2 score or its `λ`-corrected variance. `DoubleMLDIDMulti` handles **staggered treatment timing**, not Chang's Case 3 multilevel treatment *intensity* — it is not a Case 3 anchor. For Cases 2-3, validation must rest on independent equation-level fixtures (Equations 3.2/3.3 and their variance corrections) plus recovery tests on the paper's simulation DGPs (Section 4). -**Requirements checklist** (Case 1 items SHIPPED as `DMLDiD` — the -per-item conventions/deviations are the REGISTRY "DMLDiD" Notes; Case 2-3 -items remain open, tracked in ROADMAP/DEFERRED): -- [x] Neyman-orthogonal Case 1 score (3.1) implemented exactly (Abadie score + mean-zero adjustment); [ ] scores (3.2)/(3.3) — Cases 2-3 open +**Requirements checklist** (Case 1-2 items SHIPPED as `DMLDiD` — the +per-item conventions/deviations are the REGISTRY "DMLDiD" Notes; Case 3 +items remain open, tracked in DEFERRED.md): +- [x] Neyman-orthogonal Case 1 score (3.1) implemented exactly (Abadie score + mean-zero adjustment); [x] score (3.2) — SHIPPED as `chang_rcs_score` (`DMLDiD(panel=False)`); [ ] score (3.3) — Case 3 open - [x] DML2 cross-fitting: K-fold partition (D-stratified — documented deviation), nuisances fit on fold complements, never on the evaluation fold - [x] Outcome nuisance `ℓ̂` fit on the UNTREATED subsample of the auxiliary fold only (`I_kz^c`) -- [x] Scalar nuisance p̂: the global (full-sample-within-cell) convention adopted and documented (the I_k vs I_k^c printing contradiction is thereby sidestepped — see Gaps); [ ] `λ̂_k` — Case 2 open +- [x] Scalar nuisance p̂: the global (full-sample-within-cell) convention adopted and documented (the I_k vs I_k^c printing contradiction is thereby sidestepped — see Gaps); [x] `λ̂` — same global convention (`mean(T)` within cell; REGISTRY Note) - [x] Final estimator: pooled mean (equals the paper's `1/K` average at equal fold sizes — documented deviation) -- [x] Variance from the AUGMENTED score: `Ĝ_1p = -θ̃/p̂` folded in; [ ] `Ĝ_2λ (T - λ̂_k)` — Case 2 open +- [x] Variance from the AUGMENTED score: `Ĝ_1p = -θ̃/p̂` folded in; [x] `Ĝ_2λ (T - λ̂)` — explicit term in `chang_rcs_score_augmented` (`Ĝ_2λ` = sample mean of the closed-form `∂λψ₂`, `chang_rcs_lambda_slope`; the paper prints no estimator — REGISTRY Note) - [x] Strict-overlap enforcement: fitted propensities clipped to `[trim, 1-trim]` (documented deviation; paper gives no rule) -- [x] Per-fold/per-cell degenerate guards (zero treated/control, cell < K, singleton stratum) — closed skip vocabulary; [ ] Case 2 pre/post-share guards — open +- [x] Per-fold/per-cell degenerate guards (zero treated/control, cell < K, singleton stratum) — closed skip vocabulary; [x] Case 2 pre/post-share guards — four-group guard + `λ̂` extremeness warning + D×T-stratified folds (control rows in both periods per training complement by construction) - [x] Auxiliary-sample feasibility: an empty untreated training complement raises `DegenerateFoldError` (targeted, before the learner) → `cross_fit_degenerate` cell - [x] Normal-approximation inference via `safe_inference()` - [ ] Multilevel treatment (Case 3): open — DEFERRED - [ ] Case 3 guards: open — DEFERRED (see the Case 3 caution above) -- [x] Validation (Case 1): `doubleml.DoubleMLDID` (2-period) + `DoubleMLDIDBinary` (staggered per-cell, end-to-end public fit) parity spikes, doubleml==0.11.4 pinned, golden literals in-repo; [ ] Cases 2-3 equation-level fixtures — open +- [x] Validation (Case 1): `doubleml.DoubleMLDID` (2-period) + `DoubleMLDIDBinary` (staggered per-cell, end-to-end public fit) parity spikes, doubleml==0.11.4 pinned, golden literals in-repo; [x] Case 2 equation-level fixtures — SHIPPED (closed-form/oracle fixtures, derivative-identity checks, DR both directions, `DoubleMLDIDCSBinary` characterization spike — no parity oracle exists). CAVEAT: the paper's own §4 RCS simulation DGPs (pp. 17-21) are NOT replicated — the shipped recovery/coverage tests use a library-authored RCS design; replication is a tracked TODO.md row (needs the paper PDF); [ ] Case 3 fixtures — open --- @@ -200,7 +200,7 @@ items remain open, tracked in ROADMAP/DEFERRED): ### Data Structure Requirements - **Repeated outcomes (panel, primary target — SHIPPED as `DMLDiD`):** one row per unit with `(Y(0), Y(1), D, X)` — equivalently a 2-period panel converted to differences `ΔY = Y(1) - Y(0)`. The paper's design is 2-period with treatment only at `t = 1`; the shipped estimator applies it per Callaway-Sant'Anna (g, t) cell over staggered timing (a documented library extension — REGISTRY.md "DMLDiD"). -- **Repeated cross sections:** pooled rows `(Y, D, T, X)` with `T ∈ {0,1}` the post-period sampling indicator; i.i.d. within each period (Assumption 2.3). +- **Repeated cross sections (SHIPPED as `DMLDiD(panel=False)`):** pooled rows `(Y, D, T, X)` with `T ∈ {0,1}` the post-period sampling indicator; i.i.d. within each period (Assumption 2.3). The shipped lane applies Equation 3.2 per Callaway-Sant'Anna (g, t) cell over staggered timing with row-unique unit IDs (the same documented library extension as the panel lane). - **Multilevel treatment:** `(Y(0), Y(1), W, X)` with `W ∈ {0, w_1, ..., w_J}`; `W = 0` is the comparison group for every level. - Covariates `X ∈ R^d` may have `d > N` (the headline setting); no factor-structure or panel-index requirements beyond the above. diff --git a/docs/methodology/variance-conventions.md b/docs/methodology/variance-conventions.md index e77a165a8..e9b017e00 100644 --- a/docs/methodology/variance-conventions.md +++ b/docs/methodology/variance-conventions.md @@ -39,7 +39,7 @@ produced wrong figures three separate times while this inventory was drafted). | `two_stage_default` | — | None | **legitimate** | L3: Gardner two-stage variance, not the shared CR1 sandwich | | `callaway_santanna_default` | — | None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None | **legitimate** | L3: influence-function variance anchored to Stata csdid | | `dml_did` | — | None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None | **legitimate** | L3: Chang (2020) Thm 2 augmented-score plug-in variance (per-unit influence function; normal-theory safe_inference throughout, no cluster surface) | - +| `dml_did_rcs` | — | None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None | **legitimate** | L3: Chang (2020) Thm 2 lambda-corrected augmented-score plug-in variance (Case 2; per-observation influence function; normal-theory safe_inference throughout, no cluster surface) | cr1_k is the sorted multiset of K_reference counts reaching the shared clustered CR1 denominator — visible columns + the signed cluster_k_adjustment (linalg._compute_robust_vcov_numpy with @@ -132,7 +132,10 @@ output). - **L3 — CallawaySantAnna / DMLDiD / TwoStageDiD / ImputationDiD (default)**: different variance theory (influence functions / Chang (2020) augmented-score plug-in / two-stage / BJS imputation), never the shared CR1 sandwich. CS is - anchored to Stata csdid outright; DMLDiD to DoubleML. + anchored to Stata csdid outright; DMLDiD's panel lane to DoubleML at + machine precision, while its repeated-cross-section lane (panel=False) is + characterization-anchored only (DoubleML's RCS score differs and omits the + Chang lambda-correction; equation-level fixtures validate the variance). **ImputationDiD is conditional**: its pretrends=True lead regression DOES run the shared clustered CR1 — surfaced fit-time via the deprecated aggregate="event_study" or post-fit via results.aggregate('event_study') diff --git a/docs/migration-4.0.md b/docs/migration-4.0.md index 859bf13be..e7c1e5148 100644 --- a/docs/migration-4.0.md +++ b/docs/migration-4.0.md @@ -164,7 +164,7 @@ backend fail closed with a refit message. On `ImputationDiD`, `TwoStageDiD` and `ContinuousDiD`, the recompute levels still raise `NotImplementedError` when the fit used `n_bootstrap > 0` — keep the fit-time call there for now and track their open `TODO.md` rows. -`aggregate("simple")` — and, on its five adopters, `aggregate("total")` (3.10; DMLDiD joined in 3.11) — does relay, +`aggregate("simple")` — and, on its five adopters, `aggregate("total")` (3.10; DMLDiD joined in 3.11 — panel fits only, RCS fits fail `total` closed) — does relay, and `StackedDiD`, `ChaisemartinDHaultfoeuille` and `HeterogeneousAdoptionDiD` are unaffected — their `aggregate()` is a pure view over stored fields. ``` diff --git a/docs/practitioner_decision_tree.rst b/docs/practitioner_decision_tree.rst index 3c80296c8..ed9b80c44 100644 --- a/docs/practitioner_decision_tree.rst +++ b/docs/practitioner_decision_tree.rst @@ -571,8 +571,11 @@ The six scenarios above cover the most common business use cases. ATT is first-order insensitive to ML regularization bias. Covariates are REQUIRED (conditional parallel trends); aggregation is post-fit (``results.aggregate('event_study')`` feeds HonestDiD/PreTrendsPower); - set ``seed=`` for reproducible fold draws. Panel only; no - survey/cluster support (per-unit influence-function variance). + set ``seed=`` for reproducible fold draws. Panel data by default; + ``panel=False`` runs declared repeated cross sections (Chang Case 2, + λ-corrected variance). No survey/cluster support on either lane + (per-sampling-unit influence-function variance) — weighted RCS belongs + to ``CallawaySantAnna(panel=False, survey_design=...)``. For the full academic decision tree with all estimators, see :doc:`choosing_estimator`. diff --git a/docs/references.rst b/docs/references.rst index 5d284b6ee..f57aada20 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -297,7 +297,7 @@ Double/Debiased Machine Learning - **Chang, N.-C. (2020).** "Double/Debiased Machine Learning for Difference-in-Differences Models." *The Econometrics Journal*, 23(2), 177-191. https://doi.org/10.1093/ectj/utaa001 - Neyman-orthogonal DiD scores with DML2 cross-fitting for ML first stages. The Case 1 (repeated outcomes) orthogonal score and its augmented-variance companion ship as ``chang_panel_score`` / ``chang_panel_score_augmented`` in ``diff_diff/_dr_scores.py`` (consumed by the shipped ``DMLDiD`` estimator); paper review on file at ``docs/methodology/papers/chang-2020-review.md``, and the DoubleML parity anchors are committed at ``benchmarks/doubleml/chang_case1_parity.py`` (2-period) and ``benchmarks/doubleml/chang_staggered_parity.py`` (staggered per-cell). + Neyman-orthogonal DiD scores with DML2 cross-fitting for ML first stages. The Case 1 (repeated outcomes) score pair ``chang_panel_score`` / ``chang_panel_score_augmented`` and the Case 2 (repeated cross sections) family ``chang_rcs_score`` / ``chang_rcs_lambda_slope`` / ``chang_rcs_score_augmented`` ship in ``diff_diff/_dr_scores.py`` (consumed by the shipped ``DMLDiD`` estimator, ``panel=True``/``False``); paper review on file at ``docs/methodology/papers/chang-2020-review.md``. The DoubleML parity anchors are committed at ``benchmarks/doubleml/chang_case1_parity.py`` (2-period) and ``benchmarks/doubleml/chang_staggered_parity.py`` (staggered per-cell); the Case 2 lane has no parity oracle — ``benchmarks/doubleml/chang_rcs_characterization.py`` documents the divergence from DoubleML's Sant'Anna-Zhao-style RCS score. - **Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., & Robins, J. (2018).** "Double/Debiased Machine Learning for Treatment and Structural Parameters." *The Econometrics Journal*, 21(1), C1-C68. https://doi.org/10.1111/ectj.12097 diff --git a/docs/survey-roadmap.md b/docs/survey-roadmap.md index b69e13e76..3b6f67df1 100644 --- a/docs/survey-roadmap.md +++ b/docs/survey-roadmap.md @@ -271,7 +271,7 @@ the limitation and suggested alternative. | Estimator | Limitation | Alternative | |-----------|-----------|-------------| | LWDiD | Any `survey_design` / sampling weights | No weight argument exists on any path, so the failure mode is a bare `TypeError: unexpected keyword argument` rather than a descriptive error (the exception to the preamble above). The LW papers derive the transformation and exact-inference layer for unweighted panels; a weighted counterpart is DEFERRED pending user demand. Use `CallawaySantAnna` (or another survey-capable staggered estimator) when design-based variance is required. | -| DMLDiD | Any `survey_design` / sampling weights / `cluster=` | No weight or cluster argument exists on any path (bare `TypeError` on the kwarg). Chang (2020) assumes i.i.d. sampling; the per-unit influence-function variance IS unit-level clustering, and the survey/cluster extension is the tracked DEFERRED.md row (pending user demand). Use `CallawaySantAnna` when design-based variance or coarser clustering is required. | +| DMLDiD | Any `survey_design` / sampling weights / `cluster=` | No weight or cluster argument exists on any path (bare `TypeError` on the kwarg) — on the panel lane AND the declared-RCS lane (`panel=False`). Chang (2020) assumes i.i.d. sampling; the per-sampling-unit influence-function variance IS unit-level clustering, and the survey/cluster extension is the tracked DEFERRED.md row (pending user demand). NOTE the RCS lane sharpens this gap: repeated-cross-section data is typically survey data (BRFSS/ACS/CPS) — use `CallawaySantAnna(panel=False, survey_design=...)` for weighted RCS, and `CallawaySantAnna` generally when design-based variance or coarser clustering is required. | | SyntheticDiD | Replicate weights | Pre-existing limitation: no replicate-weight survey support on SDID. All three variance methods (bootstrap, placebo, jackknife) now support pweight-only and strata/PSU/FPC designs; replicate-weight designs remain rejected. | | TROP | Replicate weights | Use strata/PSU/FPC design with Rao-Wu rescaled bootstrap | | BaconDecomposition | Replicate weights | Diagnostic only, no inference | diff --git a/tests/test_dml_did.py b/tests/test_dml_did.py index 4c8e1feaf..a7fc4a93c 100644 --- a/tests/test_dml_did.py +++ b/tests/test_dml_did.py @@ -1628,3 +1628,516 @@ def test_bootstrap_summary_labels_percentile_p(self, data): s = boot.summary() assert "Boot. p" in s and "P>|z|" not in s assert "PERCENTILE" in s + + +# =========================================================================== +# Repeated cross sections (Chang Case 2, panel=False) +# =========================================================================== + + +def make_rcs_dml_data(n_rows=1500, periods=(1, 2, 3), cohorts=(0, 2, 3), seed=7, effect=2.0): + """Declared-RCS frame: one observation per row, row-unique unit IDs.""" + rng = np.random.default_rng(seed) + cohort_arr = rng.choice( + cohorts, size=n_rows, p=[0.5] + [0.5 / (len(cohorts) - 1)] * (len(cohorts) - 1) + ) + tt = rng.choice(periods, size=n_rows) + x1 = rng.normal(size=n_rows) + x2 = rng.normal(size=n_rows) + y = 1.0 + 0.5 * x1 - 0.3 * x2 + 0.2 * tt + rng.normal(scale=0.5, size=n_rows) + post = (cohort_arr > 0) & (tt >= cohort_arr) + y = y + effect * post + return pd.DataFrame( + { + "unit": np.arange(n_rows), + "time": tt, + "first_treat": cohort_arr, + "y": y, + "x1": x1, + "x2": x2, + } + ) + + +@pytest.fixture(scope="module") +def rcs_data(): + return make_rcs_dml_data() + + +@pytest.fixture(scope="module") +def rcs_fitted(rcs_data): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return DMLDiD(panel=False, seed=0).fit(rcs_data, **FIT_KW, **COV) + + +class TestRCSConstruction: + @pytest.mark.parametrize("bad", ["False", 0.0, 1, None, [True]]) + def test_panel_strict_bool(self, bad): + with pytest.raises(ValueError, match="panel must be a bool"): + DMLDiD(panel=bad) + + def test_get_set_params_roundtrip(self): + est = DMLDiD(panel=False, seed=3) + params = est.get_params() + assert params["panel"] is False + est2 = DMLDiD().set_params(**params) + assert est2.panel is False + + def test_mutated_panel_rejected_at_fit(self, rcs_data): + est = DMLDiD(panel=False, seed=0) + est.panel = "nope" + with pytest.raises(ValueError, match="panel must be a bool"): + est.fit(rcs_data, **FIT_KW, **COV) + + +class TestRCSInputValidation: + def test_duplicate_unit_ids_rejected(self, rcs_data): + df = rcs_data.copy() + df.loc[df.index[1], "unit"] = df.loc[df.index[0], "unit"] + with pytest.raises(ValueError, match="unique unit IDs"): + DMLDiD(panel=False).fit(df, **FIT_KW, **COV) + + def test_panel_frame_under_rcs_hits_unique_id_error(self, data): + # The panel duplicate-(unit,time) fixture frame has repeated unit + # IDs — under panel=False it hits the unique-ID error, not the + # duplicate-(unit,time) one. + with pytest.raises(ValueError, match="unique unit IDs"): + DMLDiD(panel=False).fit(data, **FIT_KW, **COV) + + def test_stationarity_warning_fires_only_under_rcs(self, rcs_data, data): + # The warning must state the CORRECT Assumption 2.3 interpretation: + # stable (D, X) wave composition with period-specific potential + # outcomes — NOT a stable observed-outcome distribution (trends and + # treatment effects are expected, not violations). + with pytest.warns(UserWarning, match="composition of .D, X. is stable"): + DMLDiD(panel=False, seed=0).fit(rcs_data, **FIT_KW, **COV) + with warnings.catch_warnings(record=True) as rec: + warnings.simplefilter("always") + DMLDiD(seed=0).fit(data, **FIT_KW, **COV) + assert not any("stationary cross-sectional" in str(w.message) for w in rec) + + def test_covariates_still_required(self, rcs_data): + with pytest.raises(ValueError, match="CallawaySantAnna"): + DMLDiD(panel=False).fit(rcs_data, **FIT_KW, covariates=None) + + def test_label_pipeline_witness_string_labels(self, rcs_data): + # The label pipeline is mode-independent — one witness on RCS. + df = rcs_data.copy() + df["time"] = df["time"].astype(str) + df["first_treat"] = df["first_treat"].astype(str) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = DMLDiD(panel=False, seed=0).fit(df, **FIT_KW, **COV) + assert np.isfinite(res.att) + + +class TestRCSEstimation: + def test_finite_att_se_and_recovery(self, rcs_fitted): + assert np.isfinite(rcs_fitted.att) and np.isfinite(rcs_fitted.se) + assert abs(rcs_fitted.att - 2.0) < 0.6 # DGP effect = 2.0 + + def test_determinism_same_seed(self, rcs_data): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + a = DMLDiD(panel=False, seed=5).fit(rcs_data, **FIT_KW, **COV) + b = DMLDiD(panel=False, seed=5).fit(rcs_data, **FIT_KW, **COV) + assert a.att == b.att and a.se == b.se + + def test_seed_none_entropy_surfaced(self, rcs_data): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + a = DMLDiD(panel=False).fit(rcs_data, **FIT_KW, **COV) + b = DMLDiD(panel=False).fit(rcs_data, **FIT_KW, **COV) + ent_a = {e["fold_seed"]["entropy"] for e in a.cross_fit_diagnostics.values()} + ent_b = {e["fold_seed"]["entropy"] for e in b.cross_fit_diagnostics.values()} + assert ent_a and ent_b and ent_a != ent_b + + +class TestRCSPayloadContract: + def test_if_payload_matches_cell_se(self, rcs_fitted): + kit = rcs_fitted._aggregation_kit + checked = 0 + for (g, t), entry in rcs_fitted.group_time_effects.items(): + if entry["skip_reason"] is not None or entry.get("is_reference"): + continue + ii = kit.influence[(g, t)] + se = np.sqrt(np.sum(ii["treated_inf"] ** 2) + np.sum(ii["control_inf"] ** 2)) + np.testing.assert_allclose(se, entry["se"], rtol=1e-12, atol=0) + checked += 1 + assert checked > 0 + + def test_index_arrays_disjoint_increasing_both_periods(self, rcs_fitted, rcs_data): + kit = rcs_fitted._aggregation_kit + obs_time = rcs_data["time"].to_numpy() + for (g, t), ii in kit.influence.items(): + ti, ci = ii["treated_idx"], ii["control_idx"] + if len(ti) == 0 and len(ci) == 0: + continue # universal reference cells + assert np.all(np.diff(ti) > 0) and np.all(np.diff(ci) > 0) + assert len(np.intersect1d(ti, ci)) == 0 + union_times = set(obs_time[np.concatenate([ti, ci])].tolist()) + assert len(union_times) == 2 # rows from BOTH periods + # A base-period treated row IS in treated_idx. + assert len(set(obs_time[ti].tolist())) == 2 + + def test_aggregate_replay_leaves_payload_bit_identical(self, rcs_fitted): + kit = rcs_fitted._aggregation_kit + before = { + k: {kk: np.array(vv, copy=True) for kk, vv in v.items()} + for k, v in kit.influence.items() + } + rcs_fitted.aggregate("event_study") + for k, v in kit.influence.items(): + for kk in v: + np.testing.assert_array_equal(v[kk], before[k][kk]) + + +class TestRCSDegenerateHandling: + def test_empty_four_group_zero_treated_control(self): + from tests.conftest import assert_nan_inference + + # Cohort 3 has NO control rows at its base period 2 under + # never_treated?? — construct directly: remove every control row in + # period 1 so any cell with base 1 has an empty control-base group. + df = make_rcs_dml_data(n_rows=800, seed=9) + drop = (df["first_treat"] == 0) & (df["time"] == 1) + df = df[~drop].reset_index(drop=True) + df["unit"] = np.arange(len(df)) + with warnings.catch_warnings(record=True) as rec: + warnings.simplefilter("always") + res = DMLDiD(panel=False, seed=0).fit(df, **FIT_KW, **COV) + skipped = [ + (k, e) + for k, e in res.group_time_effects.items() + if e["skip_reason"] == "zero_treated_control" + ] + assert skipped + for _k, e in skipped: + assert_nan_inference(e) + assert any("could not be estimated" in str(w.message) for w in rec) + + def test_singleton_stratum_cross_fit_degenerate(self): + from tests.conftest import assert_nan_inference + + # Exactly ONE treated row of cohort 2 in the base period 1: the + # D x T stratum is a singleton -> assign_folds raises -> skip. + df = make_rcs_dml_data(n_rows=600, seed=11) + mask = (df["first_treat"] == 2) & (df["time"] == 1) + keep_one = df.index[mask][:1] + df = df[~mask | df.index.isin(keep_one)].reset_index(drop=True) + df["unit"] = np.arange(len(df)) + with warnings.catch_warnings(record=True) as rec: + warnings.simplefilter("always") + res = DMLDiD(panel=False, seed=0).fit(df, **FIT_KW, **COV) + degen = [ + e for e in res.group_time_effects.values() if e["skip_reason"] == "cross_fit_degenerate" + ] + assert degen + for e in degen: + assert_nan_inference(e) + assert any("could not be estimated" in str(w.message) for w in rec) + + def test_learner_failure_maps_to_cross_fit_degenerate(self): + from tests.conftest import assert_nan_inference + + # A covariate CONSTANT within one cell's rows: the fail-closed + # LinearLearner raises rank-deficiency -> DegenerateFoldError -> + # cross_fit_degenerate in the RCS branch. + df = make_rcs_dml_data(n_rows=900, seed=13) + cell_rows = df["time"].isin([1, 2]) + df.loc[cell_rows, "x2"] = 1.0 + df.loc[cell_rows, "x1"] = 1.0 + with warnings.catch_warnings(record=True) as rec: + warnings.simplefilter("always") + res = DMLDiD(panel=False, seed=0).fit(df, **FIT_KW, **COV) + degen = [ + e for e in res.group_time_effects.values() if e["skip_reason"] == "cross_fit_degenerate" + ] + assert degen + for e in degen: + assert_nan_inference(e) + assert any("could not be estimated" in str(w.message) for w in rec) + + def test_injected_overflow_non_finite_score(self, rcs_data): + from tests.conftest import assert_nan_inference + + class OverflowRegressor: + def fit(self, X, y): + return self + + def predict(self, X): + return np.full(len(X), 1e308) + + with warnings.catch_warnings(record=True) as rec: + warnings.simplefilter("always") + try: + res = DMLDiD(panel=False, outcome_learner=OverflowRegressor(), seed=0).fit( + rcs_data, **FIT_KW, **COV + ) + except ValueError: + return # all cells degenerate: loud failure is acceptable + nf = [e for e in res.group_time_effects.values() if e["skip_reason"] == "non_finite_score"] + assert nf + for e in nf: + assert_nan_inference(e) + assert any("non_finite_score" in str(w.message) for w in rec) + + def test_all_degenerate_raises_before_reference_cells(self): + df = make_rcs_dml_data(n_rows=40, seed=15) + # Too few rows for any cell to cross-fit with K=5 strata. + df = df.iloc[:16].reset_index(drop=True) + df["unit"] = np.arange(len(df)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match="Could not estimate any"): + DMLDiD(panel=False, base_period="universal", seed=0).fit(df, **FIT_KW, **COV) + + def test_lam_hat_extremeness_warning(self): + # Lopsided periods: keep exactly two period-1 rows per cohort, so + # every base-1 cell has lam_hat far above 1 - pscore_trim while the + # four-group guard still passes. + df = make_rcs_dml_data(n_rows=3000, seed=17) + keep = np.zeros(len(df), dtype=bool) + at1 = df["time"] == 1 + keep[~at1.to_numpy()] = True + for cohort in (0, 2, 3): + idx = df.index[at1 & (df["first_treat"] == cohort)][:2] + keep[idx] = True + df = df[keep].reset_index(drop=True) + df["unit"] = np.arange(len(df)) + with warnings.catch_warnings(record=True) as rec: + warnings.simplefilter("always") + try: + DMLDiD(panel=False, seed=0).fit(df, **FIT_KW, **COV) + except ValueError: + pass + assert any("lam_hat" in str(w.message) for w in rec) + + +class TestRCSAggregationBootstrap: + def test_event_study_group_simple(self, rcs_fitted): + es = rcs_fitted.aggregate("event_study") + assert len(es.to_dataframe()) > 1 + gr = rcs_fitted.aggregate("group") + assert len(gr.to_dataframe()) >= 1 + si = rcs_fitted.aggregate("simple") + np.testing.assert_allclose(si.to_dataframe()["att"].iloc[0], rcs_fitted.att) + + def test_total_fails_closed(self, rcs_fitted): + with pytest.raises(NotImplementedError, match="repeated-cross-section"): + rcs_fitted.aggregate("total") + + def test_bootstrap_runs_with_per_row_weights(self, rcs_data): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + boot = DMLDiD(panel=False, seed=0, n_bootstrap=49).fit(rcs_data, **FIT_KW, **COV) + assert np.isfinite(boot.att) + assert boot.bootstrap_results is not None + es = boot.aggregate("event_study") + df_es = es.to_dataframe() + assert any("cband" in c for c in df_es.columns) or es.cband_crit_value is not None + + +class TestRCSAggWeights: + def test_weights_are_fixed_cohort_row_masses(self, rcs_fitted, rcs_data): + # aggregate('simple') equals the hand-computed agg_weight-weighted + # combination of post-treatment finite cells (WIF-consistency + # decision: fixed cohort row masses, never per-cell counts). + cohort_mass = rcs_data.groupby("first_treat").size().to_dict() + num = 0.0 + den = 0.0 + for (g, t), e in rcs_fitted.group_time_effects.items(): + if e["skip_reason"] is not None or e.get("is_reference"): + continue + if t < g - rcs_fitted.anticipation: + continue + if not np.isfinite(e["effect"]): + continue + w = e["agg_weight"] + assert w == cohort_mass[g] + num += w * e["effect"] + den += w + np.testing.assert_allclose(rcs_fitted.att, num / den, rtol=1e-12) + + def test_large_int64_cohort_labels_aggregate_and_replay(self): + # >2**53 int64 cohort labels: the float-key hazard that motivated + # leaving agg_cohort_masses unset stays off the shipped path — + # aggregation and bootstrap replay run end-to-end. + big = 2**60 + df = make_rcs_dml_data(n_rows=1200, seed=19) + tmap = {1: big + 1, 2: big + 2, 3: big + 3} + df["time"] = df["time"].map(tmap) + df["first_treat"] = df["first_treat"].map({0: 0, 2: big + 2, 3: big + 3}) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = DMLDiD(panel=False, seed=0, n_bootstrap=29).fit(df, **FIT_KW, **COV) + assert np.isfinite(res.att) + assert sorted(res.groups) == [big + 2, big + 3] + es = res.aggregate("event_study") + ets = list(es.to_dataframe()["event_time"]) + assert len(ets) == len(set(ets)) and len(ets) > 1 + + +class TestRCSResultsSurface: + def test_panel_flag_and_summary(self, rcs_fitted, fitted): + assert rcs_fitted.panel is False + s = rcs_fitted.summary() + assert "repeated cross sections" in s + assert "obs:" in s + assert "repeated cross sections" not in fitted.summary() + + def test_to_dict_json_roundtrip(self, rcs_fitted): + d = rcs_fitted.to_dict() + blob = json.dumps(d) + assert d["panel"] is False + parsed = json.loads(blob) + any_cell = next(iter(parsed["cross_fit_diagnostics"].values())) + assert "lam_hat" in any_cell and "g2_lambda" in any_cell + + def test_business_report_rcs_semantics(self, rcs_fitted, rcs_data): + from diff_diff import BusinessReport + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + rep = BusinessReport( + rcs_fitted, + data=rcs_data, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + ).full_report() + assert "observations" in rep + assert "stationary" in rep.lower() + # Correct Assumption 2.3 interpretation (stable wave composition, + # not a stable observed-Y distribution). + assert "composition of (D, X) is stable" in rep + assert "period-specific potential" in rep + + def test_target_parameter_design_aware(self, rcs_fitted, fitted): + from diff_diff._reporting_helpers import describe_target_parameter + + rcs_block = describe_target_parameter(rcs_fitted) + assert "cohort-mass-weighted" in rcs_block["name"] + panel_block = describe_target_parameter(fitted) + assert "valid-treated-count-weighted" in panel_block["name"] + + def test_diagnostic_report_bacon_skipped_on_rcs(self, rcs_fitted, rcs_data): + from diff_diff import DiagnosticReport + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + dr = DiagnosticReport( + rcs_fitted, + data=rcs_data, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + ) + rep = dr.run_all() + assert "bacon" in rep.skipped_checks + assert "requires panel data" in rep.skipped_checks["bacon"] + # Panel fits keep running the decomposition (no skip). + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + panel_df = make_staggered_dml_data() + panel_res = DMLDiD(seed=0).fit(panel_df, **FIT_KW, **COV) + panel_rep = DiagnosticReport( + panel_res, + data=panel_df, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + ).run_all() + assert "bacon" not in panel_rep.skipped_checks + + def test_practitioner_snippet_carries_panel_false(self, rcs_fitted, fitted): + from diff_diff import practitioner_next_steps + + steps = practitioner_next_steps(rcs_fitted) + text = str(steps) + assert "panel=False" in text + assert "panel=False" not in str(practitioner_next_steps(fitted)) + + def test_plot_smoke(self, rcs_fitted): + pytest.importorskip("matplotlib") + import matplotlib + + matplotlib.use("Agg") + from diff_diff import plot_event_study, plot_group_effects + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + plot_group_effects(rcs_fitted, show=False) + plot_event_study(rcs_fitted.aggregate("event_study"), show=False) + + +class TestRCSSemantics: + def test_not_yet_treated_differs(self, rcs_data): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + nt = DMLDiD(panel=False, seed=0).fit(rcs_data, **FIT_KW, **COV) + nyt = DMLDiD(panel=False, seed=0, control_group="not_yet_treated").fit( + rcs_data, **FIT_KW, **COV + ) + assert nt.att != nyt.att + + def test_anticipation_shifts_base(self, rcs_data): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + a0 = DMLDiD(panel=False, seed=0).fit(rcs_data, **FIT_KW, **COV) + a1 = DMLDiD(panel=False, seed=0, anticipation=1).fit(rcs_data, **FIT_KW, **COV) + assert set(a0.group_time_effects) != set(a1.group_time_effects) + + def test_universal_base_reference_cells(self, rcs_data): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = DMLDiD(panel=False, seed=0, base_period="universal").fit( + rcs_data, **FIT_KW, **COV + ) + refs = [e for e in res.group_time_effects.values() if e.get("is_reference")] + assert refs + for e in refs: + assert e["effect"] == 0.0 and "agg_weight" in e + assert res.reference_event_times is not None + + +class TestRCSCompleteCases: + def test_missing_outcome_row_drops_from_cell_only(self): + df = make_rcs_dml_data(n_rows=1200, seed=21) + victim = df.index[(df["first_treat"] == 2) & (df["time"] == 2)][0] + df.loc[victim, "y"] = np.nan + with pytest.warns(UserWarning, match="observation.s. were excluded"): + with warnings.catch_warnings(): + warnings.simplefilter("always") + res = DMLDiD(panel=False, seed=0).fit(df, **FIT_KW, **COV) + assert np.isfinite(res.att) + + def test_non_finite_covariate_row_excluded(self): + df = make_rcs_dml_data(n_rows=1200, seed=23) + victim = df.index[(df["first_treat"] == 0) & (df["time"] == 3)][0] + df.loc[victim, "x1"] = np.inf + with pytest.warns(UserWarning, match="observation.s. were excluded"): + with warnings.catch_warnings(): + warnings.simplefilter("always") + res = DMLDiD(panel=False, seed=0).fit(df, **FIT_KW, **COV) + assert np.isfinite(res.att) + + def test_inf_outcome_handled(self): + df = make_rcs_dml_data(n_rows=1200, seed=25) + victim = df.index[(df["first_treat"] == 3) & (df["time"] == 3)][0] + df.loc[victim, "y"] = -np.inf + with pytest.warns(UserWarning, match="observation.s. were excluded"): + with warnings.catch_warnings(): + warnings.simplefilter("always") + res = DMLDiD(panel=False, seed=0).fit(df, **FIT_KW, **COV) + assert np.isfinite(res.att) + + def test_no_warning_on_clean_input(self, rcs_data): + with warnings.catch_warnings(record=True) as rec: + warnings.simplefilter("always") + DMLDiD(panel=False, seed=0).fit(rcs_data, **FIT_KW, **COV) + assert not any("were excluded" in str(w.message) for w in rec) diff --git a/tests/test_dr_scores.py b/tests/test_dr_scores.py index 84857c212..5094b3e9b 100644 --- a/tests/test_dr_scores.py +++ b/tests/test_dr_scores.py @@ -222,3 +222,101 @@ def test_augmented_validation(self): chang_panel_score_augmented(summand, np.array([0.0, 2.0, 0.0, 1.0]), 1.0, 0.5) with pytest.raises(ValueError, match="length"): chang_panel_score_augmented(summand[:-1], D, 1.0, 0.5) + + +class TestChangRCSScoreValidation: + """Both sides of every stated Case 2 domain raise targeted errors.""" + + def _inputs(self, n=24): + rng = np.random.default_rng(5) + y = rng.normal(size=n) + D = (rng.uniform(size=n) < 0.5).astype(float) + T = (rng.uniform(size=n) < 0.5).astype(float) + m2_hat = rng.normal(size=n) + ps = np.clip(rng.uniform(size=n), 0.1, 0.9) + return y, D, T, m2_hat, ps + + def test_ps_bounds_raise(self): + from diff_diff._dr_scores import chang_rcs_score + + y, D, T, m2, ps = self._inputs() + for bad in (1.0, -0.01): + p = ps.copy() + p[0] = bad + with pytest.raises(ValueError, match=r"ps must lie in \[0, 1\)"): + chang_rcs_score(y, D, T, m2, p, 0.5, 0.5) + + @pytest.mark.parametrize("lam_hat", [0.0, 1.0, -0.2, 1.5, np.nan]) + def test_lam_hat_domain_raises(self, lam_hat): + from diff_diff._dr_scores import chang_rcs_score + + y, D, T, m2, ps = self._inputs() + with pytest.raises(ValueError, match="lam_hat"): + chang_rcs_score(y, D, T, m2, ps, 0.5, lam_hat) + + @pytest.mark.parametrize("p_hat", [0.0, 1.0, np.nan]) + def test_p_hat_domain_raises(self, p_hat): + from diff_diff._dr_scores import chang_rcs_score + + y, D, T, m2, ps = self._inputs() + with pytest.raises(ValueError, match="p_hat"): + chang_rcs_score(y, D, T, m2, ps, p_hat, 0.5) + + def test_non_binary_T_raises(self): + from diff_diff._dr_scores import chang_rcs_score + + y, D, T, m2, ps = self._inputs() + T[0] = 0.5 + with pytest.raises(ValueError, match="T must be strictly binary"): + chang_rcs_score(y, D, T, m2, ps, 0.5, 0.5) + + def test_non_binary_D_raises(self): + from diff_diff._dr_scores import chang_rcs_score + + y, D, T, m2, ps = self._inputs() + D[0] = 2.0 + with pytest.raises(ValueError, match="D must be strictly binary"): + chang_rcs_score(y, D, T, m2, ps, 0.5, 0.5) + + def test_shape_empty_scalar_nonfinite_raise(self): + from diff_diff._dr_scores import chang_rcs_score + + y, D, T, m2, ps = self._inputs() + with pytest.raises(ValueError, match="length"): + chang_rcs_score(y, D, T, m2[:-1], ps, 0.5, 0.5) + empty = np.empty(0) + with pytest.raises(ValueError, match="empty"): + chang_rcs_score(empty, empty, empty, empty, empty, 0.5, 0.5) + with pytest.raises(ValueError, match="1-dimensional"): + chang_rcs_score(1.0, 1.0, 1.0, 1.0, 0.5, 0.5, 0.5) + y_bad = y.copy() + y_bad[0] = np.inf + with pytest.raises(ValueError, match="non-finite"): + chang_rcs_score(y_bad, D, T, m2, ps, 0.5, 0.5) + + def test_lambda_slope_shares_validator(self): + from diff_diff._dr_scores import chang_rcs_lambda_slope + + y, D, T, m2, ps = self._inputs() + with pytest.raises(ValueError, match="lam_hat"): + chang_rcs_lambda_slope(y, D, T, m2, ps, 0.5, 1.0) + + def test_augmented_validation(self): + from diff_diff._dr_scores import chang_rcs_score_augmented + + y, D, T, m2, ps = self._inputs(n=6) + summand = np.zeros(6) + with pytest.raises(ValueError, match="lam_hat"): + chang_rcs_score_augmented(summand, D, T, y, m2, ps, 1.0, 0.5, 0.0) + with pytest.raises(ValueError, match="theta"): + chang_rcs_score_augmented(summand, D, T, y, m2, ps, np.nan, 0.5, 0.5) + with pytest.raises(ValueError, match="summand contains non-finite"): + bad = summand.copy() + bad[0] = np.nan + chang_rcs_score_augmented(bad, D, T, y, m2, ps, 1.0, 0.5, 0.5) + with pytest.raises(ValueError, match="summand has length"): + chang_rcs_score_augmented(summand[:-1], D, T, y, m2, ps, 1.0, 0.5, 0.5) + with pytest.raises(ValueError, match="T must be strictly binary"): + T_bad = T.copy() + T_bad[0] = 3.0 + chang_rcs_score_augmented(summand, D, T_bad, y, m2, ps, 1.0, 0.5, 0.5) diff --git a/tests/test_methodology_dml_did.py b/tests/test_methodology_dml_did.py index 0fa1f8239..fe966b9fc 100644 --- a/tests/test_methodology_dml_did.py +++ b/tests/test_methodology_dml_did.py @@ -11,7 +11,12 @@ from diff_diff import DMLDiD from diff_diff._crossfit import assign_folds, cross_fit_predict -from diff_diff._dr_scores import chang_panel_score, chang_panel_score_augmented +from diff_diff._dr_scores import ( + chang_panel_score, + chang_panel_score_augmented, + chang_rcs_score, + chang_rcs_score_augmented, +) from diff_diff._learners import LinearLearner, LogitLearner @@ -331,3 +336,358 @@ def test_misspecified_outcome_correct_propensity(self): df, **FIT_KW, covariates=["x1", "x2"] ) assert abs(res.att - theta) < 0.25 + + +# =========================================================================== +# Repeated cross sections (Chang Case 2, panel=False) +# =========================================================================== +# +# Shared fixed RCS DGP (library-authored, in the spirit of Chang Sec. 4 — +# NOT the paper's own high-dimensional Sec. 4 RCS parameterization, whose +# replication is a tracked TODO row): X ~ N(0, I_2), D ~ +# Bernoulli(sigmoid(0.5 X1 - 0.5 X2)), T ~ Bernoulli(0.5), levels +# Y = 1 + X1 + 0.5 X2 + T*(0.5 + 0.4 X1) + D + T*D*theta0 + eps. + +RCS_THETA0 = 3.0 + + +def _rcs_frame(n_rows, seed): + rng = np.random.default_rng(seed) + X = rng.standard_normal((n_rows, 2)) + g0 = 1.0 / (1.0 + np.exp(-(0.5 * X[:, 0] - 0.5 * X[:, 1]))) + D = (rng.uniform(size=n_rows) < g0).astype(int) + T = (rng.uniform(size=n_rows) < 0.5).astype(int) + trend = 0.5 + 0.4 * X[:, 0] + y = ( + 1.0 + + X[:, 0] + + 0.5 * X[:, 1] + + T * trend + + D * 1.0 + + T * D * RCS_THETA0 + + rng.normal(size=n_rows) + ) + return pd.DataFrame( + { + "unit": np.arange(n_rows), + "time": T + 1, # periods {1, 2}; cohort 2 = treated at period 2 + "first_treat": D * 2, + "y": y, + "x1": X[:, 0], + "x2": X[:, 1], + } + ) + + +class TestRCSOracleEquivalence: + def test_oracle_learners_match_closed_form(self): + # ORACLE nuisances on BOTH sides: the true DGP propensity and the + # true l20(X) = E[(T - lam0)Y | X, D=0] = 0.25*(0.4 + 0.3*x1) at + # lam0 = 0.5 (T independent of (X, D), so the x1 and eps terms are + # killed by E[T - lam0] = 0 and E[(T - lam0)T] = 0.25 survives). + # The public DMLDiD(panel=False) single-cell fit must equal the + # hand Eq 3.2 / Thm 2 pipeline at 1e-12 (algebra equivalence — + # exact for ANY fixed nuisance; using the true one keeps "oracle" + # honest). + rng = np.random.default_rng(31) + n = 600 + X = rng.standard_normal((n, 2)) + ps_true = 1 / (1 + np.exp(-X[:, 0])) + D = (rng.uniform(size=n) < ps_true).astype(float) + T = (rng.uniform(size=n) < 0.5).astype(float) + y = X[:, 0] + T * (0.4 + 0.3 * X[:, 0]) + 2.0 * T * D + rng.standard_normal(n) + df = pd.DataFrame( + { + "unit": np.arange(n), + "time": T.astype(int) + 1, + "first_treat": (D * 2).astype(int), + "y": y, + "x1": X[:, 0], + "x2": X[:, 1], + } + ) + + class OracleProp: + def fit(self, Xf, yf, sample_weight=None): + return self + + def predict_proba(self, Xf): + p = 1 / (1 + np.exp(-Xf[:, 0])) + return np.column_stack([1 - p, p]) + + class OracleReg: + def fit(self, Xf, yf, sample_weight=None): + return self + + def predict(self, Xf): + return 0.1 + 0.075 * Xf[:, 0] + + trim = 0.01 + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = DMLDiD( + propensity_learner=OracleProp(), + outcome_learner=OracleReg(), + seed=0, + pscore_trim=trim, + panel=False, + ).fit(df, **FIT_KW, covariates=["x1", "x2"]) + ps = np.clip(ps_true, trim, 1 - trim) + m2 = 0.1 + 0.075 * X[:, 0] + p_hat = float(D.mean()) + lam = float(T.mean()) + summand = chang_rcs_score(y, D, T, m2, ps, p_hat, lam) + theta = float(np.mean(summand)) + psi_bar = chang_rcs_score_augmented(summand, D, T, y, m2, ps, theta, p_hat, lam) + se = float(np.sqrt(np.mean(psi_bar**2) / n)) + np.testing.assert_allclose(res.att, theta, rtol=0, atol=1e-12) + np.testing.assert_allclose(res.se, se, rtol=0, atol=1e-12) + + def test_two_period_cell_equivalence_vs_hand_pipeline(self): + # Real (non-oracle) learners: the fit equals a hand pipeline that + # replays the SAME spawned folds + cross_fit_predict + rcs scores. + # Tolerance 1e-14 rtol, not bit-for-bit (BLAS layout lesson, B1 CI). + from diff_diff._crossfit import assign_folds, cross_fit_predict + from diff_diff._learners import LinearLearner, LogitLearner + + df = _rcs_frame(800, seed=33) + seed = 9 + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = DMLDiD(seed=seed, panel=False).fit(df, **FIT_KW, covariates=["x1", "x2"]) + D = df["first_treat"].to_numpy().astype(float) / 2.0 + T = (df["time"].to_numpy() == 2).astype(float) + y = df["y"].to_numpy() + X = df[["x1", "x2"]].to_numpy() + n = len(df) + rng = np.random.default_rng(np.random.SeedSequence(entropy=seed, spawn_key=(0, 1))) + folds = assign_folds(n, 5, rng=rng, stratify=D + 2.0 * T) + ps_res = cross_fit_predict(LogitLearner(), X, D, folds, predict_method="predict_proba") + lam = float(T.mean()) + r = (T - lam) * y + or_res = cross_fit_predict( + LinearLearner(), X, r, folds, predict_method="predict", fit_mask=(D == 0.0) + ) + ps = np.clip(ps_res.oof_predictions, 0.01, 0.99) + p_hat = float(D.mean()) + summand = chang_rcs_score(y, D, T, or_res.oof_predictions, ps, p_hat, lam) + theta = float(np.mean(summand)) + psi_bar = chang_rcs_score_augmented( + summand, D, T, y, or_res.oof_predictions, ps, theta, p_hat, lam + ) + se = float(np.sqrt(np.mean(psi_bar**2) / n)) + np.testing.assert_allclose(res.att, theta, rtol=1e-14, atol=0) + np.testing.assert_allclose(res.se, se, rtol=1e-14, atol=0) + + +class TestLambdaCorrectionRegression: + def test_reported_se_carries_the_lambda_term(self): + # Headline guard for the review's "plausible implementation bug": + # from the fit's payload + g2_lambda diagnostic, recompute the SE + # with the lambda term REMOVED and assert it differs from the + # reported SE in the pinned direction. + df = _rcs_frame(2000, seed=35) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = DMLDiD(seed=0, panel=False).fit(df, **FIT_KW, covariates=["x1", "x2"]) + (gt_key,) = [k for k, e in res.group_time_effects.items() if e["skip_reason"] is None] + entry = res.group_time_effects[gt_key] + diag = res.cross_fit_diagnostics[gt_key] + g2 = diag["g2_lambda"] + assert np.isfinite(g2) and g2 != 0.0 + kit = res._aggregation_kit + ii = kit.influence[gt_key] + idx = np.concatenate([ii["treated_idx"], ii["control_idx"]]) + phi = np.concatenate([ii["treated_inf"], ii["control_inf"]]) + n_cell = len(idx) + psi_bar = phi * n_cell # payload entries are psi_bar / n_cell + T_cell = (df["time"].to_numpy()[idx] == 2).astype(float) + lam = diag["lam_hat"] + psi_no_lambda = psi_bar - g2 * (T_cell - lam) + se_no_lambda = float(np.sqrt(np.mean(psi_no_lambda**2) / n_cell)) + np.testing.assert_allclose( + entry["se"], float(np.sqrt(np.mean(psi_bar**2) / n_cell)), rtol=1e-12 + ) + rel_gap = abs(se_no_lambda - entry["se"]) / entry["se"] + assert rel_gap > 1e-5, (se_no_lambda, entry["se"]) + + +class TestRCSATTRecovery: + def test_theta_recovery(self): + df = _rcs_frame(4000, seed=37) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = DMLDiD(seed=0, panel=False).fit(df, **FIT_KW, covariates=["x1", "x2"]) + assert abs(res.att - RCS_THETA0) < 0.3 + + def test_double_robustness_misspecified_propensity(self): + # Threshold (non-logistic) treatment rule misspecifies the logit + # propensity; correct-in-X outcome regression still recovers theta. + rng = np.random.default_rng(39) + n = 5000 + X = rng.standard_normal((n, 2)) + # Noise scale 1.5 keeps the misspecified fitted propensities away + # from the clip bounds — the RCS weight carries a 1/(lam(1-lam)) + # amplification on level outcomes, so a near-deterministic rule + # turns clipping bias into the dominant finite-sample error. + D = ((X[:, 0] + 1.5 * rng.standard_normal(n)) > 0).astype(int) + T = (rng.uniform(size=n) < 0.5).astype(int) + y = ( + 1.0 + + X[:, 0] + + 0.5 * X[:, 1] + + T * (0.5 + 0.4 * X[:, 0]) + + D * 1.0 + + T * D * RCS_THETA0 + + rng.normal(size=n) + ) + df = pd.DataFrame( + { + "unit": np.arange(n), + "time": T + 1, + "first_treat": D * 2, + "y": y, + "x1": X[:, 0], + "x2": X[:, 1], + } + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = DMLDiD(seed=0, panel=False).fit(df, **FIT_KW, covariates=["x1", "x2"]) + assert abs(res.att - RCS_THETA0) < 0.3 + + def test_double_robustness_misspecified_outcome_regression(self): + # Nonlinear-in-X outcome trend misspecifies the linear l2 learner; + # the correctly-specified (logistic) propensity still recovers theta. + rng = np.random.default_rng(41) + n = 5000 + X = rng.standard_normal((n, 2)) + g0 = 1.0 / (1.0 + np.exp(-(0.5 * X[:, 0] - 0.5 * X[:, 1]))) + D = (rng.uniform(size=n) < g0).astype(int) + T = (rng.uniform(size=n) < 0.5).astype(int) + y = ( + 1.0 + + X[:, 0] + + 0.5 * X[:, 1] ** 2 # nonlinear: linear learner misspecified + + T * (0.5 + 0.4 * X[:, 0] ** 2) + + D * 1.0 + + T * D * RCS_THETA0 + + rng.normal(size=n) + ) + df = pd.DataFrame( + { + "unit": np.arange(n), + "time": T + 1, + "first_treat": D * 2, + "y": y, + "x1": X[:, 0], + "x2": X[:, 1], + } + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = DMLDiD(seed=0, panel=False).fit(df, **FIT_KW, covariates=["x1", "x2"]) + assert abs(res.att - RCS_THETA0) < 0.3 + + +def _characterization_spike_frame(): + rng = np.random.default_rng(7) + n_rows = 4000 + cohort = rng.choice([0, 3, 4], size=n_rows, p=[0.5, 0.25, 0.25]) + tt = rng.choice([1, 2, 3, 4], size=n_rows) + X_row = rng.standard_normal((n_rows, 2)) + y = ( + 0.8 * X_row[:, 0] + - 0.4 * X_row[:, 1] + + 0.3 * tt + + 0.5 * X_row[:, 0] * tt / 4 + + rng.standard_normal(n_rows) + ) + post = (cohort > 0) & (tt >= cohort) + y = y + post * (2.0 + 0.2 * (tt - cohort)) + return pd.DataFrame( + { + "id": np.arange(n_rows), + "t": tt, + "g": cohort, + "y": y, + "x1": X_row[:, 0], + "x2": X_row[:, 1], + } + ) + + +@pytest.fixture(scope="module") +def characterization_spike_fit(): + df = _characterization_spike_frame() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return DMLDiD(seed=11, panel=False).fit( + df, outcome="y", unit="id", time="t", first_treat="g", covariates=["x1", "x2"] + ) + + +class TestRCSGoldenCharacterization: + """Golden literals from benchmarks/doubleml/chang_rcs_characterization.py + (doubleml 0.11.4, sklearn 1.9.0, SEED=7/SEED_DML=11). Part 2's CHANG-side + values (public DMLDiD(panel=False) == hand pipeline at 0.0 diff under + sklearn learners); the native reproduction swaps sklearn's lbfgs logit + for the library IRLS solver — tolerance atol 2e-4 (ATT) / 1e-5 (SE), the + B0/B1 optimizer-gap precedent. The Part 1 Chang-vs-DoubleML gaps are + additionally pinned NONZERO by recomputing them from a LIVE native fit + against the committed DoubleML ATT literals (characterization honesty: + the two scores must not coincide, and a native regression onto + DoubleML's score fails here). + """ + + GOLDEN = { + (3, 3, 2): (1.867260790882, 0.244040011648), + (3, 4, 2): (2.347579854514, 0.251882814022), + (4, 3, 2): (-0.089752999600, 0.187258740736), + (4, 4, 3): (1.822066292980, 0.260472487961), + } + # Part 1 DoubleMLDIDCSBinary ATT literals (shared folds; spike + # transcript). The gap test recomputes CHANG - DoubleML from a live + # native fit, so it fails if the native estimator ever drifts onto + # DoubleML's Sant'Anna-Zhao score. + DML_ATT = { + (3, 3, 2): 2.121300323064, + (3, 4, 2): 2.293055793734, + (4, 3, 2): -0.057963388953, + (4, 4, 3): 1.995045506760, + } + + def test_native_reproduction_matches_golden(self, characterization_spike_fit): + for (g, t_eval, _base), (att_gold, se_gold) in self.GOLDEN.items(): + e = characterization_spike_fit.group_time_effects[(g, t_eval)] + np.testing.assert_allclose(e["effect"], att_gold, rtol=0, atol=2e-4) + np.testing.assert_allclose(e["se"], se_gold, rtol=0, atol=1e-5) + + def test_characterization_gaps_are_nonzero(self, characterization_spike_fit): + # Honesty pin, recomputed LIVE: the native Chang ATT minus the + # committed DoubleML ATT must stay nonzero and bounded — a + # near-zero gap would mean the scores coincide (they must not; + # DoubleMLDIDCSBinary is NOT an oracle for Eq 3.2). Native-vs-spike + # optimizer drift is atol 2e-4, far below the 1e-3 floor; the + # smallest recorded gap is 1.73e-2. + for (g, t_eval, _base), dml_att in self.DML_ATT.items(): + gap = characterization_spike_fit.group_time_effects[(g, t_eval)]["effect"] - dml_att + assert 1e-3 < abs(gap) < 0.5 + + +@pytest.mark.slow +class TestRCSMonteCarloCoverage: + def test_coverage_sanity(self, ci_params): + n_reps = ci_params.bootstrap(200) + n = 800 + hits = 0 + for rep in range(n_reps): + df = _rcs_frame(n, seed=10_000 + rep) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = DMLDiD(seed=rep, panel=False).fit(df, **FIT_KW, covariates=["x1", "x2"]) + lo, hi = res.conf_int + hits += int(lo <= RCS_THETA0 <= hi) + coverage = hits / n_reps + lo_band, hi_band = (0.90, 0.99) if n_reps >= 100 else (0.80, 1.00) + assert lo_band <= coverage <= hi_band, coverage diff --git a/tests/test_methodology_dr_scores.py b/tests/test_methodology_dr_scores.py index 7a0f2bfa0..3c4507a9b 100644 --- a/tests/test_methodology_dr_scores.py +++ b/tests/test_methodology_dr_scores.py @@ -12,7 +12,13 @@ import numpy as np import pytest -from diff_diff._dr_scores import chang_panel_score, chang_panel_score_augmented +from diff_diff._dr_scores import ( + chang_panel_score, + chang_panel_score_augmented, + chang_rcs_lambda_slope, + chang_rcs_score, + chang_rcs_score_augmented, +) THETA0 = 2.5 @@ -184,3 +190,164 @@ def test_native_chang_estimator_matches_doubleml_goldens(self): np.testing.assert_allclose(theta, DOUBLEML_GOLDEN_ATT, atol=2e-4, rtol=0) np.testing.assert_allclose(se, DOUBLEML_GOLDEN_SE, atol=1e-5, rtol=0) + + +# --------------------------------------------------------------------------- +# Chang (2020) Case 2 (repeated cross sections) score methodology +# --------------------------------------------------------------------------- +# +# LIBRARY-AUTHORED low-dimensional RCS design in the spirit of Chang Sec. 4 +# (Gaussian X, logistic PS, theta0 = 3) — NOT the paper's own Sec. 4 RCS +# parameterization (a high-dimensional p in {100, 300} ML design, not +# extracted into the paper review; replication is a tracked TODO row). +# Rows i.i.d.: X ~ N(0, I_2); D ~ Bernoulli(sigmoid(0.5 X1 - 0.5 X2)); +# T ~ Bernoulli(lam0) independent; levels +# Y = 1 + X1 + 0.5 X2 + T*(0.5 + 0.4 X1) + D*1.0 + T*D*theta0 + eps. +# Under stationary sampling the ORACLE Case 2 outcome nuisance is +# l20(X) = E[(T - lam0) Y | X, D=0] = lam0 (1 - lam0) * (0.5 + 0.4 X1). + +RCS_THETA0 = 3.0 +RCS_LAM0 = 0.5 + + +def _rcs_dgp(n=200_000, seed=13): + rng = np.random.default_rng(seed) + X = rng.standard_normal((n, 2)) + g0 = 1.0 / (1.0 + np.exp(-(0.5 * X[:, 0] - 0.5 * X[:, 1]))) + D = (rng.uniform(size=n) < g0).astype(float) + T = (rng.uniform(size=n) < RCS_LAM0).astype(float) + trend = 0.5 + 0.4 * X[:, 0] + y = ( + 1.0 + + X[:, 0] + + 0.5 * X[:, 1] + + T * trend + + D * 1.0 + + T * D * RCS_THETA0 + + rng.normal(scale=1.0, size=n) + ) + ell20 = RCS_LAM0 * (1.0 - RCS_LAM0) * trend + return X, g0, D, T, ell20, y + + +class TestChangRCSScoreMethodology: + def test_score_mean_recovers_att_with_oracle_nuisances(self): + X, g0, D, T, ell20, y = _rcs_dgp() + ps = np.clip(g0, 1e-3, 1 - 1e-3) + summand = chang_rcs_score(y, D, T, ell20, ps, float(D.mean()), float(T.mean())) + assert abs(summand.mean() - RCS_THETA0) < 0.05 + + def test_double_robustness_misspecified_propensity(self): + # Wrong ps (constant), correct ell2 -> still recovers the ATT: the + # control-side residual E[(T-lam)Y - ell2 | X, D=0] is zero, and the + # wrong ps cancels from the treated-side weight. + X, g0, D, T, ell20, y = _rcs_dgp() + ps_wrong = np.full_like(g0, 0.3) + summand = chang_rcs_score(y, D, T, ell20, ps_wrong, float(D.mean()), float(T.mean())) + assert abs(summand.mean() - RCS_THETA0) < 0.05 + + def test_double_robustness_misspecified_outcome_regression(self): + # Wrong ell2 (zero), correct propensity -> still recovers the ATT: + # E[(D - g0)/(1 - g0) | X] = 0 kills the wrong-nuisance term. + X, g0, D, T, ell20, y = _rcs_dgp() + ps = np.clip(g0, 1e-3, 1 - 1e-3) + m_wrong = np.zeros_like(ell20) + summand = chang_rcs_score(y, D, T, m_wrong, ps, float(D.mean()), float(T.mean())) + assert abs(summand.mean() - RCS_THETA0) < 0.05 + + def test_closed_form_hand_fixture(self): + """Element-wise hand recomputation of summand, G_2lambda, psi_bar.""" + X, g0, D, T, ell20, y = _rcs_dgp(n=500, seed=6) + ps = np.clip(g0, 1e-3, 1 - 1e-3) + p_hat = float(D.mean()) + lam = float(T.mean()) + summand = chang_rcs_score(y, D, T, ell20, ps, p_hat, lam) + theta = float(summand.mean()) + psi_bar = chang_rcs_score_augmented(summand, D, T, y, ell20, ps, theta, p_hat, lam) + g2 = chang_rcs_lambda_slope(y, D, T, ell20, ps, p_hat, lam) + + # Hand formulas, written out independently (Eq 3.2 / Thm 2): + w = (D - ps) / (p_hat * lam * (1 - lam) * (1 - ps)) + hand_summand = w * ((T - lam) * y - ell20) + np.testing.assert_allclose(summand, hand_summand, atol=1e-15, rtol=0) + odds = (D - ps) / (1 - ps) + hand_g2 = float( + np.mean( + -((1 - 2 * lam) / (lam**2 * (1 - lam) ** 2)) + * (odds / p_hat) + * ((T - lam) * y - ell20) + - (y / (p_hat * lam * (1 - lam))) * odds + ) + ) + np.testing.assert_allclose(g2, hand_g2, atol=1e-12, rtol=0) + hand_psi_bar = hand_summand - D * theta / p_hat + hand_g2 * (T - lam) + np.testing.assert_allclose(psi_bar, hand_psi_bar, atol=1e-12, rtol=0) + + se = np.sqrt(np.mean(psi_bar**2) / len(D)) + hand_se = np.sqrt(np.mean(hand_psi_bar**2) / len(D)) + np.testing.assert_allclose(se, hand_se, atol=1e-15, rtol=0) + + def test_lambda_slope_first_term_algebraic_identity(self): + # term1 of d_lam psi_2 equals -((1-2*lam)/(lam*(1-lam))) * summand_i. + X, g0, D, T, ell20, y = _rcs_dgp(n=400, seed=8) + ps = np.clip(g0, 1e-3, 1 - 1e-3) + p_hat, lam = float(D.mean()), 0.4 # lam != mean(T): identity is algebraic + summand = chang_rcs_score(y, D, T, ell20, ps, p_hat, lam) + odds = (D - ps) / (1 - ps) + term1 = ( + -((1 - 2 * lam) / (lam**2 * (1 - lam) ** 2)) * (odds / p_hat) * ((T - lam) * y - ell20) + ) + np.testing.assert_allclose( + term1, -((1 - 2 * lam) / (lam * (1 - lam))) * summand, atol=1e-12, rtol=0 + ) + + def test_lambda_slope_matches_finite_difference(self): + # G_2lambda is d/d(lam) of mean(summand) holding the nuisances fixed: + # central differences of the SCORE function in lam_hat must match. + X, g0, D, T, ell20, y = _rcs_dgp(n=2_000, seed=10) + ps = np.clip(g0, 1e-3, 1 - 1e-3) + p_hat, lam = float(D.mean()), float(T.mean()) + g2 = chang_rcs_lambda_slope(y, D, T, ell20, ps, p_hat, lam) + h = 1e-6 + up = chang_rcs_score(y, D, T, ell20, ps, p_hat, lam + h).mean() + dn = chang_rcs_score(y, D, T, ell20, ps, p_hat, lam - h).mean() + np.testing.assert_allclose(g2, (up - dn) / (2 * h), rtol=1e-6, atol=1e-8) + + def test_p_derivative_identity_matches_finite_difference(self): + # d_p psi_2 = -(1/p)(psi_2 + theta) = -summand/p; check via central + # differences of mean(summand) in p_hat. + X, g0, D, T, ell20, y = _rcs_dgp(n=2_000, seed=12) + ps = np.clip(g0, 1e-3, 1 - 1e-3) + p_hat, lam = float(D.mean()), float(T.mean()) + summand = chang_rcs_score(y, D, T, ell20, ps, p_hat, lam) + expected = -float(summand.mean()) / p_hat + h = 1e-6 + up = chang_rcs_score(y, D, T, ell20, ps, p_hat + h, lam).mean() + dn = chang_rcs_score(y, D, T, ell20, ps, p_hat - h, lam).mean() + np.testing.assert_allclose(expected, (up - dn) / (2 * h), rtol=1e-5, atol=1e-8) + + def test_augmented_score_is_mean_centered_at_theta_hat(self): + # mean(psi_bar) = theta - theta*mean(D)/p + G2*(mean(T)-lam) = 0 + # exactly when p_hat = mean(D) and lam_hat = mean(T). + X, g0, D, T, ell20, y = _rcs_dgp(n=2_000, seed=14) + ps = np.clip(g0, 1e-3, 1 - 1e-3) + p_hat, lam = float(D.mean()), float(T.mean()) + summand = chang_rcs_score(y, D, T, ell20, ps, p_hat, lam) + theta = float(summand.mean()) + psi_bar = chang_rcs_score_augmented(summand, D, T, y, ell20, ps, theta, p_hat, lam) + assert abs(psi_bar.mean()) < 1e-12 + + def test_lambda_correction_changes_the_variance(self): + # The Thm 2 SE with the lambda term differs from the bare + # (summand - D*theta/p) variance — the "plausible implementation bug" + # regression at the score level. + X, g0, D, T, ell20, y = _rcs_dgp(n=5_000, seed=16) + ps = np.clip(g0, 1e-3, 1 - 1e-3) + p_hat, lam = float(D.mean()), float(T.mean()) + summand = chang_rcs_score(y, D, T, ell20, ps, p_hat, lam) + theta = float(summand.mean()) + psi_bar = chang_rcs_score_augmented(summand, D, T, y, ell20, ps, theta, p_hat, lam) + psi_no_lambda = summand - D * theta / p_hat + se_full = np.sqrt(np.mean(psi_bar**2) / len(D)) + se_no_lambda = np.sqrt(np.mean(psi_no_lambda**2) / len(D)) + assert abs(se_full - se_no_lambda) / se_full > 1e-4 diff --git a/tests/test_profile_panel.py b/tests/test_profile_panel.py index c33416555..6c76a9402 100644 --- a/tests/test_profile_panel.py +++ b/tests/test_profile_panel.py @@ -791,7 +791,8 @@ def test_guide_api_strings_resolve_against_public_api(): # Repeated-cross-section (§4.10) must not claim broad # applicability. The documented RCS-capable estimators are - # CallawaySantAnna(panel=False), TripleDifference, and + # CallawaySantAnna(panel=False), DMLDiD(panel=False), + # TripleDifference, and # StaggeredTripleDifference; EfficientDiD and # HeterogeneousAdoptionDiD explicitly reject RCS per REGISTRY.md. assert "most estimators remain applicable" not in text, ( diff --git a/tests/test_variance_conventions.py b/tests/test_variance_conventions.py index 960c5673f..28ab4ac66 100644 --- a/tests/test_variance_conventions.py +++ b/tests/test_variance_conventions.py @@ -454,6 +454,31 @@ def snapshot(self): "throughout, no cluster surface)" ), ), + dict( + key="dml_did_rcs", + # Declared repeated cross sections: every ROW becomes its own + # sampling unit (row-unique IDs), turning the shared panel fixture + # into a valid RCS frame; seed=0 pins the fold draw. + fit=lambda df: diff_diff.DMLDiD(seed=0, panel=False).fit( + df.assign( + x0=(df["unit"] % 7) * 0.1 + df["time"] * 0.01, + unit=np.arange(len(df)), + ), + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + covariates=["x0"], + ), + cr1_k=(), + tail_df=(None,) * 16, + status="legitimate", + reason=( + "L3: Chang (2020) Thm 2 lambda-corrected augmented-score " + "plug-in variance (Case 2; per-observation influence function; " + "normal-theory safe_inference throughout, no cluster surface)" + ), + ), ] _FAST_KEYS = { From 88cc8b75ba49fc14551870a0418a109e24d67dc2 Mon Sep 17 00:00:00 2001 From: igerber Date: Wed, 26 Aug 2026 12:41:14 -0400 Subject: [PATCH 2/2] docs(todo): track Case 2 lambda-slope double-computation as a perf follow-up (#794 review P3) --- TODO.md | 1 + 1 file changed, 1 insertion(+) diff --git a/TODO.md b/TODO.md index fde3148c0..5c78c9a3c 100644 --- a/TODO.md +++ b/TODO.md @@ -74,6 +74,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Reuse the demeaner's factorized codes in `absorbed_fe_rank`/`absorbed_fe_cr1_k_increment` instead of re-factorizing: at 186k rows the rank helper adds ~1.9 ms per absorbed fit (7.7% of the fastest Rust-served TWFE fit) and the K_reference increment ~3.2 ms per clustered-hc1 absorbed fit (~13%; see `docs/performance-plan.md` "Component-aware absorbed-FE rank"), and the helpers and `demean_by_groups` factorize the same group columns. Threading the codes through the call sites halves the factorize work; the `connected_components` call itself is ~1.1 ms. Deliberately not done in the correctness PRs. | `diff_diff/utils.py` | #variance-inventory | Quick | Low | | `EfficientDiD` conditional path: the largest remaining O(n) stage is the sieve/nuisance construction outside the tiled pass (~9s at 10k). (The `_ridge_solve_weights` Python-prep shave landed 2026-07-07 — the `omega_stack[rest]` fancy-index copy and tail scatter are skipped when no row is zero-masked, byte-identical outputs; the `zero_mask` abs scan itself remains, needed for correctness.) | `efficient_did_covariates.py` | CS-scaling | Mid | Low | | `_rq_fit` LP assembly is dense (`A_eq = [X, I, -I]` with dense identity blocks, rebuilt per cell fit): a `scipy.sparse` construction would cut memory and likely HiGHS time for large cells / bootstrap-heavy covariate CiC/QDiD fits. CAVEAT before doing it: a different matrix representation can change HiGHS's vertex selection at degenerate/tied QR optima - end-to-end covariate goldens are tie-selection-gated (fine), but the `qr_cases` tight coefficient matches may shift to the equal-loss branch; re-run the parity suite and re-calibrate if needed. | `diff_diff/changes_in_changes.py::_rq_fit` | covariates PR | Quick | Low | +| `_compute_dml_rcs_gt` computes the Case 2 λ-slope `Ĝ₂λ` twice per successful cell — once inside `chang_rcs_score_augmented()` and again for the `g2_lambda` diagnostic — duplicating the input validation and an O(n_cell) pass (nuisance fitting still dominates). Add a private augmented-score helper that accepts a precomputed slope, or return `(psi_bar, g2_lambda)` from an internal variant, keeping the public API unchanged | `diff_diff/_dr_scores.py`, `diff_diff/dml_did.py` | #794 | Quick | Low | | Evaluate flipping `DIFF_DIFF_SOLVE_OLS_FASTPATH` default-ON after an opt-in soak (the 2026-07 certified normal-equations Cholesky fast path, both backends). A flip needs: golden/parity-suite recapture at the tol-bounded posture (fitted ~1e-8 abs / SE ~1e-6 rel — the default today is byte-pinned in several benchmark conventions), certification-rate telemetry across real workloads (any decline is silent-correct but forfeits the speedup), and the staged default-flip protocol used for `df_convention` (v4-class change). Lifecycle tracked in docs/v4-deprecations.yaml (M-008). | `diff_diff/linalg.py::_resolve_solve_ols_fastpath`, `rust/src/linalg.rs::solve_ols_chol` | CS-scaling | Mid | Low | ### Testing / docs