diff --git a/CHANGELOG.md b/CHANGELOG.md index cfccde14..4c7bcc8a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,32 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 are unchanged; the renamed section fragments are kept backward-compatible by a client-side redirect (`docs/_static/legacy-fragment-redirect.js`) that rewrites old `#3.-Fit-Event-Study`-style deep links to their unnumbered targets. +- **Unified event-study surface: full VCV persistence for StackedDiD/TwoStageDiD + + per-row `df` provenance (4.0 program Phase 2, ledger row M-092 follow-up).** + (a) `StackedDiDResults` and `TwoStageDiDResults` gain `event_study_vcov` / + `event_study_vcov_index`: the full event-study covariance both estimators already + computed internally (pooled-regression sub-block / Gardner-GMM matrix) is now + persisted and threaded onto `EventStudyResults.vcov`, mirroring the existing + CallawaySantAnna/SunAbraham fields. StackedDiD persists in every inference mode (its + reported ES SEs are always the matrix diagonal); TwoStageDiD persists on the + analytical paths and ships `None` under bootstrap (percentile inference, no + covariance) and replicate-weight survey designs (mixed-layout replicate VCV - CS + precedent). (b) `EventStudyResults.df` is now PER-ROW (float array aligned to + `event_time`, NaN where no df governed the stored p-value/CI; scalar inputs + broadcast), and joins the pinned `to_dataframe()` schema between `n` and + `is_reference` - an unreleased-schema change, made before `aggregate()` exposes the + surface publicly. Three producers genuinely need per-row df: StackedDiD `hc2_bm` + (per-event Bell-McCaffrey Satterthwaite df), MultiPeriodDiD `hc2_bm` (per-period BM + df - its `inference_df` field stores only the post-average contrast df), and LPDiD + (per-horizon realized cluster df). (c) New `event_study_df` provenance channels + record the df each event-study row's `safe_inference` actually received: + `CallawaySantAnnaResults` (scalar; G-1 on bare-cluster fits, the conservative + min-across-horizons effective df on explicit-survey fits - `df_inference` and its + narrow HonestDiD contract are untouched), `MultiPeriodDiDResults` (per-period dict), + `StackedDiDResults` (per-event dict), `TwoStageDiDResults` (scalar survey df), + `LPDiDResults` (per-horizon dict). All channels clear under bootstrap overrides. + SunAbraham/dCDH df threading is tracked in TODO.md; their surfaces show df=NaN + ("unexposed"). - **ImputationDiD leave-one-out SE now anchored against Stata `did_imputation, leaveout` (no library behavior change).** The Borusyak-Jaravel-Spiess (2024) App. A.9 LOO variance (`leave_one_out=True`) has no runnable R reference (R `didimputation` omits LOO), so it @@ -56,6 +82,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 **Estimator equations, weighting, variance, and numerical output are unchanged.** ### Changed +- **DiagnosticReport's parallel-trends check upgrades from the Bonferroni fallback to + the joint Wald test for StackedDiD and TwoStageDiD event-study fits.** The check's + joint-Wald path runs whenever the source fit exposes `event_study_vcov`; now that + both estimators persist theirs (see Added), their pre-period tests use the full + covariance instead of per-coefficient Bonferroni. All-filtered TwoStageDiD horizons + (NaN VCV rows) are excluded from the tested family by the existing undefined- + inference collector, so the joint statistic is never computed through NaN entries. + Two accompanying guards keep the new path methodology-safe: (1) the joint-Wald + path now fail-closes on provenance - if any retained pre-period row carries a + non-finite p-value (hc2_bm Bell-McCaffrey DOF failure, collapsed replicate-survey + df), the check returns `inconclusive` instead of publishing a finite verdict + through inference the estimator refused to produce; (2) `vcov_type="hc2_bm"` + source fits skip the generic chi-square joint Wald entirely (it would discard the + CR2/BM small-sample correction; the multi-constraint AHT/HTZ CR2 test is tracked + in DEFERRED.md) and run Bonferroni over the BM-adjusted per-row p-values. See + REPORTING.md "hc2_bm parallel-trends policy". - **Diagnostic input validation in the report consumers (4.0 program Phase 2).** `BusinessReport`, `practitioner_next_steps`, and `DiagnosticReport` now route diagnostics by the new `Diagnostic` marker instead of by result-class name. Passing a diff --git a/DEFERRED.md b/DEFERRED.md index fd25dadd..1c1d43d4 100644 --- a/DEFERRED.md +++ b/DEFERRED.md @@ -48,6 +48,7 @@ exists but parity can't be verified without a local toolchain. | `StaggeredTripleDifference` per-cohort group-effect SEs include WIF (conservative vs R's `wif=NULL`); documented in REGISTRY. Could override the mixin for an exact R match (verification needs R `triplediff`). | `staggered_triple_diff.py` | #245 | Low | | **WooldridgeDiD follow-up cluster** (PR-B Stage D/E fail-closed surfaces; re-enable after R/Stata validation):
• QMLE sandwich uses `aweight` cluster adjustment `(G/(G-1))·(n-1)/(n-k)` vs Stata's `G/(G-1)` (conservative); add a `qmle` weight type if Stata goldens confirm a material difference (`wooldridge.py`, `linalg.py`).
• response-scale APE / log-link coefficient bridge for R `etwfe(family=poisson|logit)` cell-level parity — needs `emfx()` APE extraction or link-inversion with baseline-mean adjustment (`generate_wooldridge_golden.R`, `test_methodology_wooldridge.py`).
• `aggregate(weights="cohort_share")` on survey-weighted fits: `_n_g_per_cohort` uses raw `unit.nunique()`; implement design-weighted unit totals per cohort (paper W2025 §7) and lift the `ValueError` gate (`wooldridge.py`, `wooldridge_results.py`).
• unconditional inference for `cohort_share` accounting for ω̂_g sampling uncertainty (W2025 §7.5); currently NaN-closed (`wooldridge_results.py`).
• `cohort_trends=True × survey_design` and `× control_group="never_treated"` raise `NotImplementedError` (unvalidated TSL variance / trend columns spanned by placebo cell-dummies) (`wooldridge.py`).
• Stata `jwdid` golden-value `TestReferenceValues` (`tests/test_wooldridge.py`). | `wooldridge.py`, `wooldridge_results.py`, `linalg.py`, benchmarks | #216 · PR-B | Med-Low | | Extend `WooldridgeDiD` `method ∈ {logit, poisson}` with `vcov_type ∈ {classical, hc2, hc2_bm}`: composing HC2 leverage + Bell-McCaffrey DOF with the QMLE pseudo-residual sandwich needs derivation + R parity vs `clubSandwich::vcovCR(glm, type="CR2")`. Rejected at `__init__`. | `wooldridge.py` | follow-up | Medium | +| Multi-constraint CR2 parallel-trends test (AHT/HTZ) for `hc2_bm` fits: DiagnosticReport's PT check routes `vcov_type="hc2_bm"` sources to Bonferroni over the BM-adjusted per-row p-values because the generic chi-square joint Wald would discard the CR2 small-sample correction (see REPORTING.md "hc2_bm parallel-trends policy"). The proper joint test is the AHT/HTZ Wald with a Satterthwaite-style denominator df over the pre-period contrast block; needs derivation for the stacked/pooled WLS-CR2 layout + parity vs `clubSandwich::Wald_test(..., test="HTZ")`. | `diagnostic_report.py`, `linalg.py` | vcov/df round-trip PR | Low | | `PreTrendsPower` CS/SA `anticipation=1` R-parity fixture: R `pretrends` has no anticipation parameter, so the Python `_extract_pre_period_params` anticipation filter isn't R-parity-locked. Build a synthetic CS/SA result with `anticipation=1` and assert γ_p matches R's `slope_for_power()`. (Mechanism already covered by MC + full-VCV tests.) | `tests/test_methodology_pretrends.py`, `generate_pretrends_golden.R` | PR-C | Low | | Harmonize SunAbraham's HC1 within-transform finite-sample correction with `fixest::sunab()` — SA applies `n/(n-k_dm)`, fixest applies `n/(n-k_total)` (counts absorbed FE); ~1-2% SE difference, documented as a "Deviation from R" and pinned at `atol=5e-3`. Either thread `df_adjustment` or keep as an intentional, R-verified difference. | `sun_abraham.py`, `linalg.py` | follow-up | Low | | Absorbed-FE **clustered** CR1 with *non-nested* FE: for `absorb=[FE1,FE2], cluster=FE1` (e.g. `absorb=["unit","time"], cluster="unit"`), `fixest` counts the non-nested FE (time) in the CR1 `(n-1)/(n-k)` finite-sample denominator, but the clustered path uses only `k_visible`. D4 harmonized the *non-clustered* classical/hc1 full-K scale (`_absorbed_fe_vcov_scale`) and left the clustered path unchanged — correct for FE nested in the cluster, a small deviation for non-nested FE (documented in REGISTRY within-transform note). Thread a non-nested `df_adjustment` into the clustered CR1 factor; verify vs `fixest::feols(..., cluster=)`. | `linalg.py`, `estimators.py` | SE-audit D4 | Low | diff --git a/TODO.md b/TODO.md index a7fbf00c..2de0bd9f 100644 --- a/TODO.md +++ b/TODO.md @@ -23,7 +23,7 @@ Related tracking surfaces: |-------|----------|--------|--------|----------| | `SyntheticControl` conformal (CWZ 2021) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies — the residual-permutation shortcut is only valid for time-permutation-invariant proxies (SC/Lasso/DiD); an AR proxy needs innovation permutation. | `diff_diff/conformal.py`, `diff_diff/synthetic_control_results.py` | CWZ-2021 | Heavy | Low | | `ContinuousDiD` CGBS-2024 remaining extensions (earlier phases — `covariates=` reg/dr, `treatment_type="discrete"`, single-cohort `control_group="lowest_dose"` with estimand `ATT(d)−ATT(d_L)` — are already supported; see REGISTRY Note #7). Remaining (all deferred `NotImplementedError`, documented): `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate); `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF); multi-cohort **heterogeneous-support** discrete aggregation (support-aware: average each dose only over the cohorts that observe it); **multi-cohort `lowest_dose`** (within-cohort `d_L` reference + support-aware cross-cohort aggregation); and **`covariates=` × `lowest_dose`** (conditional-PT-relative-to-`d_L` estimand). Single-cohort / 2-period / shared-support multi-cohort are supported. | `continuous_did.py` | CGBS-2024 | Heavy | Low | -| `EventStudyResults` full-vcov + df round-trip for more producers (4.0 program Phase 2 follow-up). StackedDiD and TwoStageDiD compute a full event-study VCV (regression / Gardner GMM) internally but retain only marginal SEs on their result containers, so their unified surfaces carry `vcov=None`; persist the event-time submatrix + ordered horizon index on those containers and thread it in. Likewise thread producer-specific inference df into `EventStudyResults.df` (MPD `inference_df`, explicit-survey CS `survey_metadata.df_survey`, HC2-BM horizon-specific contrast df). The adapter is faithful to what containers expose today; this is additive container work, not a defect. | `diff_diff/results_base.py`, `diff_diff/stacked_did_results.py`, `diff_diff/two_stage_results.py` | #M-092 | Mid | Low | +| `event_study_df` provenance for the remaining producers (4.0 program Phase 2 follow-up; the Stacked/TwoStage vcov + CS/MPD/Stacked/TwoStage/LPDiD df round-trip landed with the per-row `EventStudyResults.df` schema). Remaining: SunAbraham (per-event BM Satterthwaite df under `hc2_bm` at `sun_abraham.py:1206`, scalar survey df at `:1256` - mechanism identical to StackedDiD's recorded-at-use dict); dCDH (its event-study inference receives a `df_inference` argument, `chaisemartin_dhaultfoeuille.py:6311`, never stored - needs tracing to the stored rows); LPDiD pooled-window df (popped and discarded at fit assembly). All purely additive: their surfaces show df=NaN ("unexposed") until threaded. | `diff_diff/sun_abraham.py`, `diff_diff/chaisemartin_dhaultfoeuille.py`, `diff_diff/lpdid.py` | #M-092 | Quick | Low | ### Performance diff --git a/diff_diff/diagnostic_report.py b/diff_diff/diagnostic_report.py index 5527595c..c33a5d6a 100644 --- a/diff_diff/diagnostic_report.py +++ b/diff_diff/diagnostic_report.py @@ -1398,24 +1398,72 @@ def _pt_event_study(self) -> Dict[str, Any]: # 2. ``event_study_vcov_index`` + ``event_study_vcov`` — the # CallawaySantAnnaResults convention, where the event-study # covariance is stored separately from the full regression vcov. + # Fail-closed inference guard for the JOINT-WALD path. The Wald + # statistic consumes only effects + covariance, so with a persisted + # vcov it would silently bypass the estimator's own fail-closed + # inference: a StackedDiD hc2_bm row whose Bell-McCaffrey DOF is + # unavailable (or a replicate-survey row whose effective df + # collapsed) carries a FINITE effect/SE and finite covariance + # entries but deliberately-NaN t/p/CI. The retained-row collector + # only drops non-finite effect/SE rows, and the nan-p guard below + # runs only on the Bonferroni fallback - so without this check a + # covariance-backed fit publishes a finite joint p-value for + # inference the estimator explicitly refused to compute. Same + # contract as the round-33/42 Bonferroni guards, applied BEFORE any + # covariance path is selected. + _n_nonfinite_p = sum( + 1 for (_, _, _, p) in pre_coefs if not (isinstance(p, (int, float)) and np.isfinite(p)) + ) + if _n_nonfinite_p > 0: + return { + "status": "ran", + "method": "inconclusive", + "joint_p_value": None, + "test_statistic": None, + "df": len(pre_coefs), + "n_pre_periods": len(pre_coefs), + "n_dropped_undefined": _n_nonfinite_p, + "per_period": per_period, + "verdict": "inconclusive", + "reason": ( + f"{_n_nonfinite_p} retained pre-period coefficient(s) " + "have non-finite per-period p-value: the source " + "estimator's inference failed closed (e.g. hc2_bm " + "Bell-McCaffrey DOF unavailable, or collapsed " + "replicate-survey df). A joint Wald over the persisted " + "covariance would silently convert that undefined " + "inference into a finite verdict; the PT test is " + "inconclusive on this fit." + ), + } vcov_for_wald: Optional[Any] = None idx_map_for_wald: Optional[Any] = None vcov_method_tag = "joint_wald" - if vcov is not None and interaction_indices is not None: - vcov_for_wald = vcov - idx_map_for_wald = interaction_indices - else: - es_vcov = getattr(r, "event_study_vcov", None) - es_vcov_index = getattr(r, "event_study_vcov_index", None) - if es_vcov is not None and es_vcov_index is not None: - vcov_for_wald = es_vcov - # ``event_study_vcov_index`` is an ordered list of relative-time - # keys; convert it into a dict mapping key -> position. - try: - idx_map_for_wald = {k: i for i, k in enumerate(es_vcov_index)} - vcov_method_tag = "joint_wald_event_study" - except TypeError: - idx_map_for_wald = None + # hc2_bm-sourced fits skip the covariance-based joint Wald entirely: + # the generic chi-square (or F(k, df_survey)) reference ignores the + # CR2 / Bell-McCaffrey small-sample correction the estimator's own + # per-row inference applies, and the proper multi-constraint CR2 + # test (AHT/HTZ, clubSandwich::Wald_test(test="HTZ")) is not + # implemented (tracked in DEFERRED.md). Bonferroni over the + # BM-adjusted per-row p-values keeps the small-sample correction. + # See REGISTRY.md StackedDiD / REPORTING.md PT notes. + if getattr(r, "vcov_type", None) != "hc2_bm": + if vcov is not None and interaction_indices is not None: + vcov_for_wald = vcov + idx_map_for_wald = interaction_indices + else: + es_vcov = getattr(r, "event_study_vcov", None) + es_vcov_index = getattr(r, "event_study_vcov_index", None) + if es_vcov is not None and es_vcov_index is not None: + vcov_for_wald = es_vcov + # ``event_study_vcov_index`` is an ordered list of + # relative-time keys; convert it into a dict mapping + # key -> position. + try: + idx_map_for_wald = {k: i for i, k in enumerate(es_vcov_index)} + vcov_method_tag = "joint_wald_event_study" + except TypeError: + idx_map_for_wald = None df_denom: Optional[float] = None if vcov_for_wald is not None and idx_map_for_wald is not None and df > 0: try: @@ -1425,7 +1473,33 @@ def _pt_event_study(self) -> Dict[str, Any]: beta_map = {k: eff for (k, eff, _, _) in pre_coefs} beta = np.array([beta_map[k] for k in keys_in_vcov], dtype=float) v_sub = np.asarray(vcov_for_wald)[np.ix_(idx, idx)] - stat = float(beta @ np.linalg.solve(v_sub, beta)) + # Rank/PSD guard. A cluster-sandwich + # covariance has rank bounded by the cluster-score + # dimension, so a singular pre-period block is a NORMAL + # edge case whenever the number of tested leads + # approaches the number of clusters. Solving through a + # singular / non-PSD block produces a garbage quadratic + # form (observed: stat ~ -1e15 -> chi2 p=1.0 -> a false + # stakeholder-facing "no violation" verdict). Require a + # finite block that is numerically full-rank and PSD + # (matrix_rank's standard tolerance: k * eps * largest + # singular value); otherwise raise into the Bonferroni + # fallback over the (finite, per Guard above) per-row + # p-values. Belt-and-braces: a Wald statistic is a PSD + # quadratic form, so any non-finite or negative result + # is invalid regardless of how it arose. + if not np.isfinite(v_sub).all(): + raise ValueError("non-finite covariance block") + v_sym = 0.5 * (v_sub + v_sub.T) + if np.linalg.matrix_rank(v_sym) < df: + raise ValueError("rank-deficient covariance block") + _eigs = np.linalg.eigvalsh(v_sym) + _scale = float(np.max(np.abs(_eigs))) if _eigs.size else 0.0 + if _scale <= 0.0 or float(_eigs.min()) < -df * np.finfo(float).eps * _scale: + raise ValueError("covariance block is not positive semi-definite") + stat = float(beta @ np.linalg.solve(v_sym, beta)) + if not np.isfinite(stat) or stat < 0.0: + raise ValueError("invalid Wald statistic") # Round-27 P1 CI review on PR #318: survey-backed # fits carry a finite ``df_survey`` on @@ -1698,8 +1772,10 @@ def _infer_cov_source(source_fit: Any) -> str: this fallback entirely. - ``"diag_fallback"`` — event-study result types with ``event_study_vcov is None`` (bootstrap or replicate-weight - CS / SA fits, plus ImputationDiD / Stacked / EfficientDiD / - TwoStageDiD / etc. which don't yet expose ``event_study_vcov``); + CS / SA / TwoStageDiD fits, plus ImputationDiD / EfficientDiD / + etc. which don't yet expose ``event_study_vcov``; StackedDiD + persists its VCV in every inference mode, so it reaches this + fallback only when no event study was requested); OR ``MultiPeriodDiDResults`` without ``interaction_indices`` (genuine diag-only path inside ``pretrends.py:_extract_pre_period_params``, no "available but unused" concern, so no downgrade applies). diff --git a/diff_diff/estimators.py b/diff_diff/estimators.py index 3e01827c..1f4043ec 100644 --- a/diff_diff/estimators.py +++ b/diff_diff/estimators.py @@ -2288,6 +2288,13 @@ def _refit_mp_absorb(w_r): period_effects = {} post_effect_values = [] post_effect_indices = [] + # Per-period df PROVENANCE for the unified event-study surface + # (row M-092): the df actually handed to each period's + # safe_inference below. Distinct from `inference_df`, which stores + # the POST-AVERAGE contrast df - under hc2_bm the per-period BM + # DOFs differ from it, so broadcasting inference_df would + # misrepresent the per-period rows. + es_df_used: Dict[Any, float] = {} assert vcov is not None for period in non_ref_periods: @@ -2305,6 +2312,11 @@ def _refit_mp_absorb(w_r): else: period_df = df t_stat, p_value, conf_int = safe_inference(effect, se, alpha=self.alpha, df=period_df) + es_df_used[period] = ( + float(period_df) + if period_df is not None and np.isfinite(period_df) and period_df > 0 + else float("nan") + ) period_effects[period] = PeriodEffect( period=period, @@ -2400,6 +2412,7 @@ def _refit_mp_absorb(w_r): if _avg_df is not None and np.isfinite(_avg_df) and _avg_df > 0 else None ), + event_study_df=es_df_used, ) self._coefficients = coefficients diff --git a/diff_diff/lpdid.py b/diff_diff/lpdid.py index 7395f4bf..eef33131 100644 --- a/diff_diff/lpdid.py +++ b/diff_diff/lpdid.py @@ -588,6 +588,7 @@ def _estimate_regression_adjustment_sample( "conf_high": np.nan, "n_obs": n_obs, "n_clusters": np.nan, + "df": np.nan, } if n_obs == 0 or sample["_entry"].nunique() < 2: return empty_result @@ -725,6 +726,9 @@ def _estimate_regression_adjustment_sample( "conf_high": conf_int[1], "n_obs": n_obs, "n_clusters": n_clusters, + # df provenance for the unified surface: the exact value handed + # to safe_inference above (NaN when None -> normal theory). + "df": float(df) if df is not None else np.nan, } def _estimate_sample( @@ -756,6 +760,7 @@ def _estimate_sample( "conf_high": np.nan, "n_obs": n_obs, "n_clusters": np.nan, + "df": np.nan, } if n_obs == 0 or sample["_entry"].nunique() < 2: @@ -848,6 +853,7 @@ def _estimate_sample( "conf_high": np.nan, "n_obs": n_obs, "n_clusters": len(pd.unique(cluster_ids)), + "df": np.nan, } use_cluster_vcov = len(pd.unique(cluster_ids)) >= 2 @@ -899,6 +905,9 @@ def _estimate_sample( "conf_high": conf_int[1], "n_obs": n_obs, "n_clusters": n_clusters, + # df provenance for the unified surface: the exact value handed + # to safe_inference above (NaN when None -> normal theory). + "df": float(df) if df is not None else np.nan, } def _estimate_survey_sample(self, sample, design, response, column_names, n_obs, survey_design): @@ -984,6 +993,10 @@ def _estimate_survey_sample(self, sample, design, response, column_names, n_obs, "conf_high": conf_int[1], "n_obs": n_obs, "n_clusters": int(resolved.n_psu), + # df provenance for the unified surface: recorded iff the value + # governed a t-reference (finite, > 0; the df<=0 sentinel and + # None both yield non-t inference). + "df": (float(df) if df is not None and np.isfinite(df) and df > 0 else np.nan), } def _estimate_horizon( @@ -1446,8 +1459,10 @@ def fit( event_study = None pooled = None + event_study_df: Optional[Dict[int, float]] = None if not only_pooled: event_rows = [] + event_study_df = {} for horizon in range(-self.pre_window, self.post_window + 1): if horizon == -1: estimate = { @@ -1472,6 +1487,12 @@ def fit( absorb=absorb, survey_design=survey_design, ) + # Pop the df provenance BEFORE building the frame row so the + # native event_study schema is unchanged; the synthetic + # horizon == -1 base row carries no df key and gets no entry. + _df_h = estimate.pop("df", None) + if _df_h is not None: + event_study_df[horizon] = float(_df_h) event_rows.append( { "horizon": horizon, @@ -1507,6 +1528,11 @@ def fit( absorb=absorb, survey_design=survey_design, ) + # Discard the pooled windows' df provenance (native pooled frame + # schema unchanged; pooled-window df threading is a tracked + # follow-up in TODO.md). + pre_estimate.pop("df", None) + post_estimate.pop("df", None) pooled = pd.DataFrame( [ { @@ -1539,6 +1565,7 @@ def fit( self.results_ = LPDiDResults( event_study=event_study, pooled=pooled, + event_study_df=event_study_df, n_obs=len(data), n_treated_units=int(treatment_by_unit.gt(0).sum()), n_control_units=int(treatment_by_unit.eq(0).sum()), diff --git a/diff_diff/lpdid_results.py b/diff_diff/lpdid_results.py index 49c1f142..7de1bd71 100644 --- a/diff_diff/lpdid_results.py +++ b/diff_diff/lpdid_results.py @@ -53,6 +53,16 @@ class LPDiDResults(BaseResults): survey_metadata: Optional[Any] = None n_strata: Optional[int] = None n_psu: Optional[int] = None + # event_study_df (spec section 5, row M-092): per-horizon df PROVENANCE - + # maps each estimated horizon to the df actually passed to its + # safe_inference (realized ``n_clusters - 1`` on the cluster/RA paths, + # the per-horizon-sample survey df ``n_PSU - n_strata`` under a survey + # design; NaN when the row used normal theory or an undefined df). Not + # derivable from the frames: the survey rule needs the per-sample + # n_strata and the cluster rule is gated on a successful vcov, neither + # of which the native tables record. The synthetic ``horizon == -1`` + # base row has no entry. None when no event study was fit. + event_study_df: Optional[Dict[int, float]] = None # ------------------------------------------------------------------ # internal helpers diff --git a/diff_diff/results.py b/diff_diff/results.py index bb39bf46..20a7993e 100644 --- a/diff_diff/results.py +++ b/diff_diff/results.py @@ -713,6 +713,15 @@ class MultiPeriodDiDResults(BaseResults): # p-value/CI without refitting. df_convention: Optional[str] = field(default=None) inference_df: Optional[float] = field(default=None) + # event_study_df (spec section 5, row M-092): PER-PERIOD df provenance - + # maps each estimated period to the df actually passed to that period's + # safe_inference (the period's OWN Bell-McCaffrey DOF under hc2_bm, the + # shared analytical df otherwise; NaN when inference used normal theory + # or an undefined df). Distinct from `inference_df` above, which stores + # the POST-AVERAGE contrast df - under hc2_bm the two genuinely differ, + # so per-period consumers must read this dict, never broadcast + # inference_df. + event_study_df: Optional[Dict[Any, float]] = field(default=None, repr=False) # --- Inference-field aliases (balance/external-adapter compatibility) --- @property diff --git a/diff_diff/results_base.py b/diff_diff/results_base.py index 6bb451f6..51e8e307 100644 --- a/diff_diff/results_base.py +++ b/diff_diff/results_base.py @@ -21,7 +21,7 @@ from __future__ import annotations from dataclasses import dataclass, field -from typing import Any, Dict, List, Optional, Tuple +from typing import Any, Dict, List, Optional, Tuple, Union, cast import numpy as np import pandas as pd @@ -125,6 +125,7 @@ def _json_safe_label(value: Any) -> Any: "cband_lower", "cband_upper", "n", + "df", "is_reference", ) @@ -188,11 +189,12 @@ class EventStudyResults(BaseResults): renumbered. vcov : np.ndarray or None Full event-study variance-covariance matrix where the RESULT - CONTAINER exposes one (today: CallawaySantAnna, SunAbraham, - MultiPeriodDiD), ordered by ``vcov_index``. StackedDiD and - TwoStageDiD compute a full VCV internally but retain only marginal - SEs on their result containers, so their surfaces carry ``vcov=None`` - until that retention is added (tracked in TODO.md). + CONTAINER exposes one (e.g. CallawaySantAnna, SunAbraham, + MultiPeriodDiD, StackedDiD, and TwoStageDiD's analytical modes), + ordered by ``vcov_index``. None when the producer records no matrix + or when the stored SEs are no longer its diagonal (bootstrap and + replicate-weight overrides clear it rather than ship an + inconsistent matrix). vcov_index : np.ndarray or None ``event_time`` labels labelling ``vcov``'s rows/columns (explicit ordering for HonestDiD / PreTrendsPower consumption). @@ -206,9 +208,18 @@ class EventStudyResults(BaseResults): Significance level of the stored intervals. source : str or None Producing results-class name (provenance). - df : float or None - Inference degrees of freedom threaded from the source, where the - producer exposes one. + df : float, np.ndarray, or None + Per-row inference degrees of freedom: ``df[i]`` is the df ACTUALLY + passed to ``safe_inference`` for row i's stored p-value/CI, threaded + from the producer. Accepts None (no df exposed -> all-NaN column), a + scalar (broadcast to every row - e.g. CallawaySantAnna, whose + explicit-survey event study applies ONE conservative df, the minimum + per-horizon effective df, to all rows), or a length-n array (per-row + producers: StackedDiD ``hc2_bm`` per-event Bell-McCaffrey df, LPDiD + per-horizon cluster df, MultiPeriodDiD ``hc2_bm`` per-period df). + NaN on any row means normal-theory inference, an undefined df, + bootstrap-overridden inference, or a producer that records none; + reference rows and rows with NaN p-values are always NaN. """ event_time: np.ndarray @@ -231,7 +242,7 @@ class EventStudyResults(BaseResults): cband_crit_value: Optional[float] = None alpha: float = 0.05 source: Optional[str] = None - df: Optional[float] = None + df: Optional[Union[float, np.ndarray]] = None _ARRAY_FIELDS = ( "att", @@ -318,6 +329,28 @@ def __post_init__(self) -> None: if arr is not None: arr[self.is_reference] = np.nan + # df is per-row PROVENANCE: the df actually passed to safe_inference + # for that row's stored p-value/CI. Normalize None (no df exposed) to + # an all-NaN column and broadcast scalars (single-df producers), then + # NaN out rows that carry no safe_inference output - reference rows + # and rows whose stored p-value is NaN (non-estimable horizons) never + # used a df. Runs AFTER reference-row normalization so p_value is + # final; explicit block rather than _ARRAY_FIELDS because df, like + # the band columns, is optional. + if self.df is None: + df_arr = np.full(n_rows, np.nan) + elif np.ndim(self.df) == 0: + df_arr = np.full(n_rows, float(cast(float, self.df))) + else: + df_arr = np.array(self.df, dtype=float) + if df_arr.shape != (n_rows,): + raise ValueError( + f"EventStudyResults field 'df' has shape {df_arr.shape}; " + f"expected ({n_rows},) to align with event_time." + ) + df_arr[~np.isfinite(self.p_value)] = np.nan + self.df = df_arr + if (self.vcov is None) != (self.vcov_index is None): raise ValueError( "EventStudyResults requires vcov and vcov_index together " @@ -375,6 +408,9 @@ def to_dataframe(self) -> pd.DataFrame: "cband_lower": self.cband_lower if self.cband_lower is not None else nan_col, "cband_upper": self.cband_upper if self.cband_upper is not None else nan_col, "n": self.n, + # __post_init__ guarantees an array; cast narrows the + # scalar-accepting field type for mypy. + "df": cast(np.ndarray, self.df), "is_reference": self.is_reference, }, columns=list(EVENT_STUDY_SCHEMA), @@ -416,7 +452,7 @@ def to_dict(self) -> Dict[str, Any]: "cband_crit_value": self.cband_crit_value, "alpha": self.alpha, "source": self.source, - "df": self.df, + "df": cast(np.ndarray, self.df).tolist(), } return out @@ -657,6 +693,23 @@ def _from_relative_dict(results: Any) -> EventStudyResults: if vcov is None or vcov_index is None: vcov = vcov_index = None + # Per-row df provenance. Primary channel: `event_study_df` (scalar or + # {event_time: df} dict of the values actually passed to safe_inference; + # producers clear it when bootstrap overrides the stored inference). + # Fallback: CallawaySantAnna's bare-cluster `df_inference` - GATED on no + # bootstrap, because bare-cluster fits populate df_inference=G-1 even + # when n_bootstrap>0 replaced the stored ES p/CIs with percentile values + # that never used that df (and bootstrap p-values are finite, so the + # container's NaN-p masking cannot catch it). + df_src: Any = getattr(results, "event_study_df", None) + if df_src is None and getattr(results, "bootstrap_results", None) is None: + df_src = getattr(results, "df_inference", None) + df_arg: Optional[Union[float, np.ndarray]] + if isinstance(df_src, dict): + df_arg = np.array([float(df_src[k]) if k in df_src else np.nan for k in keys]) + else: + df_arg = df_src + return EventStudyResults( event_time=event_time, att=att, @@ -677,7 +730,7 @@ def _from_relative_dict(results: Any) -> EventStudyResults: cband_crit_value=getattr(results, "cband_crit_value", None), alpha=getattr(results, "alpha", 0.05), source=type(results).__name__, - df=getattr(results, "df_inference", None), + df=df_arg, ) @@ -730,6 +783,16 @@ def _from_mpd(results: Any) -> EventStudyResults: vcov_sub = np.asarray(full_vcov, dtype=float)[np.ix_(idx, idx)] vcov_index_arr = np.asarray(covered) + # Per-period df: STRICTLY the `event_study_df` channel ({period: df} + # actually passed to safe_inference - per-period BM DOF under hc2_bm). + # `inference_df` is deliberately NOT a fallback here: it stores the + # POST-AVERAGE contrast df, which is the wrong provenance for the + # per-period rows under hc2_bm. + df_map = getattr(results, "event_study_df", None) + df_arg: Optional[np.ndarray] = None + if isinstance(df_map, dict): + df_arg = np.array([float(df_map[p]) if p in df_map else np.nan for p in all_periods]) + return EventStudyResults( event_time=np.asarray(all_periods), att=att, @@ -748,6 +811,7 @@ def _from_mpd(results: Any) -> EventStudyResults: vcov_index=vcov_index_arr, alpha=getattr(results, "alpha", 0.05), source=type(results).__name__, + df=df_arg, ) @@ -777,6 +841,16 @@ def _from_lpdid(results: Any) -> EventStudyResults: n = np.array(frame["n_obs"], dtype=float) is_ref = horizons == -1 + # Per-horizon df from the producer's `event_study_df` channel ({horizon: + # df actually passed to safe_inference} - cluster df or per-horizon + # survey df; NOT re-derivable from the frame, whose n_clusters column + # cannot recover the survey n_PSU - n_strata rule or the vcov-None + # guard). The synthetic horizon == -1 base row has no entry -> NaN. + df_map = getattr(results, "event_study_df", None) + df_arg: Optional[np.ndarray] = None + if isinstance(df_map, dict): + df_arg = np.array([float(df_map.get(int(h), np.nan)) for h in horizons]) + return EventStudyResults( event_time=horizons, att=att, @@ -792,6 +866,7 @@ def _from_lpdid(results: Any) -> EventStudyResults: event_time_convention="e0_first_treated", alpha=getattr(results, "alpha", 0.05), source=type(results).__name__, + df=df_arg, ) diff --git a/diff_diff/stacked_did.py b/diff_diff/stacked_did.py index fee37a32..012ae674 100644 --- a/diff_diff/stacked_did.py +++ b/diff_diff/stacked_did.py @@ -805,8 +805,12 @@ def _refit_stacked(w_r): # ---- Extract event study effects ---- event_study_effects: Optional[Dict[int, Dict[str, Any]]] = None + es_vcov: Optional[np.ndarray] = None + es_vcov_index: Optional[List[int]] = None + es_df_used: Optional[Dict[int, float]] = None if aggregate == "event_study": event_study_effects = {} + es_df_used = {} # Reference period (e = -1 - anticipation) event_study_effects[ref_period] = { "effect": 0.0, @@ -846,11 +850,21 @@ def _refit_stacked(w_r): t_stat = float("nan") p_value = float("nan") conf_int = (float("nan"), float("nan")) + # No safe_inference call happened -> no df provenance. + es_df_used[h] = float("nan") else: _df_eff = _bm_df if _bm_df is not None else _survey_df t_stat, p_value, conf_int = safe_inference( effect, se, alpha=self.alpha, df=_df_eff ) + # Record the df actually handed to safe_inference iff it + # governed a t-reference (finite, > 0); None (normal + # theory) and the df<=0 sentinels record NaN. + es_df_used[h] = ( + float(_df_eff) + if _df_eff is not None and np.isfinite(_df_eff) and _df_eff > 0 + else float("nan") + ) n_obs_h = int(np.sum((et_vals == h) & (d_vals == 1))) event_study_effects[h] = { "effect": effect, @@ -861,6 +875,17 @@ def _refit_stacked(w_r): "n_obs": n_obs_h, } + # Persist the event-time sub-block of the pooled-regression VCV + # (the reported ES SEs are exactly its diagonal in every + # inference mode - analytical sandwich, replicate refit, and + # survey TSL all reassign `vcov` before this block). The + # reference period is synthesized, never a regression column, so + # the index is the ESTIMATED event times only. + if event_times: + _delta_cols = [interaction_indices[h] for h in event_times] + es_vcov = vcov[np.ix_(_delta_cols, _delta_cols)] + es_vcov_index = [int(h) for h in event_times] + # ---- Compute overall ATT ---- # Average of post-treatment delta_h coefficients with delta-method SE # Post-treatment includes anticipation periods (h >= -anticipation) @@ -945,6 +970,9 @@ def _refit_stacked(w_r): balance=self.balance, covariates=list(covariates) if balancing else None, balance_diagnostics=balance_diagnostics, + event_study_vcov=es_vcov, + event_study_vcov_index=es_vcov_index, + event_study_df=es_df_used, ) self.is_fitted_ = True diff --git a/diff_diff/stacked_did_results.py b/diff_diff/stacked_did_results.py index 246021ba..1774668e 100644 --- a/diff_diff/stacked_did_results.py +++ b/diff_diff/stacked_did_results.py @@ -71,6 +71,27 @@ class StackedDiDResults(BaseResults): Clean control definition used. alpha : float Significance level used. + event_study_vcov : np.ndarray, optional + Full event-study variance-covariance matrix: the sub-block of the + pooled stacked-regression coefficient covariance over the estimated + ``D_sa x event-time`` interaction columns, ordered by + ``event_study_vcov_index``. The reported per-event-time SEs are + exactly ``sqrt(diag())`` of this matrix in every inference mode + (analytical hc1/hc2_bm sandwich, survey replicate refit, and survey + TSL all produce the coefficient covariance the SEs are read from). + The reference period is synthesized, never a regression column, so + it is absent from the index. None when no event study was requested. + event_study_vcov_index : list of int, optional + Event-time labels ordering ``event_study_vcov``'s rows/columns + (the estimated event times, reference excluded). + event_study_df : dict, optional + Per-event-time inference degrees of freedom PROVENANCE: maps each + estimated event time to the df actually passed to + ``safe_inference`` for its stored p-value/CI (per-event + Bell-McCaffrey Satterthwaite df under ``hc2_bm``; the scalar survey + df under survey designs), or NaN when the row used normal theory, + the df was undefined, or hc2_bm failed closed. None when no event + study was requested. """ overall_att: float @@ -118,6 +139,12 @@ class StackedDiDResults(BaseResults): balance: str = "none" covariates: Optional[List[str]] = None balance_diagnostics: Optional[Dict[Any, Dict[str, Any]]] = field(default=None) + # Unified event-study surface support (spec section 5, row M-092): the + # full ES VCV sub-block + ordered horizon index + per-event df actually + # used. See the class docstring for semantics. + event_study_vcov: Optional[np.ndarray] = field(default=None, repr=False) + event_study_vcov_index: Optional[List[int]] = field(default=None, repr=False) + event_study_df: Optional[Dict[int, float]] = field(default=None, repr=False) # --- Inference-field aliases (balance/external-adapter compatibility) --- @property diff --git a/diff_diff/staggered.py b/diff_diff/staggered.py index 6dd0c2c3..397b9012 100644 --- a/diff_diff/staggered.py +++ b/diff_diff/staggered.py @@ -1851,8 +1851,11 @@ def fit( if not (0 < self.pscore_trim < 0.5): raise ValueError(f"pscore_trim must be in (0, 0.5), got {self.pscore_trim}") - # Reset stale state from prior fit (prevents leaking event-study VCV) + # Reset stale state from prior fit (prevents leaking event-study VCV + # and its df provenance across fits / non-ES aggregates / the ES + # empty-result early return) self._event_study_vcov = None + self._event_study_df_used: Optional[float] = None # Tracker for _safe_inv lstsq fallbacks across all analytical SE # paths (PS Hessian, OR bread, event-study bread, etc.). Emit ONE @@ -2816,6 +2819,12 @@ def fit( if bootstrap_results is not None and event_study_vcov is not None: event_study_vcov = None event_study_vcov_index = None + # Same clearing rule for the ES df provenance: bootstrap replaces the + # stored ES se/p/CI with percentile values that never used the + # analytical df, so surfacing it would be false provenance. + event_study_df = getattr(self, "_event_study_df_used", None) + if bootstrap_results is not None: + event_study_df = None # Resolve canonical cluster_name + n_clusters for Results metadata. # Canonical PSU column wins when explicit: if survey_design has @@ -2881,6 +2890,7 @@ def fit( pscore_trim=self.pscore_trim, survey_metadata=survey_metadata, event_study_vcov=event_study_vcov, + event_study_df=event_study_df, event_study_vcov_index=event_study_vcov_index, panel=self.panel, allow_unbalanced_panel=self.allow_unbalanced_panel, diff --git a/diff_diff/staggered_aggregation.py b/diff_diff/staggered_aggregation.py index 7062106b..871174c5 100644 --- a/diff_diff/staggered_aggregation.py +++ b/diff_diff/staggered_aggregation.py @@ -993,6 +993,14 @@ def _aggregate_event_study( non_none_dfs = [d for d in agg_effective_dfs if d is not None] if non_none_dfs: df_survey_val = min(non_none_dfs) + # Stash the ONE df every ES row's inference is about to use (the + # conservative min under dropped replicates) as provenance for + # CallawaySantAnnaResults.event_study_df. Recorded iff it will + # govern a t-reference (finite, > 0; the df=0 replicate sentinel + # yields NaN inference, not a t-law). fit() resets this stash + # alongside _event_study_vcov. + if df_survey_val is not None and np.isfinite(df_survey_val) and df_survey_val > 0: + self._event_study_df_used = float(df_survey_val) t_stats, p_values, ci_lowers, ci_uppers = safe_inference_batch( np.array(agg_effects_list), np.array(agg_ses_list), diff --git a/diff_diff/staggered_results.py b/diff_diff/staggered_results.py index dadf0449..98357d7b 100644 --- a/diff_diff/staggered_results.py +++ b/diff_diff/staggered_results.py @@ -171,6 +171,16 @@ class CallawaySantAnnaResults(BaseResults): # Full event-study VCV matrix (Phase 7d): indexed by event_study_vcov_index event_study_vcov: Optional["np.ndarray"] = field(default=None, repr=False) event_study_vcov_index: Optional[list] = field(default=None, repr=False) + # event_study_df (spec section 5, row M-092): the ONE df actually applied + # to every stored event-study p-value/CI (safe_inference_batch in + # _aggregate_event_study). Equals G-1 on the bare-cluster-synthesize + # path; on explicit-survey fits it is the MINIMUM across per-horizon + # effective dfs (conservative by design when replicates were dropped) - + # the value shown is what was USED, not a per-horizon claim. None when + # the ES rows used normal theory or bootstrap overrode them. Distinct + # from `df_inference` below, whose narrow bare-cluster-only contract + # (HonestDiD consumer, PR #487) is unchanged. + event_study_df: Optional[float] = field(default=None, repr=False) bootstrap_results: Optional["CSBootstrapResults"] = field(default=None, repr=False) cband_crit_value: Optional[float] = None pscore_trim: float = 0.01 diff --git a/diff_diff/two_stage.py b/diff_diff/two_stage.py index 6f47c139..e3c73d88 100644 --- a/diff_diff/two_stage.py +++ b/diff_diff/two_stage.py @@ -1823,9 +1823,11 @@ def fit( # Event study and group aggregation (full-sample, for point estimates) event_study_effects = None group_effects = None + _es_vcov_ts: Optional[np.ndarray] = None + _es_vcov_index_ts: Optional[List[int]] = None if aggregate in ("event_study", "all"): - event_study_effects = self._stage2_event_study( + event_study_effects, _es_vcov_ts, _es_vcov_index_ts = self._stage2_event_study( df=df, unit=unit, time=time, @@ -1956,7 +1958,9 @@ def _refit_ts(w_r): results.append(att_r) if _sorted_es_periods_ts: - es_r = self._stage2_event_study( + # Replicate refits only need the point effects; the + # per-replicate V is irrelevant to the refit variance. + es_r, _, _ = self._stage2_event_study( df=df_tmp, unit=unit, time=time, @@ -2172,6 +2176,40 @@ def _refit_ts(w_r): _cluster_name_for_results = self.cluster if self.cluster is not None else unit _n_clusters_for_results = int(np.unique(df[cluster_var].values).size) + # Event-study VCV / df persistence, gated by inference mode (spec + # section 5, row M-092). Persist the analytical GMM matrix only when + # the stored ES SEs are actually its diagonal: + # - replicate-weight survey: reported SEs come from _vcov_rep_ts's + # [overall, ES, groups] layout (filtered/Prop-5 horizons absent) - + # ship vcov=None (CS precedent) but thread the final _survey_df + # that every recomputed ES row's safe_inference used; + # - bootstrap: percentile p/CIs replaced the analytical inference and + # no bootstrap covariance exists - clear both vcov and df; + # - otherwise: persist V + est_horizons; the scalar _survey_df (None + # on non-survey fits -> normal theory -> no df recorded) is the df + # every estimated row's safe_inference received. + _es_df_scalar: Optional[float] = ( + float(_survey_df) + if _survey_df is not None and np.isfinite(_survey_df) and _survey_df > 0 + else None + ) + if _uses_replicate_ts: + _es_vcov_final: Optional[np.ndarray] = None + _es_vcov_index_final: Optional[List[int]] = None + _es_df_final = _es_df_scalar + elif bootstrap_results is not None: + _es_vcov_final = None + _es_vcov_index_final = None + _es_df_final = None + else: + _es_vcov_final = _es_vcov_ts + _es_vcov_index_final = _es_vcov_index_ts + _es_df_final = _es_df_scalar + if event_study_effects is None: + _es_vcov_final = None + _es_vcov_index_final = None + _es_df_final = None + # Construct results self.results_ = TwoStageDiDResults( treatment_effects=treated_df, @@ -2196,6 +2234,9 @@ def _refit_ts(w_r): vcov_type=self.vcov_type, cluster_name=_cluster_name_for_results, n_clusters=_n_clusters_for_results, + event_study_vcov=_es_vcov_final, + event_study_vcov_index=_es_vcov_index_final, + event_study_df=_es_df_final, ) self.is_fitted_ = True @@ -2505,8 +2546,19 @@ def _stage2_event_study( score_pad_mask: Optional[np.ndarray] = None, cluster_ids_full: Optional[np.ndarray] = None, warn_nan: bool = True, - ) -> Dict[int, Dict[str, Any]]: - """Event study Stage 2: OLS of y_tilde on relative-time dummies.""" + ) -> Tuple[Dict[int, Dict[str, Any]], Optional[np.ndarray], Optional[List[int]]]: + """Event study Stage 2: OLS of y_tilde on relative-time dummies. + + Returns ``(effects, vcov, vcov_index)``: the per-horizon effects + dict, the full GMM variance-covariance matrix over the ESTIMATED + horizon coefficients, and the horizon labels ordering its + rows/columns. The reference period and Proposition-5 horizons are + never regression columns, so they appear in ``effects`` but not in + ``vcov_index``; all-filtered horizons (n_obs == 0) ARE columns, + with NaN-filled rows/columns from the rank guard. ``(dict, None, + None)`` on the degenerate early returns that fit no Stage-2 + regression. + """ y_tilde = df["_y_tilde"].values.copy() nan_mask = self._mask_nan_ytilde(y_tilde, warn=warn_nan) rel_times = df["_rel_time"].values @@ -2541,16 +2593,20 @@ def _stage2_event_study( UserWarning, stacklevel=2, ) - return { - ref_period: { - "effect": 0.0, - "se": 0.0, - "t_stat": np.nan, - "p_value": np.nan, - "conf_int": (0.0, 0.0), - "n_obs": 0, - } - } + return ( + { + ref_period: { + "effect": 0.0, + "se": 0.0, + "t_stat": np.nan, + "p_value": np.nan, + "conf_int": (0.0, 0.0), + "n_obs": 0, + } + }, + None, + None, + ) balance_mask = df[first_treat].isin(balanced_cohorts).values else: balance_mask = np.ones(n, dtype=bool) @@ -2591,16 +2647,20 @@ def _stage2_event_study( if len(est_horizons) == 0: # No horizons to estimate — return just reference period - return { - ref_period: { - "effect": 0.0, - "se": 0.0, - "t_stat": np.nan, - "p_value": np.nan, - "conf_int": (0.0, 0.0), - "n_obs": 0, - } - } + return ( + { + ref_period: { + "effect": 0.0, + "se": 0.0, + "t_stat": np.nan, + "p_value": np.nan, + "conf_int": (0.0, 0.0), + "n_obs": 0, + } + }, + None, + None, + ) # Build Stage 2 design: one column per horizon (no intercept) # Never-treated obs get all-zero rows (undefined relative time -> NaN) @@ -2706,7 +2766,7 @@ def _stage2_event_study( stacklevel=2, ) - return event_study_effects + return event_study_effects, V, [int(h) for h in est_horizons] def _stage2_group( self, diff --git a/diff_diff/two_stage_results.py b/diff_diff/two_stage_results.py index 422f7a3e..271f732f 100644 --- a/diff_diff/two_stage_results.py +++ b/diff_diff/two_stage_results.py @@ -152,6 +152,22 @@ class TwoStageDiDResults(BaseResults): vcov_type: str = "hc1" cluster_name: Optional[str] = None n_clusters: Optional[int] = None + # --- Unified event-study surface support (spec section 5, row M-092) --- + # event_study_vcov: the full Gardner-GMM covariance over the ESTIMATED + # horizon coefficients, ordered by event_study_vcov_index. The reference + # period and Proposition-5 horizons are never Stage-2 regression columns, + # so they appear in event_study_effects but not here; all-filtered + # horizons (n_obs == 0) ARE columns, carrying the rank guard's NaN + # rows/columns (consistent with their NaN marginal SEs). Both are None + # under bootstrap (percentile inference, no covariance) and under + # replicate-weight survey designs (the replicate VCV has a mixed + # [overall, ES, groups] layout; not persisted - CS precedent). + # event_study_df: the scalar df every estimated ES row's safe_inference + # received (the survey df; None on non-survey fits -> normal theory, and + # None under bootstrap where the stored inference never used a df). + event_study_vcov: Optional[np.ndarray] = field(default=None, repr=False) + event_study_vcov_index: Optional[List[int]] = field(default=None, repr=False) + event_study_df: Optional[float] = field(default=None, repr=False) # --- Inference-field aliases (balance/external-adapter compatibility) --- @property diff --git a/docs/api/_autosummary/diff_diff.CallawaySantAnnaResults.rst b/docs/api/_autosummary/diff_diff.CallawaySantAnnaResults.rst index 4c0b1966..8c585344 100644 --- a/docs/api/_autosummary/diff_diff.CallawaySantAnnaResults.rst +++ b/docs/api/_autosummary/diff_diff.CallawaySantAnnaResults.rst @@ -39,6 +39,7 @@ ~CallawaySantAnnaResults.df_inference ~CallawaySantAnnaResults.epv_diagnostics ~CallawaySantAnnaResults.epv_threshold + ~CallawaySantAnnaResults.event_study_df ~CallawaySantAnnaResults.event_study_effects ~CallawaySantAnnaResults.event_study_vcov ~CallawaySantAnnaResults.event_study_vcov_index diff --git a/docs/api/_autosummary/diff_diff.MeridianROIPrior.rst b/docs/api/_autosummary/diff_diff.MeridianROIPrior.rst new file mode 100644 index 00000000..4fb41831 --- /dev/null +++ b/docs/api/_autosummary/diff_diff.MeridianROIPrior.rst @@ -0,0 +1,31 @@ +diff\_diff.MeridianROIPrior +=========================== + +.. currentmodule:: diff_diff + +.. autoclass:: MeridianROIPrior + :no-members: + + + .. rubric:: Methods + + .. autosummary:: + + ~MeridianROIPrior.__init__ + ~MeridianROIPrior.to_code + ~MeridianROIPrior.to_dict + + + + + .. rubric:: Attributes + + .. autosummary:: + + ~MeridianROIPrior.parameter + ~MeridianROIPrior.per_experiment + ~MeridianROIPrior.roi_mean + ~MeridianROIPrior.roi_sd + ~MeridianROIPrior.mu + ~MeridianROIPrior.sigma + diff --git a/docs/api/_autosummary/diff_diff.MultiPeriodDiDResults.rst b/docs/api/_autosummary/diff_diff.MultiPeriodDiDResults.rst index 058bafd6..b3c2b8c4 100644 --- a/docs/api/_autosummary/diff_diff.MultiPeriodDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.MultiPeriodDiDResults.rst @@ -33,6 +33,7 @@ ~MultiPeriodDiDResults.conf_int ~MultiPeriodDiDResults.conley_lag_cutoff ~MultiPeriodDiDResults.df_convention + ~MultiPeriodDiDResults.event_study_df ~MultiPeriodDiDResults.fitted_values ~MultiPeriodDiDResults.inference_df ~MultiPeriodDiDResults.inference_method diff --git a/docs/api/_autosummary/diff_diff.StackedDiDResults.rst b/docs/api/_autosummary/diff_diff.StackedDiDResults.rst index f1362af1..33089fa4 100644 --- a/docs/api/_autosummary/diff_diff.StackedDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.StackedDiDResults.rst @@ -34,6 +34,9 @@ ~StackedDiDResults.coef_var ~StackedDiDResults.conf_int ~StackedDiDResults.covariates + ~StackedDiDResults.event_study_df + ~StackedDiDResults.event_study_vcov + ~StackedDiDResults.event_study_vcov_index ~StackedDiDResults.is_significant ~StackedDiDResults.kappa_post ~StackedDiDResults.kappa_pre diff --git a/docs/api/_autosummary/diff_diff.SunAbrahamResults.rst b/docs/api/_autosummary/diff_diff.SunAbrahamResults.rst index 527b7c8c..d55db5db 100644 --- a/docs/api/_autosummary/diff_diff.SunAbrahamResults.rst +++ b/docs/api/_autosummary/diff_diff.SunAbrahamResults.rst @@ -38,6 +38,8 @@ ~SunAbrahamResults.event_study_vcov_index ~SunAbrahamResults.is_significant ~SunAbrahamResults.p_value + ~SunAbrahamResults.reference_observed + ~SunAbrahamResults.reference_period ~SunAbrahamResults.se ~SunAbrahamResults.significance_stars ~SunAbrahamResults.survey_metadata diff --git a/docs/api/_autosummary/diff_diff.TwoStageDiDResults.rst b/docs/api/_autosummary/diff_diff.TwoStageDiDResults.rst index 6c340e7f..35aadf13 100644 --- a/docs/api/_autosummary/diff_diff.TwoStageDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.TwoStageDiDResults.rst @@ -31,6 +31,9 @@ ~TwoStageDiDResults.cluster_name ~TwoStageDiDResults.coef_var ~TwoStageDiDResults.conf_int + ~TwoStageDiDResults.event_study_df + ~TwoStageDiDResults.event_study_vcov + ~TwoStageDiDResults.event_study_vcov_index ~TwoStageDiDResults.is_significant ~TwoStageDiDResults.n_clusters ~TwoStageDiDResults.p_value diff --git a/docs/api/_autosummary/diff_diff.lpdid_results.LPDiDResults.rst b/docs/api/_autosummary/diff_diff.lpdid_results.LPDiDResults.rst index 18cfe995..ff3f4a3e 100644 --- a/docs/api/_autosummary/diff_diff.lpdid_results.LPDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.lpdid_results.LPDiDResults.rst @@ -32,6 +32,7 @@ ~LPDiDResults.covariates ~LPDiDResults.dylags ~LPDiDResults.estimand + ~LPDiDResults.event_study_df ~LPDiDResults.n_clusters ~LPDiDResults.n_psu ~LPDiDResults.n_strata diff --git a/docs/api/_autosummary/diff_diff.to_meridian_roi_prior.rst b/docs/api/_autosummary/diff_diff.to_meridian_roi_prior.rst new file mode 100644 index 00000000..71244042 --- /dev/null +++ b/docs/api/_autosummary/diff_diff.to_meridian_roi_prior.rst @@ -0,0 +1,7 @@ +diff\_diff.to\_meridian\_roi\_prior +=================================== + +.. currentmodule:: diff_diff + +.. autofunction:: to_meridian_roi_prior + :no-index: \ No newline at end of file diff --git a/docs/api/_autosummary/diff_diff.to_pymc_marketing_lift_test.rst b/docs/api/_autosummary/diff_diff.to_pymc_marketing_lift_test.rst new file mode 100644 index 00000000..e0d458a3 --- /dev/null +++ b/docs/api/_autosummary/diff_diff.to_pymc_marketing_lift_test.rst @@ -0,0 +1,7 @@ +diff\_diff.to\_pymc\_marketing\_lift\_test +========================================== + +.. currentmodule:: diff_diff + +.. autofunction:: to_pymc_marketing_lift_test + :no-index: \ No newline at end of file diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 36ff4481..85530ec0 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -1759,6 +1759,7 @@ Our implementation uses multiplier bootstrap on the GMM influence function: clus - **Note:** The GMM sandwich and bootstrap paths both use `scipy.sparse.linalg.factorized` for the Stage 1 normal-equations solve `(X'_{10} W X_{10}) gamma = X'_1 W X_2` and fall back to a **certified sparse LSMR** solve when the sparse factorization raises `RuntimeError` on a near-singular matrix (2026-07; previously dense `lstsq(toarray())`, an `O((U+T+K)²)` materialization and OOM risk on large panels — the pattern the ImputationDiD LSMR fix closed). Solver choice cannot change the output: least-squares solutions differ only by `null(X'X) = null(X_10)` components (weighted: `null(W^{1/2}X_10)`; zero-weight rows are inert in every weighted consumer), and every `gamma_hat` consumer is an `X_10`-range functional — `Psi = X_10 gamma`, the GMM score correction `c_g' gamma` with `c_g = X_{10,g}' eps_{10,g}` in `rowspace(X_10)` — so its null component annihilates. The bootstrap's `X_1_sparse @ theta_exact` consumer evaluates theta on treated rows where annihilation does not apply; parity there holds because both dense `lstsq` (SVD) and LSMR return the MIN-NORM least-squares solution, so the solvers agree on the whole vector to iterative tolerance (helper-level dense-lstsq-oracle parity test on a singular Gram + fit-level singular-design V/SE parity test). Convergence is validated fail-closed: `istop` in `{0,1,2,4,5}` certified, one uncapped-conlim retry, then `_LSMRUnconvergedError` — the analytical variance boundary reports NaN vcov/SE (raising matters: the GMM score `nan_to_num` would otherwise launder NaN scores into zeros and a finite, wrong variance) and the bootstrap boundary degenerates to the established `None` contract. All fallback sites (analytical unweighted + weighted, SpilloverDiD Wave-D meat, bootstrap) emit a `UserWarning` (silent-failure audit axis C). - **Note:** The GMM sandwich re-solves the Stage-1 unit+time fixed effects **exactly** (sparse OLS reusing the `scipy.sparse.linalg.factorized` factorization of `(X'_{10} W X_{10})` already computed for `gamma_hat`), rather than reusing the iterative alternating-projection FE (`_iterative_fe`) that produces the point estimate. The iterative solver converges only to ~1e-7 on unbalanced *untreated* panels — negligible for the ATT, but enough to perturb the variance by ~1% relative to the analytical GMM sandwich. The exact re-solve makes the analytical GMM SE match R `did2s` to ~1e-7 (`tests/test_methodology_two_stage.py::TestTwoStageDiDParityR`), mirroring ImputationDiD's exact-sparse variance path (`imputation.py` `_build_A_sparse`). It also required adding an **intercept column** to `_build_fe_design` so the first-stage column space spans the constant (the grand mean): the prior intercept-free `[unit_1.., time_1..]` layout (drop first unit + first time, no intercept) silently omitted the grand mean, which the exact residual is first-order sensitive to (the iterative point-estimate solver absorbs the grand mean into its mean-based FE, so the point estimate was unaffected). Obs whose unit or time FE are unidentified (NaN; rank-deficient / Proposition-5) fall back to the iterative residual, so those edge cases are unchanged. The reported `overall_att` still uses the iterative FE (preserving point-estimate equivalence with ImputationDiD at 1e-10); only the variance uses the exact residuals. - **Note:** `vcov_type` is permanently narrow to `{"hc1"}` (Phase 1b threading). TwoStageDiD's variance is the Gardner (2022) two-stage GMM cluster-sandwich `V = (X'_2 W X_2)^{-1} (S' S) (X'_2 W X_2)^{-1}` with the per-cluster GMM-corrected score `S_g = gamma_hat' c_g - X'_{2g} eps_{2g}`. Analytical-sandwich families `{classical, hc2, hc2_bm}` are rejected at `__init__`/`fit()`: the GMM-corrected meat folds first-stage FE estimation uncertainty into the score via the `gamma_hat' c_g` term, so there is no single hat matrix spanning both stages on which HC2 leverage or Bell-McCaffrey Satterthwaite DOF can be defined, and the Gardner first-stage correction has not been derived for the leverage-corrected or homoskedastic meat structures (no reference implementation — `clubSandwich` covers single-equation WLS/OLS CR2, not two-stage GMM; mirrors the SpilloverDiD `vcov_type="classical"` rejection). `cluster=` selects the cluster level; `cluster=None` (the default) still clusters at the `unit` column (`cluster_var = unit`), so the summary renders `"CR1 cluster-robust at , G="` rather than the generic `"HC1"` label. **Note (deviation from R):** the did2s GMM sandwich uses NO finite-sample multiplier (meat `= S' S`), so the rendered `CR1` family label carries no Stata-style `(n-1)/(n-p)` or `G/(G-1)` factor (matches R `did2s`; same FSA-free convention as ImputationDiD's Theorem 3 variance). Under bootstrap (`n_bootstrap > 0`) the analytical variance-family label is suppressed in `summary()` because `fit()` overwrites the reported SE/CI/p-value with `bootstrap_results` (mirrors `DiDResults` at `results.py:213-226`). `cluster=` combined with a replicate-weight survey design raises `NotImplementedError` (replicate-refit variance ignores `cluster=`). `vcov_type='conley'` is deferred to the TwoStageDiD Conley follow-up row in DEFERRED.md. +- **Note:** Provenance exposure (4.0 program row M-092, no numeric change): the full Gardner-GMM event-study covariance `V` is now persisted as `TwoStageDiDResults.event_study_vcov` (+ `event_study_vcov_index` = the estimated horizons, in Stage-2 column order). The reference period and Proposition-5 horizons are never Stage-2 columns and are absent from the index; an estimated horizon whose observations are all filtered keeps its column with the rank-guard's NaN row/col (consistent with its NaN marginal SE). Persisted only on the analytical paths - `None` under bootstrap (percentile inference has no covariance) and replicate-weight survey designs (the replicate VCV's [overall, ES, groups] layout is not an ES matrix; CS precedent). `event_study_df` records the scalar survey df every estimated ES row's `safe_inference` received (None on non-survey fits: normal-theory GMM inference). **Reference implementation(s):** - R: `did2s::did2s()` (Kyle Butts & John Gardner) @@ -1894,6 +1895,7 @@ For each sub-experiment `a`, entropy balancing produces nonnegative control weig The pooled estimator is `DID^{CBWSDID}_e = Σ_a (N^D_a/N^D_Ω)(Δ̄^D_{a,e} − Δ̄^{C,b}_{a,e})`, recovered by the existing Q-weighted WLS when `W_{sa}` is injected. **Estimand preservation:** because only controls are reweighted (treated cohorts and their shares `N^D_a/N^D_Ω` are unchanged), the target remains the trimmed aggregate ATT `θ^e_κ` — the refinement changes only how untreated trends are estimated. At `b_{sa}=1` this reduces to the paper's **unit-count** weighted stacked DID, which equals `StackedDiD(weighting="aggregate")` on **balanced** event windows (where unit and observation counts coincide). Validated end-to-end by `tests/test_methodology_stacked_did.py::TestCBWSDIDCovariateBalance` (closed-form `DID^{CBWSDID}_e` anchor at 1e-8), `::TestCBWSDIDEffectiveMass` (effective-mass corrector is load-bearing vs a naive `b·Q` multiply), and `::TestCBWSDIDRParity` (cross-language parity vs the R `cbwsdid` package, `refinement.method="weightit"` / `method="ebal"`, at 1e-5 — golden in `benchmarks/data/cbwsdid_golden.json`, regenerate via `benchmarks/R/generate_cbwsdid_golden.R`). - **Note:** The effective-mass `W_{sa}` is computed directly from cohort unit-counts + `Ñ^C_a` (a naive `b_{sa}·Q_aggregate` multiply is NOT equivalent — it aggregates control means with weights ∝ `(N^D_a/N^D_Ω)(Ñ^C_a/N^C_a)` instead of the required `∝ (N^D_a/N^D_Ω)`, biased unless `b` is uniform). +- **Note:** Provenance exposure (4.0 program row M-092, no numeric change): the event-time sub-block of the pooled-regression coefficient covariance is now persisted as `StackedDiDResults.event_study_vcov` (+ `event_study_vcov_index` = the estimated event times; the synthesized reference period is not a regression column and is absent). Persisted in every inference mode - analytical hc1/hc2_bm sandwich, survey replicate refit, and survey TSL all reassign the coefficient covariance the reported ES SEs are read from, so the marginal SEs are always this matrix's diagonal. `event_study_df` records, per event time, the df actually passed to `safe_inference` (the per-event Bell-McCaffrey Satterthwaite contrast df under `hc2_bm`, the scalar survey df under survey designs, NaN on the normal-theory hc1 path and on hc2_bm fail-closed rows). DiagnosticReport's parallel-trends check consumes the persisted VCV via its joint-Wald path for hc1/survey fits; `hc2_bm` fits are routed to Bonferroni over the BM-adjusted per-row p-values instead (the chi-square Wald reference would discard the CR2/BM small-sample correction; multi-constraint AHT/HTZ CR2 testing is tracked in DEFERRED.md - see REPORTING.md "hc2_bm parallel-trends policy"). - **Note:** Inference is conditional-on-the-estimated-weights cluster-robust (the existing `hc1`/`hc2_bm` path with `W_{sa}` as the WLS weights) — the paper's default. The paper's weight-re-estimating bootstrap is NOT implemented in v1 (deliberate scope; entropy balancing is smooth so the Abadie–Imbens (2008) nonsmooth-matching bootstrap caveat does not apply). `cluster` is orthogonal to `b_{sa}` (weights conditioned-on); default `unit` matches the paper. - **Note:** v1 scope — only `balance="entropy"` with `weighting="aggregate"`. `balance` + `population`/`sample_share` and `balance` + `survey_design=` raise `NotImplementedError`; matching-based balancing and the repeated `0→1/1→0` episode extension are out of scope. diff --git a/docs/methodology/REPORTING.md b/docs/methodology/REPORTING.md index e4599ad9..ead23d7d 100644 --- a/docs/methodology/REPORTING.md +++ b/docs/methodology/REPORTING.md @@ -380,9 +380,47 @@ a library setting. unused" case, so no downgrade applies). Remaining `"diag_fallback"` cases on new fits — bootstrap / - replicate-weight CS and SA, plus ImputationDiD / Stacked / - EfficientDiD / TwoStageDiD — pass through unchanged because - nothing better is available on those result types yet. + replicate-weight CS / SA, bootstrap / replicate-weight TwoStageDiD, + plus ImputationDiD / EfficientDiD — pass through unchanged because + nothing better is available on those result types in those modes. + StackedDiD and TwoStageDiD now persist `event_study_vcov` (4.0 + program row M-092): StackedDiD in EVERY inference mode (its + reported event-study SEs are always the persisted matrix's + diagonal), TwoStageDiD on the analytical paths only (bootstrap and + replicate-weight modes clear it). Where the covariance is present, + the PT check takes the joint-Wald path (subject to the hc2_bm + policy and rank guard below). Pretrends POWER is unaffected: + `compute_pretrends_power()` supports only MPD / CS / SunAbraham + fits, and Stacked/TwoStage remain outside DR's power applicability + - their new covariance is consumed by the PT check and by the + PRECOMPUTED-power provenance classifier only (a stored power result + on a fit that now exposes a full VCV is labelled + `"diag_fallback_available_full_vcov_unused"`). + +- **Note:** hc2_bm parallel-trends policy (deviation by omission, + documented). When the source fit's `vcov_type` is `"hc2_bm"`, the + PT check does NOT run the covariance-based joint Wald even when a + full event-study covariance is persisted: the generic chi-square / + F(k, df_survey) reference ignores the CR2 / Bell-McCaffrey + small-sample correction the estimator's own per-row inference + applies, and the proper multi-constraint CR2 test (AHT/HTZ, + `clubSandwich::Wald_test(test="HTZ")`) is not implemented (tracked + in DEFERRED.md). The check instead runs Bonferroni over the + BM-adjusted per-row p-values, preserving the small-sample + correction. Additionally, the joint-Wald path fail-closes on + provenance: if ANY retained pre-period row carries a non-finite + p-value (hc2_bm BM-DOF failure, collapsed replicate-survey df), the + check returns `method="inconclusive"` rather than computing a Wald + statistic through effects whose inference the estimator refused to + produce — the covariance-backed twin of the round-33/42 Bonferroni + guards. Finally, the joint-Wald solve is rank-guarded: a + cluster-sandwich covariance has rank bounded by the cluster-score + dimension, so when the number of tested leads approaches the number + of clusters the pre-period block is singular and `β′V⁻¹β` is a + garbage (possibly negative) quadratic form. The check requires a + finite, numerically full-rank, PSD block and a non-negative finite + statistic; otherwise it falls back to Bonferroni over the per-row + p-values (a documented conservative downgrade, not a defect). - **Note:** Unit-translation policy. BusinessReport does not arithmetically translate log-points to percents or level effects to diff --git a/docs/v4-deprecations.yaml b/docs/v4-deprecations.yaml index 52ec2623..5c54d1c1 100644 --- a/docs/v4-deprecations.yaml +++ b/docs/v4-deprecations.yaml @@ -931,7 +931,7 @@ rows: phase: 2 test_ref: tests/test_event_study_surface.py code_refs: [diff_diff/results_base.py, diff_diff/__init__.py] - notes: "Phase 2 unified event-study representation (spec section 5): EventStudyResults container + builders for the 14 producers (CallawaySantAnna, SunAbraham, ImputationDiD, TwoStageDiD, StackedDiD, SpilloverDiD, ContinuousDiD, EfficientDiD, WooldridgeDiD, StaggeredTripleDifference, MultiPeriodDiD, LPDiD, ChaisemartinDHaultfoeuille, HeterogeneousAdoptionDiD). Canonical quintet columns, explicit is_reference marking (successor to the retiring sentinels [M-093]), vcov+vcov_index ordering, cband columns, event_time_convention metadata. Public exposure rides aggregate(type='event_study') in Phase 2 PR (b); merged TWFE returns it in Phase 3 [M-010]. introduced_in gates the 3.9 cut, mirroring [M-091]. Born done in this introducing diff (builder is package-internal; the class is exported)." + notes: "Phase 2 unified event-study representation (spec section 5): EventStudyResults container + builders for the 14 producers (CallawaySantAnna, SunAbraham, ImputationDiD, TwoStageDiD, StackedDiD, SpilloverDiD, ContinuousDiD, EfficientDiD, WooldridgeDiD, StaggeredTripleDifference, MultiPeriodDiD, LPDiD, ChaisemartinDHaultfoeuille, HeterogeneousAdoptionDiD). Canonical quintet columns, explicit is_reference marking (successor to the retiring sentinels [M-093]), vcov+vcov_index ordering, cband columns, event_time_convention metadata. Public exposure rides aggregate(type='event_study') in Phase 2 PR (b); merged TWFE returns it in Phase 3 [M-010]. introduced_in gates the 3.9 cut, mirroring [M-091]. Born done in this introducing diff (builder is package-internal; the class is exported). Amended pre-cut (same test_ref, introduced_in 3.9 not yet released): df became PER-ROW (one entry per event time, the df each stored p/CI actually used; joins the pinned to_dataframe schema) and StackedDiD/TwoStageDiD persist their internal full ES VCVs (event_study_vcov/_index/_df container fields; mode-gated for TwoStageDiD bootstrap/replicate)." - id: M-093 kind: behavior group: results-contract diff --git a/docs/v4-design.md b/docs/v4-design.md index a025a0a3..de2ec853 100644 --- a/docs/v4-design.md +++ b/docs/v4-design.md @@ -380,7 +380,9 @@ it natively in Phase 3). This document pins the requirements it satisfies: per-event-time canonical quintet rows, explicit reference-period marking (via an `is_reference` column - no sentinel-value conventions; the n_groups==0 / n_obs==0 sentinels are retired at 4.0 [M-093]), the event-study -vcov exposed uniformly where computed, and `to_dataframe(level="event_study")` +vcov exposed uniformly where computed, per-row inference-df provenance (the +df each stored p-value/CI actually used, threaded via the producers' +`event_study_df` channels), and `to_dataframe(level="event_study")` emitting identical column schemas from every estimator. **Pickle migration.** Renamed-field classes ship `__setstate__` migration diff --git a/tests/test_diagnostic_report.py b/tests/test_diagnostic_report.py index 6a415b47..d3a15a0b 100644 --- a/tests/test_diagnostic_report.py +++ b/tests/test_diagnostic_report.py @@ -2971,3 +2971,367 @@ def test_business_report_smoke(self): def test_public_api_exports(): for name in ("DiagnosticReport", "DiagnosticReportResults", "DIAGNOSTIC_REPORT_SCHEMA_VERSION"): assert hasattr(dd, name), f"diff_diff must export {name}" + + +# --------------------------------------------------------------------------- +# ES-VCV consumers: Stacked/TwoStage now expose event_study_vcov (M-092 +# follow-up), upgrading the parallel-trends check from the Bonferroni +# fallback to the joint-Wald path. Intended behavior change, CHANGELOG'd. +# --------------------------------------------------------------------------- +class TestStackedTwoStageJointWaldUpgrade: + @staticmethod + def _panel(seed=42): + rng = np.random.default_rng(seed) + rows = [] + for u in range(60): + g = [4, 6, 0][u % 3] + for t in range(1, 11): + y = 1.0 + 0.1 * t + u * 0.01 + (1.5 if g and t >= g else 0.0) + rng.normal(0, 0.3) + rows.append({"unit": u, "time": t, "outcome": y, "first_treat": g}) + return pd.DataFrame(rows) + + def test_stacked_pt_check_uses_joint_wald(self): + from diff_diff import StackedDiD + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = StackedDiD(kappa_pre=2, kappa_post=2).fit( + self._panel(), + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + rep = DiagnosticReport(res).run_all() + pt = rep.to_dict()["parallel_trends"] + assert pt["method"] == "joint_wald_event_study", ( + f"StackedDiD persists its full ES VCV; the PT check must take the " + f"joint-Wald path, got {pt.get('method')}" + ) + assert np.isfinite(pt["joint_p_value"]) + + def test_two_stage_nan_vcov_row_excluded_from_pt_family(self): + # Induce REAL Path-B horizons (two_stage.py all-filtered branch): + # balance_e=4 drops cohort 7 (unobserved through e=4) AFTER the + # horizon grid is built from all rows, so its exclusive deep + # pre-horizons (-6, -5) stay estimated Stage-2 columns whose rows + # are all balance-filtered -> NaN effect/SE and rank-guard NaN + # rows/cols in the persisted VCV. The PT check must EXCLUDE those + # horizons from the tested family and run the joint Wald on the + # finite pre-period sub-block - never compute a statistic THROUGH + # the NaN covariance entries. + from diff_diff.two_stage import TwoStageDiD + + rng = np.random.default_rng(9) + rows = [] + for u in range(90): + g = [5, 7, 0][u % 3] + for t in range(1, 11): + y = 1.0 + 0.1 * t + u * 0.01 + (1.0 if g and t >= g else 0.0) + rng.normal(0, 0.3) + rows.append( + { + "unit": u, + "time": t, + "outcome": y, + "first_treat": g if g else np.nan, + } + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = TwoStageDiD(pretrends=True).fit( + pd.DataFrame(rows), + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + balance_e=4, + ) + # Path-B precondition: -6 is an estimated column with NaN inference + # and a NaN-filled VCV row; finite pre-periods {-4, -3, -2} remain. + assert res.event_study_effects[-6]["n_obs"] == 0 + assert not np.isfinite(res.event_study_effects[-6]["se"]) + idx = res.event_study_vcov_index + assert -6 in idx and -4 in idx + pos = idx.index(-6) + assert np.isnan(res.event_study_vcov[pos, :]).all() + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + rep = DiagnosticReport(res).run_all() + pt = rep.to_dict()["parallel_trends"] + assert pt["status"] == "ran" + # Joint Wald ran on the FINITE pre-period family only: the Path-B + # horizons (NaN se) are excluded by the collector, so the tested + # block is finite and the p-value is a real number, never + # NaN-contaminated by the excluded rows of the matrix. + assert pt["method"] == "joint_wald_event_study" + assert np.isfinite(pt["joint_p_value"]) + tested = {row["period"] for row in pt["per_period"]} + assert tested == {-4, -3, -2} + assert pt["df"] == 3 + + @staticmethod + def _hc2bm_fit(monkeypatch=None): + from diff_diff import StackedDiD + + rng = np.random.default_rng(42) + rows = [] + for u in range(60): + g = [4, 6, 0][u % 3] + for t in range(1, 11): + y = 1.0 + 0.1 * t + u * 0.01 + (1.5 if g and t >= g else 0.0) + rng.normal(0, 0.3) + rows.append({"unit": u, "time": t, "outcome": y, "first_treat": g}) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return StackedDiD(kappa_pre=2, kappa_post=2, vcov_type="hc2_bm").fit( + pd.DataFrame(rows), + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + + def test_stacked_hc2_bm_fail_closed_pt_is_inconclusive(self, monkeypatch): + # Provenance-laundering regression: when the BM contrast + # DOF is unavailable, StackedDiD fails closed (finite effect/SE, + # deliberately-NaN t/p/CI) but the persisted covariance stays + # finite. The joint-Wald path consumes only effects + covariance, + # so without the provenance guard it would publish a finite joint + # p-value for inference the estimator refused to produce. + import diff_diff.linalg as _linalg + + def _broken(*args, **kwargs): + raise ValueError("forced BM failure") + + monkeypatch.setattr(_linalg, "_compute_cr2_bm_contrast_dof", _broken) + res = self._hc2bm_fit() + # Fail-closed precondition: finite effect/SE, NaN p, finite vcov. + pre = res.event_study_effects[-2] + assert np.isfinite(pre["effect"]) and np.isfinite(pre["se"]) + assert not np.isfinite(pre["p_value"]) + assert res.event_study_vcov is not None + assert np.isfinite(np.diag(res.event_study_vcov)).all() + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + pt = DiagnosticReport(res).run_all().to_dict()["parallel_trends"] + assert pt["method"] == "inconclusive" + assert pt["verdict"] == "inconclusive" + assert pt["joint_p_value"] is None + assert pt["n_dropped_undefined"] >= 1 + + def test_stacked_hc2_bm_routes_to_bonferroni_not_chi_square_wald(self): + # Healthy hc2_bm fit: every pre-row's p-value is BM-adjusted. The + # generic chi-square joint Wald would discard that small-sample + # correction, so the PT check must use Bonferroni over the + # BM-adjusted per-row p-values despite the persisted covariance + # (REPORTING.md "hc2_bm parallel-trends policy"; AHT/HTZ tracked in + # DEFERRED.md). + res = self._hc2bm_fit() + assert res.event_study_vcov is not None + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + pt = DiagnosticReport(res).run_all().to_dict()["parallel_trends"] + assert pt["method"] == "bonferroni" + assert np.isfinite(pt["joint_p_value"]) + + def test_vcov_backed_nan_p_rows_force_inconclusive(self): + # Construction-level twin of the hc2_bm case for the collapsed + # replicate-survey df mode: retained pre-rows with finite + # effect/SE, NaN p (safe_inference df<=0 sentinel), and a persisted + # finite covariance. The joint-Wald provenance guard must return + # inconclusive, never a finite p through the covariance. + class _FakeCovBacked: + alpha = 0.05 + vcov_type = "hc1" + event_study_effects = { + -3: { + "effect": 0.10, + "se": 0.05, + "t_stat": np.nan, + "p_value": np.nan, + "conf_int": (np.nan, np.nan), + "n_obs": 40, + }, + -2: { + "effect": 0.05, + "se": 0.04, + "t_stat": 1.25, + "p_value": 0.21, + "conf_int": (-0.03, 0.13), + "n_obs": 40, + }, + -1: { + "effect": 0.0, + "se": 0.0, + "t_stat": np.nan, + "p_value": np.nan, + "conf_int": (np.nan, np.nan), + "n_obs": 0, + }, + 0: { + "effect": 0.9, + "se": 0.1, + "t_stat": 9.0, + "p_value": 0.0, + "conf_int": (0.7, 1.1), + "n_obs": 40, + }, + } + event_study_vcov = np.array( + [[0.0025, 0.0002, 0.0001], [0.0002, 0.0016, 0.0002], [0.0001, 0.0002, 0.01]] + ) + event_study_vcov_index = [-3, -2, 0] + + from diff_diff.diagnostic_report import DiagnosticReport as _DR + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + pt = _DR(_FakeCovBacked())._pt_event_study() + assert pt["method"] == "inconclusive" + assert pt["joint_p_value"] is None + assert pt["n_dropped_undefined"] == 1 + + @staticmethod + def _low_cluster_panel(): + # 12 clusters, one cohort at t=14 with a 12-lead pre window: the + # cluster-sandwich covariance has rank bounded by the cluster-score + # dimension, so the retained pre-period block is SINGULAR (rank < + # n_pre) - the normal low-G/high-lead edge case. + rng = np.random.default_rng(3) + rows = [] + for u in range(12): + g = 14 if u % 2 == 0 else 0 + for t in range(1, 18): + y = 1.0 + 0.05 * t + u * 0.05 + (1.0 if g and t >= g else 0.0) + rng.normal(0, 0.3) + rows.append({"unit": u, "time": t, "outcome": y, "first_treat": g}) + return pd.DataFrame(rows) + + def test_stacked_singular_vcov_downgrades_to_bonferroni(self): + from diff_diff import StackedDiD + + data = self._low_cluster_panel() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = StackedDiD(kappa_pre=12, kappa_post=2).fit( + data, + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + # Precondition: the retained pre-period block is genuinely singular. + V = res.event_study_vcov + idx = res.event_study_vcov_index + pre_idx = [i for i, h in enumerate(idx) if h < -1] + sub = V[np.ix_(pre_idx, pre_idx)] + sub = 0.5 * (sub + sub.T) + assert np.linalg.matrix_rank(sub) < len(pre_idx) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + pt = DiagnosticReport(res).run_all().to_dict()["parallel_trends"] + # Rank guard: never a Wald statistic through the singular block - + # the unguarded quadratic form here is ~ -2.5e14 with chi2 p=1.0. + assert pt["method"] == "bonferroni" + assert pt["test_statistic"] is None + assert np.isfinite(pt["joint_p_value"]) + + def test_two_stage_singular_vcov_downgrades_to_bonferroni(self): + from diff_diff.two_stage import TwoStageDiD + + data = self._low_cluster_panel() + data["first_treat"] = data["first_treat"].replace(0, np.nan) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = TwoStageDiD(pretrends=True).fit( + data, + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + V = res.event_study_vcov + idx = res.event_study_vcov_index + pre_idx = [i for i, h in enumerate(idx) if h < -1 and np.isfinite(V[i, i])] + sub = V[np.ix_(pre_idx, pre_idx)] + sub = 0.5 * (sub + sub.T) + assert np.linalg.matrix_rank(sub) < len(pre_idx) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + pt = DiagnosticReport(res).run_all().to_dict()["parallel_trends"] + assert pt["method"] == "bonferroni" + assert pt["test_statistic"] is None + assert np.isfinite(pt["joint_p_value"]) + + def test_joint_wald_never_negative_statistic(self): + # Belt-and-braces contract: a Wald statistic is a PSD quadratic + # form. Feed a singular-but-finite covariance directly through the + # event-study path and assert the guard rejects it rather than + # reporting a negative statistic as chi-square p=1.0. + rng = np.random.default_rng(0) + A = rng.normal(size=(3, 2)) + V = A @ A.T # PSD, rank 2 < 3 + + class _FakeSingular: + alpha = 0.05 + vcov_type = "hc1" + event_study_effects = { + -4: { + "effect": 0.10, + "se": float(np.sqrt(V[0, 0])), + "t_stat": 1.0, + "p_value": 0.3, + "conf_int": (-0.1, 0.3), + "n_obs": 40, + }, + -3: { + "effect": 0.05, + "se": float(np.sqrt(V[1, 1])), + "t_stat": 0.5, + "p_value": 0.6, + "conf_int": (-0.1, 0.2), + "n_obs": 40, + }, + -2: { + "effect": -0.07, + "se": float(np.sqrt(V[2, 2])), + "t_stat": -0.7, + "p_value": 0.5, + "conf_int": (-0.2, 0.1), + "n_obs": 40, + }, + -1: { + "effect": 0.0, + "se": 0.0, + "t_stat": np.nan, + "p_value": np.nan, + "conf_int": (np.nan, np.nan), + "n_obs": 0, + }, + 0: { + "effect": 0.9, + "se": 0.1, + "t_stat": 9.0, + "p_value": 0.0, + "conf_int": (0.7, 1.1), + "n_obs": 40, + }, + } + event_study_vcov = V + event_study_vcov_index = [-4, -3, -2] + + from diff_diff.diagnostic_report import DiagnosticReport as _DR2 + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + pt = _DR2(_FakeSingular())._pt_event_study() + assert pt["method"] == "bonferroni" + stat = pt["test_statistic"] + assert stat is None or stat >= 0.0 diff --git a/tests/test_event_study_surface.py b/tests/test_event_study_surface.py index acf3ba96..cfe0c6d0 100644 --- a/tests/test_event_study_surface.py +++ b/tests/test_event_study_surface.py @@ -1119,3 +1119,407 @@ def _native_event_map(label, native): "ci_hi": ci[1], } return out + + +# =========================================================================== +# Per-row df provenance + Stacked/TwoStage vcov persistence (M-092 follow-up) +# =========================================================================== + + +def test_df_scalar_broadcasts_per_row(): + surface = _tiny_surface(df=12.0) + assert isinstance(surface.df, np.ndarray) + assert np.isnan(surface.df[0]) # reference row: no df provenance + assert surface.df[1] == 12.0 and surface.df[2] == 12.0 + + +def test_df_none_yields_nan_column(): + surface = _tiny_surface() + assert isinstance(surface.df, np.ndarray) + assert np.isnan(surface.df).all() + + +def test_df_array_length_validated(): + with pytest.raises(ValueError, match="'df' has shape"): + _tiny_surface(df=np.array([1.0, 2.0])) + + +def test_df_nan_where_p_value_nan(): + # Row 2: finite se but NaN p (non-estimable inference) - its df must be + # stripped, because df is provenance of a STORED p-value. + surface = _tiny_surface( + p_value=np.array([np.nan, 0.0, np.nan]), + df=np.array([5.0, 12.0, 12.0]), + ) + assert np.isnan(surface.df[0]) # reference + assert surface.df[1] == 12.0 + assert np.isnan(surface.df[2]) # NaN-p row never used a df + + +def test_df_round_trips_through_to_dict(): + import json + + d = _tiny_surface(df=7.0).to_dict() + json.dumps(d) # must not raise + assert isinstance(d["df"], list) and len(d["df"]) == 3 + assert d["df"][1] == 7.0 and np.isnan(d["df"][0]) + + +def test_df_column_in_pinned_frame(): + surface = _tiny_surface(df=9.0) + frame = surface.to_dataframe() + assert "df" in frame.columns # via EVENT_STUDY_SCHEMA + np.testing.assert_array_equal(frame["df"].to_numpy(), surface.df) + + +@pytest.mark.parametrize("label", ALL_PRODUCERS) +def test_df_nan_on_reference_and_nan_p_rows(label, surfaces): + _, surface = surfaces[label] + assert isinstance(surface.df, np.ndarray) + assert surface.df.shape == surface.event_time.shape + assert np.isnan(surface.df[surface.is_reference]).all() + assert np.isnan(surface.df[~np.isfinite(surface.p_value)]).all() + + +def test_stacked_vcov_alignment(surfaces): + native, surface = surfaces["StackedDiD"] + assert surface.vcov is not None + assert surface.vcov_index is not None + idx = surface.vcov_index.tolist() + assert set(idx).issubset(set(surface.event_time.tolist())) + # The reference period is synthesized, never a regression column. + assert -1 not in idx + assert idx == native.event_study_vcov_index + # The marginal ES SEs are literally this matrix's diagonal (pure copy). + ses = {h: d["se"] for h, d in native.event_study_effects.items()} + diag = np.sqrt(np.maximum(np.diag(surface.vcov), 0.0)) + np.testing.assert_allclose(diag, [ses[h] for h in idx], rtol=1e-14) + + +def test_two_stage_vcov_alignment(surfaces): + native, surface = surfaces["TwoStageDiD"] + assert surface.vcov is not None + assert surface.vcov_index is not None + idx = surface.vcov_index.tolist() + assert set(idx).issubset(set(surface.event_time.tolist())) + assert -1 not in idx # ref_period never a Stage-2 column + assert idx == native.event_study_vcov_index + ses = {h: d["se"] for h, d in native.event_study_effects.items()} + diag = np.sqrt(np.maximum(np.diag(surface.vcov), 0.0)) + for i, h in enumerate(idx): + if np.isfinite(ses[h]): + np.testing.assert_allclose(diag[i], ses[h], rtol=1e-14) + else: + assert np.isnan(diag[i]) # rank-guard NaN row == NaN marginal SE + + +def test_two_stage_prop5_horizons_excluded_from_vcov_index(): + # No never-treated units -> late horizons are Proposition-5 unidentified + # (NaN effect, n_obs > 0). They live in the effects dict but were never + # Stage-2 columns, so the persisted vcov_index must exclude them. + from diff_diff.two_stage import TwoStageDiD + + rng = np.random.default_rng(9) + rows = [] + for u in range(60): + g = 3 if u % 2 == 0 else 5 + for t in range(1, 9): + y = 1.0 + 0.1 * t + u * 0.01 + (1.0 if t >= g else 0.0) + rng.normal(0, 0.3) + rows.append({"unit": u, "time": t, "outcome": y, "first_treat": g}) + data = pd.DataFrame(rows) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = TwoStageDiD().fit( + data, + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + prop5 = [ + h + for h, d in res.event_study_effects.items() + if d["n_obs"] > 0 and not np.isfinite(d["effect"]) + ] + assert prop5, "fixture must produce Proposition-5 horizons" + assert res.event_study_vcov_index is not None + assert not set(prop5) & set(res.event_study_vcov_index) + surface = build_event_study_surface(res) + assert surface.vcov is not None + assert set(surface.vcov_index.tolist()).issubset(set(surface.event_time.tolist())) + + +def _stacked_panel(seed=42): + rng = np.random.default_rng(seed) + rows = [] + for u in range(60): + g = [4, 6, 0][u % 3] + for t in range(1, 11): + y = 1.0 + 0.1 * t + u * 0.01 + (1.5 if g and t >= g else 0.0) + rng.normal(0, 0.3) + rows.append({"unit": u, "time": t, "outcome": y, "first_treat": g}) + return pd.DataFrame(rows) + + +def test_stacked_hc2_bm_per_row_df(): + from diff_diff import StackedDiD + + data = _stacked_panel() + kwargs = dict( + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + bm = StackedDiD(kappa_pre=2, kappa_post=2, vcov_type="hc2_bm").fit(data, **kwargs) + hc1 = StackedDiD(kappa_pre=2, kappa_post=2).fit(data, **kwargs) + + # hc2_bm: every finite-p non-reference row carries its own BM df, + # bit-equal to the producer's event_study_df provenance dict. + s_bm = build_event_study_surface(bm) + finite_p = np.isfinite(s_bm.p_value) + assert finite_p.any() + assert np.isfinite(s_bm.df[finite_p]).all() + for et, df_val in zip(s_bm.event_time.tolist(), s_bm.df.tolist()): + if et in bm.event_study_df and np.isfinite(bm.event_study_df[et]): + assert df_val == bm.event_study_df[et] + + # hc1 non-survey: normal-theory inference -> no df anywhere. + s_hc1 = build_event_study_surface(hc1) + assert np.isnan(s_hc1.df).all() + + +def test_stacked_survey_df_broadcast(): + from diff_diff import StackedDiD + from diff_diff.survey import SurveyDesign + + data = _stacked_panel() + data["w"] = 1.0 + data["strat"] = data["unit"] % 2 + data["psu_id"] = data["unit"] + design = SurveyDesign(weights="w", weight_type="pweight", strata="strat", psu="psu_id") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = StackedDiD(kappa_pre=2, kappa_post=2).fit( + data, + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + survey_design=design, + ) + assert res.survey_metadata is not None and res.survey_metadata.df_survey is not None + surface = build_event_study_surface(res) + finite_p = np.isfinite(surface.p_value) + assert finite_p.any() + vals = set(surface.df[finite_p].tolist()) + assert vals == {float(max(res.survey_metadata.df_survey, 1))} + + +def test_two_stage_survey_df_threaded(): + from diff_diff.survey import SurveyDesign + from diff_diff.two_stage import TwoStageDiD + + data = _stacked_panel() + data["w"] = 1.0 + data["strat"] = data["unit"] % 2 + data["psu_id"] = data["unit"] + design = SurveyDesign(weights="w", weight_type="pweight", strata="strat", psu="psu_id") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = TwoStageDiD().fit( + data, + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + survey_design=design, + ) + assert res.event_study_df is not None and res.event_study_df > 0 + surface = build_event_study_surface(res) + finite_p = np.isfinite(surface.p_value) + assert finite_p.any() + assert set(surface.df[finite_p].tolist()) == {float(res.event_study_df)} + + +def test_stacked_replicate_weight_vcov_diag_matches_ses(): + # Replicate-weight designs REASSIGN the coefficient covariance before + # the ES extraction loop (replicate refit), so StackedDiD persists the + # replicate VCV sub-block and the reported ES SEs remain exactly its + # diagonal - the every-inference-mode guarantee. + from diff_diff import StackedDiD + from diff_diff.survey import SurveyDesign + + data = _stacked_panel() + data["w"] = 1.0 + units = np.sort(data["unit"].unique()) + n_rep = 8 + unit_pos = {u: i for i, u in enumerate(units)} + rows = data["unit"].map(unit_pos).to_numpy() + per = max(len(units) // n_rep, 1) + rep_cols = [] + for r in range(n_rep): + w_r = np.ones(len(units)) + w_r[r * per : min((r + 1) * per, len(units))] = 0.0 + nz = w_r > 0 + w_r[nz] = w_r[nz] * n_rep / (n_rep - 1) + data[f"rep_{r}"] = w_r[rows] + rep_cols.append(f"rep_{r}") + design = SurveyDesign( + weights="w", weight_type="pweight", replicate_weights=rep_cols, replicate_method="JK1" + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = StackedDiD(kappa_pre=2, kappa_post=2).fit( + data, + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + survey_design=design, + ) + assert res.event_study_vcov is not None + ses = {h: d["se"] for h, d in res.event_study_effects.items()} + diag = np.sqrt(np.maximum(np.diag(res.event_study_vcov), 0.0)) + for i, h in enumerate(res.event_study_vcov_index): + np.testing.assert_allclose(diag[i], ses[h], rtol=1e-14) + surface = build_event_study_surface(res) + assert surface.vcov is not None + + +def _cs_cluster_panel(seed=3): + rng = np.random.default_rng(seed) + rows = [] + for u in range(80): + g = [4, 6, 0][u % 3] + for t in range(1, 9): + y = 1 + 0.1 * t + u * 0.01 + (1.0 if g and t >= g else 0.0) + rng.normal(0, 0.3) + rows.append({"unit": u, "time": t, "outcome": y, "first_treat": g, "st": u % 8}) + return pd.DataFrame(rows) + + +def test_cs_bare_cluster_df_fallback(): + from diff_diff import CallawaySantAnna + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = CallawaySantAnna(cluster="st").fit( + _cs_cluster_panel(), + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + # Bare-cluster synthesize path: the ES rows' safe_inference used G-1, + # recorded on BOTH channels (primary event_study_df; df_inference keeps + # its narrow HonestDiD contract). + assert res.event_study_df == float(res.df_inference) == 7.0 + surface = build_event_study_surface(res) + finite_p = np.isfinite(surface.p_value) + assert set(surface.df[finite_p].tolist()) == {7.0} + + +def test_cs_bare_cluster_bootstrap_df_is_nan(): + from diff_diff import CallawaySantAnna + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = CallawaySantAnna(cluster="st", n_bootstrap=20, seed=1).fit( + _cs_cluster_panel(), + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + # Bootstrap overrode the stored ES p/CIs with percentile values that + # never used the analytical df: the surface must show NO df even though + # df_inference (the HonestDiD contract field) still carries G-1 on the + # container. Bootstrap p-values are finite, so only the gated + # df_inference fallback - not the NaN-p mask - prevents the leak. + assert res.df_inference == 7 + assert res.event_study_df is None + surface = build_event_study_surface(res) + assert np.isfinite(surface.p_value).any() + assert np.isnan(surface.df).all() + + +def test_cs_survey_min_df_on_surface(): + from diff_diff import CallawaySantAnna + from diff_diff.survey import SurveyDesign + + data = _cs_cluster_panel() + data["w"] = 1.0 + design = SurveyDesign(weights="w", weight_type="pweight", strata="st", psu="unit") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = CallawaySantAnna().fit( + data, + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + survey_design=design, + ) + # Explicit-survey fits: df_inference stays None (PR #487 narrow + # HonestDiD contract - channel separation), while the surface carries + # the ONE df every ES row's safe_inference_batch actually applied. + assert res.df_inference is None + assert res.event_study_df is not None and res.event_study_df > 0 + surface = build_event_study_surface(res) + finite_p = np.isfinite(surface.p_value) + assert finite_p.any() + assert set(surface.df[finite_p].tolist()) == {float(res.event_study_df)} + + +def test_mpd_hc2_bm_per_period_df(): + from diff_diff import MultiPeriodDiD + + rng = np.random.default_rng(5) + rows = [] + for u in range(50): + tr = u % 2 + # Unbalanced panel: later periods observed for a shrinking subset, + # so the per-period BM Satterthwaite dfs genuinely differ. + t_max = 8 - (u % 4) + for t in range(1, t_max + 1): + y = 1 + 0.1 * t + 0.3 * tr + (0.8 if tr and t >= 5 else 0.0) + rng.normal(0, 0.3) + rows.append({"unit": u, "time": t, "outcome": y, "treated": tr}) + data = pd.DataFrame(rows) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = MultiPeriodDiD(vcov_type="hc2_bm", cluster="unit").fit( + data, outcome="outcome", treatment="treated", time="time", unit="unit" + ) + finite_vals = [v for v in res.event_study_df.values() if np.isfinite(v)] + assert finite_vals + # Per-period dfs vary across rows on an unbalanced panel... + assert len({round(v, 6) for v in finite_vals}) > 1 + # ...and are NOT the post-average contrast df broadcast. + assert any(v != res.inference_df for v in finite_vals) + surface = build_event_study_surface(res) + for et, df_val in zip(surface.event_time.tolist(), surface.df.tolist()): + if et in res.event_study_df and np.isfinite(res.event_study_df[et]): + assert df_val == res.event_study_df[et] + + +def test_lpdid_per_horizon_df(surfaces): + native, surface = surfaces["LPDiD"] + frame = native.event_study + ncl = dict(zip(frame["horizon"].tolist(), frame["n_clusters"].tolist())) + for et, df_val in zip(surface.event_time.tolist(), surface.df.tolist()): + if et == -1: + assert np.isnan(df_val) # synthetic base row: no provenance + elif np.isfinite(df_val): + # Non-survey cluster rule: realized per-horizon G - 1. + assert df_val == ncl[et] - 1 + # The native frame schema is UNCHANGED (df lives on event_study_df only). + assert "df" not in frame.columns diff --git a/tests/test_two_stage.py b/tests/test_two_stage.py index 8d66bacb..5d51b155 100644 --- a/tests/test_two_stage.py +++ b/tests/test_two_stage.py @@ -2635,3 +2635,111 @@ def test_main_fit_zero_weight_unit_warns_once_and_fits(self): assert len(nan_ytilde) >= 1 # main-fit warning still fires assert np.isfinite(r.overall_att) assert np.isfinite(r.overall_se) and r.overall_se > 0 + + +class TestEventStudyVcovPersistence: + """M-092 follow-up: the full Gardner-GMM ES VCV + df provenance are + persisted on TwoStageDiDResults, gated by inference mode.""" + + @staticmethod + def _panel(seed=7): + rng = np.random.default_rng(seed) + rows = [] + for u in range(90): + g = [4, 6, 0][u % 3] + for t in range(1, 11): + y = 1.0 + 0.1 * t + u * 0.01 + (1.2 if g and t >= g else 0.0) + rng.normal(0, 0.3) + rows.append( + { + "unit": u, + "time": t, + "outcome": y, + "first_treat": g if g else np.nan, + } + ) + return pd.DataFrame(rows) + + def test_analytical_fit_persists_vcov_diag_matches_ses(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = TwoStageDiD().fit( + self._panel(), + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + assert res.event_study_vcov is not None + assert res.event_study_vcov_index is not None + ses = {h: d["se"] for h, d in res.event_study_effects.items()} + diag = np.sqrt(np.maximum(np.diag(res.event_study_vcov), 0.0)) + for i, h in enumerate(res.event_study_vcov_index): + np.testing.assert_allclose(diag[i], ses[h], rtol=1e-14) + # Non-survey GMM inference is normal-theory: no df provenance. + assert res.event_study_df is None + + def test_bootstrap_clears_vcov_and_df(self): + # Bootstrap replaces the stored ES se/p/CI with percentile values: + # the analytical matrix's diagonal no longer matches the stored SEs + # and no df governed the stored inference - both must be None. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = TwoStageDiD(n_bootstrap=20, seed=1).fit( + self._panel(), + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + assert res.bootstrap_results is not None + assert res.event_study_vcov is None + assert res.event_study_vcov_index is None + assert res.event_study_df is None + + def test_replicate_weight_survey_clears_vcov_and_threads_final_df(self): + # Replicate-weight designs: reported ES SEs come from the replicate + # VCV's mixed [overall, ES, groups] layout, so the analytical GMM + # matrix must NOT be persisted; the final (possibly + # dropped-replicate-tightened) survey df that every recomputed ES + # row's safe_inference used IS threaded. + data, rep_cols = _add_survey_cols(self._panel(), n_rep=8) + design = SurveyDesign(weights="w", replicate_weights=rep_cols, replicate_method="JK1") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = TwoStageDiD().fit( + data, + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + survey_design=design, + ) + assert res.event_study_vcov is None + assert res.event_study_vcov_index is None + assert res.event_study_df is not None and res.event_study_df > 0 + assert res.event_study_df == float(res.survey_metadata.df_survey) + from diff_diff.results_base import build_event_study_surface + + surface = build_event_study_surface(res) + finite_p = np.isfinite(surface.p_value) + assert finite_p.any() + assert set(surface.df[finite_p].tolist()) == {float(res.event_study_df)} + + def test_group_only_aggregate_has_no_es_vcov(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = TwoStageDiD().fit( + self._panel(), + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="group", + ) + assert res.event_study_effects is None + assert res.event_study_vcov is None + assert res.event_study_vcov_index is None + assert res.event_study_df is None