fix(callaway): IF-based per-cell SE for no-covariate doubly-robust

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;

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_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
 

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

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
 # =============================================================================