diff --git a/CHANGELOG.md b/CHANGELOG.md index 6681fdb9..31705d37 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -209,6 +209,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 supersedes it. ### Changed +- **`CallawaySantAnna` no-covariate `estimation_method="dr"` per-cell SE is now + influence-function-based** (`sqrt(sum(phi^2))`), matching the `reg`/`ipw` branches and R's + `DRDID::drdid_panel`. It was the last per-cell SE on the ddof=1 plug-in + `sqrt(var_t/n_t + var_c/n_c)`, which deviated by O(1/n) (and `dr` is the default method, so a + plain no-covariate fit was affected). Without covariates DR reduces to difference in means, so + the per-cell SE is now bit-identical to the `reg` path. Point estimates, aggregated SEs, and + event-study SEs are unchanged — the same influence function already fed aggregation. - **`EfficientDiD` analytical covariate path rewritten as a fused unit-tiled GEMM pass** (kernel-covariance tables with `(t_pre_j, t_pre_k)` dedup, per-group kernel matrices reused across all `(g,t)` cells, batched ridge solves replacing the per-unit SVD + pseudoinverse loop; diff --git a/TODO.md b/TODO.md index 52b16fed..71119124 100644 --- a/TODO.md +++ b/TODO.md @@ -32,7 +32,6 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo | `ContinuousDiD` CGBS-2024 extensions. (a) `covariates=` kwarg — **DONE (reg/dr)**; (b) discrete-treatment saturated regression (`treatment_type="discrete"`) — **DONE**; (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0` (`control_group="lowest_dose"`) — **DONE** (discrete + continuous mass-point, single-cohort; estimand `ATT(d)−ATT(d_L)`; see REGISTRY Note #7). Remaining (all deferred `NotImplementedError`, documented): `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate); `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF); multi-cohort **heterogeneous-support** discrete aggregation (support-aware: average each dose only over the cohorts that observe it); **multi-cohort `lowest_dose`** (within-cohort `d_L` reference + support-aware cross-cohort aggregation); and **`covariates=` × `lowest_dose`** (conditional-PT-relative-to-`d_L` estimand). Single-cohort / 2-period / shared-support multi-cohort are supported. | `continuous_did.py` | CGBS-2024 | Heavy | Low | | `ImputationDiD` LOO conservative-variance refinement (BJS 2024 Supp. Appendix A.9) — a finite-sample improvement to the auxiliary-model residuals reducing overfit of `tau_tilde_g` to `epsilon`. Asymptotic Theorem-3 variance is implemented and matches R `didimputation` (which also omits LOO by default). | `imputation.py` | imputation-validation | Mid | Low | | `TwoWayFixedEffects(vcov_type in {hc2, hc2_bm})` with replicate-weight designs raises `NotImplementedError` (`twfe.py:~233`). The replicate path re-demeans per replicate, which doesn't compose with the full-dummy HC2/HC2-BM build — a correct impl needs per-replicate full-dummy refit. Workaround: `hc1` for replicate-weight CR1. | `twfe.py::fit` | follow-up | Heavy | Low | -| `CallawaySantAnna` DR no-covariate per-cell SE still uses the ddof=1 plug-in `sqrt(var_t/n_t + var_c/n_c)` rather than the IF-based `sqrt(sum(phi^2))` the reg/ipw paths switched to (O(1/n) from R's `drdid_panel`; aggregated SEs unaffected — they already consume the IF). Switch for method-uniform per-cell/aggregated consistency. | `staggered.py::_doubly_robust` (no-cov branch) | CS-scaling | Quick | Low | | `CallawaySantAnna` no-covariate ipw treats the propensity as unconditional (no correction term); R's `did` fits an intercept-only logit whose estimation effect is identically zero in the IF, so values match — decide whether to document-only (current REGISTRY note) or mirror R's code path for structural parity. | `staggered.py::_ipw_estimation` (no-cov branch) | CS-scaling | Quick | Low | | Fold the R `did` 2.5.1 ipw aggregation yardsticks (hardcoded with provenance in `test_golden_ipw_aggregation_se_vs_r_did_251`) into `csdid_golden_values.json` on the next fixture regeneration — the generator already emits the ipw `aggte` blocks; switch the test to read the JSON. | `benchmarks/R/generate_csdid_test_values.R`, `tests/test_csdid_ported.py` | CS-scaling | Quick | Low | | TWFE's HC2/HC2-BM inline full-dummy build (`twfe.py:280-315`) duplicates the dummy-construction logic in `DifferenceInDifferences(fixed_effects=...)` (`estimators.py:478-486`). Extract a shared helper, or delegate TWFE's HC2/HC2-BM path to DiD's `fixed_effects=` branch (with TWFE-specific cluster-default threading), to reduce drift risk on FE naming / survey behavior / result-surface conventions. Substantive refactor — touches both estimators. | `twfe.py::fit`, `estimators.py::DifferenceInDifferences.fit` | follow-up | Heavy | Low | diff --git a/diff_diff/staggered.py b/diff_diff/staggered.py index 225476e7..3e353e95 100644 --- a/diff_diff/staggered.py +++ b/diff_diff/staggered.py @@ -3457,17 +3457,27 @@ def _doubly_robust( else 0.0 ) else: - att = float(np.mean(treated_change) - np.mean(control_change)) - - var_t = np.var(treated_change, ddof=1) if n_t > 1 else 0.0 - var_c = np.var(control_change, ddof=1) if n_c > 1 else 0.0 + mu_t = float(np.mean(treated_change)) + mu_c = float(np.mean(control_change)) + att = mu_t - mu_c - se = float(np.sqrt(var_t / n_t + var_c / n_c)) if (n_t > 0 and n_c > 0) else 0.0 + # Influence function for the DR estimator; without covariates DR + # reduces to difference in means, so the IF matches the vectorized + # no-covariate regression path (_compute_all_att_gt_vectorized). + inf_treated = (treated_change - mu_t) / n_t + inf_control = -(control_change - mu_c) / n_c + inf_func = np.concatenate([inf_treated, inf_control]) - # Influence function for DR estimator - inf_treated = (treated_change - np.mean(treated_change)) / n_t - inf_control = (control_change - np.mean(control_change)) / n_c - inf_func = np.concatenate([inf_treated, -inf_control]) + # SE from the same IF that feeds aggregation (DRDID + # reg_did_panel/drdid_panel convention sqrt(sum(phi^2))). The + # prior ddof=1 plug-in sqrt(var_t/n_t + var_c/n_c) deviated from + # R by O(1/n) and from the reg/ipw/DR-covariate branches; the + # aggregation and bootstrap already consumed this same IF. + se = ( + float(np.sqrt(np.sum(inf_treated**2) + np.sum(inf_control**2))) + if (n_t > 0 and n_c > 0) + else 0.0 + ) return att, se, inf_func diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 75eb1bab..55f8e378 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -540,8 +540,13 @@ Aggregations: unchanged. Residual deviations kept: propensity scores are CLIPPED at `pscore_trim` (R drops at `trim.level=0.995`; differs only at extreme propensities); no-covariate ipw is treated as unconditional (R fits an intercept-only logit whose - estimation effect is identically zero in the IF, so this is presentation-only); - DR's no-covariate per-cell SE keeps its ddof=1 plug-in (O(1/n) from R; TODO row). + estimation effect is identically zero in the IF, so this is presentation-only). + DR's no-covariate per-cell SE now also uses the IF-based `sqrt(sum(phi^2))` form (it + had lagged on the ddof=1 plug-in `sqrt(var_t/n_t + var_c/n_c)`, O(1/n) from R): without + covariates DR reduces to difference in means, so its per-cell SE is now bit-identical to + the no-covariate reg path and matches R's analytical SE — point estimates and aggregated + SEs are unchanged, since the same IF already fed aggregation + (`tests/test_methodology_callaway.py::TestDRNoCovariateSEUniformity`). Side effect: reg/ipw fits with collinear covariates now route their IF breads through the rank-guarded inverse and fire the same aggregate warning as dr. Rank-0 semantics of the reg CENTERED Gram: because the intercept direction is diff --git a/tests/test_methodology_callaway.py b/tests/test_methodology_callaway.py index 64e5716d..a6c6f323 100644 --- a/tests/test_methodology_callaway.py +++ b/tests/test_methodology_callaway.py @@ -331,6 +331,98 @@ def test_estimation_methods_produce_similar_results(self): f"Estimation methods differ by {max_diff}: reg={atts[0]}, ipw={atts[1]}, dr={atts[2]}" +class TestDRNoCovariateSEUniformity: + """Without covariates, doubly-robust reduces to difference in means, so its + per-(g,t) SE must be the same influence-function form (``sqrt(sum(phi^2))``) + used by the reg/ipw branches and R's ``DRDID::drdid_panel`` — not the ddof=1 + plug-in ``sqrt(var_t/n_t + var_c/n_c)`` it historically used (O(1/n) from R). + Point estimates and aggregated SEs are unchanged because the same IF already + fed aggregation.""" + + @staticmethod + def _fit(data, method, aggregate=None): + cs = CallawaySantAnna(estimation_method=method, n_bootstrap=0) + return cs.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + aggregate=aggregate, + ) + + def test_dr_reg_per_cell_se_identical(self): + """dr and reg produce bit-identical per-cell effect AND SE without covariates.""" + data = generate_staggered_data(n_units=120, n_periods=6, never_treated_frac=0.3, seed=7) + reg = self._fit(data, "reg") + dr = self._fit(data, "dr") + + assert set(reg.group_time_effects) == set(dr.group_time_effects) + compared = 0 + for key, reg_cell in reg.group_time_effects.items(): + dr_cell = dr.group_time_effects[key] + if not np.isfinite(reg_cell["se"]): + assert not np.isfinite(dr_cell["se"]) + continue + # Both point AND SE match to machine precision (pre-fix the SE gapped + # ~1.3% via the ddof=1 plug-in; the effect always matched). + np.testing.assert_allclose(dr_cell["effect"], reg_cell["effect"], rtol=0, atol=1e-12) + np.testing.assert_allclose(dr_cell["se"], reg_cell["se"], rtol=0, atol=1e-12) + compared += 1 + assert compared >= 3, f"expected several finite cells to compare, got {compared}" + + def test_dr_no_cov_per_cell_se_hand_calc(self): + """Per-cell dr SE equals the IF form sqrt(sum(phi^2)) and is strictly + tighter than the old ddof=1 plug-in it replaced.""" + # 2-period panel: base = period 1, post = period 2; cohort g=2 vs + # never-treated (first_treat=0). One estimated cell: (g=2, t=2). + treated_y = {1: (1.0, 3.0), 2: (2.0, 5.0), 3: (0.0, 1.0), 4: (1.0, 5.0)} + control_y = {5: (0.0, 0.5), 6: (1.0, 1.0), 7: (2.0, 3.0), 8: (0.0, 1.5)} + rows = [] + for u, (y1, y2) in treated_y.items(): + rows += [(u, 1, y1, 2), (u, 2, y2, 2)] + for u, (y1, y2) in control_y.items(): + rows += [(u, 1, y1, 0), (u, 2, y2, 0)] + data = pd.DataFrame(rows, columns=["unit", "period", "outcome", "first_treat"]) + + res = self._fit(data, "dr") + cell = res.group_time_effects[(2, 2)] + + tc = np.array([y2 - y1 for (y1, y2) in treated_y.values()]) + cc = np.array([y2 - y1 for (y1, y2) in control_y.values()]) + n_t, n_c = len(tc), len(cc) + exp_att = tc.mean() - cc.mean() + exp_se = np.sqrt( + np.sum((tc - tc.mean()) ** 2) / n_t**2 + np.sum((cc - cc.mean()) ** 2) / n_c**2 + ) + old_plugin = np.sqrt(np.var(tc, ddof=1) / n_t + np.var(cc, ddof=1) / n_c) + + assert cell["effect"] == pytest.approx(exp_att, abs=1e-12) + assert cell["se"] == pytest.approx(exp_se, abs=1e-12) + # Direction sentinel: the IF-based SE is strictly smaller than the plug-in + # the fix removed (reintroducing the plug-in would trip this). + assert cell["se"] < old_plugin + + def test_dr_no_cov_aggregated_se_matches_reg(self): + """Overall and event-study SEs are identical between dr and reg (both + consume the same per-cell IF) — guards against an accidental IF change.""" + data = generate_staggered_data(n_units=120, n_periods=6, never_treated_frac=0.3, seed=11) + reg = self._fit(data, "reg", aggregate="event_study") + dr = self._fit(data, "dr", aggregate="event_study") + + np.testing.assert_allclose(dr.overall_att, reg.overall_att, rtol=0, atol=1e-12) + np.testing.assert_allclose(dr.overall_se, reg.overall_se, rtol=0, atol=1e-12) + + assert reg.event_study_effects is not None and dr.event_study_effects is not None + assert set(reg.event_study_effects) == set(dr.event_study_effects) + for e, reg_es in reg.event_study_effects.items(): + dr_es = dr.event_study_effects[e] + if np.isnan(reg_es["se"]): + assert np.isnan(dr_es["se"]) + else: + np.testing.assert_allclose(dr_es["se"], reg_es["se"], rtol=0, atol=1e-12) + + # ============================================================================= # Phase 2: R Benchmark Comparison Tests # =============================================================================