From ff3761c5dedb62694c7829d3c7e8253c54c7ae4d Mon Sep 17 00:00:00 2001 From: igerber Date: Sun, 5 Jul 2026 10:00:51 -0400 Subject: [PATCH 1/3] fix(staggered): CS panel reg/ipw influence functions + per-cell SEs to DRDID parity Two defects fixed (point estimates unchanged - the omitted terms are mean-zero, so ATTs always matched R at ~1e-11 while SEs drifted): 1. reg+covariates IF omitted DRDID::reg_did_panel's OLS estimation-effect term (asy.lin.rep.ols %*% M1): per-cell SEs 4-13% and aggregated SEs 3-20% from R, anti-conservative on some designs. The projection is evaluated in the centered basis (1/n_c + (x-xbar_c)'G^{-1}(xbar_t-xbar_c)), offset-invariant by construction, with the bread hoisted per control group on the vectorized path. 2. unweighted ipw+covariates per-cell SE used a weighted POPULATION variance never scaled by an effective sample size (~7x inflated), and the IF lacked std_ipw_did_panel's PS estimation-effect correction (aggregated ~2.4% off). Mirrors the already-correct survey branch, with the correction fused as (D-ps)*(X @ (H^-1 M2)). All panel reg/ipw per-cell SEs (incl. no-covariate) now derive from the same IF that feeds aggregation: sqrt(sum(phi^2)) (dr convention). Downstream consumers inherit: t/p/CI, aggregated SEs, event_study_vcov, cluster-override CR1, multiplier bootstrap. Validation: R-free numpy reconstructions of both DRDID IFs (incl. survey, not_yet_treated, golden-cell variants); golden per-cell + aggregated reg/ipw SE assertions enabled (~5e-12 observed); fresh R did 2.5.1 ipw aggregation yardsticks (~1e-10); behavioral tests for NaN-bread, underdetermined cells, ps fallback, universal base period, rank-guard warning (now fires for reg/ipw), uniform-weights producer equivalence. Perf guard at 100k units x 20 periods: reg within noise (+2.8/4.6%), ipw +2.8-9.5% (intrinsic per-cell PS Hessian of the new term). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01X4AzrFUMqJxcUumSH31mSr --- CHANGELOG.md | 26 ++ TODO.md | 5 +- benchmarks/R/generate_csdid_test_values.R | 14 +- diff_diff/staggered.py | 303 +++++++++++++------ docs/methodology/REGISTRY.md | 33 +- docs/survey-roadmap.md | 7 +- tests/test_csdid_ported.py | 209 ++++++++++++- tests/test_methodology_callaway.py | 351 +++++++++++++++++++++- tests/test_se_accuracy.py | 4 +- tests/test_staggered.py | 133 +++++++- tests/test_survey_r_crossvalidation.py | 7 +- 11 files changed, 970 insertions(+), 122 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index be2c23e6c..cb1a8d22d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -30,6 +30,32 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 `survey_design=SurveyDesign(psu=)`. No behavior change for unclustered fits. ### Fixed +- **`CallawaySantAnna` panel reg/ipw standard errors now match `DRDID` / + R `did` exactly** (point estimates unchanged — the omitted terms are mean-zero, + which is why ATTs always matched R at ~1e-11 while SEs drifted). Two defects: + (1) the covariate-reg influence function omitted `DRDID::reg_did_panel`'s OLS + estimation-effect term (`asy.lin.rep.ols %*% M1`) — per-cell SEs were 4-13% and + aggregated (simple/group/event-study) SEs 3-20% from R, and **anti-conservative** + on some designs (0.958x on the two-period golden fixture) despite the prior + documentation calling the plug-in "conservative"; (2) the unweighted covariate-ipw + per-cell SE used a weighted *population* variance never scaled by an effective + sample size — **~7x inflated** — and its influence function lacked + `std_ipw_did_panel`'s propensity-score estimation-effect correction + (`asy.lin.rep.ps %*% M2`; aggregated SEs ~2.4% off). The survey ipw branch already + carried both corrections (Phase 7a); the fix mirrors it method-uniformly. All panel + reg/ipw per-cell SEs (including no-covariate, where ipw's + `var_c*(1-p)/(n_c*p)` plug-in equaled R only at treated-share 0.5) are now derived + from the same influence function that feeds aggregation: `sqrt(sum(phi^2))`, the + convention dr always used. Everything downstream of the per-cell IFs inherits the + fix: t-stats/p-values/confidence intervals, aggregated SEs, `event_study_vcov` + (HonestDiD input), the bare-`cluster=` CR1 per-cell override, and + multiplier-bootstrap SEs. The reg estimation-effect projection is evaluated in the + centered basis (`1/n_c + (x-x̄_c)'G^{-1}(x̄_t-x̄_c)`), which is offset-invariant; + reg/ipw fits with collinear covariates now fire the same aggregate rank-guard + warning as dr. Post-fix parity: per-cell and aggregated reg/ipw vs the R golden + fixtures at ~5e-12, and vs fresh R `did` 2.5.1 ipw aggregations at ~1e-10 + (previously-unasserted golden SE blocks are now enabled, plus R-free numpy + reconstructions of both DRDID influence functions). - **ImputationDiD/TwoStageDiD covariate fits with zero-weight replicate designs (JK1/plain BRR) now produce finite SEs.** Replicate weights that zero out whole PSUs reach Step 1 unmasked; the previous per-estimator pandas demeaning loops divided 0/0 on zero-total-weight diff --git a/TODO.md b/TODO.md index 12c6528bf..fa6ac6489 100644 --- a/TODO.md +++ b/TODO.md @@ -32,7 +32,9 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo | `ContinuousDiD` CGBS-2024 extensions. (a) `covariates=` kwarg — **DONE (reg/dr)**; remaining: (b) discrete-treatment saturated regression (integer dose currently warned, not routed to per-level coefficients); (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0`. Also deferred from the covariate work: `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate — documented `NotImplementedError`), and `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF). | `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` reg-method aggregated SEs sit 3-20% from the R golden fixtures (`two_period` simple/group/dynamic, `dynamic_effects` dynamic) while the ATTs match at 1e-11 and the dr-method aggregated SEs match at ≤7e-6. Pre-existing (byte-identical deltas on the pre-fast-path tree); golden aggregated-SE assertions were therefore enabled for dr scenarios only (`test_golden_simple_aggregation_se`, `test_golden_event_study_aggregation_se`). Investigate the reg-path per-cell IF vs R `DRDID::reg_did_panel`'s. | `staggered.py::_outcome_regression`, `tests/test_csdid_ported.py` | CS-scaling | Mid | Medium | +| `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 | ### Performance @@ -65,6 +67,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Render `docs/methodology/REPORTING.md` and `REGISTRY.md` as in-site Sphinx pages so cross-refs can use `:doc:` instead of off-site `blob/main` URLs (stable-docs readers can otherwise land on a different revision than their package version). Two paths: (a) add `myst-parser` to `conf.py` + docs extras and link with `:doc:`, or (b) convert both to `.rst`. **Note:** REGISTRY.md is ~4.5k lines of LaTeX-heavy markdown — high risk under the `-W` (warnings-as-errors) Sphinx build; budget multiple rounds. | `docs/conf.py`, `docs/api/business_report.rst`, `docs/api/diagnostic_report.rst`, tutorials 18 & 19 | follow-up | Mid | Low | | `ImputationDiD` covariate-path variance lacks a dedicated parity anchor — only the no-covariate staggered panel is R-parity'd, though the covariate path shares the same validated projection code. Add a small dense-design **hand-calc** for the covariate projection (no external tooling), or a covariate (time-varying X) R `didimputation` golden asserting overall/ES SE parity (the golden variant needs local R). | `tests/test_methodology_imputation.py`, `benchmarks/R/generate_didimputation_golden.R` | imputation-validation | Mid | Low | | Add true half-sample BRR replicate-weight regressions per estimator family (current tests use Fay-like 0.5/1.5 perturbations; `test_survey_phase6.py` covers true BRR at the helper level). | `tests/test_replicate_weight_expansion.py` | #253 | Mid | Low | +| `test_bootstrap_se_matches_pre_pr4_baseline` (dCDH bit-identity guard) fails under FULL-SUITE order but passes standalone/as-a-file: its `pure_python` backend detection reads a state some earlier test perturbs, so it compares against the wrong per-backend baseline arm (observed failure message says `backend=pure-python` while the run produces the other arm's value). Reproduced identically on pristine `origin/main` (2026-07-05, `b56931a6`), so it pre-dates the CS reg/ipw IF fix. Make the backend detection order-robust (e.g. re-resolve from `diff_diff._backend` at call time, or isolate via monkeypatch). | `tests/test_survey_dcdh.py::TestBootstrapCellPeriod` | CS-scaling | Quick | Low | | Port the CI `` extraction into the reviewer-eval harness so `docs/tutorials/*.ipynb` cases (currently guarded out of `verify-corpus`/`run`) can be reviewed with CI-equivalent context. | `tools/reviewer-eval/adapters/ci_prompt.py` | local-review | Mid | Low | --- diff --git a/benchmarks/R/generate_csdid_test_values.R b/benchmarks/R/generate_csdid_test_values.R index 1c5e3c7a9..b25594a81 100644 --- a/benchmarks/R/generate_csdid_test_values.R +++ b/benchmarks/R/generate_csdid_test_values.R @@ -108,6 +108,13 @@ data2_ipw <- build_sim_dataset(sp2_ipw) res2_ipw <- att_gt(yname = "Y", xformla = ~X, data = data2_ipw, tname = "period", idname = "id", gname = "G", est_method = "ipw", bstrap = FALSE, cband = FALSE) +# Aggregations for the ipw scenario. Until the fixture is regenerated, these +# values are hardcoded (from did 2.5.1 / DRDID 1.3.0, 2026-07-05) in +# tests/test_csdid_ported.py::test_golden_ipw_aggregation_se_vs_r_did_251; +# on the next regeneration the test can switch to reading them from the JSON. +agg2_ipw_simple <- aggte(res2_ipw, type = "simple", bstrap = FALSE, cband = FALSE) +agg2_ipw_dynamic <- aggte(res2_ipw, type = "dynamic", bstrap = FALSE, cband = FALSE) +agg2_ipw_group <- aggte(res2_ipw, type = "group", bstrap = FALSE, cband = FALSE) scenarios$with_covariates_dr <- list( data = export_data(data2), @@ -125,7 +132,12 @@ scenarios$with_covariates_ipw <- list( data = export_data(data2_ipw), params = list(est_method = "ipw", control_group = "nevertreated", xformla = "~X", base_period = "varying"), - results = list(group_time = extract_gt(res2_ipw)) + results = list( + group_time = extract_gt(res2_ipw), + simple = list(att = agg2_ipw_simple$overall.att, se = agg2_ipw_simple$overall.se), + group = extract_agg(agg2_ipw_group), + dynamic = extract_agg(agg2_ipw_dynamic) + ) ) # --------------------------------------------------------------------------- diff --git a/diff_diff/staggered.py b/diff_diff/staggered.py index d739d8cc9..718ef1671 100644 --- a/diff_diff/staggered.py +++ b/diff_diff/staggered.py @@ -1163,15 +1163,21 @@ def _compute_all_att_gt_vectorized( ) sw_sum = float(np.sum(sw_t)) else: - att = float(np.mean(treated_change) - np.mean(control_change)) - - var_t = float(np.var(treated_change, ddof=1)) if n_t > 1 else 0.0 - var_c = float(np.var(control_change, ddof=1)) if n_c > 1 else 0.0 - se = float(np.sqrt(var_t / n_t + var_c / n_c)) if (n_t > 0 and n_c > 0) else 0.0 + mu_t = float(np.mean(treated_change)) + mu_c = float(np.mean(control_change)) + att = mu_t - mu_c # Influence function - inf_treated = (treated_change - np.mean(treated_change)) / n_t - inf_control = -(control_change - np.mean(control_change)) / n_c + inf_treated = (treated_change - mu_t) / n_t + inf_control = -(control_change - mu_c) / n_c + # SE from the same IF that feeds aggregation (DRDID convention + # sqrt(sum(phi^2)); the prior ddof=1 plug-in deviated from R's + # reg_did_panel by O(1/n)). + 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 + ) sw_sum = None # A difference-in-means ATT is finite given n_t, n_c > 0, but enforce @@ -1393,6 +1399,9 @@ def _compute_all_att_gt_covariate_reg( continue X_ctrl = None + group_bread = None + group_xbar_c = None + group_Xc_centered = None if not ctrl_has_nan: X_ctrl = np.column_stack([np.ones(n_c_base), X_ctrl_raw]) @@ -1414,6 +1423,25 @@ def _compute_all_att_gt_covariate_reg( stacklevel=2, ) + # The IF bread is shared by every task in this group when the + # control design is; hoist it (rank-guarded, so collinear + # covariates get the column-drop generalized inverse + # consistent with the DR/IPW per-cell breads). It is computed + # in the CENTERED basis — the projection + # x'(X'X)^{-1}x̄_t == 1/n_c + (x - x̄_c)'G^{-1}(x̄_t - x̄_c) + # (G = centered control Gram) is the same quantity evaluated + # offset-invariantly: the raw Gram squares the design's + # conditioning, and a large constant covariate offset would + # push the equilibrated-Gram rank check below its threshold, + # silently truncating a genuine direction. + if is_balanced and self.control_group == "never_treated": + group_xbar_c = X_ctrl_raw.mean(axis=0) + group_Xc_centered = X_ctrl_raw - group_xbar_c + group_bread = _safe_inv( + group_Xc_centered.T @ group_Xc_centered, + tracker=self._safe_inv_tracker, + ) + # Process each (g, t) pair in this group for g, t, bp_val, base_col, post_col in tasks: treated_mask = cohort_masks[g] @@ -1467,12 +1495,14 @@ def _compute_all_att_gt_covariate_reg( UserWarning, stacklevel=3, ) - att = float(np.mean(treated_change) - np.mean(control_change)) - var_t = float(np.var(treated_change, ddof=1)) if n_t > 1 else 0.0 - var_c = float(np.var(control_change, ddof=1)) if n_c > 1 else 0.0 - se = float(np.sqrt(var_t / n_t + var_c / n_c)) - inf_treated = (treated_change - np.mean(treated_change)) / n_t - inf_control = -(control_change - np.mean(control_change)) / n_c + mu_t_fb = float(np.mean(treated_change)) + mu_c_fb = float(np.mean(control_change)) + att = mu_t_fb - mu_c_fb + inf_treated = (treated_change - mu_t_fb) / n_t + inf_control = -(control_change - mu_c_fb) / n_c + # SE from the same IF that feeds aggregation (DRDID + # convention sqrt(sum(phi^2)); prior ddof=1 plug-in). + se = float(np.sqrt(np.sum(inf_treated**2) + np.sum(inf_control**2))) else: # Build per-pair X_ctrl if control_valid differs from base if is_balanced and self.control_group == "never_treated" and X_ctrl is not None: @@ -1550,11 +1580,50 @@ def _compute_all_att_gt_covariate_reg( if nan_cell: att = np.nan else: - var_t = float(np.var(treated_residuals, ddof=1)) if n_t > 1 else 0.0 - var_c = float(np.var(residuals, ddof=1)) if pair_n_c > 1 else 0.0 - se = float(np.sqrt(var_t / n_t + var_c / pair_n_c)) - inf_treated = (treated_residuals - np.mean(treated_residuals)) / n_t - inf_control = -residuals / pair_n_c + # DRDID::reg_did_panel influence function. The control + # leg carries the OLS estimation-effect term: the + # control residual projected through the rank-guarded + # bread and evaluated at the treated covariate mean + # (phi convention — every R 1/n factor cancels). + # Evaluated in the CENTERED basis, + # 1/n_c + (x - x̄_c)'G^{-1}(x̄_t - x̄_c), which is the + # same quantity as x'(X'X)^{-1}x̄_t but offset- + # invariant (see the group_bread hoist above); the + # no-covariate collapse -residuals/n_c is the leading + # 1/n_c term. Omitting this term left ATTs exact (it + # is mean-zero by the OLS normal equations) while SEs + # sat 4-13% from R per-cell. + if ( + is_balanced + and self.control_group == "never_treated" + and group_bread is not None + ): + xbar_c = group_xbar_c + Xc_centered = group_Xc_centered + bread = group_bread + else: + xbar_c = X_control_pair.mean(axis=0) + Xc_centered = X_control_pair - xbar_c + bread = _safe_inv( + Xc_centered.T @ Xc_centered, + tracker=self._safe_inv_tracker, + ) + d_tc = X_treated_pair.mean(axis=0) - xbar_c + with np.errstate(all="ignore"): + proj_c = 1.0 / pair_n_c + Xc_centered @ (bread @ d_tc) + inf_treated = (treated_residuals - att) / n_t + inf_control = -residuals * proj_c + # Per-cell SE from the same IF that feeds aggregation + # (sqrt(sum(phi^2)), the DR convention). A rank-0 + # bread yields NaN — safe_inference then NaNs the + # full inference tuple (never a silent 0.0). + var_psi = float(np.sum(inf_treated**2) + np.sum(inf_control**2)) + if not np.isfinite(var_psi): + se = float("nan") + elif var_psi > 0: + se = float(np.sqrt(var_psi)) + else: + se = 0.0 # Non-estimable (non-finite regression) cell: materialize NaN with # NO influence-function entry — uniform with the missing-period / @@ -2578,7 +2647,9 @@ def _outcome_regression( if X_treated is not None and X_control is not None and X_treated.shape[1] > 0: # Covariate-adjusted outcome regression # Fit regression on control units: E[Delta Y | X, D=0] - beta, residuals = _linear_regression( + # (residuals are recomputed below from the NaN-zeroed beta so the + # IF uses the same prediction convention as the ATT). + beta, _ = _linear_regression( X_control, control_change, rank_deficient_action=self.rank_deficient_action, @@ -2596,40 +2667,61 @@ def _outcome_regression( # ATT: survey-weighted mean of treated residuals treated_residuals = treated_change - predicted_control + # DRDID::reg_did_panel influence function (survey weights map to + # R's i.weights; only weight RATIOS enter every term, so the + # normalization constant cancels and the raw sw arrays can be + # used directly). The control leg carries the OLS estimation- + # effect term: the (weighted) control residual projected through + # the rank-guarded WLS bread and evaluated at the (weighted) + # treated covariate mean. Evaluated in the CENTERED basis, + # 1/sum(W_c) + (x - x̄_c)'G^{-1}(x̄_t - x̄_c) with G the centered + # weighted control Gram — the same quantity as x'(X'WX)^{-1}x̄_t + # but offset-invariant (the raw Gram squares the design's + # conditioning; a large constant covariate offset would push the + # equilibrated-Gram rank check below threshold and silently + # truncate a genuine direction). The no-covariate collapse is the + # leading -resid_c/n_c term. Omitting the estimation-effect term + # left ATTs exact (it is mean-zero by the WLS normal equations) + # while per-cell SEs sat 4-13% from R. phi convention: psi_R / n. + X_c_int = np.column_stack([np.ones(n_c), X_control]) + with np.errstate(all="ignore"): + resid_c = control_change - np.dot(X_c_int, beta) + if sw_treated is not None: sw_t_sum = float(np.sum(sw_treated)) - sw_c_sum = float(np.sum(sw_control)) sw_t_norm = sw_treated / sw_t_sum - sw_c_norm = sw_control / sw_c_sum att = float(np.sum(sw_t_norm * treated_residuals)) - - # Survey-weighted OR influence function. - # Mirrors unweighted: inf_treated = (resid-ATT)/n_t, - # inf_control = -resid/n_c. Survey: w_i/sum(w_group). - # WLS residuals are orthogonal to W*X by construction. - X_c_int = np.column_stack([np.ones(n_c), X_control]) - resid_c = control_change - np.dot(X_c_int, beta) - + W_c = sw_control + xbar_t = np.sum(sw_treated[:, None] * X_treated, axis=0) / sw_t_sum inf_treated = sw_t_norm * (treated_residuals - att) - inf_control = -sw_c_norm * resid_c - inf_func = np.concatenate([inf_treated, inf_control]) - - # SE: survey-weighted variance matching unweighted var_t/n_t + var_c/n_c - var_t = float(np.sum(sw_t_norm * (treated_residuals - att) ** 2)) - var_c = float(np.sum(sw_c_norm * resid_c**2)) - se = float(np.sqrt(var_t + var_c)) if (n_t > 0 and n_c > 0) else 0.0 else: att = float(np.mean(treated_residuals)) - - # Standard error using sandwich estimator - var_t = np.var(treated_residuals, ddof=1) if n_t > 1 else 0.0 - var_c = np.var(residuals, ddof=1) if n_c > 1 else 0.0 - 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 - inf_treated = (treated_residuals - np.mean(treated_residuals)) / n_t - inf_control = -residuals / n_c - inf_func = np.concatenate([inf_treated, inf_control]) + W_c = np.ones(n_c) + xbar_t = X_treated.mean(axis=0) + inf_treated = (treated_residuals - att) / n_t + + w_c_sum = float(np.sum(W_c)) + with np.errstate(all="ignore"): + xbar_c = np.sum(W_c[:, None] * X_control, axis=0) / w_c_sum + Xc_centered = X_control - xbar_c + bread = _safe_inv( + Xc_centered.T @ (W_c[:, None] * Xc_centered), + tracker=self._safe_inv_tracker, + ) + proj_c = 1.0 / w_c_sum + Xc_centered @ (bread @ (xbar_t - xbar_c)) + inf_control = -(W_c * resid_c) * proj_c + inf_func = np.concatenate([inf_treated, inf_control]) + + # Per-cell SE from the same IF that feeds aggregation + # (sqrt(sum(phi^2)), the DR convention). A rank-0 bread yields + # NaN — safe_inference then NaNs the full inference tuple. + var_psi = float(np.sum(inf_func**2)) + if not np.isfinite(var_psi): + se = float("nan") + elif var_psi > 0: + se = float(np.sqrt(var_psi)) + else: + se = 0.0 else: # Simple difference in means (no covariates) if sw_treated is not None: @@ -2651,16 +2743,19 @@ def _outcome_regression( 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 - se = float(np.sqrt(var_t / n_t + var_c / n_c)) if (n_t > 0 and n_c > 0) else 0.0 + mu_t = float(np.mean(treated_change)) + mu_c = float(np.mean(control_change)) + att = mu_t - mu_c # Influence function (for aggregation) - inf_treated = treated_change - np.mean(treated_change) - inf_control = control_change - np.mean(control_change) - inf_func = np.concatenate([inf_treated / n_t, -inf_control / n_c]) + inf_treated = (treated_change - mu_t) / n_t + inf_control = -(control_change - mu_c) / n_c + inf_func = np.concatenate([inf_treated, inf_control]) + + # SE from the same IF that feeds aggregation (DRDID convention + # sqrt(sum(phi^2)); the prior ddof=1 plug-in deviated from R's + # reg_did_panel by O(1/n)). + se = float(np.sqrt(np.sum(inf_func**2))) if (n_t > 0 and n_c > 0) else 0.0 return att, se, inf_func @@ -2698,7 +2793,6 @@ def _ipw_estimation( """ n_t = len(treated_change) n_c = len(control_change) - n_total = n_treated + n_control if X_treated is not None and X_control is not None and X_treated.shape[1] > 0: # Covariate-adjusted IPW estimation @@ -2843,27 +2937,67 @@ def _ipw_estimation( # IPW weights for control units: p(X) / (1 - p(X)) # This reweights controls to have same covariate distribution as treated weights_control = pscore_control / (1 - pscore_control) - weights_control = weights_control / np.sum(weights_control) # normalize + weights_control_norm = weights_control / np.sum(weights_control) # ATT = mean(treated) - weighted_mean(control) - att = float(np.mean(treated_change) - np.sum(weights_control * control_change)) + mu_t = float(np.mean(treated_change)) + eta_c = float(np.sum(weights_control_norm * control_change)) + att = mu_t - eta_c + + # Influence function (Hajek leading terms; mirrors the survey + # branch above with unit weights) + inf_treated = (treated_change - mu_t) / n_t + inf_control = -weights_control_norm * (control_change - eta_c) + inf_func = np.concatenate([inf_treated, inf_control]) - # Compute standard error - var_t = np.var(treated_change, ddof=1) if n_t > 1 else 0.0 + if not ps_fallback_used: + # Propensity-score estimation-effect correction + # (DRDID::std_ipw_did_panel's asy.lin.rep.ps %*% M2) — + # mirrors the survey branch above with sw=None. Its + # omission left ATTs exact (the logit score is mean-zero + # at the MLE) while aggregated SEs sat ~2.4% from R. + n_all_panel = n_t + n_c + X_all_int = np.column_stack([np.ones(n_all_panel), X_all]) + pscore_all = np.concatenate([pscore_treated, pscore_control]) - weighted_var_c = np.sum( - weights_control - * (control_change - np.sum(weights_control * control_change)) ** 2 - ) + W_ps = pscore_all * (1 - pscore_all) + # R: Hessian.ps = crossprod(X * sqrt(W)) / n + H_psi = X_all_int.T @ (W_ps[:, None] * X_all_int) / n_all_panel + H_psi_inv = _safe_inv(H_psi, tracker=self._safe_inv_tracker) - se = float(np.sqrt(var_t / n_t + weighted_var_c)) if (n_t > 0 and n_c > 0) else 0.0 + D_all = np.concatenate([np.ones(n_t), np.zeros(n_c)]) - # Influence function - inf_treated = (treated_change - np.mean(treated_change)) / n_t - inf_control = -weights_control * ( - control_change - np.sum(weights_control * control_change) - ) - inf_func = np.concatenate([inf_treated, inf_control]) + # R: M2 = colMeans(w.cont * (y - att) * X) / mean(w.cont); + # the normalized-weight subset sum matches R's full-sample + # colMeans/mean(w) after cancellation. + M2 = np.sum( + (weights_control_norm * (control_change - eta_c))[:, None] + * X_all_int[n_t:], + axis=0, + ) + + # R: asy.lin.rep.ps %*% M2 with + # asy.lin.rep.ps = score.ps %*% Hessian.ps — evaluated + # fused as (D - ps) * (X @ (H^{-1} M2)), which avoids + # materializing the (n x k) asy.lin.rep matrix (same + # algebra, one O(n*k) matvec instead of an O(n*k^2) + # matmul). psi-scale, converted to phi via /n_all_panel. + # Subtract: R adds it to inf.control, and att = treat - control. + corr = (D_all - pscore_all) * (X_all_int @ (H_psi_inv @ M2)) + inf_func = inf_func - corr / n_all_panel + + # SE from the same IF that feeds aggregation. The prior + # plug-in sqrt(var_t/n_t + weighted_var_c) used a weighted + # POPULATION variance for the control leg — never scaled by + # an effective sample size — inflating per-cell SEs ~7x + # versus R's std_ipw_did_panel. + var_psi = np.sum(inf_func**2) + if not np.isfinite(var_psi): + se = float("nan") + elif var_psi > 0: + se = float(np.sqrt(var_psi)) + else: + se = 0.0 else: # Unconditional IPW (reduces to difference in means) if sw_treated is not None: @@ -2884,24 +3018,21 @@ def _ipw_estimation( else 0.0 ) else: - p_treat = n_treated / n_total # unconditional propensity score - - 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 - - # Adjusted variance for IPW - se = float( - np.sqrt(var_t / n_t + var_c * (1 - p_treat) / (n_c * p_treat)) - if (n_t > 0 and n_c > 0 and p_treat > 0) - else 0.0 - ) + mu_t = float(np.mean(treated_change)) + mu_c = float(np.mean(control_change)) + att = mu_t - mu_c # Influence function (for aggregation) - 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]) + inf_treated = (treated_change - mu_t) / n_t + inf_control = -(control_change - mu_c) / n_c + inf_func = np.concatenate([inf_treated, inf_control]) + + # SE from the same IF that feeds aggregation. The prior + # "adjusted" plug-in var_c*(1-p_treat)/(n_c*p_treat) equals + # R's std_ipw_did_panel analytical SE only at p_treat = 0.5; + # with a constant propensity the IF sum IS R's SE (no + # estimated PS parameter, so no correction term). + se = float(np.sqrt(np.sum(inf_func**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 29ae7e3b4..3ac6e8dd8 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -460,7 +460,7 @@ methodology is introduced. - Requires never-treated units as comparison group (identified by `first_treat=0` or `never_treated=True`) - Warns if no never-treated units exist (suggests alternative comparison strategies) - Limited pre-treatment periods reduce ability to test parallel trends -- **Note (rank-guarded IF standard errors):** The analytical SE paths invert the propensity-score Hessian (`H_psi`) and outcome-regression bread (`X'WX`) across every `(g, t)` cell via `_safe_inv()`. `np.linalg.solve`/`inv` raise `LinAlgError` only on *exactly* singular matrices; a **near**-singular Gram (a constant or collinear covariate) does **not** raise, so the prior code returned a garbage inverse (entries ~1e13) that flowed straight into the SE (`overall_se` ~1e13 was reproduced for `estimation_method="dr"`). `_safe_inv()` now delegates to the shared `_rank_guarded_inv()` (`diff_diff/linalg.py`): it symmetrically equilibrates `A → D^{-1/2} A D^{-1/2}` (`D = diag(A)`), then, when rank-deficient, inverts a **column-dropped** principal submatrix — keeping the most-independent columns via pivoted QR on the **equilibrated** Gram. This is a column-drop generalized inverse in the **same family** as the point estimate / R's `lm()` (drop redundant columns, *not* a minimum-norm pseudo-inverse), but the column *selection* is computed on the equilibrated Gram and is therefore **scale-invariant** by construction — so it need not pick the *same member* of a collinear set as the point estimate's raw-pivot `_detect_rank_deficiency` under mixed-scale *exact* collinearity (see the selection caveat below); the resulting SE is well-defined and order-invariant regardless of which redundant member is dropped. The `rcond = 1e-10` relative-eigenvalue threshold on the equilibrated Gram sets the rank (chosen because a Gram squares the condition number of `X`; matches EfficientDiD's `tol/max_eigval`). The result is a **finite SE on the identified covariate subset**; an all-NaN inverse (and NaN SE) is returned **only** on true rank-0. **Column-drop equals the full-rank limit:** using a column-drop generalized inverse (rather than a minimum-norm pseudo-inverse) makes the analytical SE equal the **well-conditioned near-collinear limit** — replacing an exactly-collinear covariate with a near-collinear (full-rank) one yields the same SE to working precision (verified `se_ratio ≈ 1` across `reg`/`ipw`/`dr`), for every per-cell bread (the PS Hessian and the control- *and* treated-side OR breads). A minimum-norm pseudo-inverse would instead diverge sharply when the IF multiplier leaves `range(A)` — e.g. a control (or treated-sub-cell) bread multiplied by a mean from a cell where the covariate is *not* collinear — so it is rejected. With column-drop there is no such divergence; a covariate that is rank-deficient only within one cell still legitimately enters the other cells' full-rank fits, so the ATT and SE reflect that (poor) covariate specification, surfaced by the aggregate warning. **Selection caveat (equilibrated vs raw pivot):** because `_rank_guarded_inv` pivots the *equilibrated* Gram, the *member* it drops from a collinear set can differ from the point estimate's raw-pivot `_detect_rank_deficiency` choice **only** under mixed-scale *exact* collinearity (e.g. a covariate and its `1e8×` duplicate). This selection difference does **not** change the SE: dropping either member of an exactly-collinear pair leaves the identified subspace — hence the variance — unchanged, and the SE is order-invariant (verified for both column orders and under survey weighting). It is a documented, intentional deviation in the *generalized-inverse column choice* (equilibrated, scale-invariant) from the point-estimate's raw pivot; it is not a parity claim with the raw selection. **Caveat (point-estimate side, NOT the rank-guard):** CS/SDDD's reg/DR *point-estimate* OR fit now routes through the scale-equilibrated `solve_ols` (SVD/gelsd), so a covariate correlated with another regressor at a very large scale no longer perturbs the point-estimate ATT — or the IF SE that consistently follows it — under `reg`/`dr`; scale-invariance is pinned by tests (covariate + large constant offset → ATT(g,t) unchanged to ~1e-11). This is independent of the rank-guarded IF inverse, whose SE is invariant to which collinear member is dropped (see the selection caveat above). The well-conditioned fast path returns `np.linalg.solve(A, I)` unchanged (R-parity preserved). `fit()` emits ONE aggregate `UserWarning` (count + max condition number) reporting the dropped redundant direction(s); suppressed under `rank_deficient_action="silent"`. **`rank_deficient_action` enforcement:** `"error"` is enforced *upstream* at the point-estimate solve (`solve_ols` / `solve_logit` / the covariate-regression fit), which raises `ValueError` when the covariate **design** is rank-deficient at its `1e-7` threshold (`_detect_rank_deficiency`) **before** any influence-function SE is computed — so for a genuinely (design-)rank-deficient covariate the rank-guarded IF inverse is reached only under `"warn"`/`"silent"`. The IF guard truncates on a *stricter* relative threshold, though: `1e-10` on the **equilibrated Gram** vs `1e-7` on the design, and a Gram squares `X`'s condition number, so the IF guard drops a direction once `X`'s singular-value ratio falls below ~`1e-5` (= √`1e-10`), well above the design's `1e-7`. A cell that is *near-singular yet still full-rank by the upstream design check* (singular-value ratio between ~`1e-7` and ~`1e-5`) therefore passes the `"error"` gate without an exception and is **still** column-dropped by the IF guard — the guard does not re-raise, and (under `"error"`/`"warn"`, not `"silent"`) the aggregate warning still fires. So `"error"` blocks design-rank-deficient *covariates*; it does not promise that every near-singular IF bread/Hessian raises. Sibling of axis-A finding #17 in the Phase 2 silent-failures audit. +- **Note (rank-guarded IF standard errors):** The analytical SE paths invert the propensity-score Hessian (`H_psi`) and outcome-regression bread (`X'WX`) across every `(g, t)` cell via `_safe_inv()`. `np.linalg.solve`/`inv` raise `LinAlgError` only on *exactly* singular matrices; a **near**-singular Gram (a constant or collinear covariate) does **not** raise, so the prior code returned a garbage inverse (entries ~1e13) that flowed straight into the SE (`overall_se` ~1e13 was reproduced for `estimation_method="dr"`). `_safe_inv()` now delegates to the shared `_rank_guarded_inv()` (`diff_diff/linalg.py`): it symmetrically equilibrates `A → D^{-1/2} A D^{-1/2}` (`D = diag(A)`), then, when rank-deficient, inverts a **column-dropped** principal submatrix — keeping the most-independent columns via pivoted QR on the **equilibrated** Gram. This is a column-drop generalized inverse in the **same family** as the point estimate / R's `lm()` (drop redundant columns, *not* a minimum-norm pseudo-inverse), but the column *selection* is computed on the equilibrated Gram and is therefore **scale-invariant** by construction — so it need not pick the *same member* of a collinear set as the point estimate's raw-pivot `_detect_rank_deficiency` under mixed-scale *exact* collinearity (see the selection caveat below); the resulting SE is well-defined and order-invariant regardless of which redundant member is dropped. The `rcond = 1e-10` relative-eigenvalue threshold on the equilibrated Gram sets the rank (chosen because a Gram squares the condition number of `X`; matches EfficientDiD's `tol/max_eigval`). The result is a **finite SE on the identified covariate subset**; an all-NaN inverse (and NaN SE) is returned **only** on true rank-0. **Column-drop equals the full-rank limit:** using a column-drop generalized inverse (rather than a minimum-norm pseudo-inverse) makes the analytical SE equal the **well-conditioned near-collinear limit** — replacing an exactly-collinear covariate with a near-collinear (full-rank) one yields the same SE to working precision (verified `se_ratio ≈ 1` across `reg`/`ipw`/`dr`), for every per-cell bread (the PS Hessian — including unweighted ipw's since v3.7 — and the control- *and* treated-side OR breads, including reg's centered estimation-effect Gram `G` since v3.7). A minimum-norm pseudo-inverse would instead diverge sharply when the IF multiplier leaves `range(A)` — e.g. a control (or treated-sub-cell) bread multiplied by a mean from a cell where the covariate is *not* collinear — so it is rejected. With column-drop there is no such divergence; a covariate that is rank-deficient only within one cell still legitimately enters the other cells' full-rank fits, so the ATT and SE reflect that (poor) covariate specification, surfaced by the aggregate warning. **Selection caveat (equilibrated vs raw pivot):** because `_rank_guarded_inv` pivots the *equilibrated* Gram, the *member* it drops from a collinear set can differ from the point estimate's raw-pivot `_detect_rank_deficiency` choice **only** under mixed-scale *exact* collinearity (e.g. a covariate and its `1e8×` duplicate). This selection difference does **not** change the SE: dropping either member of an exactly-collinear pair leaves the identified subspace — hence the variance — unchanged, and the SE is order-invariant (verified for both column orders and under survey weighting). It is a documented, intentional deviation in the *generalized-inverse column choice* (equilibrated, scale-invariant) from the point-estimate's raw pivot; it is not a parity claim with the raw selection. **Caveat (point-estimate side, NOT the rank-guard):** CS/SDDD's reg/DR *point-estimate* OR fit now routes through the scale-equilibrated `solve_ols` (SVD/gelsd), so a covariate correlated with another regressor at a very large scale no longer perturbs the point-estimate ATT — or the IF SE that consistently follows it — under `reg`/`dr`; scale-invariance is pinned by tests (covariate + large constant offset → ATT(g,t) unchanged to ~1e-11). This is independent of the rank-guarded IF inverse, whose SE is invariant to which collinear member is dropped (see the selection caveat above). The well-conditioned fast path returns `np.linalg.solve(A, I)` unchanged (R-parity preserved). `fit()` emits ONE aggregate `UserWarning` (count + max condition number) reporting the dropped redundant direction(s); suppressed under `rank_deficient_action="silent"`. **`rank_deficient_action` enforcement:** `"error"` is enforced *upstream* at the point-estimate solve (`solve_ols` / `solve_logit` / the covariate-regression fit), which raises `ValueError` when the covariate **design** is rank-deficient at its `1e-7` threshold (`_detect_rank_deficiency`) **before** any influence-function SE is computed — so for a genuinely (design-)rank-deficient covariate the rank-guarded IF inverse is reached only under `"warn"`/`"silent"`. The IF guard truncates on a *stricter* relative threshold, though: `1e-10` on the **equilibrated Gram** vs `1e-7` on the design, and a Gram squares `X`'s condition number, so the IF guard drops a direction once `X`'s singular-value ratio falls below ~`1e-5` (= √`1e-10`), well above the design's `1e-7`. A cell that is *near-singular yet still full-rank by the upstream design check* (singular-value ratio between ~`1e-7` and ~`1e-5`) therefore passes the `"error"` gate without an exception and is **still** column-dropped by the IF guard — the guard does not re-raise, and (under `"error"`/`"warn"`, not `"silent"`) the aggregate warning still fires. So `"error"` blocks design-rank-deficient *covariates*; it does not promise that every near-singular IF bread/Hessian raises. Sibling of axis-A finding #17 in the Phase 2 silent-failures audit. *Variance families (`vcov_type`, IF-based):* - `hc1` (default, only accepted value) — per-unit IF variance per Callaway & Sant'Anna (2021) when `cluster=None`; cluster-robust CR1 Liang-Zeger on the IF when `cluster=X` is set (synthesizes `SurveyDesign(psu=X)` internally, threading through the same PSU machinery as explicit survey designs). When `survey_design=SurveyDesign(psu=Y)` is provided, the explicit PSU takes precedence; if `cluster=X` is also set with a different partition, emits a `UserWarning` (PSU wins). @@ -498,6 +498,31 @@ Aggregations: *Standard errors:* - Default: Analytical (influence function-based) +- **Note (v3.7 DRDID IF parity, panel reg/ipw):** Every panel reg/ipw per-cell SE is + derived from the SAME influence function that feeds aggregation, the bootstrap, the + cluster override, and `event_study_vcov`: `se = sqrt(sum(phi^2))` (DRDID convention; + dr already worked this way). Two defects were fixed to get there. (1) The reg+cov IF + omitted `DRDID::reg_did_panel`'s OLS estimation-effect term (`asy.lin.rep.ols %*% M1`) + — mean-zero by the OLS normal equations, so ATTs matched R at ~1e-11 while per-cell + SEs sat 4-13% and aggregated SEs 3-20% from the R golden fixtures (anti-conservative + on some: 0.958x on the two-period golden). (2) The unweighted ipw+cov branch computed + its per-cell SE as `sqrt(var_t/n_t + weighted_var_c)` where `weighted_var_c` was a + weighted POPULATION variance never scaled by an effective sample size (~7x inflated + per-cell SEs), and its IF lacked `std_ipw_did_panel`'s PS estimation-effect correction + (`asy.lin.rep.ps %*% M2`; aggregated SEs ~2.4% off). The survey ipw branch already + carried both (Phase 7a); the fix mirrors it. No-covariate reg/ipw per-cell SEs also + switched from ddof=1 plug-ins (ipw's `var_c*(1-p)/(n_c*p)` matched R only at p=0.5) + to the IF-based form, which is R's analytical SE exactly. Post-fix parity vs R: + per-cell and aggregated ~5e-12 on the reg goldens, ~1e-10 vs fresh did 2.5.1 ipw + aggregates (`tests/test_csdid_ported.py`, + `tests/test_methodology_callaway.py::TestDRDIDPanelIFParity`). Point estimates are + 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). + 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. - All aggregation SEs (simple, event study) include the weight influence function (WIF) adjustment, matching R's `did::aggte()`. The WIF accounts for uncertainty in estimating group-size aggregation weights. Group aggregation uses equal time weights (deterministic), @@ -636,15 +661,15 @@ The multiplier bootstrap uses random weights w_i with E[w]=0 and Var(w)=1: a base period later than `t` (matching R's `did::att_gt()`) - Does not require never-treated units: when all units are eventually treated, not-yet-treated cohorts serve as controls for each other (requires ≥2 cohorts) -- **Note:** CallawaySantAnna survey support: weights, strata, PSU, and FPC are all supported for all estimation methods (reg, ipw, dr) with or without covariates. Analytical (`n_bootstrap=0`): aggregated SEs use design-based variance via `compute_survey_if_variance()`. Bootstrap (`n_bootstrap>0`): PSU-level multiplier weights replace analytical SEs for aggregated quantities. IPW and DR with covariates use DRDID panel nuisance IF corrections (Phase 7a: PS IF correction via survey-weighted Hessian/score, OR IF correction via WLS bread and gradient; Sant'Anna & Zhao 2020, Theorem 3.1). Survey weights compose with IPW weights multiplicatively. WIF in aggregation matches R's did::wif() formula. Per-unit survey weights are extracted via `groupby(unit).first()` from the panel-normalized pweight array; on unbalanced panels the pweight normalization (`w * n_obs / sum(w)`) preserves relative unit weights since all IF/WIF formulas use weight ratios (`sw_i / sum(sw)`) where the normalization constant cancels. Scale-invariance tests pass on both balanced and unbalanced panels. +- **Note:** CallawaySantAnna survey support: weights, strata, PSU, and FPC are all supported for all estimation methods (reg, ipw, dr) with or without covariates. Analytical (`n_bootstrap=0`): aggregated SEs use design-based variance via `compute_survey_if_variance()`. Bootstrap (`n_bootstrap>0`): PSU-level multiplier weights replace analytical SEs for aggregated quantities. ALL covariate methods carry the DRDID panel nuisance IF corrections, method-uniformly across survey and unweighted paths (Phase 7a added the PS IF correction via survey-weighted Hessian/score and the DR OR IF correction via WLS bread and gradient, Sant'Anna & Zhao 2020, Theorem 3.1; v3.7 extended the PS correction to unweighted ipw and the OR estimation-effect correction to reg on both paths). Survey weights compose with IPW weights multiplicatively. WIF in aggregation matches R's did::wif() formula. Per-unit survey weights are extracted via `groupby(unit).first()` from the panel-normalized pweight array; on unbalanced panels the pweight normalization (`w * n_obs / sum(w)`) preserves relative unit weights since all IF/WIF formulas use weight ratios (`sw_i / sum(sw)`) where the normalization constant cancels. Scale-invariance tests pass on both balanced and unbalanced panels. - **Note (deviation from R):** Panel DR control augmentation is normalized by treated mass (`sw_t_sum` or `n_t`) rather than control IPW mass (`sum(w_cont)`). R's `DRDID::drdid_panel` uses `mean(w.cont)` as the control normalizer. Both are consistent asymptotically (under correct model specification, `E[w_cont] = E[D]` so the normalizers converge), but they differ in finite samples when IPW reweighting doesn't perfectly balance. The treated-mass normalization is simpler and matches the `did::att_gt` convention where ATT is defined per treated unit. Aligning to `DRDID::drdid_panel`'s exact `w.cont` normalization is deferred. - **Note:** PS nuisance IF corrections follow DRDID's M-estimation convention: `asy_lin_rep_psi` is computed on O(1) psi scale (matching R's `asy.lin.rep.ps = score %*% Hessian.ps`), then the correction `asy_lin_rep_psi @ M2` is converted to the library's O(1/n) phi convention via a single `/n` division. OR corrections use the same phi-scale pattern via `solve(X'WX)` (unnormalized Hessian). -- **Note (deviation from R):** CallawaySantAnna survey reg+covariates per-cell SE uses a conservative plug-in IF based on WLS residuals. The treated IF is `inf_treated_i = (sw_i/sum(sw_treated)) * (resid_i - ATT)` (normalized by treated weight sum, matching unweighted `(resid-ATT)/n_t`). The control IF is `inf_control_i = -(sw_i/sum(sw_control)) * wls_resid_i` (normalized by control weight sum, matching unweighted `-resid/n_c`). SE is computed as `sqrt(sum(sw_t_norm * (resid_t - ATT)^2) + sum(sw_c_norm * resid_c^2))`, the weighted analogue of the unweighted `sqrt(var_t/n_t + var_c/n_c)`. This omits the semiparametrically efficient nuisance correction from DRDID's `reg_did_panel` — WLS residuals are orthogonal to the weighted design matrix by construction, so the first-order IF term is asymptotically valid but may be conservative. SEs pass weight-scale-invariance tests. The efficient DRDID correction is deferred to future work. +- **Note:** CallawaySantAnna panel reg+covariates per-cell IF (survey AND unweighted) carries `DRDID::reg_did_panel`'s full OLS estimation-effect term (v3.7). The treated IF is unchanged: `inf_treated_i = (sw_i/sum(sw_treated)) * (resid_i - ATT)` (unweighted collapse `(resid-ATT)/n_t`). The control IF is `inf_control_i = -(sw_i * wls_resid_i) * proj_i` where `proj_i = 1/sum(sw_control) + (x_i - x̄_c)' G^{-1} (x̄_t - x̄_c)` with `G` the centered (weighted) control Gram and `x̄_t`/`x̄_c` the (weighted) treated/control covariate means — algebraically identical to R's `asy.lin.rep.ols %*% M1 / mean(w.treat)` but evaluated in the centered basis, which is offset-invariant (the raw Gram squares the design's conditioning; a large constant covariate offset would otherwise push the equilibrated rank check below threshold and silently truncate a genuine direction). The intercept-only collapse is exactly `-resid_c/n_c`. HISTORY: the pre-v3.7 IF used the raw-residual plug-in (`-sw_c_norm * wls_resid_i`), documented at the time as "asymptotically valid but may be conservative" — that rationale was WRONG: the omitted term is mean-zero (WLS orthogonality), which is why point estimates always matched R at ~1e-11, but its variance contribution is first-order, leaving per-cell SEs 4-13% and aggregated SEs 3-20% from R, and ANTI-conservative on some fixtures (0.958x on the two-period golden). Per-cell parity vs R goldens is now ~5e-12 (`tests/test_csdid_ported.py`, `tests/test_methodology_callaway.py::TestDRDIDPanelIFParity`). SEs pass weight-scale-invariance tests (only weight ratios enter every term). - **Note (deviation from R):** Per-cell ATT(g,t) SEs under survey weights use influence-function-based variance (matching R's `did::att_gt` analytical SE path) rather than full Taylor-series linearization. When strata/PSU/FPC are present, analytical aggregated SEs (`n_bootstrap=0`) use `compute_survey_if_variance()` on the combined IF/WIF; bootstrap aggregated SEs (`n_bootstrap>0`) use PSU-level multiplier weights. - **Note:** Repeated cross-sections (`panel=False`, Phase 7b): supports surveys like BRFSS, ACS annual, and CPS monthly where units are not followed over time. Uses cross-sectional DRDID (Sant'Anna & Zhao 2020, Section 4): `reg` matches `DRDID::reg_did_rc` (Eq 2.2), `dr` matches `DRDID::drdid_rc` (locally efficient, Eq 3.3+3.4 with 4 OLS fits), `ipw` matches `DRDID::std_ipw_did_rc`. Per-observation influence functions instead of per-unit. All three estimation methods support covariates and survey weights. - **Note:** Panel and RCS influence functions use the library-wide `phi_i = psi_i / n` convention (SE = `sqrt(sum(phi^2))`, algebraically equivalent to R's `sd(psi)*sqrt(n-1)/n`). Leading IF terms are computed on psi scale and divided by n; PS nuisance corrections are computed on psi scale (`score @ solve(Hessian)`) with a single `/n` conversion to phi. -- **Note:** Non-survey DR path also includes nuisance IF corrections (PS + OR), matching the survey path structure (Phase 7a). Previously used plug-in IF only. +- **Note:** Non-survey DR path also includes nuisance IF corrections (PS + OR), matching the survey path structure (Phase 7a). Previously used plug-in IF only. As of v3.7 the non-survey reg and ipw paths carry their corrections too (OR estimation-effect / PS score), so the nuisance-IF treatment is method-uniform. **Reference implementation(s):** - R: `did::att_gt()` (Callaway & Sant'Anna's official package) diff --git a/docs/survey-roadmap.md b/docs/survey-roadmap.md index 06d4f1d4c..afc40702f 100644 --- a/docs/survey-roadmap.md +++ b/docs/survey-roadmap.md @@ -273,6 +273,7 @@ the limitation and suggested alternative. **Warning/fallback (no error):** MultiPeriodDiD with `wild_bootstrap` + `survey_design` warns and falls back to analytical inference. -**Conservative approach (no error):** CallawaySantAnna `reg`+covariates -uses conservative plug-in IF rather than efficient DRDID nuisance IF -correction (see REGISTRY.md). +**Resolved (2026-07):** CallawaySantAnna `reg`+covariates (survey and +unweighted) now carries the full `DRDID::reg_did_panel` estimation-effect +IF correction; the former "conservative plug-in IF" deviation is removed +(see REGISTRY.md, CallawaySantAnna standard-error notes). diff --git a/tests/test_csdid_ported.py b/tests/test_csdid_ported.py index 36fc2c3f1..129a00b4f 100644 --- a/tests/test_csdid_ported.py +++ b/tests/test_csdid_ported.py @@ -966,7 +966,14 @@ def test_golden_default_dgp_dr(self, golden_values): ), f"ATT(g={g}, t={t}): Python={py_att:.6f}, R={r_att:.6f}" def test_golden_default_dgp_reg(self, golden_values): - """Regression method with covariates matches R.""" + """Regression method with covariates matches R - ATT AND per-cell SE. + + The SE assertion pins the ``DRDID::reg_did_panel`` influence function + (including the OLS estimation-effect term ``asy.lin.rep.ols @ M1``). + Before that term was added the per-cell SEs sat 4-13% from these + fixtures; observed agreement is now ~5e-12 (deterministic OLS), so + the 1e-8 tolerance is pure platform headroom. + """ if "with_covariates_reg" not in golden_values: pytest.skip("Scenario not in golden values") results, expected = self._run_scenario( @@ -978,14 +985,26 @@ def test_golden_default_dgp_reg(self, golden_values): for i, (g, t) in enumerate(zip(r_gt["group"], r_gt["time"])): g, t = int(g), int(t) if (g, t) in results.group_time_effects: - py_att = results.group_time_effects[(g, t)]["effect"] + eff = results.group_time_effects[(g, t)] r_att = r_gt["att"][i] assert ( - abs(py_att - r_att) < 0.02 - ), f"REG ATT(g={g}, t={t}): Py={py_att:.6f}, R={r_att:.6f}" + abs(eff["effect"] - r_att) < 0.02 + ), f"REG ATT(g={g}, t={t}): Py={eff['effect']:.6f}, R={r_att:.6f}" + assert ( + abs(eff["se"] - r_gt["se"][i]) < 1e-8 + ), f"REG SE(g={g}, t={t}): Py={eff['se']:.10f}, R={r_gt['se'][i]:.10f}" def test_golden_default_dgp_ipw(self, golden_values): - """IPW method matches R.""" + """IPW method matches R - ATT AND per-cell SE. + + The SE assertion pins the ``DRDID::std_ipw_did_panel`` influence + function (Hajek terms + the PS estimation-effect correction + ``asy.lin.rep.ps @ M2``) and the IF-based per-cell SE. Before the + fix the per-cell SEs were ~7x these fixtures (a weighted population + variance was never scaled by an effective sample size). Observed + agreement is now ~5e-12; 1e-6 leaves headroom for IRLS-vs-glm logit + convergence differences across platforms/BLAS. + """ if "with_covariates_ipw" not in golden_values: pytest.skip("Scenario not in golden values") results, expected = self._run_scenario( @@ -997,11 +1016,14 @@ def test_golden_default_dgp_ipw(self, golden_values): for i, (g, t) in enumerate(zip(r_gt["group"], r_gt["time"])): g, t = int(g), int(t) if (g, t) in results.group_time_effects: - py_att = results.group_time_effects[(g, t)]["effect"] + eff = results.group_time_effects[(g, t)] r_att = r_gt["att"][i] assert ( - abs(py_att - r_att) < 0.05 - ), f"IPW ATT(g={g}, t={t}): Py={py_att:.6f}, R={r_att:.6f}" + abs(eff["effect"] - r_att) < 0.05 + ), f"IPW ATT(g={g}, t={t}): Py={eff['effect']:.6f}, R={r_att:.6f}" + assert ( + abs(eff["se"] - r_gt["se"][i]) < 1e-6 + ), f"IPW SE(g={g}, t={t}): Py={eff['se']:.10f}, R={r_gt['se'][i]:.10f}" def test_golden_default_dgp_dr_covariates(self, golden_values): """DR method WITH covariates matches R `did` att AND se. @@ -1013,9 +1035,8 @@ def test_golden_default_dgp_dr_covariates(self, golden_values): TripleDifference + R's ``lm()``/QR); this pins both the ATT(g,t) and the analytical IF SE against R. Tolerances are empirical: DR ATT(g,t) agrees with R to ~1e-3 (DR small-sample numerics) and the DR SE to ~1e-4. - (``reg``/``ipw`` SE are NOT pinned here — they have a pre-existing - divergence from R's ``did`` unrelated to the OR solver: ipw SE depends on - the propensity IF, not the OR fit.) + (``reg``/``ipw`` per-cell SEs are pinned in their own golden tests + above, now that both carry the DRDID estimation-effect IF terms.) """ if "with_covariates_dr" not in golden_values: pytest.skip("Scenario not in golden values") @@ -1296,11 +1317,10 @@ def test_golden_event_study_aggregation_se(self, golden_values): this pins parity on the estimated horizons (observed agreement ~7e-6). - NOT enabled (pre-existing, unrelated to the IF-assembly fast path - - verified byte-identical deltas on the pre-rewrite tree): the - ``two_period`` and ``dynamic_effects`` aggregated SEs, where the reg - method's aggregated SEs sit 3-20% from these R fixtures while the - ATTs match at 1e-11. Tracked in TODO.md. + The ``two_period`` and ``dynamic_effects`` (reg) aggregated SEs - + previously 3-20% off because the reg IF lacked the OLS + estimation-effect term - are now enabled in + ``test_golden_reg_aggregation_se`` below. """ if "varying_vs_universal" not in golden_values: pytest.skip("Scenario not in golden values") @@ -1333,3 +1353,160 @@ def test_golden_event_study_aggregation_se(self, golden_values): ), f"{bp} dyn SE(e={e}): Py={es['se']:.6f}, R={r_se}" n_checked += 1 assert n_checked >= 8, f"{bp}: too few comparable horizons" + + def test_golden_reg_group_time_se(self, golden_values): + """Per-cell reg SEs match R across the remaining reg golden scenarios. + + Complements ``test_golden_default_dgp_reg``: ``two_period`` (single + 2x2 cell) and ``dynamic_effects`` (3 cohorts x 3 periods, including + PRE-treatment cells whose varying base period is t-1). These pin the + ``DRDID::reg_did_panel`` IF on cell shapes the default scenario does + not cover. Observed agreement ~5e-12 (deterministic OLS); 1e-8 is + platform headroom. + """ + for name in ("two_period", "dynamic_effects"): + if name not in golden_values: + pytest.skip("Scenario not in golden values") + scenario = golden_values[name] + data = _golden_to_df(scenario["data"]) + cs = CallawaySantAnna(estimation_method="reg") + results = cs.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X"], + ) + r_gt = scenario["results"]["group_time"] + groups = np.atleast_1d(r_gt["group"]) + times = np.atleast_1d(r_gt["time"]) + atts = np.atleast_1d(r_gt["att"]) + ses = np.atleast_1d(r_gt["se"]) + for g, t, r_att, r_se in zip(groups, times, atts, ses): + eff = results.group_time_effects[(int(g), int(t))] + assert ( + abs(eff["effect"] - r_att) < 1e-8 + ), f"{name} ATT(g={g},t={t}): Py={eff['effect']:.10f}, R={r_att:.10f}" + assert ( + abs(eff["se"] - r_se) < 1e-8 + ), f"{name} SE(g={g},t={t}): Py={eff['se']:.10f}, R={r_se:.10f}" + + def test_golden_reg_aggregation_se(self, golden_values): + """Aggregated reg ATTs and SEs match R ``aggte`` (simple/group/dynamic). + + Enables the previously-unasserted aggregation blocks of the two reg + golden scenarios: ``two_period`` (simple + group + dynamic of a + single cell) and ``dynamic_effects`` (dynamic, all horizons incl. + pre-treatment). Before the reg IF carried the OLS estimation-effect + term these SEs sat 3-20% from the fixtures while the ATTs matched at + 1e-11. Observed post-fix agreement is ~5e-12 (deterministic OLS + + the aggregation machinery already validated at 1e-12 on the dr + scenarios); 1e-8 is platform headroom. + """ + if "two_period" not in golden_values or "dynamic_effects" not in golden_values: + pytest.skip("Scenario not in golden values") + + # two_period: one (2,2) cell, so simple == group == dynamic(e=0) + scenario = golden_values["two_period"] + data = _golden_to_df(scenario["data"]) + cs = CallawaySantAnna(estimation_method="reg") + results = cs.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X"], + aggregate="all", + ) + r = scenario["results"] + assert abs(results.overall_att - r["simple"]["att"]) < 1e-8 + assert abs(results.overall_se - r["simple"]["se"]) < 1e-8 + g_att = np.atleast_1d(r["group"]["att"])[0] + g_se = np.atleast_1d(r["group"]["se"])[0] + grp = results.group_effects[2] + assert abs(grp["effect"] - g_att) < 1e-8 + assert abs(grp["se"] - g_se) < 1e-8 + d_att = np.atleast_1d(r["dynamic"]["att"])[0] + d_se = np.atleast_1d(r["dynamic"]["se"])[0] + es0 = results.event_study_effects[0] + assert abs(es0["effect"] - d_att) < 1e-8 + assert abs(es0["se"] - d_se) < 1e-8 + + # dynamic_effects: full event study, pre and post horizons + scenario = golden_values["dynamic_effects"] + data = _golden_to_df(scenario["data"]) + cs = CallawaySantAnna(estimation_method="reg") + results = cs.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X"], + aggregate="event_study", + ) + r_dyn = scenario["results"]["dynamic"] + for i, e in enumerate(r_dyn["egt"]): + es = results.event_study_effects.get(int(e)) + assert es is not None, f"missing event time e={e}" + assert ( + abs(es["effect"] - r_dyn["att"][i]) < 1e-8 + ), f"dyn ATT(e={e}): Py={es['effect']:.10f}, R={r_dyn['att'][i]:.10f}" + assert ( + abs(es["se"] - r_dyn["se"][i]) < 1e-8 + ), f"dyn SE(e={e}): Py={es['se']:.10f}, R={r_dyn['se'][i]:.10f}" + + def test_golden_ipw_aggregation_se_vs_r_did_251(self, golden_values): + """Aggregated ipw SEs match R did 2.5.1 ``aggte`` on the golden data. + + The golden JSON carries no aggregation blocks for the ipw scenario, + so the expected values are hardcoded from a fresh R run (did 2.5.1, + DRDID 1.3.0, ``bstrap=FALSE``, computed 2026-07-05 on the EXACT + ``with_covariates_ipw`` fixture data; see + benchmarks/R/generate_csdid_test_values.R which will fold these into + the JSON on its next regeneration). Discriminating: the pre-fix + Python simple SE was 0.32125722 (missing PS estimation-effect + correction, ~2.4% off). Observed post-fix agreement ~1e-10; 1e-6 + covers the fixtures' 8-decimal quantization + IRLS-vs-glm headroom. + """ + if "with_covariates_ipw" not in golden_values: + pytest.skip("Scenario not in golden values") + scenario = golden_values["with_covariates_ipw"] + data = _golden_to_df(scenario["data"]) + cs = CallawaySantAnna(estimation_method="ipw") + results = cs.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X"], + aggregate="all", + ) + # R did 2.5.1: aggte(type="simple", bstrap=FALSE, cband=FALSE) + assert abs(results.overall_att - 0.68857256) < 1e-6 + assert abs(results.overall_se - 0.31382230) < 1e-6 + # aggte(type="dynamic"): (att.egt, se.egt) per event time + r_dynamic = { + -2: (-0.58336476, 0.27944371), + -1: (-0.40369753, 0.22022753), + 0: (0.68317989, 0.20153325), + 1: (0.69816671, 0.38903747), + 2: (0.68521672, 0.64844936), + } + for e, (r_att, r_se) in r_dynamic.items(): + es = results.event_study_effects[e] + assert abs(es["effect"] - r_att) < 1e-6, f"ipw dyn ATT(e={e})" + assert abs(es["se"] - r_se) < 1e-6, f"ipw dyn SE(e={e})" + # aggte(type="group"): (att.egt, se.egt) per cohort + r_group = { + 2: (0.91085897, 0.44399703), + 3: (0.60976619, 0.33782836), + 4: (0.24708101, 0.30682305), + } + for g, (r_att, r_se) in r_group.items(): + grp = results.group_effects[g] + assert abs(grp["effect"] - r_att) < 1e-6, f"ipw group ATT(g={g})" + assert abs(grp["se"] - r_se) < 1e-6, f"ipw group SE(g={g})" diff --git a/tests/test_methodology_callaway.py b/tests/test_methodology_callaway.py index 33c75b0de..2d12dc68d 100644 --- a/tests/test_methodology_callaway.py +++ b/tests/test_methodology_callaway.py @@ -10,16 +10,19 @@ with multiple time periods. Journal of Econometrics, 225(2), 200-230. """ +import json import subprocess import unittest.mock import warnings -from typing import Any, Dict, Tuple +from pathlib import Path +from typing import Any, Dict, Optional, Tuple import numpy as np import pandas as pd import pytest -from diff_diff import CallawaySantAnna +from diff_diff import CallawaySantAnna, SurveyDesign +from diff_diff.linalg import solve_logit from diff_diff.prep import generate_staggered_data from diff_diff.bootstrap_utils import ( generate_bootstrap_weights_batch as _generate_bootstrap_weights_batch, @@ -1985,3 +1988,347 @@ def test_or_fit_offset_invariant(self): assert eff["effect"] == pytest.approx( shifted.group_time_effects[gt]["effect"], abs=1e-7 ), f"ATT{gt} not offset-invariant" + + +# ============================================================================= +# DRDID panel influence-function parity (numpy cross-checks) +# ============================================================================= + +_GOLDEN_VALUES_PATH = ( + Path(__file__).parents[1] / "benchmarks" / "data" / "csdid_golden_values.json" +) + + +def _reg_did_panel_ref( + dY: np.ndarray, + D: np.ndarray, + X: np.ndarray, + weights: Optional[np.ndarray] = None, +) -> Tuple[float, float]: + """Independent numpy reconstruction of ``DRDID::reg_did_panel``. + + Returns ``(att, se)`` where ``se`` is the analytical influence-function + standard error ``sd(psi) * sqrt(n-1) / n``. The control-side psi carries + the OLS estimation-effect term ``asy.lin.rep.ols @ M1`` (the WLS residual + projected through the bread and evaluated at the treated covariate mean); + its intercept-only collapse is ``-resid_c / n_c``. + """ + n = len(D) + w = np.ones(n) if weights is None else weights / np.mean(weights) + Xi = np.column_stack([np.ones(n), X]) + c = D == 0 + Wc = w[c] + beta, *_ = np.linalg.lstsq( + np.sqrt(Wc)[:, None] * Xi[c], np.sqrt(Wc) * dY[c], rcond=None + ) + pred = Xi @ beta + w_treat = w * D + mw_t = w_treat.mean() + att = float(np.sum(w_treat * (dY - pred)) / np.sum(w_treat)) + resid_c = dY[c] - pred[c] + XpX_inv = np.linalg.inv(Xi[c].T @ (Wc[:, None] * Xi[c]) / n) + asy_lin_rep_ols = (Wc * resid_c)[:, None] * Xi[c] @ XpX_inv + M1 = (w_treat[:, None] * Xi).sum(axis=0) / n + psi = np.zeros(n) + psi[~c] = w_treat[~c] * (dY[~c] - pred[~c] - att) / mw_t + psi[c] = -(asy_lin_rep_ols @ M1) / mw_t + se = float(psi.std(ddof=1) * np.sqrt(n - 1) / n) + return att, se + + +def _std_ipw_did_panel_ref( + dY: np.ndarray, + D: np.ndarray, + X: np.ndarray, + ps: np.ndarray, + weights: Optional[np.ndarray] = None, +) -> Tuple[float, float]: + """Independent numpy reconstruction of ``DRDID::std_ipw_did_panel``. + + Hajek IPW with the propensity-score estimation-effect correction + ``asy.lin.rep.ps @ M2``. Takes the fitted propensity scores as input so + the check isolates the IF algebra from the logit solver. + """ + n = len(D) + w = np.ones(n) if weights is None else weights / np.mean(weights) + Xi = np.column_stack([np.ones(n), X]) + w_t = w * D + w_c = w * ps * (1 - D) / (1 - ps) + mw_t, mw_c = w_t.mean(), w_c.mean() + eta_t = (w_t * dY).mean() / mw_t + eta_c = (w_c * dY).mean() / mw_c + att = float(eta_t - eta_c) + W_ps = ps * (1 - ps) * w + H_inv = np.linalg.inv(Xi.T @ (W_ps[:, None] * Xi)) * n + score = (w * (D - ps))[:, None] * Xi + asy_lin_rep_ps = score @ H_inv + inf_t = (w_t * dY - w_t * eta_t) / mw_t + inf_c1 = w_c * dY - w_c * eta_c + M2 = ((w_c * (dY - eta_c))[:, None] * Xi).sum(axis=0) / n + psi = inf_t - (inf_c1 + asy_lin_rep_ps @ M2) / mw_c + se = float(psi.std(ddof=1) * np.sqrt(n - 1) / n) + return att, se + + +def _make_cov_panel( + seed: int, + n_units: int = 400, + n_periods: int = 2, + cohorts: Tuple[int, ...] = (2,), +) -> pd.DataFrame: + """Two-covariate staggered panel with X-dependent selection and trends. + + The selection index is bounded so every cell's fitted propensity stays + well inside (0.05, 0.95) - the estimator's clip at ``pscore_trim`` and + R's ``trim.level`` drop are both no-ops, keeping the reconstruction and + the estimator on identical weights. + """ + rng = np.random.default_rng(seed) + X1 = rng.normal(size=n_units) + X2 = rng.uniform(-1.0, 1.0, size=n_units) + p_treat = 1.0 / (1.0 + np.exp(-(0.6 * X1 - 0.5 * X2))) + u = rng.uniform(size=n_units) + first_treat = np.zeros(n_units, dtype=int) + # Split the treated mass evenly across cohorts. + share = 0.55 + for j, g in enumerate(cohorts): + lo = share * p_treat * j / len(cohorts) + hi = share * p_treat * (j + 1) / len(cohorts) + first_treat[(u >= lo) & (u < hi)] = g + alpha = rng.normal(size=n_units) + 0.5 * X1 + weight = rng.uniform(0.5, 2.0, size=n_units) + rows = [] + for t in range(1, n_periods + 1): + eps = rng.normal(scale=0.5, size=n_units) + trend = 0.5 * t + (0.3 * X1 - 0.2 * X2) * t + effect = np.where( + (first_treat > 0) & (t >= first_treat), 1.0 + 0.5 * (t - first_treat), 0.0 + ) + y = alpha + trend + effect + eps + rows.append( + pd.DataFrame( + { + "unit": np.arange(n_units), + "period": t, + "first_treat": first_treat, + "outcome": y, + "X1": X1, + "X2": X2, + "weight": weight, + } + ) + ) + return pd.concat(rows, ignore_index=True) + + +def _cell_arrays( + df: pd.DataFrame, + g: int, + t: int, + base: int, + control_cohorts: Tuple[int, ...], +) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Unit-level (dY, D, X, weights) arrays for one (g, t) cell.""" + keep = list(control_cohorts) + [g] + sub = df[df["first_treat"].isin(keep)] + wide = sub.pivot(index="unit", columns="period", values="outcome") + meta = sub.drop_duplicates(subset="unit").set_index("unit") + units = wide[[base, t]].dropna().index + dY = (wide.loc[units, t] - wide.loc[units, base]).to_numpy() + D = (meta.loc[units, "first_treat"] == g).to_numpy(dtype=float) + X = meta.loc[units, ["X1", "X2"]].to_numpy() + w = meta.loc[units, "weight"].to_numpy() + return dY, D, X, w + + +class TestDRDIDPanelIFParity: + """Numpy cross-checks of CS panel reg/ipw per-cell IF standard errors + against independent reconstructions of ``DRDID::reg_did_panel`` and + ``DRDID::std_ipw_did_panel`` (the estimators R's ``did::att_gt`` + delegates to for panel data). + + These pin the estimation-effect (nuisance) IF terms: the OLS + ``asy.lin.rep.ols @ M1`` term for reg and the PS score correction + ``asy.lin.rep.ps @ M2`` for ipw - the terms whose omission left ATTs + exact (both are mean-zero) while per-cell/aggregated SEs diverged from + R by 4-13% / 3-20% (reg) and ~7x per-cell (ipw). + """ + + def _fit( + self, df: pd.DataFrame, method: str, survey_design: Any = None, **cs_kwargs + ) -> Any: + fit_kwargs = {"survey_design": survey_design} if survey_design is not None else {} + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return CallawaySantAnna(estimation_method=method, **cs_kwargs).fit( + df, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X1", "X2"], + **fit_kwargs, + ) + + @staticmethod + def _assert_ps_interior(X: np.ndarray, D: np.ndarray) -> np.ndarray: + """Fit the cell logit and assert the clip/trim guards are no-ops.""" + _, ps = solve_logit(X, D) + assert ps.min() > 0.02 and ps.max() < 0.98, "DGP must keep ps interior" + return ps + + def test_reg_if_crosscheck(self): + """Unweighted reg+cov single cell: fitted per-cell/overall SE must + equal the reg_did_panel reconstruction (vectorized covariate path).""" + df = _make_cov_panel(seed=90210) + res = self._fit(df, "reg") + dY, D, X, _ = _cell_arrays(df, g=2, t=2, base=1, control_cohorts=(0,)) + att, se = _reg_did_panel_ref(dY, D, X) + cell = res.group_time_effects[(2, 2)] + assert cell["effect"] == pytest.approx(att, rel=1e-9) + assert cell["se"] == pytest.approx(se, rel=1e-9) + assert res.overall_se == pytest.approx(se, rel=1e-9) + + def test_reg_if_crosscheck_not_yet_treated(self): + """not_yet_treated forces the per-pair (non-hoisted) bread branch of + the vectorized producer; both post cells must match reg_did_panel.""" + df = _make_cov_panel(seed=90211, n_periods=3, cohorts=(2, 3)) + res = self._fit(df, "reg", control_group="not_yet_treated") + # (2,2): controls = never-treated + cohort 3 (not yet treated at t=2) + dY, D, X, _ = _cell_arrays(df, g=2, t=2, base=1, control_cohorts=(0, 3)) + att, se = _reg_did_panel_ref(dY, D, X) + cell = res.group_time_effects[(2, 2)] + assert cell["effect"] == pytest.approx(att, rel=1e-9) + assert cell["se"] == pytest.approx(se, rel=1e-9) + # (3,3): only never-treated remain + dY, D, X, _ = _cell_arrays(df, g=3, t=3, base=2, control_cohorts=(0,)) + att, se = _reg_did_panel_ref(dY, D, X) + cell = res.group_time_effects[(3, 3)] + assert cell["effect"] == pytest.approx(att, rel=1e-9) + assert cell["se"] == pytest.approx(se, rel=1e-9) + + def test_reg_if_crosscheck_survey(self): + """Survey-weighted reg+cov (general path): fitted per-cell SE must + equal the reg_did_panel reconstruction with i.weights.""" + df = _make_cov_panel(seed=90212) + res = self._fit(df, "reg", survey_design=SurveyDesign(weights="weight")) + dY, D, X, w = _cell_arrays(df, g=2, t=2, base=1, control_cohorts=(0,)) + att, se = _reg_did_panel_ref(dY, D, X, weights=w) + cell = res.group_time_effects[(2, 2)] + assert cell["effect"] == pytest.approx(att, rel=1e-9) + assert cell["se"] == pytest.approx(se, rel=1e-9) + + def test_ipw_if_crosscheck(self): + """Unweighted ipw+cov: fitted per-cell/overall SE must equal the + std_ipw_did_panel reconstruction (PS estimation effect included).""" + df = _make_cov_panel(seed=90213) + res = self._fit(df, "ipw") + dY, D, X, _ = _cell_arrays(df, g=2, t=2, base=1, control_cohorts=(0,)) + ps = self._assert_ps_interior(X, D) + att, se = _std_ipw_did_panel_ref(dY, D, X, ps) + cell = res.group_time_effects[(2, 2)] + assert cell["effect"] == pytest.approx(att, rel=1e-8) + assert cell["se"] == pytest.approx(se, rel=1e-8) + assert res.overall_se == pytest.approx(se, rel=1e-8) + + def test_ipw_if_crosscheck_survey(self): + """Survey-weighted ipw+cov: the survey branch already carried the PS + correction (Phase 7a) - this pins it against the i.weights + reconstruction.""" + df = _make_cov_panel(seed=90214) + res = self._fit(df, "ipw", survey_design=SurveyDesign(weights="weight")) + dY, D, X, w = _cell_arrays(df, g=2, t=2, base=1, control_cohorts=(0,)) + _, ps = solve_logit(X, D, weights=w) + assert ps.min() > 0.02 and ps.max() < 0.98 + att, se = _std_ipw_did_panel_ref(dY, D, X, ps, weights=w) + cell = res.group_time_effects[(2, 2)] + assert cell["effect"] == pytest.approx(att, rel=1e-8) + assert cell["se"] == pytest.approx(se, rel=1e-8) + + +@pytest.fixture(scope="module") +def golden(): + """Golden R fixture scenarios; skip when the file is absent.""" + if not _GOLDEN_VALUES_PATH.exists(): + pytest.skip( + "Golden values file not found; " + "run: Rscript benchmarks/R/generate_csdid_test_values.R" + ) + with open(_GOLDEN_VALUES_PATH) as f: + return json.load(f)["scenarios"] + + +class TestDRDIDGoldenIFReconstruction: + """Reconstruct the DRDID per-cell IF SEs in numpy from the golden-fixture + data and pin BOTH directions per cell: + + 1. reconstruction == R's golden SEs (the ground truth; this held before + the fix and guards the reconstruction itself), and + 2. fitted == reconstruction (the fix contract). + + Covers pre-treatment cells and varying base periods (base = t-1 for + t < g, else g-1), which the DGP-based cross-checks above do not. + """ + + @staticmethod + def _df(scenario: dict) -> pd.DataFrame: + d = scenario["data"] + return pd.DataFrame( + { + "unit": d["unit"], + "period": d["period"], + "first_treat": d["first_treat"], + "outcome": d["outcome"], + "X": d["X"], + } + ) + + @staticmethod + def _cell(df: pd.DataFrame, g: int, t: int): + base = t - 1 if t < g else g - 1 + sub = df[df["first_treat"].isin([0, g])] + wide = sub.pivot(index="unit", columns="period", values="outcome") + meta = sub.drop_duplicates(subset="unit").set_index("unit") + units = wide[[base, t]].dropna().index + dY = (wide.loc[units, t] - wide.loc[units, base]).to_numpy() + D = (meta.loc[units, "first_treat"] == g).to_numpy(dtype=float) + X = meta.loc[units, "X"].to_numpy()[:, None] + return dY, D, X + + def _run(self, golden, name: str, method: str, ref): + scenario = golden[name] + df = self._df(scenario) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = CallawaySantAnna(estimation_method=method).fit( + df, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X"], + ) + gt = scenario["results"]["group_time"] + for g, t, r_att, r_se in zip(gt["group"], gt["time"], gt["att"], gt["se"]): + dY, D, X = self._cell(df, int(g), int(t)) + if method == "ipw": + _, ps = solve_logit(X, D) + att, se = ref(dY, D, X, ps) + else: + att, se = ref(dY, D, X) + # Direction 1: reconstruction reproduces R's golden values. + assert att == pytest.approx(r_att, rel=1e-6), f"recon att ({g},{t})" + assert se == pytest.approx(r_se, rel=1e-6), f"recon se ({g},{t})" + # Direction 2: the fitted estimator matches the reconstruction. + cell = res.group_time_effects[(int(g), int(t))] + assert cell["effect"] == pytest.approx(att, rel=1e-7), f"fit att ({g},{t})" + assert cell["se"] == pytest.approx(se, rel=1e-7), f"fit se ({g},{t})" + + def test_reg_golden_cells(self, golden): + self._run(golden, "with_covariates_reg", "reg", _reg_did_panel_ref) + + def test_reg_golden_cells_multi_cohort(self, golden): + self._run(golden, "dynamic_effects", "reg", _reg_did_panel_ref) + + def test_ipw_golden_cells(self, golden): + self._run(golden, "with_covariates_ipw", "ipw", _std_ipw_did_panel_ref) diff --git a/tests/test_se_accuracy.py b/tests/test_se_accuracy.py index c56ca1503..e56876585 100644 --- a/tests/test_se_accuracy.py +++ b/tests/test_se_accuracy.py @@ -148,7 +148,9 @@ def test_individual_att_gt_se_formula(self, benchmark_data): """ Test that individual ATT(g,t) SE follows expected formula. - For simple diff-in-means: SE = sqrt(var_t/n_t + var_c/n_c) + For simple diff-in-means the per-cell SE is derived from the + influence function, sqrt(sum(IF^2)) (== the ddof=0 form of + sqrt(var_t/n_t + var_c/n_c)), matching DRDID::reg_did_panel. """ cs = CallawaySantAnna( estimation_method="reg", # outcome regression (simpler) diff --git a/tests/test_staggered.py b/tests/test_staggered.py index 1bed85d59..d4b0b34e5 100644 --- a/tests/test_staggered.py +++ b/tests/test_staggered.py @@ -1611,12 +1611,16 @@ def test_constant_covariate_finite_se_matches_drop_one(self, method): with_const.overall_att, drop_one.overall_att, rtol=1e-9 ) - def test_constant_covariate_emits_single_rank_guard_warning(self): + @pytest.mark.parametrize("method", ["reg", "ipw", "dr"]) + def test_constant_covariate_emits_single_rank_guard_warning(self, method): + # reg/ipw now route their IF breads (OLS bread / PS Hessian) through + # _safe_inv like dr, so the aggregate warning fires for ALL methods + # on a collinear design (previously dr/survey-ipw only). data = generate_staggered_data_with_covariates(seed=789) data["xc"] = 5.0 with warnings.catch_warnings(record=True) as caught: warnings.simplefilter("always") - CallawaySantAnna(estimation_method="dr").fit( + CallawaySantAnna(estimation_method=method).fit( data, "outcome", "unit", "time", "first_treat", covariates=["x1", "xc"] ) rank_guard = [ @@ -1651,11 +1655,16 @@ def test_clustered_constant_covariate_finite_se(self): assert np.isfinite(res.overall_se) assert res.overall_se < 1.0 - def test_rank0_bread_propagates_nan_not_zero(self, monkeypatch): + @pytest.mark.parametrize("method", ["reg", "ipw", "dr"]) + def test_rank0_bread_propagates_nan_not_zero(self, monkeypatch, method): # rank-0 is unreachable through covariates alone (the always-present # intercept guarantees rank >= 1), so simulate an all-NaN bread to # exercise the NaN-masking fix: var_psi becomes NaN and must yield a NaN - # SE, NOT 0.0 via the old ``var_psi > 0 else 0.0`` guard. + # SE, NOT 0.0 via the old ``var_psi > 0 else 0.0`` guard. reg/ipw now + # derive their per-cell SE from _safe_inv-based IF terms too (OLS + # estimation-effect bread / PS Hessian), so all three methods share + # the contract. The point estimate does NOT depend on the bread, so + # it stays finite (NaN inference on an estimable cell, not _nan_gt_entry). from tests.conftest import assert_nan_inference import diff_diff.staggered as staggered_mod @@ -1667,10 +1676,11 @@ def _all_nan_inv(A, tracker=None): data = generate_staggered_data_with_covariates(seed=7) with warnings.catch_warnings(): warnings.simplefilter("ignore") - res = CallawaySantAnna(estimation_method="dr").fit( + res = CallawaySantAnna(estimation_method=method).fit( data, "outcome", "unit", "time", "first_treat", covariates=["x1", "x2"] ) for cell in res.group_time_effects.values(): + assert np.isfinite(cell["effect"]), "point estimate is bread-independent" assert np.isnan(cell["se"]), "rank-0 bread must give NaN SE, not 0.0" assert_nan_inference(cell) assert np.isnan(res.overall_se) @@ -1923,6 +1933,119 @@ def test_exact_duplicate_covariate_survey_weighted(self, method): np.testing.assert_allclose(big_ab.overall_se, big_ba.overall_se, rtol=1e-7) +class TestRegIpwIFBehavior: + """Behavioral contracts of the DRDID-parity reg/ipw per-cell IF/SE fix + (estimation-effect terms + per-cell SE = sqrt(sum(IF^2))).""" + + def test_ipw_pscore_fallback_uses_uncorrected_if_se(self, monkeypatch): + """When the per-cell logit fails and pscore_fallback="unconditional" + kicks in, the PS estimation-effect correction is SKIPPED (a constant + propensity has no estimated parameter) and the cell collapses to the + difference-in-means IF - i.e. per-cell effect AND se must equal the + no-covariate ipw fit on the same data.""" + import diff_diff.staggered as staggered_mod + + def _failing_logit(*args, **kwargs): + raise ValueError("forced logit failure for fallback test") + + data = generate_staggered_data_with_covariates(seed=31) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + no_cov = CallawaySantAnna(estimation_method="ipw").fit( + data, "outcome", "unit", "time", "first_treat" + ) + monkeypatch.setattr(staggered_mod, "solve_logit", _failing_logit) + fallback = CallawaySantAnna( + estimation_method="ipw", pscore_fallback="unconditional" + ).fit( + data, "outcome", "unit", "time", "first_treat", covariates=["x1", "x2"] + ) + assert set(fallback.group_time_effects) == set(no_cov.group_time_effects) + for key, cell in fallback.group_time_effects.items(): + ref = no_cov.group_time_effects[key] + np.testing.assert_allclose(cell["effect"], ref["effect"], rtol=1e-12) + np.testing.assert_allclose(cell["se"], ref["se"], rtol=1e-12) + + def test_underdetermined_control_cell_reg_no_crash(self): + """Cells with fewer controls than covariate columns (n_c < k+1) fit a + reduced design; the rank-guarded IF bread column-drops and the + interpolating fit leaves ~zero control residuals, so the SE is finite + (treated-side variation only), never a crash or a silent 0.""" + rng = np.random.default_rng(11) + n_units, k = 60, 5 + rows = [] + for u in range(n_units): + ft = 0 if u < 3 else 2 # only 3 never-treated controls + x = rng.normal(size=k) + base = rng.normal() + for t in (1, 2): + y = base + 0.4 * t + ( + 1.0 if (ft == 2 and t >= 2) else 0.0 + ) + rng.normal(0, 0.3) + rows.append( + {"unit": u, "time": t, "first_treat": ft, "outcome": y, + **{f"x{j}": x[j] for j in range(k)}} + ) + data = pd.DataFrame(rows) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = CallawaySantAnna( + estimation_method="reg", rank_deficient_action="silent" + ).fit( + data, "outcome", "unit", "time", "first_treat", + covariates=[f"x{j}" for j in range(k)], + ) + cell = res.group_time_effects[(2, 2)] + assert np.isfinite(cell["effect"]) + assert np.isfinite(cell["se"]) and cell["se"] > 0 + + def test_universal_base_period_anticipation_reg_smoke(self): + """reg+cov under base_period="universal" + anticipation=1: every + produced cell has finite inference and the reference period + t = g-1-anticipation is excluded from the grid (no degenerate + zero-IF cell exists to divide by).""" + data = generate_staggered_data_with_covariates(seed=97) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = CallawaySantAnna( + estimation_method="reg", base_period="universal", anticipation=1 + ).fit( + data, "outcome", "unit", "time", "first_treat", + covariates=["x1", "x2"], + ) + assert len(res.group_time_effects) > 0 + for (g, t), cell in res.group_time_effects.items(): + assert t != g - 2, "universal reference period must be excluded" + if cell["skip_reason"] is None: + assert np.isfinite(cell["se"]), f"cell ({g},{t})" + assert np.isfinite(res.overall_se) + + def test_uniform_survey_weights_match_unweighted_per_cell_se(self): + """Uniform survey weights route reg+cov through the general + (survey-branch) producer while the unweighted fit takes the + vectorized producer; both now share the same DRDID IF algebra, so + per-cell effects AND SEs must agree.""" + from diff_diff.survey import SurveyDesign + + data = generate_staggered_data_with_covariates(seed=53) + data["w_ones"] = 1.0 + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + unweighted = CallawaySantAnna(estimation_method="reg").fit( + data, "outcome", "unit", "time", "first_treat", + covariates=["x1", "x2"], + ) + uniform = CallawaySantAnna(estimation_method="reg").fit( + data, "outcome", "unit", "time", "first_treat", + covariates=["x1", "x2"], + survey_design=SurveyDesign(weights="w_ones"), + ) + for key, cell in unweighted.group_time_effects.items(): + ref = uniform.group_time_effects[key] + np.testing.assert_allclose(cell["effect"], ref["effect"], rtol=1e-9) + np.testing.assert_allclose(cell["se"], ref["se"], rtol=1e-9) + + class TestCallawaySantAnnaRankDeficiencyPaths: """Tests for rank-deficiency handling in DR and reg not_yet_treated paths.""" diff --git a/tests/test_survey_r_crossvalidation.py b/tests/test_survey_r_crossvalidation.py index 99b04b001..f978c78ce 100644 --- a/tests/test_survey_r_crossvalidation.py +++ b/tests/test_survey_r_crossvalidation.py @@ -258,9 +258,10 @@ class TestCSWeighted: diff-diff computes design-based SEs via compute_survey_if_variance(). Therefore we compare ATT only and verify that diff-diff's design-based SE is finite, positive, and differs from the naive (unweighted) SE. - For cs_weighted_cov, the reg+covariates SE additionally differs due to - the documented deviation in REGISTRY.md (conservative plug-in IF vs - semiparametrically efficient IF). + (The survey reg+covariates per-cell IF now carries the full DRDID + estimation-effect term — the former "conservative plug-in IF" deviation + is removed — but the design-based aggregated variance still differs + from R's survey-naive SEs by construction.) """ def _fit_scenario(self, r_results, key): From 9cb5f26a2a67cdd82c3cd4121efd90d6ce65ede5 Mon Sep 17 00:00:00 2001 From: igerber Date: Sun, 5 Jul 2026 10:40:56 -0400 Subject: [PATCH 2/3] fix(staggered): rank-0 centered reg bread collapses to zero correction, not NaN CI review P1: with a constant as the ONLY reg covariate the centered estimation-effect Gram is rank-0, and _safe_inv's all-NaN sentinel NaN-poisoned per-cell/aggregated SEs for cells whose intercept-only projection is fully identified (pre-fix raw-Gram column-drop and the no-covariate fit both give finite SEs). New _centered_or_bread helper (shared by all three reg sites) maps the rank-0 sentinel to a zero correction - proj collapses to 1/sum(W_c), equal to the no-covariate fit - while partial deficiency and the aggregate rank-guard warning are unchanged. ipw/dr breads are over [1, X] and keep NaN propagation (true pathology); the rank-0 monkeypatch test is re-scoped to ipw/dr with constant-only reg equivalence tests (unweighted + survey) added. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01X4AzrFUMqJxcUumSH31mSr --- diff_diff/staggered.py | 37 ++++++++++------ docs/methodology/REGISTRY.md | 9 ++++ tests/test_staggered.py | 83 +++++++++++++++++++++++++++++++++--- 3 files changed, 111 insertions(+), 18 deletions(-) diff --git a/diff_diff/staggered.py b/diff_diff/staggered.py index 718ef1671..1662184fe 100644 --- a/diff_diff/staggered.py +++ b/diff_diff/staggered.py @@ -1437,10 +1437,7 @@ def _compute_all_att_gt_covariate_reg( if is_balanced and self.control_group == "never_treated": group_xbar_c = X_ctrl_raw.mean(axis=0) group_Xc_centered = X_ctrl_raw - group_xbar_c - group_bread = _safe_inv( - group_Xc_centered.T @ group_Xc_centered, - tracker=self._safe_inv_tracker, - ) + group_bread = self._centered_or_bread(group_Xc_centered.T @ group_Xc_centered) # Process each (g, t) pair in this group for g, t, bp_val, base_col, post_col in tasks: @@ -1604,10 +1601,7 @@ def _compute_all_att_gt_covariate_reg( else: xbar_c = X_control_pair.mean(axis=0) Xc_centered = X_control_pair - xbar_c - bread = _safe_inv( - Xc_centered.T @ Xc_centered, - tracker=self._safe_inv_tracker, - ) + bread = self._centered_or_bread(Xc_centered.T @ Xc_centered) d_tc = X_treated_pair.mean(axis=0) - xbar_c with np.errstate(all="ignore"): proj_c = 1.0 / pair_n_c + Xc_centered @ (bread @ d_tc) @@ -2616,6 +2610,28 @@ def fit( self.is_fitted_ = True return self.results_ + def _centered_or_bread(self, gram: np.ndarray) -> np.ndarray: + """Rank-guarded inverse of the CENTERED control covariate Gram. + + Unlike the ``[1, X]`` breads elsewhere, the intercept direction is + handled analytically here (the ``1/sum(W_c)`` leading term of the + reg estimation-effect projection), so a rank-0 centered Gram — zero + within-control covariate variation, e.g. a constant covariate as + the only regressor — is benign: the correction on the identified + (intercept-only) subset is exactly zero, and the projection + collapses to ``1/sum(W_c)``. Map ``_safe_inv``'s all-NaN rank-0 + sentinel to a zero matrix so the cell keeps the finite SE on the + identified subset (rank-guard REGISTRY contract; equals both the + no-covariate fit and the pre-centered raw-Gram column-drop) instead + of NaN-poisoning the IF. Partial deficiency is untouched — + ``_safe_inv`` already returns a zero-filled column-dropped inverse, + and the aggregate rank-guard warning still fires via the tracker. + """ + bread = _safe_inv(gram, tracker=self._safe_inv_tracker) + if bread.size and not np.any(np.isfinite(bread)): + return np.zeros_like(bread) + return bread + def _outcome_regression( self, treated_change: np.ndarray, @@ -2704,10 +2720,7 @@ def _outcome_regression( with np.errstate(all="ignore"): xbar_c = np.sum(W_c[:, None] * X_control, axis=0) / w_c_sum Xc_centered = X_control - xbar_c - bread = _safe_inv( - Xc_centered.T @ (W_c[:, None] * Xc_centered), - tracker=self._safe_inv_tracker, - ) + bread = self._centered_or_bread(Xc_centered.T @ (W_c[:, None] * Xc_centered)) proj_c = 1.0 / w_c_sum + Xc_centered @ (bread @ (xbar_t - xbar_c)) inf_control = -(W_c * resid_c) * proj_c inf_func = np.concatenate([inf_treated, inf_control]) diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 3ac6e8dd8..c1026a91f 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -523,6 +523,15 @@ Aggregations: DR's no-covariate per-cell SE keeps its ddof=1 plug-in (O(1/n) from R; TODO row). 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 + handled analytically by the `1/sum(W_c)` leading term, a rank-0 centered Gram + (zero within-control covariate variation, e.g. a constant as the only + covariate) is the benign identified-subset case — the estimation-effect + correction is exactly zero and the fit collapses to the no-covariate reg fit + with finite SEs (`_centered_or_bread` maps `_safe_inv`'s all-NaN rank-0 + sentinel to a zero correction; the aggregate rank-guard warning still fires). + This differs from the `[1, X]` breads (ipw PS Hessian, dr), where all-NaN is + a true pathology that NaN-propagates. - All aggregation SEs (simple, event study) include the weight influence function (WIF) adjustment, matching R's `did::aggte()`. The WIF accounts for uncertainty in estimating group-size aggregation weights. Group aggregation uses equal time weights (deterministic), diff --git a/tests/test_staggered.py b/tests/test_staggered.py index d4b0b34e5..26f8c72a9 100644 --- a/tests/test_staggered.py +++ b/tests/test_staggered.py @@ -1655,16 +1655,20 @@ def test_clustered_constant_covariate_finite_se(self): assert np.isfinite(res.overall_se) assert res.overall_se < 1.0 - @pytest.mark.parametrize("method", ["reg", "ipw", "dr"]) + @pytest.mark.parametrize("method", ["ipw", "dr"]) def test_rank0_bread_propagates_nan_not_zero(self, monkeypatch, method): # rank-0 is unreachable through covariates alone (the always-present # intercept guarantees rank >= 1), so simulate an all-NaN bread to # exercise the NaN-masking fix: var_psi becomes NaN and must yield a NaN - # SE, NOT 0.0 via the old ``var_psi > 0 else 0.0`` guard. reg/ipw now - # derive their per-cell SE from _safe_inv-based IF terms too (OLS - # estimation-effect bread / PS Hessian), so all three methods share - # the contract. The point estimate does NOT depend on the bread, so - # it stays finite (NaN inference on an estimable cell, not _nan_gt_entry). + # SE, NOT 0.0 via the old ``var_psi > 0 else 0.0`` guard. ipw's PS + # Hessian and dr's breads are over [1, X], so all-NaN there is a true + # pathology that must propagate. reg is deliberately EXCLUDED: its + # estimation-effect bread is over the CENTERED covariate Gram (the + # intercept is handled analytically), where rank-0 is the benign + # constant-covariate case mapped to a zero correction — see + # test_reg_constant_only_covariate_matches_no_covariate below. + # The point estimate does NOT depend on the bread, so it stays + # finite (NaN inference on an estimable cell, not _nan_gt_entry). from tests.conftest import assert_nan_inference import diff_diff.staggered as staggered_mod @@ -2020,6 +2024,73 @@ def test_universal_base_period_anticipation_reg_smoke(self): assert np.isfinite(cell["se"]), f"cell ({g},{t})" assert np.isfinite(res.overall_se) + def test_reg_constant_only_covariate_matches_no_covariate(self): + """A constant as the ONLY reg covariate makes the CENTERED + estimation-effect Gram rank-0 (all-zero). The correction on the + identified (intercept-only) subset is exactly zero, so effects AND + SEs must equal the no-covariate fit - finite, never NaN (the rank-0 + centered bread maps to a zero correction, not an all-NaN inverse).""" + data = generate_staggered_data_with_covariates(seed=789) + data["xc"] = 5.0 + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + no_cov = CallawaySantAnna(estimation_method="reg").fit( + data, "outcome", "unit", "time", "first_treat" + ) + const_only = CallawaySantAnna(estimation_method="reg").fit( + data, "outcome", "unit", "time", "first_treat", covariates=["xc"] + ) + assert np.isfinite(const_only.overall_se) + np.testing.assert_allclose( + const_only.overall_att, no_cov.overall_att, rtol=1e-12 + ) + np.testing.assert_allclose( + const_only.overall_se, no_cov.overall_se, rtol=1e-9 + ) + for key, cell in const_only.group_time_effects.items(): + ref = no_cov.group_time_effects[key] + np.testing.assert_allclose(cell["effect"], ref["effect"], rtol=1e-12) + np.testing.assert_allclose(cell["se"], ref["se"], rtol=1e-9) + + def test_reg_constant_only_covariate_matches_no_covariate_survey(self): + """Survey-weighted twin of the constant-only-covariate case: the + weighted centered Gram is also rank-0, and the general + (survey-branch) producer must likewise collapse to the + no-covariate survey fit with finite SEs.""" + from diff_diff.survey import SurveyDesign + + rng = np.random.default_rng(17) + data = generate_staggered_data_with_covariates(seed=789) + data["xc"] = 5.0 + weights = pd.DataFrame( + { + "unit": data["unit"].unique(), + "weight": rng.uniform(0.5, 2.0, size=data["unit"].nunique()), + } + ) + data = data.merge(weights, on="unit") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + no_cov = CallawaySantAnna(estimation_method="reg").fit( + data, "outcome", "unit", "time", "first_treat", + survey_design=SurveyDesign(weights="weight"), + ) + const_only = CallawaySantAnna(estimation_method="reg").fit( + data, "outcome", "unit", "time", "first_treat", + covariates=["xc"], + survey_design=SurveyDesign(weights="weight"), + ) + assert np.isfinite(const_only.overall_se) + np.testing.assert_allclose( + const_only.overall_att, no_cov.overall_att, rtol=1e-12 + ) + np.testing.assert_allclose( + const_only.overall_se, no_cov.overall_se, rtol=1e-9 + ) + for key, cell in const_only.group_time_effects.items(): + ref = no_cov.group_time_effects[key] + np.testing.assert_allclose(cell["se"], ref["se"], rtol=1e-9) + def test_uniform_survey_weights_match_unweighted_per_cell_se(self): """Uniform survey weights route reg+cov through the general (survey-branch) producer while the unweighted fit takes the From e226c0f35b012079c0da3d45f964dc4a1afbd2ce Mon Sep 17 00:00:00 2001 From: igerber Date: Sun, 5 Jul 2026 10:45:56 -0400 Subject: [PATCH 3/3] docs(staggered): update rank-0 comments for centered reg bread semantics CI review P3: two comments still said "a rank-0 bread yields NaN" on the reg paths, contradicting _centered_or_bread's zero-correction collapse. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01X4AzrFUMqJxcUumSH31mSr --- diff_diff/staggered.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/diff_diff/staggered.py b/diff_diff/staggered.py index 1662184fe..1d4ef2476 100644 --- a/diff_diff/staggered.py +++ b/diff_diff/staggered.py @@ -1609,8 +1609,10 @@ def _compute_all_att_gt_covariate_reg( inf_control = -residuals * proj_c # Per-cell SE from the same IF that feeds aggregation # (sqrt(sum(phi^2)), the DR convention). A rank-0 - # bread yields NaN — safe_inference then NaNs the - # full inference tuple (never a silent 0.0). + # CENTERED Gram is benign — _centered_or_bread maps it + # to a zero correction (finite SE on the identified + # intercept subset); the NaN guard below covers only + # non-finite inputs, never a silent 0.0. var_psi = float(np.sum(inf_treated**2) + np.sum(inf_control**2)) if not np.isfinite(var_psi): se = float("nan") @@ -2726,8 +2728,10 @@ def _outcome_regression( inf_func = np.concatenate([inf_treated, inf_control]) # Per-cell SE from the same IF that feeds aggregation - # (sqrt(sum(phi^2)), the DR convention). A rank-0 bread yields - # NaN — safe_inference then NaNs the full inference tuple. + # (sqrt(sum(phi^2)), the DR convention). A rank-0 CENTERED Gram + # is benign — _centered_or_bread maps it to a zero correction + # (finite SE on the identified intercept subset); the NaN guard + # below covers only non-finite inputs. var_psi = float(np.sum(inf_func**2)) if not np.isfinite(var_psi): se = float("nan")