diff --git a/CHANGELOG.md b/CHANGELOG.md index fa8b97873..26aa08f49 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,99 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- **Event-study container consumers + dCDH post-fit `aggregate()`** (v4 + program 2(b) PR-1; ledger row [M-026] flips to shimmed, rows [M-092] / + [M-093] amended). + - `compute_honest_did`, `compute_pretrends_power` and + `plot_event_study` now accept the unified `EventStudyResults` + container from `CallawaySantAnnaResults.aggregate('event_study')` + directly, and `plot_honest_event_study` renders an + `HonestDiDResults` whose `original_results` is such a container - + the post-fit route no longer dead-ends at the consumers. HonestDiD + bounds are identical across routes (route-parity gated at 1e-14); + PreTrendsPower's extraction is bit-exact (its Monte-Carlo power + carries scipy's inherent MVN-CDF jitter either way). Admission in + HonestDiD/PreTrendsPower is SOURCE-SCOPED to CS-produced containers: + dCDH l1 containers are rejected BY DESIGN (their placebo semantics + need HonestDiD's native dCDH branch), other producers arrive with + their own `aggregate()` migrations; the plotters take no source + guard. Inputs whose coefficients mix reference normalizations fail + closed in HonestDiD and PreTrendsPower on BOTH routes - a container + with multiple marked reference rows, and any universal-base fit + whose cohorts carry more than one positional-base event time (the + `reference_event_times` guard below). + - `EventStudyResults` gains four optional provenance fields (declared + last - positional compatibility preserved): `base_period`, + `anticipation`, `df_survey` (the fit's resolved scalar inference + df; `0.0` = replicate design with undefined df, which the per-row + `df` column cannot encode), and `reference_event_times` (the + common-reference provenance described below). Threaded by the + builders, including the + requested-but-empty path; `to_dict()` gains the four keys. The + `df_survey` typing widened `Optional[int]` -> `Optional[float]` + through the HonestDiD call path. + - `ChaisemartinDHaultfoeuille.fit(aggregate=)` is DEPRECATED + ([M-026]): it never computed anything (every non-None value raised + "reserved for Phase 3"), so supplying it warns and a non-None value + raises `ValueError` pointing at the new post-fit route. + `ChaisemartinDHaultfoeuilleResults.aggregate()` ships as a pure + VIEW: `'event_study'` returns the unified container (Phase-1 fits + the 2-row l=1 view), `'simple'` a one-row `AggregationResult` + relaying the overall estimand bit-exactly with estimand-aware + `target`/`n`/`n_kind` (DID_M/switcher-cells, DID_1/groups, + delta/no-scalar-count) and a `df` column resolved from the ACTUAL + inference path (the `L_max>=2` delta keeps its analytical survey df + even under bootstrap, per the REGISTRY delta-SE note). Because + nothing is recomputed, bootstrap fits are permitted - unlike + CallawaySantAnna's kit-based `aggregate()`, which stays fail-closed. + - Consumer hardening that landed with the wave: `plot_event_study` + renders MULTI-reference containers via row-aligned hollow anchors + (never dropped, never filled estimates; explicit `reference_period=` + re-basing of such a surface around a non-anchor fails closed), and + warns when a container's stored intervals sit at a different level + than the requested `alpha` (stored bootstrap/t intervals cannot be + re-leveled from the SE) - `plot_honest_event_study` names the same + mismatch against the HonestDiD level. HonestDiD's MPD and dCDH + branches now drop `se == 0` rows exactly like the CS/container + branches (undefined inference is never laundered into finite + bounds). PreTrendsPower warns on CS `base_period='varying'` input on + BOTH routes: its `linear` violation benchmark assumes level + coefficients against a common reference, while varying-base + pre-treatment effects are consecutive-period comparisons (REGISTRY + PreTrendsPower note; violation-vector transformation tracked in + TODO.md). + - Common-reference guard: CS `base_period='universal'` fits (and + their containers) record `reference_event_times` - the distinct + per-cohort positional-base event times. On gapped time grids the + bases land at different event times, and a cohort's base can + OVERLAP another cohort's estimated horizon where no reference-only + row marks it; HonestDiD and PreTrendsPower now FAIL CLOSED on both + input routes instead of silently returning bounds/power over + coefficients normalized against different bases (REGISTRY HonestDiD + common-reference-guard note). Regular (ungapped) grids carry a + single entry and are unaffected; varying-base fits carry None. + - HonestDiD's no-reference (varying-base) pre/post split is now + anticipation-aware on both routes: with `anticipation=k` the window + `[e=-k, -1]` carries anticipated treatment effects (the REGISTRY + contract PreTrendsPower already applied), so the clean pre-trend + set is `e < -k` and `beta_post` starts at `-k` - previously the + split at 0 misclassified anticipated effects as pre-trend + coefficients. The pretrends varying-base warning is scoped to the + `linear` violation benchmark it concerns. + - PreTrendsPower's `linear` violation weights are REFERENCE-ANCHORED + on the CS-universal, SunAbraham and container routes (`|t - t_ref|`; + Roth's violation is normalized like the coefficients and vanishes at + the omitted period - MPD already anchored via its reference helper): + raw treatment-relative labels overstated each pre-period's + hypothesized violation by the reference offset and correspondingly + understated the γ-unit MDV. Universal-base CS/SA linear + pretrends power and MDV numbers CHANGE to the Roth-correct values + (e.g. weights `[3, 2]` -> `[2, 1]` for pre `[-3, -2]` around + `t_ref = -1`); varying-base input keeps raw labels under its + existing warning (REGISTRY linear-convention note records the + correction). + ### Deprecated - **Semantic rename wave** (v4 program 2(c)-ii PR-B, ledger rows [M-030] [M-031] [M-043] [M-044] [M-045]..[M-047] [M-084] [M-086] [M-087] [M-095] diff --git a/TODO.md b/TODO.md index 01436c8ea..301e30d29 100644 --- a/TODO.md +++ b/TODO.md @@ -21,10 +21,10 @@ Related tracking surfaces: | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| +| PreTrendsPower `violation='linear'` on CS `base_period='varying'` input targets the wrong alternative: `δ_pre = M · \|t\|` assumes level coefficients against a common reference, but varying-base pre-treatment effects are consecutive-period comparisons (constant increments under a linear trend). Both CS-sourced routes now WARN (REGISTRY PreTrendsPower Note), and universal-base GAPPED grids fail closed via the `reference_event_times` common-reference guard; what remains is the varying-base resolution - either transforming the violation vector through each coefficient's actual base mapping (needs per-horizon base provenance) or requiring `base_period='universal'` for the linear benchmark - a per-estimator methodology decision with a hand-calculated linear-violation gate | `diff_diff/pretrends.py` | 2(b) PR-1 R5 | Mid | Medium | | `run_placebo_test`/`run_all_placebo_tests` expose ONE overloaded `time` param (calendar for timing/group tests, 0/1 post dummy forwarded into the M-137/M-138 callees - works only because a two-period 0/1 calendar column is both); redesign the wrapper signatures (dual `time=`/`post=` params or a split) so each callee gets an honestly-named argument | `diff_diff/diagnostics.py` | 2(c)-ii PR-B | Mid | Medium | | `absorbed_fe_rank` N>=3 general rank: the helper keeps `sum(levels-1)` for 3+ absorbed dims, exact for independent connected dims but an over-count for duplicated/nested triples (measured `a(5),b(4),c==b(4)`: true 7 vs formula 10) and for disconnected N-way graphs. Two-way is component-exact. Deriving general N-way FE rank is a hypergraph problem; do it with a reference (fixest's `fixef.rm`/reghdfe df_a) rather than a guess. See `docs/methodology/variance-conventions.md` D3. | `diff_diff/utils.py` | #variance-inventory | Mid | Low | | `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 | -| Make the post-fit `results.aggregate("event_study")` container consumable downstream. `EventStudyResults` is rejected by all THREE consumers that read a CS event study — `compute_honest_did` (`honest_did.py`, dispatches on `CallawaySantAnnaResults` and raises `TypeError`), `compute_pretrends_power` (`pretrends.py`, same), and `plot_event_study` (`visualization`, same) — so `fit(aggregate="event_study")` is still the only route for them and their error messages say so explicitly. Needs an `EventStudyResults` branch in each extraction path (consuming `event_time` / `is_reference` / `vcov` / `vcov_index` / per-row `df`) PLUS `base_period` and `anticipation` provenance, which the unified container does not carry and HonestDiD needs for its universal-base-period warning and pre-period classification. Gate with end-to-end tests: `compute_honest_did(res.aggregate("event_study"))` at `base_period="universal"`, and `compute_pretrends_power(...)` at `anticipation=1`. | `diff_diff/honest_did.py`, `diff_diff/pretrends.py`, `diff_diff/results_base.py` | #726 | Mid | Medium | | Re-run the R-dependent benchmark refresh so `docs/benchmarks.rst`'s TWFE "SE Rel Diff 0.1%" cell reflects the 3.9 K_reference convergence (expected 0.1% -> 0.0%; the table is generated, never hand-edited — the movement is noted in the CHANGELOG entry). | `docs/benchmarks.rst`, `benchmarks/R/` | #variance-inventory | Quick | Low | | `SunAbraham`: a cohort not observed at its own reference relative period (`e = -1 - anticipation`) makes that cohort's block collinear, so QR drops an unnamed column (`dropping 1 of 12 columns (column 9)`) and `overall_att` comes back **NaN**. Found by auditing the sibling estimator while fixing the ETWFE analogue (#724); PRE-EXISTING, not introduced there. Lower severity than #724 — that returned a silently WRONG finite number, this returns NaN with a rank warning — but the event-study surface still looks complete, so a user may not notice the loss. SA already omits its reference explicitly and tracks `_reference_observed`, so the fix is per-cohort support for that flag rather than the ETWFE-style redesign. | `diff_diff/sun_abraham.py` | #724-audit | Mid | Low | | Define `N_g` (W2025 Eqs. 7.4/7.6) for UNBALANCED panels where comparison-support filtering removes every observation of some units in an estimated cohort, then replace the fail-closed guard with the defined behavior. `_n_g_per_cohort` is read off the final sample, so those units vanish from the cohort-share weights; measured on a cohort supplied with 100 units of which 90 appear only at a dropped period, `aggregate(weights="cohort_share")` moves 1.8078 -> 3.8157. The paper assumes a balanced panel and does not say whether `N_g` counts the supplied cohort or the surviving units, and the two disagree materially, so `aggregate` currently raises naming the cohorts and counts ([M-125]); `weights="cell"` is unaffected and balanced panels never trip it. Settle the estimand (likely: count the supplied cohort, since ATT(g,t) is a cohort-level quantity, but that weights units with no retained observation) and gate with a test computing Eq. 7.4 by hand on unequal cohort sizes. | `diff_diff/wooldridge_results.py`, `diff_diff/wooldridge.py` | #729-followup | Mid | Medium | diff --git a/diff_diff/aggregation.py b/diff_diff/aggregation.py index 6340c9950..34bcf6a55 100644 --- a/diff_diff/aggregation.py +++ b/diff_diff/aggregation.py @@ -25,7 +25,7 @@ """ from dataclasses import dataclass, field -from typing import Any, Dict, Optional, Tuple +from typing import Any, ClassVar, Dict, Optional, Tuple import numpy as np import pandas as pd @@ -134,9 +134,13 @@ class AggregationResult(BaseResults): aligned estimands over the same labels (ContinuousDiD's ATT(d) and ACRT(d) become 2N rows). ``"att"`` where an estimator has one. att, se, t_stat, p_value : np.ndarray - The canonical quintet, per row. On a bootstrapped fit ``t_stat`` / - ``p_value`` / the interval are the producer's percentile-bootstrap - statistics carried through unchanged - NOT recomputed analytically. + The canonical quintet, per row, carrying WHATEVER inference the fit + stored - never recomputed. On a bootstrapped fit that usually means + the producer's percentile-bootstrap statistics carried through + unchanged; view-relay producers can mix regimes per row where the + fit itself did (dCDH's ``L_max >= 2`` cost-benefit delta keeps + analytical ``safe_inference`` even under ``n_bootstrap > 0`` - see + the REGISTRY Phase 2 cost-benefit delta SE note). conf_int_lower, conf_int_upper : np.ndarray Interval bounds at the fit's ``alpha``. n : np.ndarray @@ -156,8 +160,10 @@ class AggregationResult(BaseResults): one would be a fabricated number. df : np.ndarray Per-row inference degrees of freedom, NaN where none governed the - stored p-value. Entirely NaN on bootstrap fits, whose percentile - inference uses no df. + stored p-value. NaN on percentile-bootstrap rows (no df governs + them); a bootstrapped fit's rows can still carry a finite df where + the fit kept analytical inference for that row (the dCDH delta + case above). alpha : float Significance level the interval was computed at. estimator : str or None @@ -460,12 +466,14 @@ class AggregationMixin: #: Aggregation types this results class implements. A subset of #: :data:`AGGREGATION_VOCABULARY` plus any documented per-estimator extra. - _AGGREGATE_SUPPORTED: Tuple[str, ...] = () + #: ClassVar so that dataclass results classes overriding these hooks with + #: an annotation do not grow a spurious ``__init__`` field. + _AGGREGATE_SUPPORTED: ClassVar[Tuple[str, ...]] = () #: Types for which ``balance_e`` is meaningful. CallawaySantAnna threads it #: only through event-study aggregation, so accepting it elsewhere would #: silently ignore a user's argument. - _AGGREGATE_BALANCE_E_TYPES: Tuple[str, ...] = ("event_study",) + _AGGREGATE_BALANCE_E_TYPES: ClassVar[Tuple[str, ...]] = ("event_study",) def aggregate( self, @@ -505,10 +513,14 @@ def aggregate( ) raise ValueError(f"Unsupported aggregation type {type!r}. Supported: {known}.{extra}") if balance_e is not None and type not in self._AGGREGATE_BALANCE_E_TYPES: - usable = ", ".join(repr(t) for t in self._AGGREGATE_BALANCE_E_TYPES) + if self._AGGREGATE_BALANCE_E_TYPES: + usable = ", ".join(repr(t) for t in self._AGGREGATE_BALANCE_E_TYPES) + applies = f"It applies to: {usable}." + else: + applies = "It applies to no aggregation type on this estimator." raise ValueError( f"balance_e is not used by aggregate(type={type!r}) and would be " - f"silently ignored. It applies to: {usable}." + f"silently ignored. {applies}" ) self._aggregate_validate_weights(weights) return self._aggregate_compute(type, weights=weights, balance_e=balance_e) diff --git a/diff_diff/chaisemartin_dhaultfoeuille.py b/diff_diff/chaisemartin_dhaultfoeuille.py index 690afee0c..6d119c7e8 100644 --- a/diff_diff/chaisemartin_dhaultfoeuille.py +++ b/diff_diff/chaisemartin_dhaultfoeuille.py @@ -37,7 +37,12 @@ import pandas as pd from diff_diff._base import BaseEstimator -from diff_diff._deprecation import NOT_SUPPLIED, require_arg, resolve_renamed_kwarg +from diff_diff._deprecation import ( + NOT_SUPPLIED, + require_arg, + resolve_renamed_kwarg, + warn_deprecated_kwarg, +) from diff_diff.chaisemartin_dhaultfoeuille_bootstrap import ( ChaisemartinDHaultfoeuilleBootstrapMixin, ) @@ -401,7 +406,9 @@ class ChaisemartinDHaultfoeuille(ChaisemartinDHaultfoeuilleBootstrapMixin, BaseE variance (BRR/Fay/JK1/JKn/SDR) - TWFE decomposition diagnostic from Theorem 1 of AER 2020 - Only ``aggregate`` on :meth:`fit` still raises ``NotImplementedError``. + Aggregation is a post-fit step: ``results.aggregate('event_study')`` / + ``results.aggregate('simple')`` (the deprecated fit-time ``aggregate=`` + never computed anything and now warns; row M-026). Parameters ---------- @@ -876,7 +883,7 @@ def fit( time: Any = NOT_SUPPLIED, treatment: Any = NOT_SUPPLIED, # ---------- forward-compat parameters ---------- - aggregate: Optional[str] = None, + aggregate: Any = NOT_SUPPLIED, L_max: Optional[int] = None, covariates: Any = NOT_SUPPLIED, trends_linear: Optional[bool] = None, @@ -916,8 +923,14 @@ def fit( non-binary (ordinal or continuous) treatment are supported. Non-binary treatment requires ``L_max >= 1``. aggregate : str, optional - **Reserved for Phase 3.** Must be ``None``; any other value - raises ``NotImplementedError``. + DEPRECATED since 3.9, removed in 4.0 (ledger row M-026). + Passing it emits a ``FutureWarning``. The parameter never + computed aggregations here (it was reserved and any non-None + value has always raised); aggregation is a POST-FIT step: + ``results.aggregate('event_study')`` for the multi-horizon + surface, ``results.aggregate('simple')`` for the overall + estimand. A non-None value now raises ``ValueError`` pointing + at that post-fit route. L_max : int, optional Maximum event-study horizon. When set, computes ``DID_l`` for ``l = 1, ..., L_max`` using the per-group building block @@ -1083,8 +1096,29 @@ def fit( # ------------------------------------------------------------------ # Step 2: Forward-compat gates # ------------------------------------------------------------------ + # M-026: fit(aggregate=) is deprecated - the aggregation entry point + # is post-fit (results.aggregate()). The param never computed + # anything here (any non-None value has always raised), so the shim + # warns on supply and rejects non-None values with a pointer at the + # post-fit route. + if aggregate is not NOT_SUPPLIED: + warn_deprecated_kwarg( + "ChaisemartinDHaultfoeuille.fit", + "aggregate", + "aggregation is a post-fit step: " + "results.aggregate('event_study') / results.aggregate('simple')", + ) + if aggregate is not None: + raise ValueError( + f"fit(aggregate={aggregate!r}) never computed aggregations " + "on ChaisemartinDHaultfoeuille (the parameter was reserved " + "and always raised). Aggregate post-fit instead: " + "results.aggregate('event_study') for the multi-horizon " + "surface (Phase-1 fits without L_max return the 2-row " + "l=1 view), or results.aggregate('simple') for the " + "overall estimand." + ) _check_forward_compat_gates( - aggregate=aggregate, L_max=L_max, controls=controls, trends_linear=trends_linear, @@ -4385,7 +4419,6 @@ def fit( def _check_forward_compat_gates( - aggregate: Optional[str], L_max: Optional[int], controls: Optional[List[str]], trends_linear: Optional[bool], @@ -4395,15 +4428,11 @@ def _check_forward_compat_gates( """Raise ``NotImplementedError`` for any non-default Phase 3 parameter. Phase 2 parameters (``L_max``) are validated inline in ``fit()`` - after period detection. The ``aggregate`` parameter is still - reserved for Phase 3. + after period detection. The deprecated ``aggregate`` parameter is + resolved in ``fit()`` itself (row M-026: FutureWarning on supply, + ValueError on any non-None value, pointing at the post-fit + ``results.aggregate()`` route). """ - if aggregate is not None: - raise NotImplementedError( - f"aggregate={aggregate!r} is reserved for Phase 3 of dCDH. " - "Multi-horizon event study effects are computed automatically " - "when L_max is set. See ROADMAP.md Phase 3." - ) # L_max is validated inline in fit() after period detection (needs # the period count). Not gated here. # controls gate lifted — DID^X covariate residualization implemented. diff --git a/diff_diff/chaisemartin_dhaultfoeuille_results.py b/diff_diff/chaisemartin_dhaultfoeuille_results.py index 443131723..ff1712c45 100644 --- a/diff_diff/chaisemartin_dhaultfoeuille_results.py +++ b/diff_diff/chaisemartin_dhaultfoeuille_results.py @@ -23,7 +23,7 @@ from __future__ import annotations from dataclasses import dataclass, field -from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple +from typing import TYPE_CHECKING, Any, ClassVar, Dict, List, Optional, Tuple if TYPE_CHECKING: from diff_diff.honest_did import HonestDiDResults @@ -32,8 +32,9 @@ import pandas as pd from diff_diff._deprecation import deprecated_field_property +from diff_diff.aggregation import AggregationMixin, AggregationResult from diff_diff.results import _get_significance_stars -from diff_diff.results_base import BaseResults +from diff_diff.results_base import BaseResults, build_event_study_surface __all__ = [ "ChaisemartinDHaultfoeuilleResults", @@ -182,7 +183,7 @@ class DCDHBootstrapResults: @dataclass -class ChaisemartinDHaultfoeuilleResults(BaseResults): +class ChaisemartinDHaultfoeuilleResults(BaseResults, AggregationMixin): """ Results from de Chaisemartin-D'Haultfoeuille (dCDH) Phase 1 estimation. @@ -700,6 +701,87 @@ def __setstate__(self, state: Dict[str, Any]) -> None: state["units"] = state.pop("groups") self.__dict__.update(state) + # ------------------------------------------------------------------ + # Post-fit aggregation (row M-026, on the M-122 contract): dCDH's + # aggregate() is a pure VIEW over stored fields - nothing is + # recomputed - so bootstrap fits are PERMITTED (unlike + # CallawaySantAnna's kit-based recompute, which fails closed): each + # row relays whatever inference the fit stored. + # ------------------------------------------------------------------ + # ClassVar: on a dataclass a bare annotation would turn this routing + # configuration into an ``__init__`` field. + _AGGREGATE_SUPPORTED: ClassVar[Tuple[str, ...]] = ("simple", "event_study") + # dCDH has no balance_e machinery on any aggregation level. + _AGGREGATE_BALANCE_E_TYPES: ClassVar[Tuple[str, ...]] = () + + def _overall_inference_df(self) -> float: + """The df the STORED overall p-value/CI actually used (NaN = none). + + Deliberately NOT the ``event_study_df`` channel alone: that + provenance is cleared under bootstrap, while the ``L_max >= 2`` + cost-benefit delta keeps analytical ``safe_inference`` with a + (possibly finite survey) df even when ``n_bootstrap > 0`` - + REGISTRY ``Note (Phase 2 cost-benefit delta SE)``. Resolution: + + - percentile-bootstrap ``DID_M`` / ``DID_1`` -> NaN (no df); + - analytical paths (including the delta under bootstrap) -> + ``event_study_df`` when recorded, else the post-fit-refreshed + ``survey_metadata.df_survey``, else NaN (z-inference). + """ + is_delta = self.L_max is not None and self.L_max >= 2 + if self.bootstrap_results is not None and not is_delta: + return float("nan") + if self.event_study_df is not None: + return float(self.event_study_df) + sm = self.survey_metadata + if sm is not None and getattr(sm, "df_survey", None) is not None: + return float(sm.df_survey) + return float("nan") + + def _aggregate_compute( + self, level: str, *, weights: Optional[str], balance_e: Optional[int] + ) -> Any: + if level == "event_study": + # The unified container over ``event_study_effects``: Phase-1 + # fits (``L_max=None``) return the 2-row l=1 view; ``L_max >= 1`` + # fits the multi-horizon surface (``l1_first_switch`` + # convention, ``n_kind`` per ``results_base._from_dcdh``). + return build_event_study_surface(self) + + # level == "simple": one-row view of the stored overall estimand - + # DID_M (L_max=None), DID_1 (L_max==1), or the cost-benefit delta + # (L_max>=2; all-NaN by design under trends_linear, where the + # estimand label points at linear_trends_effects instead). + if self.L_max is None: + # N_S: switching (g, t) CELLS (one group may contribute several). + n_val, n_kind = float(self.n_switcher_cells), "switcher_cells" + elif self.L_max == 1: + # n_switcher_cells is repurposed to the eligible-switcher GROUP + # count N_1 once L_max >= 1 (class docstring). + n_val, n_kind = float(self.n_switcher_cells), "groups" + else: + # The delta is weighted over horizon-specific N_l - no truthful + # scalar count exists. + n_val, n_kind = float("nan"), None + ci = self.overall_conf_int if self.overall_conf_int is not None else (np.nan, np.nan) + return AggregationResult( + level="simple", + label=np.array(["overall"], dtype=object), + target=np.array([self._estimand_label()], dtype=object), + att=np.array([self.overall_att], dtype=float), + se=np.array([self.overall_se], dtype=float), + t_stat=np.array([self.overall_t_stat], dtype=float), + p_value=np.array([self.overall_p_value], dtype=float), + conf_int_lower=np.array([ci[0]], dtype=float), + conf_int_upper=np.array([ci[1]], dtype=float), + n=np.array([n_val], dtype=float), + df=np.array([self._overall_inference_df()], dtype=float), + alpha=self.alpha, + n_kind=n_kind, + weight=np.array([1.0], dtype=float), + estimator=type(self).__name__.replace("Results", ""), + ) + # ------------------------------------------------------------------ # Repr / properties # ------------------------------------------------------------------ diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index f89639d3b..93a748600 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -240,17 +240,21 @@ by_cohort = results.aggregate('group') # -> AggregationResult ``` Fit-time `aggregate=` / `balance_e=` are DEPRECATED since 3.9 (removed in 4.0, -ledger rows M-020 / M-117) and emit a `FutureWarning`. They remain the required -path in two cases: - -1. **Bootstrap inference.** `aggregate()` raises on a bootstrapped fit rather - than substituting analytical inference for percentile-bootstrap statistics. -2. **Downstream consumers.** `plot_event_study`, `compute_honest_did` and - `compute_pretrends_power` still read the fit-time surface and do not yet - accept the post-fit container. +ledger rows M-020 / M-117) and emit a `FutureWarning`. The downstream consumers +all accept the post-fit container directly - `plot_event_study`, +`compute_honest_did` and `compute_pretrends_power` each take +`results.aggregate('event_study')` - so the only case still requiring the +fit-time path is **bootstrap inference**: CallawaySantAnna's `aggregate()` +raises on a bootstrapped fit rather than substituting analytical inference for +percentile-bootstrap statistics. ```python -# Both cases: aggregate at fit time. +# Post-fit route (recommended): aggregate once, feed any consumer. +es = results.aggregate('event_study') +plot_event_study(es) +honest = compute_honest_did(es, M=1.0) + +# Bootstrap fits only: aggregate at fit time. boot = CallawaySantAnna(estimation_method="dr", n_bootstrap=999, seed=42) plotted = boot.fit(data, outcome='outcome', unit='unit', time='period', first_treat='first_treat', aggregate='event_study') @@ -291,7 +295,7 @@ est.fit( # ---- multi-horizon ---- L_max: int | None = None, # Max horizon; None = l=1 only # ---- covariates and extensions ---- - aggregate: str | None = None, # Reserved; raises NotImplementedError + aggregate: Any = NOT_SUPPLIED, # DEPRECATED (M-026): warns; non-None raises - use results.aggregate() post-fit covariates: list[str] | None = None, # DID^X residualization-style covariates trends_linear: bool | None = None, # DID^{fd} group-specific linear trends trends_nonparam: Any | None = None, # DID^s state-set-specific trends @@ -301,7 +305,7 @@ est.fit( ) -> ChaisemartinDHaultfoeuilleResults ``` -`L_max` controls multi-horizon computation. `covariates`, `trends_linear`, `trends_nonparam`, `honest_did`, `heterogeneity`, `design2`, and `survey_design` are all supported; only `aggregate` still raises `NotImplementedError`. The 3.8 names `group=` and `controls=` remain accepted as deprecated aliases (FutureWarning; removed in 4.0). +`L_max` controls multi-horizon computation. `covariates`, `trends_linear`, `trends_nonparam`, `honest_did`, `heterogeneity`, `design2`, and `survey_design` are all supported. Fit-time `aggregate=` is DEPRECATED (row M-026: it never computed anything here; supplying it warns and any non-None value raises) - aggregation is post-fit: `results.aggregate('event_study')` returns the unified `EventStudyResults` view (Phase-1 fits the 2-row l=1 view, `L_max>=1` the multi-horizon surface) and `results.aggregate('simple')` a one-row `AggregationResult` relaying the overall estimand (DID_M / DID_1 / delta) bit-exactly; both are pure views, so bootstrap fits are permitted. The dCDH container is deliberately NOT accepted by `compute_honest_did`/`compute_pretrends_power` (their l1 placebo semantics need the native dCDH branch - pass the results object itself). The 3.8 names `group=` and `controls=` remain accepted as deprecated aliases (FutureWarning; removed in 4.0). **Usage:** diff --git a/diff_diff/guides/llms-practitioner.txt b/diff_diff/guides/llms-practitioner.txt index e158f4684..2650ef568 100644 --- a/diff_diff/guides/llms-practitioner.txt +++ b/diff_diff/guides/llms-practitioner.txt @@ -344,8 +344,11 @@ Bounds on the treatment effect under violations of parallel trends. Works with MultiPeriodDiD, CallawaySantAnna, and ChaisemartinDHaultfoeuille (dCDH) results. -- For CS: requires `aggregate='event_study'` or `aggregate='all'` so that - event study effects are available. +- For CS: pass the post-fit container - `compute_honest_did( + results.aggregate('event_study'))` - no refit needed (the fit-time + `aggregate=` route is deprecated since 3.9). EXCEPTION: on a + BOOTSTRAPPED fit (`n_bootstrap > 0`) post-fit `aggregate()` raises; + use the retained fit-time `aggregate='event_study'` there. - For dCDH: requires `L_max >= 1` (multi-horizon mode). Bounds use placebo estimates `DID^{pl}_l` as pre-period coefficients rather than standard event-study pre-treatment coefficients, and use diagonal variance (no diff --git a/diff_diff/guides/llms.txt b/diff_diff/guides/llms.txt index 42731e296..9fc1cfe1d 100644 --- a/diff_diff/guides/llms.txt +++ b/diff_diff/guides/llms.txt @@ -21,7 +21,7 @@ diagnostic steps produces unreliable results. 4. **Choose estimator** — staggered adoption → CS/SA/BJS (NOT plain TWFE); few treated units → SDiD; factor confounding → TROP; simple 2x2 → DiD. Run `BaconDecomposition` to diagnose TWFE bias. 5. **Estimate** — `estimator.fit(data, ...)`. Always print the cluster count first and choose inference method based on the result (cluster-robust if >= 50 clusters, wild bootstrap if fewer). 6. **Sensitivity analysis** — `compute_honest_did(results)` for bounds under PT violations (MultiPeriodDiD, CS, or dCDH), `run_all_placebo_tests()` for 2x2 falsification, specification comparisons for staggered designs. -7. **Heterogeneity** — CS: `results.aggregate('group')`/`.aggregate('event_study')` post-fit, no refit (fit-time `aggregate=`/`balance_e=` are deprecated since 3.9, removed in 4.0; they remain only for `compute_honest_did` / `compute_pretrends_power` / `plot_event_study`, which still read the fit-time surface); SA: `results.event_study_effects`/`to_dataframe(level='cohort')`; subgroup re-estimation. +7. **Heterogeneity** — CS: `results.aggregate('group')`/`.aggregate('event_study')` post-fit, no refit (fit-time `aggregate=`/`balance_e=` are deprecated since 3.9, removed in 4.0; `compute_honest_did` / `compute_pretrends_power` / `plot_event_study` all accept the post-fit `results.aggregate('event_study')` container directly; EXCEPTION: a BOOTSTRAPPED CS fit's `aggregate()` raises — use the fit-time aggregation there); dCDH: `results.aggregate('event_study')`/`.aggregate('simple')` post-fit views (bootstrap fits included — pure views); SA: `results.event_study_effects`/`to_dataframe(level='cohort')`; subgroup re-estimation. 8. **Robustness** — compare 2-3 estimators (CS vs SA vs BJS), MUST report with and without covariates (shows whether conditioning drives identification), present pre-trends and sensitivity bounds. After estimation, call `practitioner_next_steps(results)` for context-aware diff --git a/diff_diff/honest_did.py b/diff_diff/honest_did.py index e5aa61afb..0f5cae5a4 100644 --- a/diff_diff/honest_did.py +++ b/diff_diff/honest_did.py @@ -29,7 +29,7 @@ from diff_diff.results import ( MultiPeriodDiDResults, ) -from diff_diff.results_base import Diagnostic +from diff_diff.results_base import Diagnostic, _validate_vcov_subblock from diff_diff.utils import _get_critical_value # ============================================================================= @@ -209,7 +209,7 @@ class HonestDiDResults(Diagnostic): event_study_bounds: Optional[Dict[Any, Dict[str, float]]] = field(default=None, repr=False) # Survey design metadata (Phase 7d) survey_metadata: Optional[Any] = field(default=None, repr=False) - df_survey: Optional[int] = field(default=None, repr=False) + df_survey: Optional[float] = field(default=None, repr=False) def _ci_is_finite(self) -> bool: """Check if CI endpoints are finite (not NaN/inf).""" @@ -577,16 +577,258 @@ def plot( # ============================================================================= +def _extract_container_params( + surface: Any, +) -> Tuple[np.ndarray, np.ndarray, int, int, List[Any], List[Any], Optional[float]]: + """Container branch of ``_extract_event_study_params``. + + Consumes the unified ``EventStudyResults`` surface produced by + ``CallawaySantAnnaResults.aggregate('event_study')``. Admission is + SOURCE-SCOPED: containers from other producers are rejected rather + than silently admitted - widening is a per-estimator methodology + decision in each estimator's own ``aggregate()`` migration (row + M-093), not a side effect of the container existing. dCDH containers + in particular are rejected BY DESIGN: their ``l1_first_switch`` + placebo rows need the native dCDH branch's mandatory reinterpretation + warning and consecutive-horizon trimming. + """ + import warnings + + if surface.source != "CallawaySantAnnaResults": + raise TypeError( + "HonestDiD accepts EventStudyResults containers produced by " + "CallawaySantAnnaResults.aggregate('event_study') only " + f"(got source={surface.source!r}). For other estimators pass " + "the native results object where supported " + "(MultiPeriodDiDResults, CallawaySantAnnaResults, or " + "ChaisemartinDHaultfoeuilleResults); container admission for " + "further producers arrives with their own aggregate() " + "migrations." + ) + if surface.time_scale != "relative": + raise TypeError( + "HonestDiD requires a relative-time event-study container; " + f"got time_scale={surface.time_scale!r}." + ) + + # Common-reference guard: reference_event_times is the producer's + # cohort-level normalization-base provenance. More than one entry + # means the coefficients were normalized against DIFFERENT bases + # (CS universal on a gapped grid) - including the layout where one + # cohort's base OVERLAPS another cohort's estimated horizon, which + # no is_reference row marks. Mirrors the native CS branch. A + # universal container WITHOUT the field (hand-built - CS-produced + # containers always record it) cannot be verified: warn fail-safe + # rather than fail open silently. + _ref_e = surface.reference_event_times + if _ref_e is None and surface.base_period == "universal": + warnings.warn( + "This CallawaySantAnna event-study container carries no " + "reference_event_times provenance, so a common reference " + "period cannot be verified (universal base on a gapped time " + "grid may mix cohort-specific bases). CS-produced containers " + "record it - re-aggregate from the fitted results object.", + UserWarning, + stacklevel=4, + ) + if _ref_e is not None and len(set(_ref_e)) > 1: + raise ValueError( + "HonestDiD requires event-study coefficients normalized " + "against one common reference period, but this " + "CallawaySantAnna base_period='universal' fit selected " + "cohort-specific positional bases at event times " + f"{sorted(set(_ref_e))} (gapped time grid). On such grids a " + "cohort's base can overlap another cohort's estimated " + "horizon, so the coefficients are normalized against " + "different bases and are not jointly interpretable under " + "Rambachan-Roth's delta_0 = 0 normalization. Re-estimate on " + "a consecutive (ungapped) time grid so every cohort's base " + "falls at the same event time." + ) + + # Universal-base interpretation warning. Fail-safe: a container with + # base_period=None (unknown provenance, e.g. hand-built) cannot be + # verified, so it warns too. CS-produced containers always carry the + # real value. + if surface.base_period != "universal": + provenance_clause = ( + "This container carries no base_period provenance, so the " + "base-period regime cannot be verified. " + if surface.base_period is None + else "With base_period='varying', pre-treatment coefficients " + "use consecutive comparisons (not a common reference period), " + "which changes the meaning of the parallel trends restriction. " + ) + warnings.warn( + "HonestDiD sensitivity analysis on CallawaySantAnna results " + "requires base_period='universal' for valid interpretation. " + + provenance_clause + + "Re-run with CallawaySantAnna(base_period='universal') for " + "methodologically valid HonestDiD bounds.", + UserWarning, + stacklevel=4, + ) + + # Row filter: is_reference supersedes the fit-time n_groups==0 sniff; + # non-finite-SE rows (genuinely non-estimable horizons) drop exactly as + # the fit-time branch drops them. se == 0 rows drop too: safe_inference + # treats a zero SE as undefined inference (NaN t/p/CI), and admitting + # such a row here would launder undefined source inference into finite + # sensitivity bounds (mirrors the native CS branch and pretrends). + # Containers are publicly constructible: duplicate event-time labels + # would make every label-keyed subset below ambiguous. + _all_labels = surface.event_time.tolist() + if len(set(_all_labels)) != len(_all_labels): + raise ValueError( + "The event-study container carries duplicate event_time " + f"labels ({_all_labels}); each horizon must appear exactly " + "once." + ) + + keep = (~surface.is_reference) & np.isfinite(surface.se) & (surface.se > 0) + rel_times = [t for t, k in zip(surface.event_time.tolist(), keep) if k] + + ref_rows = surface.event_time[surface.is_reference].tolist() + if len(ref_rows) > 1: + # No route recommendation here: the native-results path picks the + # first n_groups==0 marker and then fails its own consecutive-grid + # validation on the same gapped layout, so pointing users there + # would be a dead end. + raise ValueError( + "HonestDiD cannot consume an event-study container with " + f"multiple reference rows ({sorted(ref_rows)}): its " + "consecutive-grid contract is defined around a single omitted " + "reference. Multiple references arise for CallawaySantAnna " + "base_period='universal' on a gapped time grid, where each " + "cohort's positional base is its own reference-only horizon. " + "Re-estimate on a consecutive (ungapped) time grid so a " + "single common reference is materialized." + ) + ref_period = ref_rows[0] if ref_rows else None + + if ref_period is not None: + pre_times = sorted(t for t in rel_times if t < ref_period) + post_times = sorted(t for t in rel_times if t > ref_period) + else: + # No reference row (varying base): the anticipation window + # [e = -k, -1] carries anticipated TREATMENT effects (REGISTRY + # anticipation contract; pretrends applies the same cutoff), so + # the clean pre-period set is e < -k and beta_post starts at -k. + # Splitting at 0 would misclassify anticipated effects as + # pre-trend coefficients. Mirrors the native CS branch. + _post_start = -int(surface.anticipation or 0) + pre_times = sorted(t for t in rel_times if t < _post_start) + post_times = sorted(t for t in rel_times if t >= _post_start) + + if len(pre_times) == 0: + raise ValueError( + "No pre-period effects with finite estimates found in the " + "event-study container. HonestDiD requires at least one " + "identified pre-period coefficient." + ) + + # Consecutive-grid validation, mirroring the fit-time branch: for a + # universal base pre[-1]+1 = ref and ref+1 = post[0] (gap of 2); for a + # varying base pre ends at -1 and post starts at 0 (gap of 1). + if pre_times and post_times: + ref_gap = post_times[0] - pre_times[-1] + has_gap = ref_gap != (2 if ref_period is not None else 1) + else: + has_gap = False + for block in [pre_times, post_times]: + if len(block) >= 2: + for i in range(len(block) - 1): + if block[i + 1] - block[i] != 1: + has_gap = True + break + if has_gap: + raise ValueError( + "HonestDiD requires a consecutive event-time grid " + "around the omitted reference period. Retained " + f"pre-periods {pre_times} and post-periods " + f"{post_times} have gaps. This can happen when " + "some event-study horizons have non-finite SEs. " + "Ensure all event-study periods have valid estimates, " + "or use balance_e to restrict to a balanced subset." + ) + + # beta_hat/sigma are subset in EXPLICIT [sorted pre; sorted post] + # order: the fit-side Rambachan-Roth split takes the first num_pre + # entries as beta_pre, and a hand-built container's rows need not + # arrive sorted - row-order subsetting would silently misalign the + # coefficient blocks and the covariance. + _row_of = {t: i for i, t in enumerate(surface.event_time.tolist())} + ordered_times = list(pre_times) + list(post_times) + idx_ordered = np.asarray([_row_of[t] for t in ordered_times], dtype=int) + beta_hat = np.asarray(surface.att[idx_ordered], dtype=float) + ses = np.asarray(surface.se[idx_ordered], dtype=float) + + if surface.vcov is not None and surface.vcov_index is not None: + vcov_labels = list(surface.vcov_index.tolist()) + if len(set(vcov_labels)) != len(vcov_labels): + raise ValueError( + "The event-study container's vcov_index carries duplicate " + f"labels ({vcov_labels}); the covariance sub-block is " + "ambiguous." + ) + missing = [t for t in ordered_times if t not in vcov_labels] + if missing: + # A SUPPLIED covariance whose index omits a retained horizon is + # inconsistent - fail loud rather than silently degrading to a + # diagonal approximation (the pretrends helper's convention; + # diagonal fallback is reserved for vcov is None). + raise ValueError( + f"The event-study container's vcov_index is missing " + f"retained horizon(s) {missing}; cannot extract the " + f"covariance sub-block. Available index: {vcov_labels}." + ) + idx = [vcov_labels.index(t) for t in ordered_times] + # allow_singular=False: Rambachan-Roth inference assumes covariance + # eigenvalues bounded away from zero (PreTrendsPower keeps its + # documented singular handling and passes the default True). + sigma = _validate_vcov_subblock( + surface.vcov[np.ix_(idx, idx)], ses, "HonestDiD", allow_singular=False + ) + else: + # Container-specific message: a vcov-less container has several + # documented causes (bootstrap/replicate overrides, producers that + # record no matrix) - do not attribute it to bootstrap. + warnings.warn( + "Event-study container carries no full covariance matrix; " + "using a diagonal approximation from the stored standard " + "errors. Cross-event-time covariance is unavailable on this " + "surface.", + UserWarning, + stacklevel=4, + ) + sigma = np.diag(ses**2) + + # Scalar df provenance threaded by the builder - exact fit-time parity, + # including the replicate-undefined 0.0 sentinel. + return ( + beta_hat, + sigma, + len(pre_times), + len(post_times), + pre_times, + post_times, + surface.df_survey, + ) + + def _extract_event_study_params( results: Union[MultiPeriodDiDResults, Any], -) -> Tuple[np.ndarray, np.ndarray, int, int, List[Any], List[Any], Optional[int]]: +) -> Tuple[np.ndarray, np.ndarray, int, int, List[Any], List[Any], Optional[float]]: """ Extract event study parameters from results objects. Parameters ---------- - results : MultiPeriodDiDResults, CallawaySantAnnaResults, or ChaisemartinDHaultfoeuilleResults - Estimation results with event study structure. + results : MultiPeriodDiDResults, CallawaySantAnnaResults, ChaisemartinDHaultfoeuilleResults, or EventStudyResults + Estimation results with event study structure, or the unified + event-study container produced by + ``CallawaySantAnnaResults.aggregate('event_study')`` (CS-sourced + containers only; see ``_extract_container_params``). Returns ------- @@ -602,24 +844,66 @@ def _extract_event_study_params( Pre-period identifiers. post_periods : list Post-period identifiers. - df_survey : int or None - Survey degrees of freedom for t-distribution inference. + df_survey : float or None + Survey degrees of freedom for t-distribution inference + (``0.0`` = replicate design with undefined df). """ + from diff_diff.results_base import EventStudyResults + + if isinstance(results, EventStudyResults): + return _extract_container_params(results) + if isinstance(results, MultiPeriodDiDResults): # Extract from MultiPeriodDiD pre_periods = results.pre_periods post_periods = results.post_periods - # Filter periods with finite effects/SEs, maintaining pre-then-post order + # Filter periods with finite effects/SEs, maintaining pre-then-post + # order. se <= 0 drops too: safe_inference treats it as undefined + # inference, and admitting such a row would launder NaN source + # inference into finite sensitivity bounds (mirrors the CS and + # container branches and pretrends). finite_periods = { p for p in results.period_effects.keys() if np.isfinite(results.period_effects[p].effect) and np.isfinite(results.period_effects[p].se) + and results.period_effects[p].se > 0 } pre_estimated = [p for p in pre_periods if p in finite_periods] post_estimated = [p for p in post_periods if p in finite_periods] + + # The RR constraint builders treat retained coefficients as + # CONSECUTIVE around the reference (they index by position, not + # label), so a dropped interior or reference-adjacent horizon + # would silently change the smoothness / first-difference + # geometry and return wrong bounds. Require the retained + # pre-periods to be a contiguous SUFFIX of the estimable pre grid + # (ending adjacent to the reference) and the retained + # post-periods a contiguous PREFIX (starting immediately after + # it). Leading-pre / trailing-post drops keep valid geometry; + # anything else fails closed (mirrors the CS branch's + # consecutive-grid validation). + _pre_grid = [p for p in pre_periods if p in results.period_effects] + _post_grid = [p for p in post_periods if p in results.period_effects] + _pre_ok = pre_estimated == _pre_grid[len(_pre_grid) - len(pre_estimated) :] + _post_ok = post_estimated == _post_grid[: len(post_estimated)] + if not (_pre_ok and _post_ok): + _dropped_pre = [p for p in _pre_grid if p not in finite_periods] + _dropped_post = [p for p in _post_grid if p not in finite_periods] + raise ValueError( + "HonestDiD requires consecutive estimated horizons around " + "the reference period: retained pre-periods must end " + "immediately before it and retained post-periods must " + "start immediately after it, with no interior gaps (the " + "Rambachan-Roth restrictions are built positionally). " + "Horizons with undefined inference (non-finite or zero " + f"SE) break that grid here: dropped pre {_dropped_pre}, " + f"dropped post {_dropped_post}. Only leading pre-periods " + "and trailing post-periods can be dropped safely." + ) + all_estimated = pre_estimated + post_estimated if not all_estimated: @@ -674,13 +958,19 @@ def _extract_event_study_params( if df_inference is not None: df_survey = int(df_inference) + # Return the ESTIMATED label lists - the ones beta_hat/sigma were + # built from - not the declared results.pre_periods/post_periods + # (which include the reference period and any dropped + # zero/non-finite-SE horizons). HonestDiDResults.pre_periods_used / + # post_periods_used relay these verbatim, so misaligned labels + # would claim excluded horizons entered the estimand. return ( beta_hat, sigma, num_pre_periods, num_post_periods, - pre_periods, - post_periods, + pre_estimated, + post_estimated, df_survey, ) @@ -692,10 +982,57 @@ def _extract_event_study_params( if isinstance(results, CallawaySantAnnaResults): if results.event_study_effects is None: raise ValueError( - "CallawaySantAnnaResults must have event_study_effects for HonestDiD. " - "Fit with aggregate='event_study'. HonestDiD reads the fit-time " - "surface and does not yet accept the post-fit container returned " - "by results.aggregate('event_study')." + "CallawaySantAnnaResults must have event_study_effects for " + "HonestDiD. Either pass the post-fit container " + "(compute_honest_did(results.aggregate('event_study'))) or " + "fit with the deprecated aggregate='event_study' to " + "populate the fit-time surface." + ) + + # Common-reference guard, mirroring the container branch: + # multiple distinct cohort base event times = coefficients + # normalized against different bases (gapped universal + # grid, including the overlapped-base layout no + # reference-only row marks). Provenance-less universal + # results (pre-3.9 pickles, replace()-stripped copies) + # must not FAIL OPEN: derive the bases from the + # materialized reference cells; warn fail-safe when even + # those are absent. + _ref_e = getattr(results, "reference_event_times", None) + if _ref_e is None and getattr(results, "base_period", None) == "universal": + _gte = getattr(results, "group_time_effects", None) or {} + _derived = {t - g for (g, t), _d in _gte.items() if _d.get("is_reference")} + if _derived: + _ref_e = tuple(sorted(_derived)) + else: + import warnings + + warnings.warn( + "This CallawaySantAnna base_period='universal' " + "result carries no reference_event_times " + "provenance and no materialized reference " + "cells, so a common reference period cannot " + "be verified (a gapped time grid may mix " + "cohort-specific bases). Re-fit with the " + "current version to record provenance.", + UserWarning, + stacklevel=4, + ) + if _ref_e is not None and len(set(_ref_e)) > 1: + raise ValueError( + "HonestDiD requires event-study coefficients " + "normalized against one common reference period, " + "but this CallawaySantAnna base_period='universal' " + "fit selected cohort-specific positional bases at " + f"event times {sorted(set(_ref_e))} (gapped time " + "grid). On such grids a cohort's base can overlap " + "another cohort's estimated horizon, so the " + "coefficients are normalized against different " + "bases and are not jointly interpretable under " + "Rambachan-Roth's delta_0 = 0 normalization. " + "Re-estimate on a consecutive (ungapped) time grid " + "so every cohort's base falls at the same event " + "time." ) # Warn if not using universal base period (R's HonestDiD requires it) @@ -715,11 +1052,16 @@ def _extract_event_study_params( ) # Extract event study effects by relative time - # Filter out normalization constraints (n_groups=0) and non-finite SEs + # Filter out normalization constraints (n_groups=0), non-finite + # SEs, and zero SEs (safe_inference treats se <= 0 as undefined + # inference - admitting such a row would launder NaN source + # inference into finite sensitivity bounds; mirrors pretrends). event_effects = { t: data for t, data in results.event_study_effects.items() - if data.get("n_groups", 1) > 0 and np.isfinite(data.get("se", np.nan)) + if data.get("n_groups", 1) > 0 + and np.isfinite(data.get("se", np.nan)) + and float(data.get("se", 0.0)) > 0 } rel_times = sorted(event_effects.keys()) @@ -742,9 +1084,16 @@ def _extract_event_study_params( pre_times = [t for t in rel_times if t < ref_period] post_times = [t for t in rel_times if t > ref_period] else: - # Varying base or no reference marker: split at t < 0 / t >= 0 - pre_times = [t for t in rel_times if t < 0] - post_times = [t for t in rel_times if t >= 0] + # Varying base or no reference marker: the anticipation + # window [e = -k, -1] carries anticipated TREATMENT + # effects (REGISTRY anticipation contract; pretrends + # applies the same cutoff), so the clean pre-period set + # is e < -k and beta_post starts at -k. Splitting at 0 + # would misclassify anticipated effects as pre-trend + # coefficients. Mirrors the container branch. + _post_start = -int(getattr(results, "anticipation", 0) or 0) + pre_times = [t for t in rel_times if t < _post_start] + post_times = [t for t in rel_times if t >= _post_start] if len(pre_times) == 0: raise ValueError( @@ -901,16 +1250,20 @@ def _extract_event_study_params( "for HonestDiD." ) - # Filter for finite SEs in both surfaces + # Filter for finite, strictly positive SEs in both surfaces + # (se <= 0 = undefined inference per safe_inference; a zero + # row would enter Sigma with zero variance and launder NaN + # source inference into finite bounds - mirrors the MPD, CS + # and container branches and pretrends). placebo_finite = { h: data for h, data in results.placebo_event_study.items() - if np.isfinite(data.get("se", np.nan)) + if np.isfinite(data.get("se", np.nan)) and float(data.get("se", 0.0)) > 0 } effects_finite = { h: data for h, data in results.event_study_effects.items() - if np.isfinite(data.get("se", np.nan)) + if np.isfinite(data.get("se", np.nan)) and float(data.get("se", 0.0)) > 0 } pre_times = sorted(placebo_finite.keys()) # -P, ..., -1 @@ -1039,7 +1392,8 @@ def _largest_consecutive_block(times, boundary_val): raise TypeError( f"Unsupported results type: {type(results)}. " "Expected MultiPeriodDiDResults, CallawaySantAnnaResults, " - "or ChaisemartinDHaultfoeuilleResults." + "ChaisemartinDHaultfoeuilleResults, or an EventStudyResults " + "container from CallawaySantAnnaResults.aggregate('event_study')." ) @@ -1440,7 +1794,7 @@ def _compute_flci( ub: float, se: float, alpha: float = 0.05, - df: Optional[int] = None, + df: Optional[float] = None, ) -> Tuple[float, float]: """ Compute Fixed Length Confidence Interval (FLCI). @@ -1485,7 +1839,7 @@ def _compute_flci( return ci_lb, ci_ub -def _cv_alpha(t: float, alpha: float, df: Optional[int] = None) -> float: +def _cv_alpha(t: float, alpha: float, df: Optional[float] = None) -> float: """ Compute the (1-alpha) quantile of the folded distribution |X|. @@ -1512,6 +1866,14 @@ def _cv_alpha(t: float, alpha: float, df: Optional[int] = None) -> float: target = 1 - alpha t = abs(t) + if df is not None and (not np.isfinite(df) or df <= 0): + # A PROVIDED nonpositive/non-finite df is the replicate-undefined + # sentinel (df_survey=0): inference is undefined and must fail + # closed to NaN - matching utils._get_critical_value's + # t.ppf(., 0) = NaN on the naive path - never silently fall + # through to normal theory (that path is reserved for df=None). + return float("nan") + if df is not None and df > 0: # Folded non-central t: P(|nct(df,t)| <= x) = F(x;df,t) - F(-x;df,t) from scipy.stats import nct as nct_dist @@ -1845,7 +2207,7 @@ def _flci_optimal_h( h0: float, M: float, alpha: float, - df: Optional[int], + df: Optional[float], levels=(40, 120, 120), ) -> float: """Argmin over the estimator SD ``h`` of the half-length @@ -1894,7 +2256,7 @@ def _flci_solve( num_post: int, M: float, alpha: float = 0.05, - df: Optional[int] = None, + df: Optional[float] = None, ) -> Tuple[float, float, Optional[np.ndarray]]: """ Compute the optimal Fixed Length Confidence Interval for Delta^SD. @@ -1996,7 +2358,7 @@ def _compute_optimal_flci( num_post: int, M: float, alpha: float = 0.05, - df: Optional[int] = None, + df: Optional[float] = None, ) -> Tuple[float, float]: """Optimal Delta^SD FLCI ``(ci_lb, ci_ub)`` (Rambachan & Roth 2023 §4.1). Thin wrapper over :func:`_flci_solve` (which also returns the optimal affine @@ -2458,8 +2820,16 @@ def fit( Parameters ---------- - results : MultiPeriodDiDResults, CallawaySantAnnaResults, or ChaisemartinDHaultfoeuilleResults - Results from event study estimation. + results : MultiPeriodDiDResults, CallawaySantAnnaResults, ChaisemartinDHaultfoeuilleResults, or EventStudyResults + Results from event study estimation, or the unified event-study + container from + ``CallawaySantAnnaResults.aggregate('event_study')``. On the + container route the scalar inference df arrives via the + container's ``df_survey`` provenance field (bounds and CIs match + the native route exactly), while the stored + ``HonestDiDResults.survey_metadata`` is None - the container + carries no survey-metadata object. That field's only inferential + consumer is the df extraction, so no number diverges. M : float, optional Override the M parameter for this fit. @@ -2612,7 +2982,7 @@ def _compute_smoothness_bounds( num_pre: int, num_post: int, M: float, - df: Optional[int] = None, + df: Optional[float] = None, ) -> Tuple[float, float, float, float]: """Compute bounds under smoothness restriction (Delta^SD). @@ -2674,7 +3044,7 @@ def _compute_rm_bounds( Mbar: float, pre_periods: List, results: Any, - df: Optional[int] = None, + df: Optional[float] = None, ) -> Tuple[float, float, float, float]: """Compute bounds under relative magnitudes restriction (Delta^RM). @@ -2718,7 +3088,7 @@ def _compute_combined_bounds( M: float, pre_periods: List, results: Any, - df: Optional[int] = None, + df: Optional[float] = None, ) -> Tuple[float, float, float, float]: """Compute bounds under combined smoothness + RM restriction.""" import warnings @@ -2830,7 +3200,7 @@ def sensitivity_analysis( Parameters ---------- - results : MultiPeriodDiDResults, CallawaySantAnnaResults, or ChaisemartinDHaultfoeuilleResults + results : MultiPeriodDiDResults, CallawaySantAnnaResults, ChaisemartinDHaultfoeuilleResults, or EventStudyResults Results from event study estimation. M_grid : list of float, optional Grid of M values to evaluate. If None, uses default grid @@ -2932,7 +3302,7 @@ def breakdown_value( Parameters ---------- - results : MultiPeriodDiDResults, CallawaySantAnnaResults, or ChaisemartinDHaultfoeuilleResults + results : MultiPeriodDiDResults, CallawaySantAnnaResults, ChaisemartinDHaultfoeuilleResults, or EventStudyResults Results from event study estimation. tol : float Tolerance for binary search. @@ -2989,7 +3359,7 @@ def compute_honest_did( Parameters ---------- - results : MultiPeriodDiDResults, CallawaySantAnnaResults, or ChaisemartinDHaultfoeuilleResults + results : MultiPeriodDiDResults, CallawaySantAnnaResults, ChaisemartinDHaultfoeuilleResults, or EventStudyResults Results from event study estimation. method : str Type of restriction ("smoothness", "relative_magnitude", "combined"). @@ -3031,7 +3401,7 @@ def sensitivity_plot( Parameters ---------- - results : MultiPeriodDiDResults, CallawaySantAnnaResults, or ChaisemartinDHaultfoeuilleResults + results : MultiPeriodDiDResults, CallawaySantAnnaResults, ChaisemartinDHaultfoeuilleResults, or EventStudyResults Results from event study estimation. method : str Type of restriction. diff --git a/diff_diff/practitioner.py b/diff_diff/practitioner.py index 9ca7a7d6a..b0cf41a49 100644 --- a/diff_diff/practitioner.py +++ b/diff_diff/practitioner.py @@ -378,23 +378,57 @@ def _handle_multi_period(results: Any): def _handle_cs(results: Any): + # Post-fit aggregate() RAISES on a bootstrapped fit (percentile + # statistics are not retained for re-aggregation), so the guidance + # must route those fits through the retained fit-time aggregation + # instead of advice that cannot run. + is_bootstrap = getattr(results, "bootstrap_results", None) is not None + if is_bootstrap: + sensitivity_why = ( + "Bounds the treatment effect under plausible violations of " + "parallel trends. This fit is BOOTSTRAPPED, and post-fit " + "aggregate() raises on bootstrap fits - refit with the " + "fit-time aggregation to populate the event-study surface." + ) + sensitivity_code = ( + "from diff_diff import compute_honest_did\n" + "# Bootstrap fit: post-fit aggregate() raises - use the\n" + "# fit-time aggregation:\n" + "results = cs.fit(data, ..., aggregate='event_study')\n" + "honest = compute_honest_did(results, method='relative_magnitude', M=1.0)\n" + "print(honest.summary())" + ) + heterogeneity_code = ( + "# Bootstrap fit: aggregate at fit time:\n" + "results = cs.fit(data, ..., aggregate='all')\n" + "print(results.group_effects) # Per-cohort ATTs\n" + "print(results.event_study_effects) # Dynamic effects" + ) + else: + sensitivity_why = ( + "Bounds the treatment effect under plausible violations of " + "parallel trends. Aggregate the event study post-fit — no " + "refit needed." + ) + sensitivity_code = ( + "from diff_diff import compute_honest_did\n" + "# Aggregate post-fit; the container feeds HonestDiD directly:\n" + "es = results.aggregate('event_study')\n" + "honest = compute_honest_did(es, method='relative_magnitude', M=1.0)\n" + "print(honest.summary())" + ) + heterogeneity_code = ( + "# Aggregate post-fit - no refit needed:\n" + "print(results.aggregate('group').to_dataframe()) # Per-cohort ATTs\n" + "print(results.aggregate('event_study').to_dataframe()) # Dynamic effects" + ) steps = [ _parallel_trends_step(staggered=True), _step( baker_step=6, label="Run HonestDiD sensitivity analysis", - why=( - "Bounds the treatment effect under plausible violations of " - "parallel trends. Requires event study effects — refit with " - "aggregate='event_study' or 'all' if not already done." - ), - code=( - "from diff_diff import compute_honest_did\n" - "# CS results must have event_study_effects:\n" - "results = cs.fit(data, ..., aggregate='event_study')\n" - "honest = compute_honest_did(results, method='relative_magnitude', M=1.0)\n" - "print(honest.summary())" - ), + why=sensitivity_why, + code=sensitivity_code, step_name="sensitivity", ), _step( @@ -405,12 +439,7 @@ def _handle_cs(results: Any): "dynamic effects over time. Inspect group and event study " "aggregations." ), - code=( - "# Re-fit with aggregate='all' to get all aggregations:\n" - "results = cs.fit(data, ..., aggregate='all')\n" - "print(results.group_effects) # Per-cohort ATTs\n" - "print(results.event_study_effects) # Dynamic effects" - ), + code=heterogeneity_code, step_name="heterogeneity", ), _robustness_compare_step("SA, BJS, or Gardner"), diff --git a/diff_diff/pretrends.py b/diff_diff/pretrends.py index 6e06caabb..8f39d0af0 100644 --- a/diff_diff/pretrends.py +++ b/diff_diff/pretrends.py @@ -35,7 +35,7 @@ from diff_diff._base import BaseEstimator from diff_diff.results import MultiPeriodDiDResults -from diff_diff.results_base import Diagnostic +from diff_diff.results_base import Diagnostic, _validate_vcov_subblock def _compute_nis_acceptance_prob( @@ -239,6 +239,44 @@ def _extract_event_study_vcov_subblock( return np.asarray(es_vcov)[np.ix_(indices, indices)], "full_pre_period_vcov" +def _extract_container_vcov_subblock( + surface: Any, + pre_periods: List[int], + ses: np.ndarray, +) -> Tuple[np.ndarray, str]: + """Container sibling of ``_extract_event_study_vcov_subblock``. + + Reads the unified ``EventStudyResults`` container's ``vcov`` / + ``vcov_index`` fields instead of the fit-time ``event_study_vcov`` + channel. Deliberately a SEPARATE helper: the fit-time one is shared by + both the CS and SA branches, and generalizing it would put those paths + in the blast radius of a container-only change. + """ + if surface.vcov is None or surface.vcov_index is None: + return np.diag(ses**2), "diag_fallback" + + index_labels = list(surface.vcov_index.tolist()) + if len(set(index_labels)) != len(index_labels): + raise ValueError( + "The event-study container's vcov_index carries duplicate " + f"labels ({index_labels}); the covariance sub-block is " + "ambiguous." + ) + try: + indices = [index_labels.index(t) for t in pre_periods] + except ValueError as e: + raise ValueError( + f"The event-study container's vcov_index is missing one of the " + f"pre-period labels {pre_periods}; cannot extract sub-block. " + f"Available index: {index_labels}. Original error: {e}" + ) from e + + sub = _validate_vcov_subblock( + np.asarray(surface.vcov)[np.ix_(indices, indices)], ses, "PreTrendsPower" + ) + return sub, "full_pre_period_vcov" + + # ============================================================================= # Results Classes # ============================================================================= @@ -1006,8 +1044,11 @@ def _extract_pre_period_params( Parameters ---------- - results : MultiPeriodDiDResults or similar - Results object from event study estimation. + results : MultiPeriodDiDResults, CallawaySantAnnaResults, SunAbrahamResults, or EventStudyResults + Results object from event study estimation, or the unified + event-study container from + ``CallawaySantAnnaResults.aggregate('event_study')`` + (CS-sourced containers only). pre_periods : list of int, optional Explicit list of pre-treatment periods. If None, uses results.pre_periods. @@ -1043,6 +1084,11 @@ def _extract_pre_period_params( result type (which would diverge from the actual extraction path the moment the routing changes — see PR-B Step 3). """ + from diff_diff.results_base import EventStudyResults + + if isinstance(results, EventStudyResults): + return self._extract_container_pre_period_params(results, pre_periods) + if isinstance(results, MultiPeriodDiDResults): # Get pre-period information - use explicit pre_periods if provided if pre_periods is not None: @@ -1123,9 +1169,82 @@ def _extract_pre_period_params( if results.event_study_effects is None: raise ValueError( "CallawaySantAnnaResults must have event_study_effects. " - "Fit with aggregate='event_study'. PreTrendsPower reads the " - "fit-time surface and does not yet accept the post-fit container " - "returned by results.aggregate('event_study')." + "Either pass the post-fit container " + "(compute_pretrends_power(results.aggregate(" + "'event_study'), ...)) or fit with the deprecated " + "aggregate='event_study' to populate the fit-time " + "surface." + ) + + # Common-reference guard, mirroring the container branch + # and HonestDiD: multiple distinct cohort base event + # times = pre-period coefficients normalized against + # DIFFERENT bases (gapped universal grid, including the + # overlapped-base layout no reference-only row marks). + # Provenance-less universal results (pre-3.9 pickles) + # must not FAIL OPEN: derive from the materialized + # reference cells; warn fail-safe when absent. + _ref_e = getattr(results, "reference_event_times", None) + if _ref_e is None and getattr(results, "base_period", None) == "universal": + _gte = getattr(results, "group_time_effects", None) or {} + _derived = {t - g for (g, t), _d in _gte.items() if _d.get("is_reference")} + if _derived: + _ref_e = tuple(sorted(_derived)) + else: + warnings.warn( + "This CallawaySantAnna base_period='universal' " + "result carries no reference_event_times " + "provenance and no materialized reference " + "cells, so a common reference period cannot " + "be verified (a gapped time grid may mix " + "cohort-specific bases). Re-fit with the " + "current version to record provenance.", + UserWarning, + stacklevel=4, + ) + if _ref_e is not None and len(set(_ref_e)) > 1: + raise ValueError( + "PreTrendsPower requires pre-period coefficients " + "normalized against one common reference period, " + "but this CallawaySantAnna base_period='universal' " + "fit selected cohort-specific positional bases at " + f"event times {sorted(set(_ref_e))} (gapped time " + "grid). The hypothesized violation delta is " + "defined relative to one reference, so power/MDV " + "over mixed-base coefficients target an " + "ill-defined alternative. Re-estimate on a " + "consecutive (ungapped) time grid so every " + "cohort's base falls at the same event time." + ) + + # Varying-base interpretation warning (twin of HonestDiD's + # universal-base warning): the built-in ``linear`` + # violation constructs delta as a slope on RELATIVE TIME, + # which assumes pre-treatment coefficients are LEVELS + # against one common reference. With base_period="varying" + # (the CS default), pre-treatment effects are + # consecutive-period comparisons - under a linear trend + # they are constant increments, not values proportional to + # |event_time| - so linear power/MDV target a different + # violation shape. Scoped to the LINEAR benchmark: the + # constant/last_period/custom vectors are user-specified + # in coefficient space. + if ( + getattr(results, "base_period", None) != "universal" + and getattr(self, "violation_type", "linear") == "linear" + ): + warnings.warn( + "PreTrendsPower on CallawaySantAnna results with " + "base_period='varying': pre-treatment coefficients " + "are consecutive-period comparisons, but the " + "'linear' violation benchmark assumes level " + "coefficients against a common reference period, " + "so linear power/MDV target a different violation " + "shape. Re-run with CallawaySantAnna(" + "base_period='universal') for Roth-faithful linear " + "benchmarks.", + UserWarning, + stacklevel=4, ) # Get pre-period effects. Anticipation-aware cutoff per @@ -1157,6 +1276,25 @@ def _extract_pre_period_params( if not pre_effects: raise ValueError("No pre-treatment periods found in event study.") + # An explicitly requested pre-period subset is honored, + # never silently ignored (mirrors the MPD branch and the + # container branch). + if pre_periods is not None: + _missing = [t for t in sorted(pre_periods) if t not in pre_effects] + if _missing: + raise ValueError( + f"Requested pre_periods {_missing} are not " + f"eligible pre-treatment periods in this event " + f"study (eligible: {sorted(pre_effects)})." + ) + pre_effects = {t: pre_effects[t] for t in sorted(pre_periods)} + if not pre_effects: + raise ValueError( + "Pre-trends power analysis requires at least " + "one pre-period; the requested pre_periods " + "subset is empty." + ) + pre_periods = sorted(pre_effects.keys()) n_pre = len(pre_periods) @@ -1171,7 +1309,23 @@ def _extract_pre_period_params( results, pre_periods, ses ) - relative_times = np.asarray(pre_periods, dtype=float) + # Roth's linear violation is defined relative to the + # NORMALIZED ZERO (his omitted period is labeled t=0): + # delta must vanish at the omitted reference, so relative + # times are measured FROM it - |t - t_ref| - not from + # treatment (the MPD branch's + # _coerce_relative_times_from_reference precedent). + # Universal fits carry the anchor as the singleton + # reference_event_times; varying fits have no common + # reference (raw labels retained under the varying-base + # warning; transformation tracked in TODO.md). + if _ref_e is not None and len(set(_ref_e)) == 1: + _anchor = float(next(iter(set(_ref_e)))) + relative_times = np.asarray( + [float(t) - _anchor for t in pre_periods], dtype=float + ) + else: + relative_times = np.asarray(pre_periods, dtype=float) return effects, ses, vcov, n_pre, relative_times, covariance_source except ImportError: pass @@ -1202,6 +1356,24 @@ def _extract_pre_period_params( if not pre_effects: raise ValueError("No pre-treatment periods found in event study.") + # Honor an explicitly requested pre-period subset (mirrors + # the CS branch above). + if pre_periods is not None: + _missing = [t for t in sorted(pre_periods) if t not in pre_effects] + if _missing: + raise ValueError( + f"Requested pre_periods {_missing} are not " + f"eligible pre-treatment periods in this event " + f"study (eligible: {sorted(pre_effects)})." + ) + pre_effects = {t: pre_effects[t] for t in sorted(pre_periods)} + if not pre_effects: + raise ValueError( + "Pre-trends power analysis requires at least " + "one pre-period; the requested pre_periods " + "subset is empty." + ) + pre_periods = sorted(pre_effects.keys()) n_pre = len(pre_periods) @@ -1217,16 +1389,204 @@ def _extract_pre_period_params( results, pre_periods, ses ) - relative_times = np.asarray(pre_periods, dtype=float) + # SunAbraham omits e = -1 - anticipation as its COMMON + # reference (saturated-regression omitted category): + # anchor Roth's linear violation there (see the CS + # branch for the normalization rationale). + _sa_ref = getattr(results, "reference_period", None) + if _sa_ref is None: + _sa_ref = -1 - int(getattr(results, "anticipation", 0) or 0) + relative_times = np.asarray( + [float(t) - float(_sa_ref) for t in pre_periods], dtype=float + ) return effects, ses, vcov, n_pre, relative_times, covariance_source except ImportError: pass raise TypeError( f"Unsupported results type: {type(results)}. " - "Expected MultiPeriodDiDResults, CallawaySantAnnaResults, or SunAbrahamResults." + "Expected MultiPeriodDiDResults, CallawaySantAnnaResults, " + "SunAbrahamResults, or an EventStudyResults container from " + "CallawaySantAnnaResults.aggregate('event_study')." ) + def _extract_container_pre_period_params( + self, + surface: Any, + pre_periods: Optional[List[int]] = None, + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, int, Optional[np.ndarray], str]: + """Container branch of ``_extract_pre_period_params``. + + Consumes the unified ``EventStudyResults`` surface produced by + ``CallawaySantAnnaResults.aggregate('event_study')``. Admission is + SOURCE-SCOPED: containers from other producers are rejected rather + than silently admitted - widening is a per-estimator methodology + decision in each estimator's own ``aggregate()`` migration (row + M-093). dCDH containers in particular are rejected BY DESIGN: + their ``l1_first_switch`` placebo rows are a different estimand + from pre-trend coefficients, and PreTrendsPower has no native dCDH + branch either. + """ + if surface.source != "CallawaySantAnnaResults": + raise TypeError( + "PreTrendsPower accepts EventStudyResults containers " + "produced by CallawaySantAnnaResults.aggregate(" + "'event_study') only " + f"(got source={surface.source!r}). For other estimators " + "pass the native results object where supported " + "(MultiPeriodDiDResults, CallawaySantAnnaResults, or " + "SunAbrahamResults); container admission for further " + "producers arrives with their own aggregate() migrations." + ) + if surface.time_scale != "relative": + raise TypeError( + "PreTrendsPower requires a relative-time event-study " + f"container; got time_scale={surface.time_scale!r}." + ) + + # Containers are publicly constructible: duplicate event-time + # labels would make every label-keyed subset below ambiguous. + _all_labels = surface.event_time.tolist() + if len(set(_all_labels)) != len(_all_labels): + raise ValueError( + "The event-study container carries duplicate event_time " + f"labels ({_all_labels}); each horizon must appear " + "exactly once." + ) + + # Common-reference guard: cohort-level normalization-base + # provenance (mirrors the fit-time CS branch and HonestDiD). A + # universal container WITHOUT the field (hand-built - CS-produced + # containers always record it) cannot be verified: warn + # fail-safe rather than fail open silently. + _ref_e = surface.reference_event_times + if _ref_e is None and surface.base_period == "universal": + warnings.warn( + "This CallawaySantAnna event-study container carries no " + "reference_event_times provenance, so a common reference " + "period cannot be verified (universal base on a gapped " + "time grid may mix cohort-specific bases). CS-produced " + "containers record it - re-aggregate from the fitted " + "results object.", + UserWarning, + stacklevel=4, + ) + if _ref_e is not None and len(set(_ref_e)) > 1: + raise ValueError( + "PreTrendsPower requires pre-period coefficients " + "normalized against one common reference period, but " + "this CallawaySantAnna base_period='universal' fit " + "selected cohort-specific positional bases at event " + f"times {sorted(set(_ref_e))} (gapped time grid). The " + "hypothesized violation delta is defined relative to " + "one reference, so power/MDV over mixed-base " + "coefficients target an ill-defined alternative. " + "Re-estimate on a consecutive (ungapped) time grid so " + "every cohort's base falls at the same event time." + ) + + # Varying-base interpretation warning, mirroring the fit-time CS + # branch (and HonestDiD's container fail-safe): base_period=None + # (unknown provenance, e.g. hand-built) cannot be verified, so it + # warns too. CS-produced containers always carry the real value. + # Scoped to the LINEAR benchmark (the built-in construction the + # warning is about); other violation vectors are user-specified + # in coefficient space. + if ( + surface.base_period != "universal" + and getattr(self, "violation_type", "linear") == "linear" + ): + provenance_clause = ( + "This container carries no base_period provenance, so the " + "base-period regime cannot be verified. " + if surface.base_period is None + else "With base_period='varying', pre-treatment " + "coefficients are consecutive-period comparisons, but the " + "'linear' violation benchmark assumes level coefficients " + "against a common reference period, so linear power/MDV " + "target a different violation shape. " + ) + warnings.warn( + "PreTrendsPower on a CallawaySantAnna event-study " + "container without base_period='universal'. " + + provenance_clause + + "Re-run with CallawaySantAnna(base_period='universal') " + "for Roth-faithful linear benchmarks.", + UserWarning, + stacklevel=4, + ) + + # Anticipation-aware cutoff, mirroring the fit-time CS branch: with + # anticipation=k, true pre-periods are t < -k. The container's + # anticipation provenance field carries the fit's value (None on a + # hand-built container -> 0). + _ant = int(surface.anticipation or 0) + _pre_cutoff = -_ant + keep = ( + (~surface.is_reference) + & np.isfinite(surface.se) + & (surface.se > 0) + & (surface.event_time < _pre_cutoff) + ) + + if not keep.any(): + raise ValueError("No pre-treatment periods found in event study.") + + # CHRONOLOGICAL order, not row order: positional violation + # patterns (last_period assigns weights[-1] to the final entry; + # user-supplied custom weights are positional too) silently + # mis-map on a hand-built container whose rows arrive unsorted. + keep_idx = np.where(keep)[0] + keep_idx = keep_idx[np.argsort(surface.event_time[keep_idx])] + pre_list = list(surface.event_time[keep_idx].tolist()) + effects = np.asarray(surface.att[keep_idx], dtype=float) + ses = np.asarray(surface.se[keep_idx], dtype=float) + + # An explicitly requested pre-period subset is honored, never + # silently ignored: every requested label must be an eligible + # pre-period on this surface, and effects/SEs/VCV subset to + # exactly that (sorted) set. + if pre_periods is not None: + requested = sorted(pre_periods) + missing = [t for t in requested if t not in pre_list] + if missing: + raise ValueError( + f"Requested pre_periods {missing} are not eligible " + f"pre-treatment periods on this event-study surface " + f"(eligible: {pre_list})." + ) + sel = [i for i, t in enumerate(pre_list) if t in set(requested)] + pre_list = [pre_list[i] for i in sel] + effects = effects[sel] + ses = ses[sel] + if not pre_list: + raise ValueError( + "Pre-trends power analysis requires at least one " + "pre-period; the requested pre_periods subset is empty." + ) + + vcov, covariance_source = _extract_container_vcov_subblock(surface, pre_list, ses) + + # Anchor Roth's linear violation at the omitted reference (see the + # native CS branch for the normalization rationale): the single + # marked reference row, else the singleton reference_event_times + # provenance (a universal fit whose base row was trimmed), else no + # common reference exists (varying) and raw labels are retained + # under the varying-base warning. + _ref_rows = surface.event_time[surface.is_reference].tolist() + _anchor: Optional[float] + if len(_ref_rows) == 1: + _anchor = float(_ref_rows[0]) + elif _ref_e is not None and len(set(_ref_e)) == 1: + _anchor = float(next(iter(set(_ref_e)))) + else: + _anchor = None + if _anchor is not None: + relative_times = np.asarray([float(t) - _anchor for t in pre_list], dtype=float) + else: + relative_times = np.asarray(pre_list, dtype=float) + return effects, ses, np.asarray(vcov), len(pre_list), relative_times, covariance_source + def _compute_power( self, M: float, @@ -1530,7 +1890,7 @@ def fit( Parameters ---------- - results : MultiPeriodDiDResults, CallawaySantAnnaResults, or SunAbrahamResults + results : MultiPeriodDiDResults, CallawaySantAnnaResults, SunAbrahamResults, or EventStudyResults Results from an event study estimation. M : float, optional Specific violation magnitude to evaluate. If None, evaluates at diff --git a/diff_diff/results_base.py b/diff_diff/results_base.py index 3935fb888..36895777f 100644 --- a/diff_diff/results_base.py +++ b/diff_diff/results_base.py @@ -240,6 +240,40 @@ class EventStudyResults(BaseResults): 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. + base_period : str or None + Producer provenance: the fit's base-period regime where the + producer has one (CallawaySantAnna vocabulary: ``"varying"`` or + ``"universal"``). None when the producer has no such notion. + HonestDiD reads this for its universal-base interpretation + warning. + anticipation : int or None + Producer provenance: the fit's anticipation window in periods, + where the producer has one. None when the producer has no such + notion. PreTrendsPower reads this to exclude anticipation-window + rows (``event_time >= -anticipation``) from the pre-trend set. + df_survey : float or None + Producer provenance: the fit's resolved SCALAR inference df, + with the established semantics of the fit-time consumers - + ``survey_metadata.df_survey`` where present (``0.0`` = replicate + design with an undefined df, which fails closed to NaN critical + values downstream), else ``df_inference`` (the bare-``cluster=`` + carrier), else None (no scalar df notion). Exists beside the + per-row ``df`` column because that column CANNOT encode the + replicate-undefined sentinel: ``__post_init__`` forces per-row + df to NaN wherever the p-value is non-finite. + reference_event_times : tuple or None + Producer provenance: the DISTINCT per-cohort normalization-base + event times (CallawaySantAnna ``base_period="universal"``: each + cohort's positional base period minus its cohort, deduplicated, + sorted). More than one entry means the coefficients were + normalized against DIFFERENT bases (gapped time grid) - and on + such grids a cohort's base can OVERLAP another cohort's + estimated horizon, where NO reference-only row exists to mark + it, so this field (not ``is_reference``) is the authoritative + common-reference signal. HonestDiD and PreTrendsPower fail + closed when it carries more than one entry. None when the + producer records no such notion (varying base, non-CS + producers, hand-built surfaces). """ event_time: np.ndarray @@ -263,6 +297,22 @@ class EventStudyResults(BaseResults): alpha: float = 0.05 source: Optional[str] = None df: Optional[Union[float, np.ndarray]] = None + # Provenance fields declared LAST so every pre-existing field keeps its + # positional index in the generated __init__ (the constructor signature + # is public API). + base_period: Optional[str] = None + anticipation: Optional[int] = None + df_survey: Optional[float] = None + # Distinct per-cohort normalization-base EVENT TIMES (CallawaySantAnna + # base_period="universal": each cohort's positional base minus its + # cohort, deduplicated and sorted). NOT the same thing as the + # ``reference_periods`` property (the is_reference-marked rows): on a + # gapped grid a cohort's base can OVERLAP another cohort's estimated + # horizon, where no reference-only row exists to mark it - this field + # is the consumer-facing signal that coefficients were normalized + # against more than one base. None when the producer records no such + # notion (varying base, non-CS producers, hand-built surfaces). + reference_event_times: Optional[Tuple[Any, ...]] = None _ARRAY_FIELDS = ( "att", @@ -482,6 +532,16 @@ def to_dict(self) -> Dict[str, Any]: "alpha": self.alpha, "source": self.source, "df": cast(np.ndarray, self.df).tolist(), + "base_period": self.base_period, + "anticipation": self.anticipation, + "df_survey": self.df_survey, + "reference_event_times": ( + # _json_safe_label per element: CS period arithmetic yields + # numpy scalars, which json.dumps cannot serialize. + [_json_safe_label(v) for v in self.reference_event_times] + if self.reference_event_times is not None + else None + ), } return out @@ -563,6 +623,115 @@ def _absent(results: Any) -> ValueError: return ValueError(f"{name} carries no event-study surface - {hint}.") +def _validate_vcov_subblock( + sigma: np.ndarray, + ses: np.ndarray, + consumer: str, + *, + allow_singular: bool = True, +) -> np.ndarray: + """Integrity checks for a consumer-bound covariance sub-block. + + Containers are publicly constructible, so the consumer boundary + validates what the producers guarantee by construction: finite + entries, symmetry, a diagonal equal to the stored ``se**2`` (the + container contract clears ``vcov`` rather than ship a matrix whose + diagonal disagrees with the stored SEs), and no material + indefiniteness. Every tolerance is RELATIVE to the matrix scale - an + absolute floor would wave through a materially indefinite low-scale + matrix (e.g. diagonal 1e-10 with -1e-10 eigenvalues). + + ``allow_singular=False`` additionally rejects singular/near-singular + sub-blocks: Rambachan-Roth's inference assumes eigenvalues bounded + away from zero, so HonestDiD passes False; PreTrendsPower keeps its + documented singular-covariance handling with the default True. + """ + sigma = np.asarray(sigma, dtype=float) + if not np.all(np.isfinite(sigma)): + raise ValueError( + f"{consumer}: the event-study container's covariance " + "sub-block contains non-finite entries." + ) + _scale = float(np.max(np.abs(sigma))) if sigma.size else 0.0 + if not np.allclose(sigma, sigma.T, rtol=1e-8, atol=1e-12 * _scale): + raise ValueError( + f"{consumer}: the event-study container's covariance " "sub-block is not symmetric." + ) + ses_arr = np.asarray(ses, dtype=float) + # Pure relative: retained rows have se > 0, so the diagonal target is + # strictly positive and an absolute atol would mask low-scale + # inconsistencies. + if not np.allclose(np.diag(sigma), ses_arr**2, rtol=1e-6, atol=0.0): + raise ValueError( + f"{consumer}: the event-study container's covariance diagonal " + "is inconsistent with the stored standard errors (the " + "container contract clears vcov rather than ship a matrix " + "whose diagonal disagrees with se**2)." + ) + if sigma.size: + eigs = np.linalg.eigvalsh((sigma + sigma.T) / 2.0) + _eig_scale = float(np.max(np.abs(eigs))) + if _eig_scale > 0.0 and float(eigs.min()) < -1e-8 * _eig_scale: + raise ValueError( + f"{consumer}: the event-study container's covariance " + "sub-block is indefinite (most negative eigenvalue " + f"{float(eigs.min()):.3e} at scale {_eig_scale:.3e}). " + "Positive semi-definiteness is required." + ) + if not allow_singular and _eig_scale > 0.0 and float(eigs.min()) < 1e-10 * _eig_scale: + raise ValueError( + f"{consumer}: the event-study container's covariance " + "sub-block is singular or near-singular (smallest " + f"eigenvalue {float(eigs.min()):.3e} at scale " + f"{_eig_scale:.3e}). Rambachan-Roth inference assumes " + "eigenvalues bounded away from zero; drop collinear " + "horizons or re-estimate before running HonestDiD." + ) + return sigma + + +def _resolve_scalar_df_survey(results: Any) -> Optional[float]: + """Resolve the producer's SCALAR inference df for container provenance. + + Mirrors the fit-time consumers' preference order (honest_did): + ``survey_metadata.df_survey`` where present, with a replicate design + whose df is undefined mapping to the ``0.0`` sentinel (fails closed to + NaN critical values downstream); else the bare-``cluster=`` + ``df_inference`` carrier; else None. Deliberate local sibling of + ``aggregation.resolve_inference_df``: this module cannot import it + (aggregation.py imports results_base - the dependency is one-way); + folding the copies together is tracked in TODO.md (df-resolution / + adapter-naming consolidation row). + """ + sm = getattr(results, "survey_metadata", None) + if sm is not None: + df_survey = getattr(sm, "df_survey", None) + if df_survey is not None: + return float(df_survey) + if getattr(sm, "replicate_method", None) is not None: + return 0.0 + df_inference = getattr(results, "df_inference", None) + if df_inference is not None: + return float(df_inference) + return None + + +def _provenance_kwargs(results: Any) -> Dict[str, Any]: + """Producer-provenance fields threaded onto the container. + + ``getattr`` with a None default: producers that declare no + ``base_period``/``anticipation`` notion yield None - values are never + invented. + """ + ref_e = getattr(results, "reference_event_times", None) + return { + "base_period": getattr(results, "base_period", None), + "anticipation": getattr(results, "anticipation", None), + "df_survey": _resolve_scalar_df_survey(results), + "reference_event_times": tuple(ref_e) if ref_e is not None else None, + } + + def _empty_surface(results: Any) -> EventStudyResults: """Zero-row surface for a requested-but-empty event study.""" empty_f = np.empty(0, dtype=float) @@ -579,6 +748,7 @@ def _empty_surface(results: Any) -> EventStudyResults: time_scale="relative", alpha=getattr(results, "alpha", 0.05), source=type(results).__name__, + **_provenance_kwargs(results), ) @@ -761,6 +931,7 @@ def _from_relative_dict(results: Any) -> EventStudyResults: alpha=getattr(results, "alpha", 0.05), source=type(results).__name__, df=df_arg, + **_provenance_kwargs(results), ) @@ -842,6 +1013,7 @@ def _from_mpd(results: Any) -> EventStudyResults: alpha=getattr(results, "alpha", 0.05), source=type(results).__name__, df=df_arg, + **_provenance_kwargs(results), ) @@ -897,6 +1069,7 @@ def _from_lpdid(results: Any) -> EventStudyResults: alpha=getattr(results, "alpha", 0.05), source=type(results).__name__, df=df_arg, + **_provenance_kwargs(results), ) @@ -985,6 +1158,7 @@ def _from_dcdh(results: Any) -> EventStudyResults: # so a scalar is faithful - the container broadcasts it and NaN-masks # the synthesized reference row at 0 plus any NaN-p rows. df=getattr(results, "event_study_df", None), + **_provenance_kwargs(results), ) @@ -1022,6 +1196,7 @@ def _from_had(results: Any) -> EventStudyResults: cband_crit_value=getattr(results, "cband_crit_value", None), alpha=getattr(results, "alpha", 0.05), source=type(results).__name__, + **_provenance_kwargs(results), ) diff --git a/diff_diff/staggered.py b/diff_diff/staggered.py index 8cc315284..0370a2c07 100644 --- a/diff_diff/staggered.py +++ b/diff_diff/staggered.py @@ -1860,10 +1860,15 @@ def fit( - "group": Aggregate by treatment cohort - "all": Compute all aggregations - Fit-time aggregation remains ONLY as the temporary compatibility - surface for consumers that still read it (``compute_honest_did``, - ``compute_pretrends_power`` and ``plot_event_study``, tracked in - ``TODO.md``); it is not the recommended path for new code. + ``compute_honest_did``, ``compute_pretrends_power`` and + ``plot_event_study`` all accept the post-fit container from + ``results.aggregate('event_study')`` directly, so no consumer + requires the fit-time surface anymore - EXCEPT on bootstrapped + fits (``n_bootstrap > 0``), where post-fit ``aggregate()`` + raises (percentile inference cannot be reproduced from the + retained analytical state) and fit-time aggregation remains + the supported route. Otherwise it remains only as the + deprecated compatibility path through 3.9. balance_e : int, optional DEPRECATED since 3.9, removed in 4.0 (ledger row M-117). Passing it emits a ``FutureWarning``; it moves onto the post-fit call as @@ -2687,6 +2692,25 @@ def fit( "control_inf": np.array([]), } + # Common-reference provenance (universal base only): the DISTINCT + # per-cohort base EVENT TIMES (base - g). On gapped grids the + # positional bases land at different event times - and a cohort's + # base can OVERLAP another cohort's estimated horizon, where NO + # reference-only row marks it - so consumers that require one + # common reference (HonestDiD / PreTrendsPower) read this field, + # not the marker rows. Computed for every cohort with a valid + # base, independent of the materialization conditions above. + # Varying base has no constant per-cohort reference -> None. + reference_event_times: Optional[Tuple[Any, ...]] = None + if self.base_period == "universal": + _ref_event_times = set() + for _g in treatment_groups: + _b = self._select_base_period(_g, _g, precomputed["observed_sorted"]) + if _b is not None: + _ref_event_times.add(_b - _g) + if _ref_event_times: + reference_event_times = tuple(sorted(_ref_event_times)) + # Compute overall ATT (simple aggregation) overall_att, overall_se, overall_effective_df = self._aggregate_simple( group_time_effects, influence_func_info, df, unit, precomputed @@ -2732,6 +2756,14 @@ def fit( precomputed, ) event_study_effects = es_aggregation.effects + # The stored fit-time surface is THIS aggregation's (possibly + # balance_e-restricted): its common-reference provenance must + # describe the RETAINED cohorts, not the fit-wide set, or the + # native consumer route would disagree with the equivalent + # post-fit container (route parity). The fit-wide tuple + # computed above stands only when no event-study surface was + # produced. + reference_event_times = es_aggregation.reference_event_times if aggregate in ["group", "all"]: group_effects = self._aggregate_by_group( @@ -2981,6 +3013,7 @@ def fit( cluster_name=cluster_name_for_results, n_clusters=n_clusters_for_results, df_inference=df_inference_for_results, + reference_event_times=reference_event_times, ) # Attach the post-fit aggregation kit (spec section 6, rows M-020/M-117). diff --git a/diff_diff/staggered_aggregation.py b/diff_diff/staggered_aggregation.py index e37fa4baf..68d2d1a7a 100644 --- a/diff_diff/staggered_aggregation.py +++ b/diff_diff/staggered_aggregation.py @@ -50,6 +50,14 @@ class EventStudyAggregation: df_used: Optional[float] = None vcov: Optional[np.ndarray] = None vcov_index: Optional[List[Any]] = None + #: Distinct base EVENT TIMES of the cohorts RETAINED by this + #: aggregation (derived from the materialized universal-base + #: is_reference cells; None when there are none, e.g. varying base). + #: Surface-faithful, unlike the fit-level fit-wide tuple: balance_e + #: can drop the cohort responsible for a second base, and the + #: container's common-reference guard must reflect the cohorts that + #: actually entered the reported estimand. + reference_event_times: Optional[Tuple[Any, ...]] = None def fixed_cohort_agg_weights( @@ -956,6 +964,23 @@ def _aggregate_event_study( ) effects_by_e = balanced_effects + # Common-reference provenance for THIS aggregation's surface: the + # distinct base event times of the RETAINED cohorts, read off the + # materialized universal-base is_reference cells (varying-base + # fits have none -> None). Surface-faithful by construction: + # balance_e can drop the cohort responsible for a second base, in + # which case the fit-level fit-wide tuple would over-restrict the + # balanced container. + _retained_cohorts = {g for cells in effects_by_e.values() for (g, _t), _eff, _w in cells} + _ref_es = { + t - g + for (g, t), data in group_time_effects.items() + if data.get("is_reference") and g in _retained_cohorts + } + es_reference_event_times: Optional[Tuple[Any, ...]] = ( + tuple(sorted(_ref_es)) if _ref_es else None + ) + # Universal base period: each cohort's positional base is materialized in # `group_time_effects` / `influence_func_info` (with a zero effect and a # zero influence function) by `fit()` before aggregation, so it is already @@ -1202,6 +1227,7 @@ def _aggregate_event_study( df_used=es_df_used, vcov=event_study_vcov, vcov_index=event_study_vcov_index, + reference_event_times=es_reference_event_times, ) def _aggregate_by_group( diff --git a/diff_diff/staggered_results.py b/diff_diff/staggered_results.py index 66a523fad..3a8e96d54 100644 --- a/diff_diff/staggered_results.py +++ b/diff_diff/staggered_results.py @@ -239,6 +239,15 @@ class CallawaySantAnnaResults(BaseResults, AggregationMixin): # here rather than set dynamically so it is a typed part of the contract. # Excluded from repr and equality: it is internal bookkeeping, not a # reportable result, and its arrays would make `==` raise. + # Distinct per-cohort normalization-base EVENT TIMES under + # base_period="universal" (each cohort's positional base minus its + # cohort, deduplicated, sorted): the common-reference provenance + # HonestDiD / PreTrendsPower need on gapped grids, where a cohort's + # base can overlap another cohort's estimated horizon and no + # reference-only event-study row marks it. None on varying-base fits + # (no constant per-cohort reference exists). Declared last among the + # public fields (positional-compat convention). + reference_event_times: Optional[Tuple[Any, ...]] = None _aggregation_kit: Optional[Any] = field(default=None, repr=False, compare=False) # --- Inference-field aliases (balance/external-adapter compatibility) --- @@ -352,6 +361,11 @@ def _aggregate_compute( event_study_vcov=es.vcov, event_study_vcov_index=es.vcov_index, event_study_df=es.df_used, + # Surface-faithful common-reference provenance: the aggregation + # recomputes it over the RETAINED cohorts (balance_e can drop + # the cohort responsible for a second base; the fit-level + # fit-wide tuple would over-restrict the balanced container). + reference_event_times=es.reference_event_times, ) return build_event_study_surface(carrier) diff --git a/diff_diff/utils.py b/diff_diff/utils.py index 2fa236b37..cf729f89a 100644 --- a/diff_diff/utils.py +++ b/diff_diff/utils.py @@ -33,10 +33,10 @@ _NUMERICAL_EPS = 1e-10 # Small constant to prevent division by zero # Cache for critical values to avoid repeated scipy calls -_critical_value_cache: Dict[Tuple[float, Optional[int]], float] = {} +_critical_value_cache: Dict[Tuple[float, Optional[float]], float] = {} -def _get_critical_value(alpha: float, df: Optional[int] = None) -> float: +def _get_critical_value(alpha: float, df: Optional[float] = None) -> float: """Return cached critical value for (alpha, df) pair.""" key = (alpha, df) if key not in _critical_value_cache: diff --git a/diff_diff/visualization/_event_study.py b/diff_diff/visualization/_event_study.py index 09a7fb00d..9beec3fdf 100644 --- a/diff_diff/visualization/_event_study.py +++ b/diff_diff/visualization/_event_study.py @@ -1,6 +1,7 @@ """Event study visualization functions.""" -from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple, Union +import warnings +from typing import TYPE_CHECKING, Any, Dict, List, Optional, Set, Tuple, Union import numpy as np import pandas as pd @@ -12,6 +13,7 @@ from diff_diff.honest_did import HonestDiDResults from diff_diff.imputation import ImputationDiDResults from diff_diff.results import MultiPeriodDiDResults + from diff_diff.results_base import EventStudyResults from diff_diff.stacked_did import StackedDiDResults from diff_diff.staggered import CallawaySantAnnaResults from diff_diff.sun_abraham import SunAbrahamResults @@ -26,6 +28,7 @@ "TwoStageDiDResults", "StackedDiDResults", "ChaisemartinDHaultfoeuilleResults", + "EventStudyResults", pd.DataFrame, ] @@ -67,9 +70,12 @@ def plot_event_study( Parameters ---------- - results : MultiPeriodDiDResults, CallawaySantAnnaResults, or DataFrame, optional - Results object from MultiPeriodDiD, CallawaySantAnna, or a DataFrame - with columns 'period', 'effect', 'se' (and optionally 'conf_int_lower', + results : MultiPeriodDiDResults, CallawaySantAnnaResults, EventStudyResults, or DataFrame, optional + Results object from MultiPeriodDiD, CallawaySantAnna, a unified + :class:`~diff_diff.results_base.EventStudyResults` container (from + any producer's ``aggregate('event_study')`` or + ``build_event_study_surface``), or a DataFrame with columns + 'period', 'effect', 'se' (and optionally 'conf_int_lower', 'conf_int_upper'). If None, must provide effects and se directly. effects : dict, optional Dictionary mapping periods to effect estimates. Used if results is None. @@ -87,7 +93,13 @@ def plot_event_study( post_periods : list, optional List of post-treatment periods. Used for shading. alpha : float, default=0.05 - Significance level for confidence intervals. + Significance level for confidence intervals. Applies to intervals + recomputed from the SE (dict/DataFrame-without-cband inputs and + explicit-normalization replots). An ``EventStudyResults`` + container's STORED intervals are drawn at the fit's own level - + bootstrap/t-based intervals cannot be re-leveled from the SE - and + a ``UserWarning`` names the mismatch when it differs from + ``alpha``. figsize : tuple, default=(10, 6) Figure size (width, height) in inches. title : str, default="Event Study" @@ -183,6 +195,7 @@ def plot_event_study( # Extract data from results if provided ci_lower_override = None ci_upper_override = None + reference_marks: Optional[Set[Any]] = None if results is not None: ( effects, @@ -194,14 +207,22 @@ def plot_event_study( reference_inferred, ci_lower_override, ci_upper_override, + pw_lower_override, + pw_upper_override, + reference_marks, ) = _extract_plot_data(results, periods, pre_periods, post_periods, reference_period) # If reference was inferred from results, it was NOT explicitly provided if reference_inferred: reference_period_explicit = False - # Suppress simultaneous confidence band overrides when user opts out + # Channel selection: the band channel when the user wants + # simultaneous bands; otherwise the POINTWISE channel - the + # producer's stored intervals (NaN preserved), for input routes + # that carry one. Routes without a pointwise channel keep the + # legacy behavior (overrides cleared -> pointwise recomputation + # from the SE downstream). if not use_cband: - ci_lower_override = None - ci_upper_override = None + ci_lower_override = pw_lower_override + ci_upper_override = pw_upper_override elif effects is None or se is None: raise ValueError("Must provide either 'results' or both 'effects' and 'se'") @@ -229,6 +250,39 @@ def plot_event_study( # Set reference SE to NaN (it's now a constraint, not an estimate) # This follows fixest convention where the omitted category has no SE/CI se = {p: (np.nan if p == reference_period else s) for p, s in se.items()} + # REGISTRY (Event Study Plotting): after explicit normalization, + # CIs are RECOMPUTED from the normalized effects and original + # SEs, and the reference CI becomes (NaN, NaN). Any simultaneous + # bands were computed around the UN-normalized effects - keeping + # them would draw intervals centered on stale estimates (and a + # finite interval on the constraint row) - so discard them and + # fall through to the pointwise recomputation below. + ci_lower_override = None + ci_upper_override = None + + # An EventStudyResults container's stored intervals are at the FIT's + # level - bootstrap-percentile / survey-t / Bell-McCaffrey intervals + # cannot be reconstructed from the SE at another level - so ``alpha`` + # does not apply to them. Warn instead of silently relabeling coverage. + # (The explicit-normalization path above discards the overrides and + # recomputes at the requested ``alpha``, so it never reaches here with + # overrides active.) + if results is not None and ci_lower_override is not None: + from diff_diff.results_base import EventStudyResults + + if isinstance(results, EventStudyResults): + stored_alpha = getattr(results, "alpha", None) + if stored_alpha is not None and not np.isclose(float(stored_alpha), alpha): + warnings.warn( + f"plot_event_study(alpha={alpha}) does not apply to an " + "EventStudyResults container: the stored intervals " + f"drawn here are at the fit's alpha={stored_alpha} " + f"({(1 - float(stored_alpha)) * 100:g}% coverage). " + "Re-aggregate from a fit at the desired level to " + "change the plotted coverage.", + UserWarning, + stacklevel=2, + ) plot_data = [] for period in periods: @@ -258,7 +312,11 @@ def plot_event_study( "se": std_err, "ci_lower": ci_lower, "ci_upper": ci_upper, - "is_reference": period == reference_period, + # Row-aligned marks (multi-reference containers) hollow + # every normalization anchor; the scalar covers the + # single-reference routes. + "is_reference": period == reference_period + or (reference_marks is not None and period in reference_marks), } ) @@ -359,19 +417,21 @@ def _render_event_study_mpl( ref_x = period_to_x[reference_period] ax.axvline(x=ref_x, color="gray", linestyle=":", linewidth=1, zorder=1) - # Plot error bars (only for entries with finite CI) + # Plot error bars (only for entries with finite CI). Bars are drawn + # ENDPOINT-BASED - anchored at the interval midpoint with symmetric + # yerr, so the segment is exactly [ci_lower, ci_upper]. Centering on + # the estimate breaks on stored percentile/bootstrap intervals, which + # need not contain it (negative yerr -> matplotlib ValueError). has_ci = df["ci_lower"].notna() & df["ci_upper"].notna() if has_ci.any(): df_with_ci = df[has_ci] x_with_ci = [period_to_x[p] for p in df_with_ci["period"]] - yerr = [ - df_with_ci["effect"] - df_with_ci["ci_lower"], - df_with_ci["ci_upper"] - df_with_ci["effect"], - ] + ci_mid = (df_with_ci["ci_lower"] + df_with_ci["ci_upper"]) / 2.0 + ci_half = (df_with_ci["ci_upper"] - df_with_ci["ci_lower"]).abs() / 2.0 ax.errorbar( x_with_ci, - df_with_ci["effect"], - yerr=yerr, + ci_mid, + yerr=ci_half, fmt="none", color=color, capsize=capsize, @@ -563,7 +623,20 @@ def _extract_plot_data( pre_periods: Optional[List[Any]], post_periods: Optional[List[Any]], reference_period: Optional[Any], -) -> Tuple[Dict, Dict, List, List, List, Any, bool, Optional[Dict], Optional[Dict]]: +) -> Tuple[ + Dict, + Dict, + List, + List, + List, + Any, + bool, + Optional[Dict], + Optional[Dict], + Optional[Dict], + Optional[Dict], + Optional[Set[Any]], +]: """ Extract plotting data from various result types. @@ -585,10 +658,142 @@ def _extract_plot_data( True if reference_period was auto-detected from results rather than explicitly provided by the user. ci_lower_override : dict or None - Simultaneous confidence band lower bounds, if available. + BAND-channel lower bounds (simultaneous bands, falling back + per-row to a stored pointwise interval where the producer carries + one), if available. Selected by the caller when ``use_cband=True``. ci_upper_override : dict or None - Simultaneous confidence band upper bounds, if available. + Band-channel upper bounds. + pw_lower_override : dict or None + POINTWISE-channel lower bounds: the producer's STORED intervals + (NaN preserved - a stored NaN means undefined inference and must + not be recomputed from the SE). Selected when ``use_cband=False``. + None for input routes that carry no stored intervals. + pw_upper_override : dict or None + Pointwise-channel upper bounds. + reference_marks : set or None + Row-aligned reference marking for surfaces with MULTIPLE + reference rows (the CS gapped-grid universal case): every listed + period renders hollow even though the scalar ``reference_period`` + cannot name them all. None on every single/zero-reference route. """ + # Unified EventStudyResults container (any producer - unlike HonestDiD / + # PreTrendsPower there is NO source scoping here: plotting is + # label-faithful for every producer, convention and time scale). + from diff_diff.results_base import EventStudyResults + + if isinstance(results, EventStudyResults): + surface = results + keys = list(surface.event_time.tolist()) + effects = {k: float(a) for k, a in zip(keys, surface.att)} + se = {k: float(s) for k, s in zip(keys, surface.se)} + + # Interval channels. The container carries the producer's ACTUAL + # inference (bootstrap-percentile, survey-t, Bell-McCaffrey + # intervals differ from att +/- z*se), so both channels cover + # EVERY row and preserve stored NaN (undefined inference must not + # be recomputed from a finite SE - e.g. a zero-SE row): + # - band channel: the simultaneous band where finite, else the + # stored pointwise interval (use_cband=True selection); + # - pointwise channel: the stored intervals verbatim + # (use_cband=False selection). + # (Explicit reference_period= normalization still discards the + # selected overrides in the renderer and recomputes around the + # normalized effects, per the REGISTRY plotting contract.) + band_lo: Dict[Any, float] = {} + band_hi: Dict[Any, float] = {} + pw_lo: Dict[Any, float] = {} + pw_hi: Dict[Any, float] = {} + for i, k in enumerate(keys): + ci_lo = float(surface.conf_int_lower[i]) + ci_hi = float(surface.conf_int_upper[i]) + pw_lo[k], pw_hi[k] = ci_lo, ci_hi + cb_lo = float(surface.cband_lower[i]) if surface.cband_lower is not None else np.nan + cb_hi = float(surface.cband_upper[i]) if surface.cband_upper is not None else np.nan + if np.isfinite(cb_lo) and np.isfinite(cb_hi): + band_lo[k], band_hi[k] = cb_lo, cb_hi + else: + band_lo[k], band_hi[k] = ci_lo, ci_hi + ci_lower_override: Optional[Dict[Any, float]] = band_lo + ci_upper_override: Optional[Dict[Any, float]] = band_hi + pw_lower_override: Optional[Dict[Any, float]] = pw_lo + pw_upper_override: Optional[Dict[Any, float]] = pw_hi + + if periods is None: + periods = keys + + # The renderer hollows rows by a per-row is_reference flag; the + # scalar reference_period is the single-reference convenience (it + # also drives the optional vertical reference line and explicit + # normalization). Exactly one is_reference row -> that label. + # MULTIPLE reference rows (legal on the container - the CS + # gapped-grid universal case carries one per cohort's positional + # base) are carried ROW-ALIGNED through ``reference_marks`` so + # every normalization anchor renders hollow - never silently + # dropped, never presented as a filled estimate. Explicitly + # renormalizing such a surface around a period that is NOT one of + # its reference rows fails closed: each anchor is a constraint + # under its own cohort base, so no single shift represents them + # faithfully. + reference_inferred = False + reference_marks: Optional[Set[Any]] = None + ref_rows = surface.event_time[surface.is_reference].tolist() + if len(ref_rows) > 1: + if reference_period is not None and reference_period not in ref_rows: + raise ValueError( + "Cannot renormalize an event-study container with " + f"multiple reference rows ({sorted(ref_rows)}) around " + f"reference_period={reference_period!r}: each reference " + "row is a normalization constraint under its own base, " + "so re-basing to a different period is not defined. " + "Plot without reference_period=, or re-estimate with a " + "single common reference." + ) + reference_marks = set(ref_rows) + elif reference_period is None and len(ref_rows) == 1: + reference_period = ref_rows[0] + reference_inferred = True + + if pre_periods is None or post_periods is None: + if surface.time_scale == "relative": + # Anticipation-aware split: with anticipation=k the window + # [e=-k, -1] carries anticipated TREATMENT effects + # (REGISTRY contract; HonestDiD/pretrends apply the same + # boundary), so pre-shading covers e < -k only. + _post_start = -int(surface.anticipation or 0) + derived_pre = [p for p in periods if p < _post_start] + derived_post = [p for p in periods if p >= _post_start] + else: + # Calendar labels (possibly str/Timestamp): numeric `p < 0` + # is undefined - split POSITIONALLY around the reference + # row (rows before it in event_time order are pre); with + # no reference row, all rows are post. + if reference_period is not None and reference_period in keys: + ref_pos = keys.index(reference_period) + derived_pre = [p for p in periods if p in keys and keys.index(p) < ref_pos] + derived_post = [p for p in periods if p in keys and keys.index(p) > ref_pos] + else: + derived_pre = [] + derived_post = [p for p in periods] + if pre_periods is None: + pre_periods = derived_pre + if post_periods is None: + post_periods = derived_post + + return ( + effects, + se, + periods, + pre_periods, + post_periods, + reference_period, + reference_inferred, + ci_lower_override, + ci_upper_override, + pw_lower_override, + pw_upper_override, + reference_marks, + ) + # Handle DataFrame input if isinstance(results, pd.DataFrame): if "period" not in results.columns: @@ -625,6 +830,9 @@ def _extract_plot_data( False, ci_lower_override, ci_upper_override, + None, + None, + None, ) # Handle MultiPeriodDiDResults @@ -665,6 +873,9 @@ def _extract_plot_data( ref_inferred, None, None, + None, + None, + None, ) # Handle ChaisemartinDHaultfoeuilleResults (dCDH event study) @@ -719,6 +930,9 @@ def _extract_plot_data( True, # inferred ci_lower_override if has_cband else None, ci_upper_override if has_cband else None, + None, + None, + None, ) # Handle CallawaySantAnnaResults (event study aggregation) @@ -768,11 +982,15 @@ def _extract_plot_data( if reference_period is None: reference_period = -1 + # Anticipation-aware split, mirroring the container branch: with + # anticipation=k the window [e=-k, -1] carries anticipated + # treatment effects, so pre-shading covers e < -k only. + _post_start = -int(getattr(results, "anticipation", 0) or 0) if pre_periods is None: - pre_periods = [p for p in periods if p < 0] + pre_periods = [p for p in periods if p < _post_start] if post_periods is None: - post_periods = [p for p in periods if p >= 0] + post_periods = [p for p in periods if p >= _post_start] return ( effects, @@ -784,13 +1002,16 @@ def _extract_plot_data( reference_inferred, ci_lower_override if has_cband else None, ci_upper_override if has_cband else None, + None, + None, + None, ) raise TypeError( f"Cannot extract plot data from {type(results).__name__}. " "Expected MultiPeriodDiDResults, CallawaySantAnnaResults, " "SunAbrahamResults, ImputationDiDResults, " - "ChaisemartinDHaultfoeuilleResults, or DataFrame." + "ChaisemartinDHaultfoeuilleResults, EventStudyResults, or DataFrame." ) @@ -870,7 +1091,54 @@ def plot_honest_event_study( raise ValueError("HonestDiDResults must have original_results to plot event study") # Extract data from original results - if hasattr(original_results, "period_effects"): + from diff_diff.results_base import EventStudyResults as _ESR + + stored_ci_lower: Optional[Dict[Any, float]] = None + stored_ci_upper: Optional[Dict[Any, float]] = None + if isinstance(original_results, _ESR): + # Unified container route (HonestDiD.fit stored the container the + # user passed): non-reference rows carry the estimates, and the + # container's STORED intervals are the producer's actual inference + # (bootstrap-percentile / survey-t / Bell-McCaffrey - NOT + # att +/- z*se); NaN stored bounds mean undefined inference and + # are preserved rather than recomputed from a finite SE. + keys = list(original_results.event_time.tolist()) + # Rows plotted: the rows HonestDiD retained (non-reference, finite + # positive SE) PLUS reference rows, which render as + # normalization-only anchors (stored att=0.0, NaN inference - no + # standard or honest interval). Undefined NON-reference rows were + # absent from beta_hat, so painting them with the aggregate honest + # interval would show sensitivity inference that was never + # computed, and per-period bounds are keyed on the retained set. + _retained = [ + i + for i, (r, s) in enumerate(zip(original_results.is_reference, original_results.se)) + if r or (np.isfinite(s) and float(s) > 0) + ] + effects_dict = {keys[i]: float(original_results.att[i]) for i in _retained} + se_dict = {keys[i]: float(original_results.se[i]) for i in _retained} + stored_ci_lower = {keys[i]: float(original_results.conf_int_lower[i]) for i in _retained} + stored_ci_upper = {keys[i]: float(original_results.conf_int_upper[i]) for i in _retained} + # Infer the container's single reference for the hollow-marker + # contract when the caller passed none (multi-reference containers + # cannot reach HonestDiD - the common-reference guard rejects + # them - so the single-scalar convention suffices here). + _ref_labels = [keys[i] for i, r in enumerate(original_results.is_reference) if r] + if reference_period is None and len(_ref_labels) == 1: + reference_period = _ref_labels[0] + if periods is None: + periods = sorted(effects_dict.keys()) + else: + _missing = [p for p in periods if p not in effects_dict] + if _missing: + raise ValueError( + f"Requested periods {_missing} are not retained by " + "HonestDiD on this event-study container: rows with " + "undefined inference (non-finite or zero SE) are " + "excluded from the sensitivity analysis and carry no " + "honest interval." + ) + elif hasattr(original_results, "period_effects"): # MultiPeriodDiDResults effects_dict = {p: pe.effect for p, pe in original_results.period_effects.items()} se_dict = {p: pe.se for p, pe in original_results.period_effects.items()} @@ -887,22 +1155,65 @@ def plot_honest_event_study( else: raise TypeError("Cannot extract event study data from original_results") - # Compute CIs + # Original CIs: the container's STORED intervals where available + # (NaN preserved); z-reconstruction from the SE only for input routes + # that carry no stored intervals (MPD / fit-time CS dict surfaces). alpha_val = honest_results.alpha + if stored_ci_lower is not None: + container_alpha = getattr(original_results, "alpha", None) + if container_alpha is not None and not np.isclose(float(container_alpha), alpha_val): + # Stored intervals cannot be re-leveled from the SE; label the + # mismatch instead of silently mixing coverage levels. + warnings.warn( + "plot_honest_event_study: the original confidence intervals " + "drawn are the container's stored intervals at " + f"alpha={container_alpha}, while the HonestDiD intervals " + f"are at alpha={alpha_val}. Re-aggregate the container at " + "the HonestDiD level for matching coverage.", + UserWarning, + stacklevel=2, + ) z = scipy_stats.norm.ppf(1 - alpha_val / 2) effects = [effects_dict[p] for p in periods] - original_ci_lower = [effects_dict[p] - z * se_dict[p] for p in periods] - original_ci_upper = [effects_dict[p] + z * se_dict[p] for p in periods] + if stored_ci_lower is not None and stored_ci_upper is not None: + original_ci_lower = [stored_ci_lower[p] for p in periods] + original_ci_upper = [stored_ci_upper[p] for p in periods] + else: + original_ci_lower = [effects_dict[p] - z * se_dict[p] for p in periods] + original_ci_upper = [effects_dict[p] + z * se_dict[p] for p in periods] # Get honest bounds if available for each period + _nan_bounds = {"ci_lb": np.nan, "ci_ub": np.nan} if honest_results.event_study_bounds: - honest_ci_lower = [honest_results.event_study_bounds[p]["ci_lb"] for p in periods] - honest_ci_upper = [honest_results.event_study_bounds[p]["ci_ub"] for p in periods] + # The reference row is a normalization constraint with no honest + # interval: tolerate its absence from per-period bounds (other + # rows stay strict - a missing real row is a caller error). + honest_ci_lower = [ + ( + honest_results.event_study_bounds.get(p, _nan_bounds) + if p == reference_period + else honest_results.event_study_bounds[p] + )["ci_lb"] + for p in periods + ] + honest_ci_upper = [ + ( + honest_results.event_study_bounds.get(p, _nan_bounds) + if p == reference_period + else honest_results.event_study_bounds[p] + )["ci_ub"] + for p in periods + ] else: - # Use scalar bounds applied to all periods - honest_ci_lower = [honest_results.ci_lb] * len(periods) - honest_ci_upper = [honest_results.ci_ub] * len(periods) + # Scalar bounds apply to every ESTIMATED period; the reference is + # a constraint, never painted with the aggregate honest interval. + honest_ci_lower = [ + np.nan if p == reference_period else honest_results.ci_lb for p in periods + ] + honest_ci_upper = [ + np.nan if p == reference_period else honest_results.ci_ub for p in periods + ] if backend == "plotly": return _render_honest_event_study_plotly( @@ -985,15 +1296,18 @@ def _render_honest_event_study_mpl( # Zero line ax.axhline(y=0, color="gray", linestyle="--", linewidth=1, alpha=0.5) + # Interval bars are drawn ENDPOINT-BASED (midpoint anchor, symmetric + # yerr) so each segment is exactly [lower, upper]: stored percentile/ + # bootstrap intervals - and honest bounds - need not contain the point + # estimate, and estimate-centered yerr goes negative there. + # Plot original CIs (thinner, background) - yerr_orig = [ - [e - lower for e, lower in zip(effects, original_ci_lower)], - [u - e for e, u in zip(effects, original_ci_upper)], - ] + mid_orig = [(lo + hi) / 2.0 for lo, hi in zip(original_ci_lower, original_ci_upper)] + half_orig = [abs(hi - lo) / 2.0 for lo, hi in zip(original_ci_lower, original_ci_upper)] ax.errorbar( x_vals, - effects, - yerr=yerr_orig, + mid_orig, + yerr=half_orig, fmt="none", color=original_color, capsize=capsize - 1, @@ -1003,14 +1317,12 @@ def _render_honest_event_study_mpl( ) # Plot honest CIs (thicker, foreground) - yerr_honest = [ - [e - lower for e, lower in zip(effects, honest_ci_lower)], - [u - e for e, u in zip(effects, honest_ci_upper)], - ] + mid_honest = [(lo + hi) / 2.0 for lo, hi in zip(honest_ci_lower, honest_ci_upper)] + half_honest = [abs(hi - lo) / 2.0 for lo, hi in zip(honest_ci_lower, honest_ci_upper)] ax.errorbar( x_vals, - effects, - yerr=yerr_honest, + mid_honest, + yerr=half_honest, fmt="none", color=honest_color, capsize=capsize, diff --git a/docs/api/_autosummary/diff_diff.CallawaySantAnnaResults.rst b/docs/api/_autosummary/diff_diff.CallawaySantAnnaResults.rst index 8a850cfc1..223336b92 100644 --- a/docs/api/_autosummary/diff_diff.CallawaySantAnnaResults.rst +++ b/docs/api/_autosummary/diff_diff.CallawaySantAnnaResults.rst @@ -52,6 +52,7 @@ ~CallawaySantAnnaResults.panel ~CallawaySantAnnaResults.pscore_fallback ~CallawaySantAnnaResults.pscore_trim + ~CallawaySantAnnaResults.reference_event_times ~CallawaySantAnnaResults.se ~CallawaySantAnnaResults.significance_stars ~CallawaySantAnnaResults.survey_metadata diff --git a/docs/api/_autosummary/diff_diff.ChaisemartinDHaultfoeuilleResults.rst b/docs/api/_autosummary/diff_diff.ChaisemartinDHaultfoeuilleResults.rst index 50fc42056..9c89a435c 100644 --- a/docs/api/_autosummary/diff_diff.ChaisemartinDHaultfoeuilleResults.rst +++ b/docs/api/_autosummary/diff_diff.ChaisemartinDHaultfoeuilleResults.rst @@ -12,6 +12,7 @@ .. autosummary:: ~ChaisemartinDHaultfoeuilleResults.__init__ + ~ChaisemartinDHaultfoeuilleResults.aggregate ~ChaisemartinDHaultfoeuilleResults.print_summary ~ChaisemartinDHaultfoeuilleResults.summary ~ChaisemartinDHaultfoeuilleResults.to_dataframe diff --git a/docs/api/_autosummary/diff_diff.EventStudyResults.rst b/docs/api/_autosummary/diff_diff.EventStudyResults.rst index e820809be..0e64dfcf3 100644 --- a/docs/api/_autosummary/diff_diff.EventStudyResults.rst +++ b/docs/api/_autosummary/diff_diff.EventStudyResults.rst @@ -24,12 +24,16 @@ .. autosummary:: ~EventStudyResults.alpha + ~EventStudyResults.anticipation + ~EventStudyResults.base_period ~EventStudyResults.cband_crit_value ~EventStudyResults.cband_lower ~EventStudyResults.cband_upper ~EventStudyResults.df + ~EventStudyResults.df_survey ~EventStudyResults.event_time_convention ~EventStudyResults.n_kind + ~EventStudyResults.reference_event_times ~EventStudyResults.reference_period ~EventStudyResults.reference_periods ~EventStudyResults.source diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index cc15d77f4..2ca0d3b70 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -1061,6 +1061,8 @@ Dynamic placebos `DID^{pl}_l` look backward from each group's reference period, - **Note (Phase 2 cost-benefit delta SE):** When `L_max >= 2`, `overall_att` holds the cost-benefit `delta`. Its SE is computed via the delta method from per-horizon SEs: `SE(delta) = sqrt(sum w_l^2 * SE(DID_l)^2)`, treating horizons as independent (conservative under Assumption 8). When bootstrap is enabled, per-horizon bootstrap SEs flow through the delta-method formula, so `overall_se` reflects bootstrap-derived per-horizon uncertainty but the delta aggregation itself uses normal-theory (not bootstrap percentile). This is an intentional exception to the general bootstrap-inference-surface contract: `overall_p_value` and `overall_conf_int` for `delta` use `safe_inference(delta, delta_se)`, not percentile bootstrap, because the delta is a derived aggregate rather than a directly bootstrapped estimand. +- **Note (post-fit `aggregate()` is a view - row M-026):** `results.aggregate('event_study')` and `results.aggregate('simple')` are pure VIEWS over stored fields - nothing is recomputed. `'event_study'` returns the unified `EventStudyResults` container (Phase-1 `L_max=None` fits: the 2-row l=1 view; `L_max >= 1`: the multi-horizon `l1_first_switch` surface); `'simple'` a one-row `AggregationResult` relaying `overall_att/se/t/p/CI` bit-exactly with the estimand-aware `target` label (DID_M / DID_1 / delta, or the trends-linear first-difference label whose overall row is all-NaN by design). Because it is a view, BOOTSTRAP FITS ARE PERMITTED - deviating from CallawaySantAnna's fail-closed `aggregate()`, whose kit-based recompute could silently substitute analytical inference; here each row relays exactly the inference the fit stored, including the delta's analytical-with-df numbers under bootstrap per the `Note (Phase 2 cost-benefit delta SE)` above (the view's `df` column resolves from the actual inference path, not the bootstrap-cleared `event_study_df` channel). The dCDH event-study CONTAINER is deliberately rejected by `compute_honest_did`/`compute_pretrends_power` - its l1 placebo semantics need HonestDiD's native dCDH branch (mandatory reinterpretation warning + horizon trimming); pass the results object itself. + - **Note (dynamic placebo SE - library extension):** Dynamic placebos `DID^{pl}_l` (negative horizons in `placebo_event_study`) now have analytical SE and bootstrap SE when `L_max >= 1`. The placebo IF uses the same cohort-recentered structure as positive horizons, applied to backward outcome differences `Y_{g, F_g-1-l} - Y_{g, F_g-1}` with the dual-eligibility control pool (forward + backward observation required). The paper's Theorem 1 variance result is stated for `DID_l`, not `DID^{pl}_l` - this extension applies the same IF/variance structure to the placebo estimand as a library enhancement. The single-period placebo `DID_M^pl` (`L_max=None`) retains NaN SE because the per-period aggregation path has no IF derivation. *Standard errors (Web Appendix Section 3.7.3 of the dynamic companion paper):* @@ -4646,6 +4648,7 @@ CRITICAL: δ_pre = β_pre pins pre-treatment violations to observed coefficients - **Note:** Phase 7d: survey variance support. When input results carry `survey_metadata` with `df_survey`, Delta^SD smoothness uses folded non-central t critical values (`scipy.stats.nct`); Delta^RM and naive FLCI paths use `_get_critical_value(alpha, df)` (standard t-distribution). `df_survey=0` → NaN inference. CallawaySantAnnaResults stores `event_study_vcov` (full cross-event-time VCV from IF vectors), which HonestDiD uses instead of the diagonal fallback. For replicate-weight designs, the event-study VCV falls back to diagonal (multivariate replicate VCV deferred). - **Note (deviation from R):** When HonestDiD receives bootstrap-fitted CallawaySantAnna results (`n_bootstrap > 0`), the full event-study covariance is unavailable (cleared to prevent mixing analytical VCV with bootstrap SEs). HonestDiD falls back to `diag(se^2)` from the bootstrap SEs with a UserWarning. R's `honest_did.AGGTEobj` computes a full covariance from the influence function matrix; implementing bootstrap event-study covariance is deferred. For full covariance structure in HonestDiD, use analytical SEs (`n_bootstrap=0`). - **Note (deviation from R):** When CallawaySantAnna results are passed to HonestDiD, `base_period != "universal"` emits a warning but does not error. R's `honest_did::honest_did.AGGTEobj` requires universal base period. Our implementation warns because the varying-base pre-treatment coefficients use consecutive comparisons (not a common reference), which changes the parallel-trends restriction interpretation. +- **Note (common-reference guard — gapped universal grids):** CallawaySantAnna `base_period="universal"` fits record `reference_event_times` — the distinct per-cohort positional-base event times (also carried on the `EventStudyResults` container). On a gapped time grid the positional bases land at different event times, and a cohort's base can OVERLAP another cohort's estimated horizon, where the aggregated row is a real estimate and NO reference-only row marks the anchor — `is_reference` cannot see it. When the field carries more than one entry, HonestDiD raises on BOTH input routes (native results and container) rather than returning bounds over coefficients normalized against different bases: Rambachan-Roth's `δ₀ = 0` normalization is defined around one common reference. PreTrendsPower applies the same guard (its hypothesized violation `δ` is likewise defined relative to one reference). Varying-base fits carry `reference_event_times=None` (no constant per-cohort reference exists) and are covered by the varying-base warnings instead. Universal-base inputs WITHOUT the field are never failed open: the native route re-derives the bases from the materialized reference cells (covers pre-3.9 pickles), and a provenance-less universal container (hand-built; CS-produced containers always record the field) warns fail-safe that the common reference cannot be verified. **Reference implementation(s):** - R: `HonestDiD` package (Rambachan & Roth's official package) @@ -4711,7 +4714,8 @@ Violation types: - **Note (paper-supported alternative — Wald pretest form):** the library retains the Wald noncentral-χ² form as `pretest_form='wald'`. NIS is the paper's primary analysis convention (used for all 12 surveyed papers' empirical exercises in Section I), but the Wald form is also a paper-supported alternative: Roth's Propositions 1, 3, and 4 apply to any (measurable) acceptance region for the conditional moments (Props 1+3) and to any convex acceptance region for the variance-reduction guarantee (Prop 4). The Wald ellipsoid is convex, so all four propositions apply. Wald is faster (no MVN CDF call) and matches the pre-PR-B shipped numerical baseline. Use Wald for backwards-compat / speed; use NIS for canonical paper alignment and R `pretrends` parity. -- **Note (convention — `linear` violation pattern, γ-unit MDV):** `_get_violation_weights('linear')` consumes actual pre-period relative-time labels threaded through `fit()` (PR-B 2026-05-17 resolution of the PR-A linear-pattern deviation). When `relative_times` is provided (e.g., `[-3, -2, -1]` for a regular grid or `[-5, -3, -1]` for an irregular grid), weights = `|t|` directly with NO L2 normalization, so `δ_pre = M · |t|` reflects Roth's `δ_t = γ · t` convention and the reported MDV equals γ. Callers that bypass `fit()` and supply only `n_pre` retain the previous count-based, L2-normalized `[n_pre-1, ..., 0]` direction (preserves shipped Wald numerical baselines for unit tests). **MPD period-label coverage:** for `MultiPeriodDiDResults`, the relative-time derivation in `_extract_pre_period_params` supports numeric labels (`int` / `float` / `np.int64`) and `pandas.Period` / `pandas.Timestamp` / `np.datetime64` (via Period or Timedelta arithmetic with units of frequency / days respectively). For genuinely non-numeric or unordered labels (string period IDs, unranked categoricals), the helper emits an explicit `UserWarning` and falls back to the legacy count-based normalized direction — the reported MDV is then NOT in Roth's γ units. Users on string period IDs who need γ-unit MDV should re-fit with numeric labels. +- **Note (convention — `linear` violation pattern, γ-unit MDV):** `_get_violation_weights('linear')` consumes actual pre-period relative-time labels threaded through `fit()` (PR-B 2026-05-17 resolution of the PR-A linear-pattern deviation). When `relative_times` is provided, weights = `|t|` directly with NO L2 normalization, so `δ_pre = M · |t|` reflects Roth's `δ_t = γ · t` convention and the reported MDV equals γ. **The relative times are REFERENCE-ANCHORED on common-reference routes (2026-08 correction):** Roth's violation is normalized like the coefficients — it must vanish at the omitted period (his labeling puts the reference at `t = 0`) — so the threaded values are `t - t_ref`, not raw treatment-relative event labels: MPD derives them from `reference_period` (`_coerce_relative_times_from_reference`), the CS universal route anchors at the singleton `reference_event_times`, SunAbraham at its omitted `e = -1 - anticipation`, and the `EventStudyResults` container at its single marked reference row (else the singleton provenance). Raw labels are retained ONLY where no common reference exists (CS `base_period="varying"` — covered by the varying-base warning and the TODO.md transformation row). Before the correction the CS/SA/container routes passed raw treatment-relative labels, overstating the violation at each pre-period by the reference offset (e.g. weights `[3, 2]` instead of `[2, 1]` for pre `[-3, -2]` around `t_ref = -1`) and understating γ-unit MDV correspondingly. Callers that bypass `fit()` and supply only `n_pre` retain the previous count-based, L2-normalized `[n_pre-1, ..., 0]` direction (preserves shipped Wald numerical baselines for unit tests). **MPD period-label coverage:** for `MultiPeriodDiDResults`, the relative-time derivation in `_extract_pre_period_params` supports numeric labels (`int` / `float` / `np.int64`) and `pandas.Period` / `pandas.Timestamp` / `np.datetime64` (via Period or Timedelta arithmetic with units of frequency / days respectively). For genuinely non-numeric or unordered labels (string period IDs, unranked categoricals), the helper emits an explicit `UserWarning` and falls back to the legacy count-based normalized direction — the reported MDV is then NOT in Roth's γ units. Users on string period IDs who need γ-unit MDV should re-fit with numeric labels. +- **Note (CS varying-base input — warned, not rejected):** the `linear` construction `δ_pre = M · |t|` assumes pre-period coefficients are LEVELS against one common reference period. CallawaySantAnna `base_period="varying"` (the CS default) pre-treatment effects are consecutive-period comparisons — under a linear underlying trend they are constant increments, not values proportional to `|t|` — so linear power/MDV computed on them target a different violation shape. Both CS-sourced input routes (the native results object and the `EventStudyResults` container, which also warns fail-safe on `base_period=None` provenance) emit a `UserWarning` recommending `base_period='universal'`, mirroring HonestDiD's universal-base warning. Universal-base fits on GAPPED grids do not merely warn — they fail closed via the common-reference guard (`reference_event_times` provenance; see the HonestDiD note), because cohort-specific positional bases make the violation space ill-defined regardless of the pattern. Transforming the varying-base violation vector through each coefficient's actual base mapping (or requiring universal base) is tracked in TODO.md. *Standard errors:* - Power calculations are exact (no sampling variability — power is computed against a hypothesized population trend, not estimated) @@ -4892,6 +4896,13 @@ should be a deliberate user choice. - Reference period CI becomes (NaN, NaN) after normalization (explicit only) - Reference period is plotted with hollow marker (both explicit and auto-inferred) - Reference period error bars: removed for explicit, retained for auto-inferred +- Multiple reference rows (an `EventStudyResults` container from CallawaySantAnna + `base_period="universal"` on a gapped grid): every marked row plots hollow at 0 with no + error bar via row-aligned marking - never dropped, never shown as a filled estimate; the + scalar reference (and its vertical line) stays unset. Explicit `reference_period=` + normalization is accepted only when the chosen period is one of the marked rows (a no-op + shift); any other period raises `ValueError`, because each anchor is a constraint under + its own cohort base and no single shift represents them faithfully **Reference implementation(s):** - R: `fixest::coefplot()` with reference category shown at 0 with no CI diff --git a/docs/v4-deprecations.yaml b/docs/v4-deprecations.yaml index d47d8a9af..837808632 100644 --- a/docs/v4-deprecations.yaml +++ b/docs/v4-deprecations.yaml @@ -300,11 +300,12 @@ rows: introduced_in: "3.9" deprecated_in: "3.9" removed_in: "4.0" - status: planned - phase: 2 + status: shimmed + phase: 5 warning: FutureWarning - code_refs: [diff_diff/chaisemartin_dhaultfoeuille.py] - notes: "dCDH placebo/normalized surfaces stay estimator-native; only the aggregation entry point moves." + test_ref: tests/test_aggregate_contract.py + code_refs: [diff_diff/chaisemartin_dhaultfoeuille.py, diff_diff/chaisemartin_dhaultfoeuille_results.py, diff_diff/aggregation.py, diff_diff/results_base.py] + notes: "Shimmed in 3.9: fit(aggregate=) warns via the shared NOT_SUPPLIED sentinel (a plain fit() never warns). The param NEVER computed aggregations on dCDH - every non-None value raised 'reserved for Phase 3' NotImplementedError since introduction - so there is no legacy surface to keep populating: a non-None value now raises ValueError pointing at the post-fit route, and the kwargs-splitting chaisemartin_dhaultfoeuille() wrapper forwards the shim unchanged. The successor is a pure VIEW, not a kit recompute: aggregate('event_study') returns build_event_study_surface(self) (Phase-1 L_max=None fits return the 2-row l=1 view - NOT an error; L_max>=1 the multi-horizon l1_first_switch surface), aggregate('simple') a one-row AggregationResult relaying overall_att/se/t/p/CI bit-exact with target = _estimand_label() (DID_M / DID_1 / delta / the trends-linear first-difference label whose overall row is all-NaN by design) and estimand-aware n/n_kind (N_S switcher_cells at L_max=None; N_1 groups at L_max==1; NaN/None for the delta, which averages horizon-specific N_l). SUPPORTED SUBSET simple|event_study - group/calendar fail closed via the mixin message; balance_e applies to no level here (empty _AGGREGATE_BALANCE_E_TYPES); weights= rejected. BOOTSTRAP FITS ARE PERMITTED, deviating from CallawaySantAnna's fail-closed rule, because nothing is recomputed - each row relays the STORED inference: percentile-bootstrap DID_M/DID_1 rows carry NaN df, while the L_max>=2 delta keeps analytical safe_inference with a possibly-finite survey df even under bootstrap (REGISTRY Note, Phase 2 cost-benefit delta SE), resolved by _overall_inference_df from event_study_df / the post-fit-refreshed survey_metadata.df_survey. dCDH placebo/normalized surfaces stay estimator-native; only the aggregation entry point moves. NOTE: the dCDH event-study CONTAINER is deliberately rejected by compute_honest_did/compute_pretrends_power (source-scoped admission; the l1 placebo semantics need honest_did's native dCDH branch) - see M-093." - id: M-027 kind: param group: aggregate-postfit @@ -983,7 +984,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). 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). Completed pre-cut by the remaining producer channels: SunAbraham (per-event dict) and de Chaisemartin-D'Haultfoeuille (scalar) event_study_df, plus LPDiD pooled_df for the headline pre/post windows - every producer whose inference records a df now exposes it." + 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). Completed pre-cut by the remaining producer channels: SunAbraham (per-event dict) and de Chaisemartin-D'Haultfoeuille (scalar) event_study_df, plus LPDiD pooled_df for the headline pre/post windows - every producer whose inference records a df now exposes it. Amended pre-cut a second time (2(b) PR-1, with M-026): three optional PROVENANCE fields appended last - base_period, anticipation, and df_survey (the fit's resolved SCALAR inference df beside the per-row channel: survey_metadata.df_survey with replicate-undefined mapping to the 0.0 fail-closed sentinel, else df_inference, else None - the per-row df column cannot encode that sentinel because __post_init__ NaNs it wherever p is non-finite) - threaded by the builders (the _empty_surface early return included) so the container consumers (HonestDiD's universal-base check, PreTrendsPower's anticipation cutoff, the honest df extraction) read fit-faithful values instead of dropping them. Amended pre-cut a third time (same PR): a fourth provenance field reference_event_times (also new on CallawaySantAnnaResults, computed at fit under base_period=universal) - the DISTINCT per-cohort positional-base event times, the common-reference signal is_reference cannot carry on gapped grids where a cohort's base overlaps another cohort's estimated horizon; HonestDiD and PreTrendsPower fail closed on more than one entry, on BOTH input routes (REGISTRY HonestDiD common-reference-guard Note)." - id: M-093 kind: behavior group: results-contract @@ -993,8 +994,8 @@ rows: removed_in: null status: planned phase: 5 - code_refs: [diff_diff/staggered_results.py, diff_diff/sun_abraham.py, diff_diff/imputation_results.py, diff_diff/two_stage_results.py, diff_diff/stacked_did_results.py, diff_diff/efficient_did_results.py, diff_diff/continuous_did_results.py, diff_diff/wooldridge_results.py, diff_diff/chaisemartin_dhaultfoeuille_results.py, diff_diff/lpdid_results.py, diff_diff/staggered_triple_diff_results.py, diff_diff/results.py, diff_diff/had.py, diff_diff/visualization/_event_study.py] - notes: "4.0 sentinel retirement + schema enforcement (spec section 5): the n_groups==0 / n_obs==0 reference-row sentinels retire; every estimator's to_dataframe(level='event_study') emits the [M-092] column schema; the plotter / HonestDiD / PreTrendsPower consume the unified surface. Thirteen ES-carrying source modules enumerated file-by-file (results.py covers MultiPeriodDiD + SpilloverDiD) plus the plotter. behavior-at-done requires test_ref." + code_refs: [diff_diff/staggered_results.py, diff_diff/sun_abraham.py, diff_diff/imputation_results.py, diff_diff/two_stage_results.py, diff_diff/stacked_did_results.py, diff_diff/efficient_did_results.py, diff_diff/continuous_did_results.py, diff_diff/wooldridge_results.py, diff_diff/chaisemartin_dhaultfoeuille_results.py, diff_diff/lpdid_results.py, diff_diff/staggered_triple_diff_results.py, diff_diff/results.py, diff_diff/had.py, diff_diff/visualization/_event_study.py, diff_diff/honest_did.py, diff_diff/pretrends.py] + notes: "4.0 sentinel retirement + schema enforcement (spec section 5): the n_groups==0 / n_obs==0 reference-row sentinels retire; every estimator's to_dataframe(level='event_study') emits the [M-092] column schema; the plotter / HonestDiD / PreTrendsPower consume the unified surface. Thirteen ES-carrying source modules enumerated file-by-file (results.py covers MultiPeriodDiD + SpilloverDiD) plus the plotter. behavior-at-done requires test_ref. PARTIAL PRE-CUT DELIVERY (2(b) PR-1, with M-026): the consumer half shipped for CS-SOURCED containers - compute_honest_did, compute_pretrends_power (both with SOURCE-SCOPED admission: source == CallawaySantAnnaResults only, rejecting dCDH l1 containers BY DESIGN since their placebo semantics need honest_did's native branch, and every other producer pending its own aggregate() migration) and plot_event_study / plot_honest_event_study (no source scoping - plotting is label-faithful). Admission widening is each later shim PR's methodology decision, not automatic. This row stays planned: its transition remains the 4.0 sentinel retirement + full-producer schema enforcement." # ---- Behavior policies (schema-tracked, spec-governed; no reality probe) - - id: M-080 diff --git a/docs/v4-design.md b/docs/v4-design.md index d96d5c8d5..c6bef7cd4 100644 --- a/docs/v4-design.md +++ b/docs/v4-design.md @@ -407,7 +407,14 @@ 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, 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")` +`event_study_df` channels), producer-provenance fields for the container +consumers (`base_period`, `anticipation`, the scalar `df_survey` +channel beside the per-row one - the per-row column cannot encode the +replicate-undefined 0.0 sentinel - and `reference_event_times`, the +distinct per-cohort positional-base event times whose multi-entry case +fails the HonestDiD/PreTrendsPower common-reference guard on gapped +universal grids; amended into [M-092] pre-cut with the +2(b) consumer delivery), 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_aggregate_contract.py b/tests/test_aggregate_contract.py index 1594bf666..249f89770 100644 --- a/tests/test_aggregate_contract.py +++ b/tests/test_aggregate_contract.py @@ -835,3 +835,277 @@ def test_staggered_triple_diff_overall_att_es_still_works(): ) assert res.overall_att_es is not None assert np.isfinite(res.overall_att_es) + + +# --------------------------------------------------------------------------- # +# dCDH (row M-026): fit(aggregate=) shim + the VIEW-based aggregate() +# --------------------------------------------------------------------------- # + +DCDH_KW = dict(outcome="outcome", unit="unit", time="period", treatment="treat") + + +def _dcdh_panel(seed=5, n_units=40, n_periods=6, switch_t=4): + rng = np.random.RandomState(seed) + rows = [] + for u in range(n_units): + s_t = switch_t if u < n_units // 2 else 10**6 + for t in range(1, n_periods + 1): + d = 1 if t >= s_t else 0 + rows.append( + { + "unit": u, + "period": t, + "outcome": u / 10 + 0.2 * t + 1.5 * d + rng.randn() * 0.3, + "treat": d, + } + ) + return pd.DataFrame(rows) + + +def _dcdh_survey_panel(seed=7, n_units=60, n_periods=6, switch_t=4): + df = _dcdh_panel(seed=seed, n_units=n_units, n_periods=n_periods, switch_t=switch_t) + df["survey_weights"] = 1.0 + 0.1 * (df["unit"] % 5) + df["strata"] = df["unit"] % 4 + df["psu"] = df["unit"] + return df + + +def _fit_dcdh(data, *, est_kw=None, **fit_kw): + from diff_diff.chaisemartin_dhaultfoeuille import ChaisemartinDHaultfoeuille + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return ChaisemartinDHaultfoeuille(**(est_kw or {})).fit(data, **DCDH_KW, **fit_kw) + + +@pytest.fixture(scope="module") +def dcdh_panel(): + return _dcdh_panel() + + +@pytest.fixture(scope="module") +def dcdh_fitted(dcdh_panel): + """Phase-1 fit (L_max=None).""" + return _fit_dcdh(dcdh_panel) + + +class TestDcdhShim: + def test_plain_fit_does_not_warn(self, dcdh_panel): + from diff_diff.chaisemartin_dhaultfoeuille import ChaisemartinDHaultfoeuille + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + ChaisemartinDHaultfoeuille().fit(dcdh_panel, **DCDH_KW) + assert [w for w in caught if issubclass(w.category, FutureWarning)] == [] + + def test_aggregate_kwarg_warns_even_at_none(self, dcdh_panel): + from diff_diff.chaisemartin_dhaultfoeuille import ChaisemartinDHaultfoeuille + + with pytest.warns(FutureWarning, match=r"fit\(aggregate=\) is deprecated"): + ChaisemartinDHaultfoeuille().fit(dcdh_panel, aggregate=None, **DCDH_KW) + + def test_non_none_value_warns_then_raises(self, dcdh_panel): + from diff_diff.chaisemartin_dhaultfoeuille import ChaisemartinDHaultfoeuille + + with pytest.warns(FutureWarning, match="aggregate"): + with pytest.raises(ValueError, match=r"results\.aggregate"): + ChaisemartinDHaultfoeuille().fit(dcdh_panel, aggregate="event_study", **DCDH_KW) + + def test_wrapper_forwarded_aggregate_warns(self, dcdh_panel): + # chaisemartin_dhaultfoeuille() splits **kwargs by signature and + # forwards non-__init__ names into fit(), so the shim is reachable + # through the wrapper too. + from diff_diff.chaisemartin_dhaultfoeuille import chaisemartin_dhaultfoeuille + + with pytest.warns(FutureWarning, match=r"fit\(aggregate=\) is deprecated"): + chaisemartin_dhaultfoeuille( + dcdh_panel, + outcome="outcome", + group="unit", + time="period", + treatment="treat", + aggregate=None, + ) + + +class TestDcdhAggregate: + def _assert_surface_matches_builder(self, res): + from diff_diff.results_base import build_event_study_surface + + es = res.aggregate("event_study") + assert isinstance(es, EventStudyResults) + built = build_event_study_surface(res) + # The dataclass's generated == raises on ndarray fields; compare + # to_dataframe rows per the file's precedent. + a, b = es.to_dataframe(), built.to_dataframe() + assert list(a.columns) == list(b.columns) + assert a.shape == b.shape + for col in a.columns: + av, bv = a[col].to_numpy(), b[col].to_numpy() + if av.dtype.kind in "fc": + assert np.allclose(av.astype(float), bv.astype(float), equal_nan=True) + else: + assert list(av) == list(bv) + return es + + def test_simple_view_bit_exact_phase1(self, dcdh_fitted): + agg = dcdh_fitted.aggregate("simple") + assert isinstance(agg, AggregationResult) + assert agg.level == "simple" + assert list(agg.label) == ["overall"] + assert list(agg.target) == ["DID_M"] + assert float(agg.att[0]) == dcdh_fitted.overall_att + assert float(agg.se[0]) == dcdh_fitted.overall_se + assert float(agg.t_stat[0]) == dcdh_fitted.overall_t_stat + assert float(agg.p_value[0]) == dcdh_fitted.overall_p_value + assert float(agg.conf_int_lower[0]) == dcdh_fitted.overall_conf_int[0] + assert float(agg.conf_int_upper[0]) == dcdh_fitted.overall_conf_int[1] + assert float(agg.n[0]) == float(dcdh_fitted.n_switcher_cells) + assert agg.n_kind == "switcher_cells" + # Non-survey analytical inference is z-based: no df. + assert np.isnan(agg.df[0]) + assert agg.estimator == "ChaisemartinDHaultfoeuille" + + def test_simple_view_lmax1_groups(self, dcdh_panel): + res = _fit_dcdh(dcdh_panel, L_max=1) + agg = res.aggregate("simple") + assert list(agg.target) == ["DID_1"] + assert agg.n_kind == "groups" + assert float(agg.n[0]) == float(res.n_switcher_cells) + + def test_simple_view_lmax2_delta(self, dcdh_panel): + res = _fit_dcdh(dcdh_panel, L_max=2) + agg = res.aggregate("simple") + assert list(agg.target) == ["delta"] + # The delta averages horizon-specific N_l: no truthful scalar count. + assert np.isnan(agg.n[0]) + assert agg.n_kind is None + assert float(agg.att[0]) == res.overall_att + + def test_simple_view_trends_linear_all_nan_relay(self, dcdh_panel): + # trends_linear + L_max>=2 suppresses the delta by design: every + # overall_* field is NaN and the estimand label points at + # linear_trends_effects. The view relays the all-NaN row honestly. + res = _fit_dcdh(dcdh_panel, L_max=2, trends_linear=True) + agg = res.aggregate("simple") + assert "fd" in str(agg.target[0]) + assert np.isnan(agg.att[0]) + assert np.isnan(agg.se[0]) + assert np.isnan(agg.p_value[0]) + assert np.isnan(agg.conf_int_lower[0]) and np.isnan(agg.conf_int_upper[0]) + assert np.isnan(agg.df[0]) + + def test_simple_view_bootstrap_percentile_relay(self, dcdh_panel): + # Bootstrap fits are PERMITTED (pure view): the row relays the + # stored percentile-bootstrap inference; df is NaN (no df used). + res = _fit_dcdh(dcdh_panel, est_kw=dict(n_bootstrap=49, seed=3)) + agg = res.aggregate("simple") + assert float(agg.att[0]) == res.overall_att + assert float(agg.se[0]) == res.overall_se + assert np.isnan(agg.df[0]) + + def test_event_study_container_threads_survey_df(self): + # CQ1 (local review R3): the dCDH builder threads the scalar + # df_survey provenance too - a survey fit's container must carry + # survey_metadata.df_survey, not None. + from diff_diff.survey import SurveyDesign + + df = _dcdh_survey_panel() + sd = SurveyDesign(weights="survey_weights", strata="strata", psu="psu") + res = _fit_dcdh(df, L_max=2, survey_design=sd) + surface = res.aggregate("event_study") + assert res.survey_metadata is not None + assert surface.df_survey == float(res.survey_metadata.df_survey) + + def test_simple_view_survey_analytical_df(self): + # Analytical survey fit: the stored p/CI used the survey df; the + # view relays it (event_study_df carries it here). + from diff_diff.survey import SurveyDesign + + df = _dcdh_survey_panel() + sd = SurveyDesign(weights="survey_weights", strata="strata", psu="psu") + res = _fit_dcdh(df, L_max=2, survey_design=sd) + agg = res.aggregate("simple") + assert res.survey_metadata is not None + expected = res.survey_metadata.df_survey + assert expected is not None and np.isfinite(agg.df[0]) + assert float(agg.df[0]) == float(expected) + + def test_simple_view_lmax2_survey_bootstrap_finite_df(self): + # THE df-provenance pin: under n_bootstrap>0 the event_study_df + # channel is cleared, but the L_max>=2 delta's stored p/CI still + # came from analytical safe_inference with the survey df (REGISTRY + # Note, Phase 2 cost-benefit delta SE). The view must report that + # finite df, not NaN. + from diff_diff.survey import SurveyDesign + + df = _dcdh_survey_panel() + sd = SurveyDesign(weights="survey_weights", strata="strata", psu="psu") + res = _fit_dcdh(df, L_max=2, survey_design=sd, est_kw=dict(n_bootstrap=49, seed=3)) + assert res.event_study_df is None # cleared under bootstrap + assert np.isfinite(res.overall_p_value) # delta stayed analytical + agg = res.aggregate("simple") + assert res.survey_metadata is not None + assert float(agg.df[0]) == float(res.survey_metadata.df_survey) + + def test_event_study_view_phase1_two_rows(self, dcdh_fitted): + es = self._assert_surface_matches_builder(dcdh_fitted) + # Phase-1 (L_max=None): the 2-row l=1 view, l1 convention - NOT an + # error (fit populates event_study_effects={1: ...} on this path). + assert es.event_time.tolist() == [0, 1] + assert es.event_time_convention == "l1_first_switch" + assert es.n_kind == "switcher_cells" + + def test_event_study_view_multi_horizon(self, dcdh_panel): + res = _fit_dcdh(dcdh_panel, L_max=2) + es = self._assert_surface_matches_builder(res) + assert 2 in es.event_time.tolist() + assert es.n_kind == "groups" + + def test_balance_e_rejected_empty_vocabulary(self, dcdh_fitted): + with pytest.raises(ValueError, match="no aggregation type on this estimator"): + dcdh_fitted.aggregate("event_study", balance_e=1) + + def test_weights_rejected(self, dcdh_fitted): + with pytest.raises(ValueError, match="does not accept a weights selector"): + dcdh_fitted.aggregate("simple", weights="cell") + + @pytest.mark.parametrize("bad", ["group", "calendar", "all", "nonsense"]) + def test_unsupported_types_fail_closed(self, dcdh_fitted, bad): + with pytest.raises(ValueError, match="Unsupported aggregation type"): + dcdh_fitted.aggregate(bad) + + def test_mixin_hooks_are_not_dataclass_fields(self): + # Regression: on a dataclass results class, annotating the mixin + # routing hooks without ClassVar turns them into __init__ fields, + # widening the public constructor/repr/equality surface. Enforced + # dynamically for EVERY dataclass that mixes AggregationMixin in, + # so later 2(b) waves are enrolled automatically. + import dataclasses + import inspect + + import diff_diff + from diff_diff.aggregation import AggregationMixin + + hooks = ("_AGGREGATE_SUPPORTED", "_AGGREGATE_BALANCE_E_TYPES") + checked = [] + for name in dir(diff_diff): + obj = getattr(diff_diff, name) + if ( + inspect.isclass(obj) + and issubclass(obj, AggregationMixin) + and obj is not AggregationMixin + and dataclasses.is_dataclass(obj) + ): + checked.append(name) + # dataclasses.fields() (not __dataclass_fields__, which + # also lists ClassVar pseudo-fields) = the real + # init/repr/eq surface. + fields = {f.name for f in dataclasses.fields(obj)} + params = inspect.signature(obj.__init__).parameters + for hook in hooks: + assert hook not in fields, f"{name}.{hook} leaked into fields" + assert hook not in params, f"{name}.{hook} leaked into __init__" + # The roster must at least cover the two shipped mixin adopters. + assert "CallawaySantAnnaResults" in checked + assert "ChaisemartinDHaultfoeuilleResults" in checked diff --git a/tests/test_business_report.py b/tests/test_business_report.py index ba9618739..df1a65d0b 100644 --- a/tests/test_business_report.py +++ b/tests/test_business_report.py @@ -2605,9 +2605,11 @@ def test_full_vcov_path_no_downgrade_on_real_cs_fit(self, cs_fit): assert block.get("status") == "ran", "pretrends_power should run on cs_fit" # Deterministic fixture pins (cs_fit at seed=7, treatment_effect=1.5): - # cov_source = full_pre_period_vcov; max_abs_pre_violation ≈ 0.375 - # (γ * max(|t|) where pre-periods are [-4, -3, -2]); |att| ≈ 1.779; - # mdv_share_of_att ≈ 0.211, well under 0.25 → tier = well_powered. + # cov_source = full_pre_period_vcov; max_abs_pre_violation ≈ 0.401 + # (γ * max(|t - t_ref|): the linear violation is REFERENCE-ANCHORED + # per Roth's normalization, so pre-periods [-4, -3, -2] around the + # universal e=-1 reference carry offsets [-3, -2, -1]); |att| ≈ + # 1.779; mdv_share_of_att ≈ 0.225, under 0.25 → tier = well_powered. # Codex R12 P1: this ratio is now `max_abs_pre_violation / |att|`, # the level-scale max pre-period violation under the MDV (post-PR-B # Step 4 linear MDV is in Roth's γ units, a slope; the level-scale @@ -2616,11 +2618,12 @@ def test_full_vcov_path_no_downgrade_on_real_cs_fit(self, cs_fit): "cs_fit is analytical CS with event_study_vcov populated — " "PR-B routing must report full_pre_period_vcov" ) - # max_abs_pre_violation = mdv * max(|t|) = 0.0937 * 4 ≈ 0.375 + # max_abs_pre_violation = mdv * max(|t - t_ref|) = 0.1337 * 3 ≈ 0.401 assert block.get("max_abs_pre_violation") is not None - assert 0.35 < block["max_abs_pre_violation"] < 0.40, ( + assert 0.38 < block["max_abs_pre_violation"] < 0.42, ( f"cs_fit max_abs_pre_violation={block['max_abs_pre_violation']} " - "should be ≈ 0.375 (γ ≈ 0.094 × max|t|=4)" + "should be ≈ 0.401 (γ ≈ 0.134 × max|t - t_ref|=3, " + "reference-anchored)" ) ratio = block["mdv_share_of_att"] assert ratio is not None and ratio < 0.25, ( diff --git a/tests/test_chaisemartin_dhaultfoeuille.py b/tests/test_chaisemartin_dhaultfoeuille.py index 26582b5c4..5892e01a1 100644 --- a/tests/test_chaisemartin_dhaultfoeuille.py +++ b/tests/test_chaisemartin_dhaultfoeuille.py @@ -306,28 +306,32 @@ def data(self): def _est(self): return ChaisemartinDHaultfoeuille() - def test_aggregate_simple_raises_not_implemented(self, data): - # aggregate is reserved for Phase 3; require aggregate=None exactly - with pytest.raises(NotImplementedError, match="Phase 3"): - self._est().fit( - data, - outcome="outcome", - unit="group", - time="period", - treatment="treatment", - aggregate="simple", - ) + def test_aggregate_simple_raises_value_error(self, data): + # M-026: fit(aggregate=) is deprecated and never computed anything; + # a non-None value warns then raises ValueError pointing at the + # post-fit results.aggregate() route. + with pytest.warns(FutureWarning, match="aggregate"): + with pytest.raises(ValueError, match="results.aggregate"): + self._est().fit( + data, + outcome="outcome", + unit="group", + time="period", + treatment="treatment", + aggregate="simple", + ) - def test_aggregate_event_study_raises_not_implemented(self, data): - with pytest.raises(NotImplementedError, match="Phase 3"): - self._est().fit( - data, - outcome="outcome", - unit="group", - time="period", - treatment="treatment", - aggregate="event_study", - ) + def test_aggregate_event_study_raises_value_error(self, data): + with pytest.warns(FutureWarning, match="aggregate"): + with pytest.raises(ValueError, match="results.aggregate"): + self._est().fit( + data, + outcome="outcome", + unit="group", + time="period", + treatment="treatment", + aggregate="event_study", + ) def test_L_max_validation(self, data): """L_max is now a Phase 2 feature: positive int or None accepted, diff --git a/tests/test_event_study_consumers.py b/tests/test_event_study_consumers.py new file mode 100644 index 000000000..b626b369c --- /dev/null +++ b/tests/test_event_study_consumers.py @@ -0,0 +1,1539 @@ +"""EventStudyResults consumer gates (2(b) PR-1, rows M-092/M-093 pre-cut half). + +``compute_honest_did``, ``compute_pretrends_power`` and ``plot_event_study`` / +``plot_honest_event_study`` accept the unified event-study container produced +by ``CallawaySantAnnaResults.aggregate('event_study')``. The gates: + +- END-TO-END (the TODO row's acceptance criteria): HonestDiD on a + ``base_period='universal'`` container; PreTrendsPower on an + ``anticipation=1`` container. +- ROUTE PARITY: for the same fit, the container route reproduces the + native route. HonestDiD outputs are deterministic and compared at + equality; PreTrendsPower's extraction tuple is compared bit-exactly, + while its end-to-end power gets a STOCHASTIC tolerance - scipy's Genz + multivariate-normal CDF is internally randomized (two native calls on + identical inputs differ at ~1e-5), so power equality at 1e-14 is not a + property even of the native route. +- SOURCE-SCOPED ADMISSION: honest/pretrends accept CS-sourced containers + only; dCDH l1 containers, calendar containers, non-CS/e0 containers and + hand-built source=None containers fail closed. The plotters take no + source guard (label-faithful rendering). +""" + +import warnings + +import numpy as np +import pandas as pd +import pytest + +from diff_diff import CallawaySantAnna, compute_honest_did, compute_pretrends_power +from diff_diff.pretrends import PreTrendsPower +from diff_diff.results_base import EventStudyResults + +FIT_KW = dict(outcome="y", unit="unit", time="time", first_treat="first_treat") + + +def _panel(seed=11, n_units=80, n_periods=8): + rng = np.random.RandomState(seed) + rows = [] + for u in range(n_units): + g = 4 if u < n_units // 3 else (6 if u < 2 * n_units // 3 else 0) + ui = rng.randn() * 2 + for t in range(1, n_periods + 1): + post = 1 if (g > 0 and t >= g) else 0 + rows.append( + { + "unit": u, + "time": t, + "first_treat": g, + "y": ui + 0.3 * t + 2.0 * post + rng.randn() * 0.5, + "cluster_col": u % 20, + "survey_weights": 1.0 + 0.1 * (u % 5), + "strata": u % 4, + "psu": u, + } + ) + return pd.DataFrame(rows) + + +def _fit_cs(data, **cs_kw): + # The deprecated fit-time aggregate= populates the NATIVE surface the + # route-parity tests compare against; the container side re-aggregates + # from the kit either way. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return CallawaySantAnna(**cs_kw).fit(data, aggregate="event_study", **FIT_KW) + + +@pytest.fixture(scope="module") +def panel(): + return _panel() + + +@pytest.fixture(scope="module") +def cs_universal(panel): + return _fit_cs(panel, base_period="universal") + + +@pytest.fixture(scope="module") +def cs_varying(panel): + return _fit_cs(panel) + + +def _tiny_container(**overrides): + """Hand-built 4-row relative container (one reference row at -1).""" + kwargs = dict( + event_time=np.array([-2, -1, 0, 1]), + att=np.array([0.1, 0.0, 1.9, 2.1]), + se=np.array([0.1, np.nan, 0.12, 0.13]), + t_stat=np.array([1.0, np.nan, 15.8, 16.2]), + p_value=np.array([0.3, np.nan, 0.0, 0.0]), + conf_int_lower=np.array([-0.1, np.nan, 1.66, 1.85]), + conf_int_upper=np.array([0.3, np.nan, 2.14, 2.35]), + is_reference=np.array([False, True, False, False]), + n=np.array([10.0, np.nan, 10.0, 10.0]), + source="CallawaySantAnnaResults", + ) + kwargs.update(overrides) + return EventStudyResults(**kwargs) + + +# --------------------------------------------------------------------------- # +# End-to-end acceptance (the TODO row's gates) +# --------------------------------------------------------------------------- # + + +class TestEndToEnd: + def test_honest_did_on_universal_container(self, cs_universal): + surface = cs_universal.aggregate("event_study") + h = compute_honest_did(surface, M=0.5) + assert np.isfinite(h.lb) and np.isfinite(h.ub) + assert np.isfinite(h.ci_lb) and np.isfinite(h.ci_ub) + + def test_pretrends_power_on_anticipation_container(self, panel): + res = _fit_cs(panel, anticipation=1) + surface = res.aggregate("event_study") + assert surface.anticipation == 1 + p = compute_pretrends_power(surface, M=0.1) + assert np.isfinite(p.power) + + +# --------------------------------------------------------------------------- # +# Route parity: container route == native route on the same fit +# --------------------------------------------------------------------------- # + +_HONEST_FIELDS = ("lb", "ub", "ci_lb", "ci_ub", "original_estimate", "original_se", "df_survey") + + +def _assert_honest_parity(res): + surface = res.aggregate("event_study") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + h_native = compute_honest_did(res, M=0.5) + h_container = compute_honest_did(surface, M=0.5) + for attr in _HONEST_FIELDS: + a, b = getattr(h_native, attr), getattr(h_container, attr) + if a is None or b is None: + assert a is b, attr + else: + np.testing.assert_allclose( + np.asarray(a, dtype=float), + np.asarray(b, dtype=float), + atol=1e-14, + rtol=1e-14, + equal_nan=True, + err_msg=attr, + ) + # Documented divergence: the container carries no survey-metadata + # object, so the stored field is None on the container route (its only + # inferential consumer is the df extraction, replaced by df_survey). + assert h_container.survey_metadata is None + return h_native, h_container + + +class TestHonestRouteParity: + def test_universal(self, cs_universal): + _assert_honest_parity(cs_universal) + + def test_varying(self, cs_varying): + _assert_honest_parity(cs_varying) + + def test_anticipation(self, panel): + _assert_honest_parity(_fit_cs(panel, anticipation=1)) + + def test_anticipation_window_is_post_not_pre(self, panel): + # REGISTRY anticipation contract: with anticipation=k the window + # [e=-k, -1] carries anticipated TREATMENT effects, so the clean + # pre-trend set is e < -k and beta_post starts at -k. Splitting + # at 0 misclassified e=-1 as a pre-trend coefficient IDENTICALLY + # on both routes - the parity gate alone could not catch it, so + # this pins the semantics directly. + res = _fit_cs(panel, anticipation=1) + surface = res.aggregate("event_study") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + h_native = compute_honest_did(res, M=0.5) + h_container = compute_honest_did(surface, M=0.5) + for h in (h_native, h_container): + assert -1 not in h.pre_periods_used + assert h.post_periods_used[0] == -1 + + def test_bare_cluster_df_threads(self, panel): + res = _fit_cs(panel, cluster="cluster_col") + h_native, h_container = _assert_honest_parity(res) + # bare-cluster fits carry df_inference -> finite scalar df on BOTH routes + assert h_container.df_survey is not None + assert np.isfinite(float(h_container.df_survey)) + + def test_survey_df_threads(self, panel): + from diff_diff.survey import SurveyDesign + + sd = SurveyDesign(weights="survey_weights", strata="strata", psu="psu") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = CallawaySantAnna().fit( + _panel(), survey_design=sd, aggregate="event_study", **FIT_KW + ) + h_native, h_container = _assert_honest_parity(res) + assert h_container.df_survey is not None + assert np.isfinite(float(h_container.df_survey)) + + def test_zero_se_rows_dropped_on_both_routes(self, panel): + # A zero-SE row carries undefined inference (safe_inference NaNs + # its t/p/CI); admitting it would launder that into finite honest + # bounds. Both routes drop it identically. Container side: + surface = _tiny_container( + base_period="universal", + se=np.array([0.0, np.nan, 0.12, 0.13]), + t_stat=np.array([np.nan, np.nan, 15.8, 16.2]), + p_value=np.array([np.nan, np.nan, 0.0, 0.0]), + conf_int_lower=np.array([np.nan, np.nan, 1.66, 1.85]), + conf_int_upper=np.array([np.nan, np.nan, 2.14, 2.35]), + ) + # The only pre-period row has se == 0 -> dropped -> no pre periods. + with pytest.raises(ValueError, match="pre-period"): + compute_honest_did(surface, M=0.5) + # Native side: inject a zero-SE pre row into a real fit's surface + # and assert the same drop (the row disappears from the retained + # pre set rather than entering beta with sigma 0). + res = _fit_cs(panel, base_period="universal") + from diff_diff.honest_did import _extract_event_study_params + + pre_key = min(k for k in res.event_study_effects if k < -1) + res.event_study_effects[pre_key] = dict( + res.event_study_effects[pre_key], se=0.0, p_value=np.nan + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + try: + out = _extract_event_study_params(res) + except ValueError: + out = None # grid gap after the drop - also a valid fail-closed + if out is not None: + assert pre_key not in out[4] + + def test_misaligned_container_vcov_raises(self): + # A SUPPLIED covariance whose index omits a retained horizon is + # inconsistent: fail loud, never silently degrade to diagonal + # (diag fallback is reserved for vcov=None). + surface = _tiny_container( + base_period="universal", + vcov=np.eye(2), + vcov_index=np.array([-2, 0]), # omits retained horizon 1 + ) + with pytest.raises(ValueError, match="vcov_index is missing"): + compute_honest_did(surface, M=0.5) + + def test_replicate_undefined_sentinel_relays(self): + # The container's df_survey=0.0 sentinel (replicate design with an + # undefined df) passes through to HonestDiDResults.df_survey exactly + # as the fit-time branch's sentinel does. + surface = _tiny_container(base_period="universal", df_survey=0.0) + h = compute_honest_did(surface, M=0.5) + assert h.df_survey == 0.0 + + @pytest.mark.parametrize("method", ["smoothness", "relative_magnitude"]) + def test_replicate_undefined_df_fails_closed_to_nan_ci(self, method): + # df_survey=0.0 means UNDEFINED inference: every FLCI path must + # yield NaN CI endpoints (never a silent normal-theory fallback) - + # incl. the optimal smoothness-FLCI with a full covariance, whose + # _cv_alpha/_flci_solve guards fail closed on a provided df <= 0. + surface = _tiny_container( + base_period="universal", + df_survey=0.0, + vcov=np.diag([0.01, 0.0144, 0.0169]), + vcov_index=np.array([-2, 0, 1]), + ) + h = compute_honest_did(surface, method=method, M=0.5) + assert np.isnan(h.ci_lb) and np.isnan(h.ci_ub) + + +class TestPretrendsRouteParity: + def test_extraction_bit_exact(self, panel): + res = _fit_cs(panel, anticipation=1) + surface = res.aggregate("event_study") + pt = PreTrendsPower() + e1, s1, v1, n1, r1, src1 = pt._extract_pre_period_params(res) + e2, s2, v2, n2, r2, src2 = pt._extract_pre_period_params(surface) + assert src1 == src2 == "full_pre_period_vcov" + assert n1 == n2 + np.testing.assert_array_equal(r1, r2) + # bit-exact: the container relays stored values verbatim + assert np.array_equal(e1, e2) + assert np.array_equal(s1, s2) + assert np.array_equal(v1, v2) + + def test_power_within_stochastic_tolerance(self, cs_varying): + surface = cs_varying.aggregate("event_study") + p_native = compute_pretrends_power(cs_varying, M=0.1) + p_container = compute_pretrends_power(surface, M=0.1) + # scipy's MVN CDF is internally randomized: identical inputs differ + # at ~1e-5 across calls, so this is a smoke bound, not 1e-14. + assert abs(p_native.power - p_container.power) < 1e-3 + + def test_explicit_pre_periods_honored_on_both_routes(self, panel): + # An explicitly requested pre-period subset must subset effects/ + # SEs/VCV on BOTH routes - never be silently ignored. + res = _fit_cs(panel) + surface = res.aggregate("event_study") + pt = PreTrendsPower() + full = pt._extract_pre_period_params(surface) + subset_labels = [int(t) for t in full[4][:2]] + e1, s1, v1, n1, r1, src1 = pt._extract_pre_period_params(res, subset_labels) + e2, s2, v2, n2, r2, src2 = pt._extract_pre_period_params(surface, subset_labels) + assert n1 == n2 == len(subset_labels) + np.testing.assert_array_equal(r1, r2) + assert np.array_equal(e1, e2) and np.array_equal(s1, s2) and np.array_equal(v1, v2) + assert v2.shape == (len(subset_labels), len(subset_labels)) + + def test_invalid_explicit_pre_periods_raise(self, panel): + res = _fit_cs(panel) + surface = res.aggregate("event_study") + pt = PreTrendsPower() + with pytest.raises(ValueError, match="not eligible"): + pt._extract_pre_period_params(surface, [999]) + with pytest.raises(ValueError, match="not eligible"): + pt._extract_pre_period_params(res, [999]) + + def test_empty_explicit_pre_periods_raise(self, panel): + # pre_periods=[] passes the per-label eligibility check vacuously; + # without a post-subset guard it reached zero-dimensional matrix + # logic downstream (opaque reshape error). All three explicit- + # subset paths reject it with the user-facing message. + from diff_diff import SunAbraham + + res = _fit_cs(panel) + surface = res.aggregate("event_study") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + sa = SunAbraham().fit(panel, **FIT_KW) + pt = PreTrendsPower() + for target in (res, surface, sa): + with pytest.raises(ValueError, match="at least one pre-period"): + pt._extract_pre_period_params(target, []) + + +# --------------------------------------------------------------------------- # +# Universal-base warning (fail-safe on missing provenance) +# --------------------------------------------------------------------------- # + + +class TestUniversalBaseWarning: + def test_varying_container_warns(self, cs_varying): + surface = cs_varying.aggregate("event_study") + assert surface.base_period == "varying" + with pytest.warns(UserWarning, match="base_period='universal'"): + compute_honest_did(surface, M=0.5) + + def test_missing_provenance_warns(self): + surface = _tiny_container(base_period=None) + with pytest.warns(UserWarning, match="no base_period provenance"): + compute_honest_did(surface, M=0.5) + + def test_universal_container_silent(self, cs_universal): + surface = cs_universal.aggregate("event_study") + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + compute_honest_did(surface, M=0.5) + assert not [w for w in caught if "base_period" in str(w.message)] + + +class TestCommonReferenceGuard: + """Cohort-level normalization-base provenance (reference_event_times). + + CS base_period='universal' on a GAPPED grid selects cohort-specific + positional bases. In the OVERLAP layout ({1,2,3,5}, cohorts {2,3,5}) + cohort 5's base (period 3, e=-2) coincides with cohort 3's estimated + pre-trend horizon, so the aggregated e=-2 row is a real estimate and + NO reference-only row marks that anchor - is_reference-based guards + cannot see it. The reference_event_times provenance field is the + authoritative signal: more than one distinct entry means the + coefficients were normalized against different bases, and HonestDiD / + PreTrendsPower fail closed on BOTH routes. + """ + + @staticmethod + def _fit_gapped(periods, cohorts, **fit_kw): + rng = np.random.RandomState(7) + rows = [] + coh = list(cohorts) + [0] + for u in range(160): + g = coh[u % len(coh)] + ufe = rng.randn() * 2 + for t in periods: + post = 1 if (g > 0 and t >= g) else 0 + rows.append( + { + "unit": u, + "time": t, + "first_treat": g, + "y": ufe + 0.3 * t + 2.0 * post + rng.randn() * 0.5, + } + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return CallawaySantAnna(n_bootstrap=0, base_period="universal").fit( + pd.DataFrame(rows), aggregate="event_study", **fit_kw, **FIT_KW + ) + + def test_overlap_layout_provenance(self): + res = self._fit_gapped((1, 2, 3, 5), (2, 3, 5)) + # Cohorts 2,3 base at e=-1; cohort 5's base (period 3) at e=-2. + assert tuple(int(e) for e in res.reference_event_times) == (-2, -1) + surface = res.aggregate("event_study") + assert tuple(int(e) for e in surface.reference_event_times) == (-2, -1) + # The overlapped anchor is INVISIBLE to is_reference: only e=-1 is + # a reference-only row; e=-2 aggregates cohort 3's real estimate. + marked = sorted(int(t) for t in surface.event_time[surface.is_reference]) + assert marked == [-1] + assert surface.to_dict()["reference_event_times"] == [-2, -1] + + def test_overlap_fails_closed_on_all_four_routes(self): + res = self._fit_gapped((1, 2, 3, 5), (2, 3, 5)) + surface = res.aggregate("event_study") + for consumer, target in ( + (compute_honest_did, res), + (compute_honest_did, surface), + (compute_pretrends_power, res), + (compute_pretrends_power, surface), + ): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match="common reference"): + consumer(target, M=0.5) + + def test_non_overlap_gapped_fails_closed(self): + # The {1,3,6} layout materializes every anchor as its own + # reference-only row; the provenance guard still fires first with + # the actionable common-reference message on both routes. + res = self._fit_gapped((1, 3, 6), (3, 6)) + assert tuple(int(e) for e in res.reference_event_times) == (-3, -2) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match="common reference"): + compute_honest_did(res, M=0.5) + with pytest.raises(ValueError, match="common reference"): + compute_honest_did(res.aggregate("event_study"), M=0.5) + + def test_regular_universal_singleton_passes(self, cs_universal): + assert tuple(int(e) for e in cs_universal.reference_event_times) == (-1,) + surface = cs_universal.aggregate("event_study") + assert tuple(int(e) for e in surface.reference_event_times) == (-1,) + h = compute_honest_did(surface, M=0.5) + assert np.isfinite(h.lb) and np.isfinite(h.ub) + + def test_varying_fit_carries_none(self, cs_varying): + # Varying base has no constant per-cohort reference: the field is + # None (unknown/NA), never invented - the varying-base WARNINGS + # cover that regime instead. + assert cs_varying.reference_event_times is None + assert cs_varying.aggregate("event_study").reference_event_times is None + + def test_balance_e_recomputes_provenance_over_retained_cohorts(self): + # The FIT-level tuple is fit-wide; the CONTAINER's must reflect + # the cohorts the aggregation actually retained. balance_e can + # drop the cohort responsible for the second base - a stale tuple + # would reject a balanced surface whose remaining cohorts share + # one reference. + res = self._fit_gapped((1, 2, 3, 4, 6), (3, 6)) + # Cohort 3 base at e=-1 (period 2); cohort 6 base at e=-2 (period + # 4, the positional neighbor on the gapped grid). + assert tuple(int(e) for e in res.reference_event_times) == (-2, -1) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + unbalanced = res.aggregate("event_study") + balanced = res.aggregate("event_study", balance_e=1) + + # Unbalanced surface: both cohorts retained -> both bases -> the + # common-reference guard fires. + assert tuple(int(e) for e in unbalanced.reference_event_times) == (-2, -1) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match="common reference"): + compute_pretrends_power(unbalanced, M=0.1) + + # balance_e=1 retains only cohort 3 (the only cohort with an + # effect at e=1): the surface-faithful provenance is the single + # remaining base, and the consumer accepts the surface. + assert tuple(int(e) for e in balanced.reference_event_times) == (-1,) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + p = compute_pretrends_power(balanced, M=0.1) + assert np.isfinite(p.power) + + def test_missing_native_provenance_derives_from_reference_cells(self): + # A provenance-less universal result (pre-3.9 pickle or + # replace()-stripped copy) must not FAIL OPEN: the cohort bases + # are re-derived from the materialized reference cells, so the + # mixed-base layout still fails closed on both consumers. + import dataclasses + + res = self._fit_gapped((1, 2, 3, 5), (2, 3, 5)) + stripped = dataclasses.replace(res, reference_event_times=None) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match="common reference"): + compute_honest_did(stripped, M=0.5) + with pytest.raises(ValueError, match="common reference"): + compute_pretrends_power(stripped, M=0.1) + + def test_missing_container_provenance_warns(self): + # A hand-built universal container without the field cannot be + # verified (no cells to derive from): warn fail-safe, never fail + # open silently. CS-produced containers always record the field. + surface = _tiny_container(base_period="universal") + assert surface.reference_event_times is None + with pytest.warns(UserWarning, match="no reference_event_times provenance"): + compute_honest_did(surface, M=0.5) + with pytest.warns(UserWarning, match="no reference_event_times provenance"): + compute_pretrends_power(surface, M=0.1) + + def test_fit_time_balance_e_provenance_matches_surface(self): + # The deprecated fit-time aggregate="event_study" + balance_e + # stores a RESTRICTED surface: its provenance must describe the + # retained cohorts too, or the native route would reject a fit + # whose equivalent post-fit container is accepted (route parity). + res = self._fit_gapped((1, 2, 3, 4, 6), (3, 6), balance_e=1) + assert tuple(int(e) for e in res.reference_event_times) == (-1,) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + p = compute_pretrends_power(res, M=0.1) + assert np.isfinite(p.power) + + def test_to_dict_reference_event_times_json_safe(self): + # CS period arithmetic yields numpy scalars; to_dict must emit + # JSON-serializable labels. + import json + + res = self._fit_gapped((1, 3, 6), (3, 6)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + surface = res.aggregate("event_study") + d = surface.to_dict() + assert d["reference_event_times"] == [-3, -2] + json.dumps(d) # must not raise on numpy-labeled provenance + + +class TestContainerIntegrity: + """Hand-built containers with malformed rows/covariance fail closed. + + Containers are publicly constructible: consumers subset by explicit + [sorted pre; sorted post] label order (row order is not trusted) and + validate covariance integrity at the boundary. + """ + + @staticmethod + def _container(order, vcov=None, vcov_index=None, se_override=None): + data = { + -3: (0.10, 0.10, -0.10, 0.30), + -2: (0.05, 0.10, -0.15, 0.25), + -1: (0.0, np.nan, np.nan, np.nan), + 0: (1.9, 0.12, 1.66, 2.14), + 1: (2.1, 0.13, 1.85, 2.35), + } + rows = [data[t] for t in order] + se = np.array([r[1] for r in rows]) + if se_override is not None: + se = se_override + return EventStudyResults( + event_time=np.array(order), + att=np.array([r[0] for r in rows]), + se=se, + t_stat=np.array([np.nan] * len(order)), + p_value=np.array([np.nan] * len(order)), + conf_int_lower=np.array([r[2] for r in rows]), + conf_int_upper=np.array([r[3] for r in rows]), + is_reference=np.array([t == -1 for t in order]), + n=np.array([10.0] * len(order)), + source="CallawaySantAnnaResults", + base_period="universal", + reference_event_times=(-1,), + vcov=vcov, + vcov_index=vcov_index, + ) + + def test_permuted_rows_produce_identical_bounds(self): + # Interleaved rows must yield the SAME bounds as the sorted + # container: beta_hat/sigma are subset in [sorted pre; sorted + # post] order, never row order (the fit-side split takes the + # first num_pre entries as beta_pre). + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + h_sorted = compute_honest_did(self._container([-3, -2, -1, 0, 1]), M=0.5) + h_perm = compute_honest_did(self._container([-3, 0, -2, -1, 1]), M=0.5) + assert h_perm.lb == h_sorted.lb and h_perm.ub == h_sorted.ub + assert h_perm.pre_periods_used == h_sorted.pre_periods_used == [-3, -2] + assert h_perm.post_periods_used == h_sorted.post_periods_used == [0, 1] + # Pretrends is label-aligned elementwise: power invariant too. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + p_sorted = compute_pretrends_power(self._container([-3, -2, -1, 0, 1]), M=0.1) + p_perm = compute_pretrends_power(self._container([-3, 0, -2, -1, 1]), M=0.1) + assert abs(p_sorted.power - p_perm.power) < 1e-3 # MVN-CDF jitter + + def test_reversed_rows_last_period_violation_invariant(self): + # Positional violation patterns (last_period assigns weights[-1] + # to the FINAL entry) require chronological pre-period order, not + # row order: a reversed hand-built container must produce the + # same power as the sorted one. + pt_kwargs = dict(M=0.5, violation_type="last_period") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + p_sorted = compute_pretrends_power(self._container([-3, -2, -1, 0, 1]), **pt_kwargs) + p_rev = compute_pretrends_power(self._container([-2, -3, -1, 0, 1]), **pt_kwargs) + assert abs(p_sorted.power - p_rev.power) < 1e-3 # MVN-CDF jitter + # Extraction-level exactness: chronological labels either way. + pt = PreTrendsPower() + rel_sorted = pt._extract_pre_period_params(self._container([-3, -2, -1, 0, 1]))[4] + rel_rev = pt._extract_pre_period_params(self._container([-2, -3, -1, 0, 1]))[4] + np.testing.assert_array_equal(rel_sorted, rel_rev) + + def test_duplicate_event_time_labels_rejected(self): + surface = self._container([-3, -2, -1, 0, 0]) + for consumer in (compute_honest_did, compute_pretrends_power): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match="duplicate event_time"): + consumer(surface, M=0.5) + + @pytest.mark.parametrize( + "corruption, match", + [ + ("nonfinite", "non-finite"), + ("asymmetric", "not symmetric"), + ("indefinite", "indefinite"), + ("diag_mismatch", "inconsistent with the stored standard errors"), + ("dup_index", "duplicate\\s+labels|carries duplicate"), + ], + ) + def test_malformed_covariance_rejected(self, corruption, match): + order = [-3, -2, -1, 0, 1] + ses = np.array([0.10, 0.10, 0.12, 0.13]) # retained rows, sorted + vcov = np.diag(ses**2) + vcov_index = np.array([-3, -2, 0, 1]) + if corruption == "nonfinite": + vcov = vcov.copy() + vcov[0, 1] = np.nan + vcov[1, 0] = np.nan + elif corruption == "asymmetric": + vcov = vcov.copy() + vcov[0, 1] = 0.005 # not mirrored + elif corruption == "indefinite": + vcov = vcov.copy() + # off-diagonal larger than the diagonal product -> negative eig + vcov[0, 1] = vcov[1, 0] = 0.02 + elif corruption == "diag_mismatch": + vcov = vcov.copy() + vcov[0, 0] = 0.5 # != se**2 + elif corruption == "dup_index": + vcov_index = np.array([-3, -3, 0, 1]) + surface = self._container(order, vcov=vcov, vcov_index=vcov_index) + for consumer in (compute_honest_did, compute_pretrends_power): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match=match): + consumer(surface, M=0.5) + + def test_low_scale_indefinite_rejected(self): + # Tolerances are RELATIVE to the matrix scale: a uniformly tiny + # indefinite matrix (diag 1e-10, eigenvalues [-1e-10, 3e-10]) + # must not slip under an absolute floor. + vcov = np.diag(np.full(4, 1e-10)) + vcov[0, 1] = vcov[1, 0] = 2e-10 + se_override = np.array([1e-5, 1e-5, np.nan, 1e-5, 1e-5]) + surface = self._container( + [-3, -2, -1, 0, 1], + vcov=vcov, + vcov_index=np.array([-3, -2, 0, 1]), + se_override=se_override, + ) + for consumer in (compute_honest_did, compute_pretrends_power): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match="indefinite"): + consumer(surface, M=0.5) + + def test_singular_covariance_honest_rejects_pretrends_accepts(self): + # Perfectly-correlated pre-rows: PSD but SINGULAR. HonestDiD + # rejects (Rambachan-Roth assumes covariance eigenvalues bounded + # away from zero); PreTrendsPower keeps its documented + # singular-covariance handling. + ses = np.array([0.1, 0.1, 0.12, 0.13]) + vcov = np.diag(ses**2) + vcov[0, 1] = vcov[1, 0] = 0.01 # corr = 1 between the pre rows + surface = self._container( + [-3, -2, -1, 0, 1], vcov=vcov, vcov_index=np.array([-3, -2, 0, 1]) + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match="singular"): + compute_honest_did(surface, M=0.5) + p = compute_pretrends_power(surface, M=0.1) + assert p is not None + + def test_valid_covariance_still_accepted(self): + ses = np.array([0.10, 0.10, 0.12, 0.13]) + vcov = np.diag(ses**2) + surface = self._container( + [-3, -2, -1, 0, 1], vcov=vcov, vcov_index=np.array([-3, -2, 0, 1]) + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + h = compute_honest_did(surface, M=0.5) + p = compute_pretrends_power(surface, M=0.1) + assert np.isfinite(h.lb) and np.isfinite(h.ub) + assert np.isfinite(p.power) + + +class TestLinearViolationAnchoring: + """Roth's linear violation is anchored at the omitted reference. + + Roth labels the omitted period t=0, so delta = gamma*t vanishes there + by construction; translated to estimator-native labels the threaded + relative times are t - t_ref, NOT raw treatment-relative labels + (which overstate the violation by the reference offset). MPD already + anchored via _coerce_relative_times_from_reference; these pin the + CS-universal, SunAbraham and container routes. + """ + + def test_universal_weights_are_reference_relative(self, panel): + res = _fit_cs(panel, base_period="universal") + surface = res.aggregate("event_study") + keep = ( + (~surface.is_reference) + & np.isfinite(surface.se) + & (surface.se > 0) + & (surface.event_time < 0) + ) + labels = surface.event_time[keep].astype(float) + pt = PreTrendsPower() + rel_native = pt._extract_pre_period_params(res)[4] + rel_container = pt._extract_pre_period_params(surface)[4] + expected = labels - (-1.0) # anchored at the e=-1 reference + np.testing.assert_array_equal(rel_native, expected) + np.testing.assert_array_equal(rel_container, expected) + + def test_universal_anticipation_anchor_and_weights(self, panel): + # anticipation=1: reference at e=-2; pre labels t < -1 anchor + # there, and the hand-calculated weight vector is |t - t_ref|. + res = _fit_cs(panel, base_period="universal", anticipation=1) + surface = res.aggregate("event_study") + assert surface.reference_period == -2 + keep = ( + (~surface.is_reference) + & np.isfinite(surface.se) + & (surface.se > 0) + & (surface.event_time < -1) + ) + expected = surface.event_time[keep].astype(float) + 2.0 + pt = PreTrendsPower() + rel = pt._extract_pre_period_params(surface)[4] + np.testing.assert_array_equal(rel, expected) + w = pt._get_violation_weights(len(rel), relative_times=rel) + np.testing.assert_array_equal(w, np.abs(expected)) + + def test_anchor_invariant_to_label_origin(self): + # Containers identical up to a label SHIFT (reference at -1 vs 0) + # extract IDENTICAL reference-relative times - the violation is + # anchored at the reference, not at the label origin. + common = dict( + att=np.array([0.1, 0.05, 0.0, 1.9]), + se=np.array([0.1, 0.1, np.nan, 0.12]), + t_stat=np.array([1.0, 0.5, np.nan, 15.8]), + p_value=np.array([0.3, 0.6, np.nan, 0.0]), + conf_int_lower=np.array([-0.1, -0.15, np.nan, 1.66]), + conf_int_upper=np.array([0.3, 0.25, np.nan, 2.14]), + is_reference=np.array([False, False, True, False]), + n=np.array([10.0, 10.0, np.nan, 10.0]), + source="CallawaySantAnnaResults", + base_period="universal", + ) + a = EventStudyResults( + event_time=np.array([-3, -2, -1, 0]), reference_event_times=(-1,), **common + ) + b = EventStudyResults( + event_time=np.array([-2, -1, 0, 1]), reference_event_times=(0,), **common + ) + pt = PreTrendsPower() + rel_a = pt._extract_pre_period_params(a)[4] + rel_b = pt._extract_pre_period_params(b)[4] + np.testing.assert_array_equal(rel_a, rel_b) + np.testing.assert_array_equal(rel_a, np.array([-2.0, -1.0])) + + def test_sun_abraham_anchor(self, panel): + from diff_diff import SunAbraham + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + sa = SunAbraham().fit(panel, **FIT_KW) + assert sa.reference_period == -1 + pt = PreTrendsPower() + effects, _, _, n_pre, rel, _ = pt._extract_pre_period_params(sa) + # Anchored at the omitted e = -1 - anticipation = -1. + labels = sorted( + t + for t, d in sa.event_study_effects.items() + if t < 0 and np.isfinite(d.get("se", np.nan)) and float(d.get("se", 0.0)) > 0 + ) + np.testing.assert_array_equal(rel, np.asarray(labels, dtype=float) + 1.0) + + +class TestVaryingBasePretrendsWarning: + """Twin of HonestDiD's universal-base warning, on both pretrends routes. + + The built-in ``linear`` violation constructs delta as a slope on + relative time (level coefficients against one common reference); CS + varying-base pre-treatment effects are consecutive-period comparisons, + so linear power/MDV target a different violation shape (REGISTRY + PreTrendsPower Note; full fix tracked in TODO.md). + """ + + def test_varying_native_warns(self, cs_varying): + with pytest.warns(UserWarning, match="base_period='universal'"): + compute_pretrends_power(cs_varying, M=0.1) + + def test_varying_container_warns(self, cs_varying): + surface = cs_varying.aggregate("event_study") + with pytest.warns(UserWarning, match="base_period='universal'"): + compute_pretrends_power(surface, M=0.1) + + def test_missing_provenance_warns(self): + surface = _tiny_container(base_period=None) + with pytest.warns(UserWarning, match="no base_period provenance"): + compute_pretrends_power(surface, M=0.1) + + def test_non_linear_violation_does_not_warn(self, cs_varying): + # The warning concerns the built-in LINEAR construction only; + # constant/last_period/custom vectors are user-specified in + # coefficient space. + surface = cs_varying.aggregate("event_study") + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + compute_pretrends_power(cs_varying, M=0.1, violation_type="constant") + compute_pretrends_power(surface, M=0.1, violation_type="constant") + assert not [w for w in caught if "base_period" in str(w.message)] + + def test_universal_silent_on_both_routes(self, panel): + res = _fit_cs(panel, base_period="universal") + surface = res.aggregate("event_study") + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + compute_pretrends_power(res, M=0.1) + compute_pretrends_power(surface, M=0.1) + assert not [w for w in caught if "base_period" in str(w.message)] + + +# --------------------------------------------------------------------------- # +# Provenance threading (incl. the requested-but-empty path) +# --------------------------------------------------------------------------- # + + +class TestProvenanceThreading: + def test_container_carries_fit_provenance(self, panel): + res = _fit_cs(panel, base_period="universal", anticipation=1) + surface = res.aggregate("event_study") + assert surface.base_period == "universal" + assert surface.anticipation == 1 + assert surface.df_survey is None # no survey design, no cluster df + + def test_bare_cluster_container_df(self, panel): + res = _fit_cs(panel, cluster="cluster_col") + surface = res.aggregate("event_study") + assert surface.df_survey is not None and np.isfinite(surface.df_survey) + + +# --------------------------------------------------------------------------- # +# Source-scoped admission (fail-closed) +# --------------------------------------------------------------------------- # + + +def _dcdh_container(): + from diff_diff.chaisemartin_dhaultfoeuille import ChaisemartinDHaultfoeuille + + rng = np.random.RandomState(5) + rows = [] + for u in range(30): + s_t = 4 if u < 15 else 10**6 + for t in range(1, 7): + d = 1 if t >= s_t else 0 + rows.append( + { + "unit": u, + "period": t, + "outcome": u / 10 + 0.2 * t + 1.5 * d + rng.randn() * 0.3, + "treat": d, + } + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = ChaisemartinDHaultfoeuille().fit( + pd.DataFrame(rows), outcome="outcome", unit="unit", time="period", treatment="treat" + ) + return res.aggregate("event_study") + + +class TestSourceScopedAdmission: + def test_dcdh_container_rejected_by_honest(self): + surface = _dcdh_container() + assert surface.event_time_convention == "l1_first_switch" + with pytest.raises(TypeError, match="CallawaySantAnnaResults.aggregate"): + compute_honest_did(surface, M=0.5) + + def test_dcdh_container_rejected_by_pretrends_without_dead_route(self): + surface = _dcdh_container() + with pytest.raises(TypeError) as exc_info: + compute_pretrends_power(surface, M=0.1) + msg = str(exc_info.value) + # pretrends' native accepted set has NO dCDH branch - the message + # must name ITS OWN natives, never point dCDH at a dead route. + assert "SunAbrahamResults" in msg + assert "natively" not in msg or "ChaisemartinDHaultfoeuille" not in msg + + def test_hand_built_source_none_rejected(self): + surface = _tiny_container(source=None) + with pytest.raises(TypeError, match="source=None"): + compute_honest_did(surface, M=0.5) + with pytest.raises(TypeError, match="source=None"): + compute_pretrends_power(surface, M=0.1) + + def test_non_cs_e0_source_rejected(self): + surface = _tiny_container(source="ImputationDiDResults") + with pytest.raises(TypeError, match="ImputationDiDResults"): + compute_honest_did(surface, M=0.5) + with pytest.raises(TypeError, match="ImputationDiDResults"): + compute_pretrends_power(surface, M=0.1) + + def test_calendar_scale_rejected(self): + # Belt-and-suspenders: even a CS-sourced container is rejected on a + # non-relative time scale (CS never emits calendar). + surface = _tiny_container( + event_time=np.array(["2018", "2019", "2020", "2021"], dtype=object), + time_scale="calendar", + ) + with pytest.raises(TypeError, match="relative"): + compute_honest_did(surface, M=0.5) + with pytest.raises(TypeError, match="relative"): + compute_pretrends_power(surface, M=0.1) + + def test_multiple_reference_rows_fail_closed_in_honest(self): + # DELIBERATE deviation from the fit-time branch, which silently + # splits around the FIRST n_groups==0 marker in dict order: the + # container branch refuses - the consecutive-grid contract is + # defined around a single omitted reference. + surface = _tiny_container( + is_reference=np.array([True, True, False, False]), + att=np.array([0.0, 0.0, 1.9, 2.1]), + se=np.array([np.nan, np.nan, 0.12, 0.13]), + t_stat=np.array([np.nan, np.nan, 15.8, 16.2]), + p_value=np.array([np.nan, np.nan, 0.0, 0.0]), + conf_int_lower=np.array([np.nan, np.nan, 1.66, 1.85]), + conf_int_upper=np.array([np.nan, np.nan, 2.14, 2.35]), + n=np.array([np.nan, np.nan, 10.0, 10.0]), + base_period="universal", + ) + with pytest.raises(ValueError, match="multiple reference rows") as exc_info: + compute_honest_did(surface, M=0.5) + # Message-level pin: the native-results route fails its own + # consecutive-grid validation on the same gapped layout, so the + # error must NOT recommend it - it recommends re-estimation on a + # consecutive grid instead. + msg = str(exc_info.value) + assert "native" not in msg + assert "consecutive" in msg and "Re-estimate" in msg + + +# --------------------------------------------------------------------------- # +# Plotting (no source guard: label-faithful for any producer) +# --------------------------------------------------------------------------- # + + +class TestPlotting: + @pytest.fixture(autouse=True) + def _agg_backend(self): + matplotlib = pytest.importorskip("matplotlib") + matplotlib.use("Agg") + yield + import matplotlib.pyplot as plt + + plt.close("all") + + def test_cs_container_plots(self, cs_universal): + from diff_diff.visualization import plot_event_study + + surface = cs_universal.aggregate("event_study") + ax = plot_event_study(surface, show=False) + assert ax is not None + + def test_dcdh_l1_container_plots(self): + from diff_diff.visualization import plot_event_study + + surface = _dcdh_container() + ax = plot_event_study(surface, show=False) + assert ax is not None + + @staticmethod + def _multi_ref_container(): + return _tiny_container( + is_reference=np.array([True, True, False, False]), + att=np.array([0.0, 0.0, 1.9, 2.1]), + se=np.array([np.nan, np.nan, 0.12, 0.13]), + t_stat=np.array([np.nan, np.nan, 15.8, 16.2]), + p_value=np.array([np.nan, np.nan, 0.0, 0.0]), + conf_int_lower=np.array([np.nan, np.nan, 1.66, 1.85]), + conf_int_upper=np.array([np.nan, np.nan, 2.14, 2.35]), + n=np.array([np.nan, np.nan, 10.0, 10.0]), + ) + + @staticmethod + def _hollow_marker_xs(ax): + """Positional x of single-point markers drawn hollow (white face).""" + out = set() + for line in ax.lines: + xs = line.get_xdata() + if len(xs) == 1 and line.get_markerfacecolor() == "white": + out.add(float(xs[0])) + return out + + def test_multi_reference_rows_render_hollow(self): + # Several is_reference rows (the CS gapped-grid universal case) are + # carried ROW-ALIGNED via reference_marks: every normalization + # anchor renders hollow at 0 - never silently dropped, never + # presented as a filled estimate. The scalar reference stays None. + from diff_diff.visualization import plot_event_study + from diff_diff.visualization._event_study import _extract_plot_data + + surface = self._multi_ref_container() + out = _extract_plot_data(surface, None, None, None, None) + effects, _, periods, _, _, ref, ref_inferred, *_rest, marks = out + assert ref is None and ref_inferred is False + assert periods == [-2, -1, 0, 1] + assert effects[-2] == 0.0 and effects[-1] == 0.0 + assert marks == {-2, -1} + + ax = plot_event_study(surface, show=False) + # Periods -2, -1 sit at positional x 0, 1. + assert self._hollow_marker_xs(ax) == {0.0, 1.0} + + def test_multi_reference_explicit_renormalization(self): + # Re-basing a multi-reference surface around a period that is NOT + # one of its anchors is undefined (each anchor constrains its own + # cohort base) -> fail closed. Choosing one of the anchors is a + # no-op shift and keeps every anchor hollow. + from diff_diff.visualization import plot_event_study + + surface = self._multi_ref_container() + with pytest.raises(ValueError, match="multiple reference rows"): + plot_event_study(surface, reference_period=1, show=False) + ax = plot_event_study(surface, reference_period=-1, show=False) + assert self._hollow_marker_xs(ax) == {0.0, 1.0} + + def test_anticipation_window_shaded_as_post(self, panel): + # REGISTRY anticipation contract, mirrored from the HonestDiD/ + # pretrends boundary: with anticipation=k the window [e=-k, -1] + # carries anticipated TREATMENT effects, so the pre-treatment + # shading covers e < -k only - on the container route AND the + # native fit-time dict route. + from diff_diff.visualization._event_study import _extract_plot_data + + res = _fit_cs(panel, anticipation=1) + surface = res.aggregate("event_study") + for target in (surface, res): + out = _extract_plot_data(target, None, None, None, None) + pre, post = out[3], out[4] + assert -1 in post and -1 not in pre + assert all(p < -1 for p in pre) + + def test_plotly_off_center_intervals_render_endpoints(self): + # Twin of the matplotlib endpoint-based pin: the plotly renderers + # draw the stored interval endpoints verbatim (filled band / + # segment traces, no estimate-centered arithmetic), so intervals + # wholly above the estimate render exactly and never raise. + plotly = pytest.importorskip("plotly") + assert plotly is not None + from diff_diff.visualization import plot_event_study, plot_honest_event_study + + surface = _tiny_container( + # wholly ABOVE the estimates + conf_int_lower=np.array([0.3, np.nan, 2.2, 2.5]), + conf_int_upper=np.array([0.5, np.nan, 2.6, 2.9]), + base_period="universal", + reference_event_times=(-1,), + ) + fig = plot_event_study(surface, show=False, backend="plotly") + band_y = [round(float(v), 6) for v in fig.data[0].y] + # upper endpoints forward, lower endpoints reversed - verbatim. + assert band_y == [0.5, 2.6, 2.9, 2.5, 2.2, 0.3] + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + h = compute_honest_did(surface, M=0.5) + h.event_study_bounds = { + int(t): {"ci_lb": float(a) + 0.2, "ci_ub": float(a) + 0.6} + for t, a, r in zip(surface.event_time, surface.att, surface.is_reference) + if not r + } + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + fig2 = plot_honest_event_study(h, show=False, backend="plotly") + assert fig2 is not None + + def test_off_center_stored_intervals_render(self): + # Percentile/bootstrap intervals need not contain the point + # estimate; estimate-centered yerr went negative and crashed + # matplotlib. Endpoint-based bars draw exactly [lower, upper]. + from diff_diff.visualization import plot_event_study + + surface = _tiny_container( + # interval wholly ABOVE the estimate at -2, wholly BELOW at 1 + conf_int_lower=np.array([0.3, np.nan, 1.66, 1.5]), + conf_int_upper=np.array([0.5, np.nan, 2.14, 1.9]), + ) + ax = plot_event_study(surface, show=False) + assert self._interval_set(ax) == {(0.3, 0.5), (1.66, 2.14), (1.5, 1.9)} + + def test_container_all_post_rows_undefined_fails_closed(self): + # Pin: an all-invalid post block never reaches zero-dimensional + # optimization - HonestDiD.fit's centralized num_post check + # rejects it with a clear message (verified for the container + # route; the MPD twin is pinned in test_honest_did.py). + surface = _tiny_container( + se=np.array([0.1, np.nan, 0.0, np.nan]), + t_stat=np.array([1.0, np.nan, np.nan, np.nan]), + p_value=np.array([0.3, np.nan, np.nan, np.nan]), + conf_int_lower=np.array([-0.1, np.nan, np.nan, np.nan]), + conf_int_upper=np.array([0.3, np.nan, np.nan, np.nan]), + base_period="universal", + reference_event_times=(-1,), + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match="No post-period effects"): + compute_honest_did(surface, M=0.5) + + def test_honest_plot_skips_rows_excluded_by_honest(self): + # A zero-SE TRAILING row survives the consecutive-grid check but + # is excluded from beta_hat; the honest plotter must not paint it + # with the aggregate honest interval (scalar-bounds path), must + # accept per-period bounds keyed on the RETAINED set, and must + # reject an explicit request for the excluded row. + from diff_diff.visualization import plot_honest_event_study + + surface = _tiny_container( + se=np.array([0.1, np.nan, 0.12, 0.0]), + t_stat=np.array([1.0, np.nan, 15.8, np.nan]), + p_value=np.array([0.3, np.nan, 0.0, np.nan]), + conf_int_lower=np.array([-0.1, np.nan, 1.66, np.nan]), + conf_int_upper=np.array([0.3, np.nan, 2.14, np.nan]), + base_period="universal", + reference_event_times=(-1,), + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + h = compute_honest_did(surface, M=0.5) + assert h.post_periods_used == [0] # row 1 excluded (zero SE) + + # Scalar-bounds path: only retained rows are plotted. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + ax = plot_honest_event_study(h, show=False) + tick_labels = [t.get_text() for t in ax.get_xticklabels()] + assert "1" not in tick_labels and "-2" in tick_labels and "0" in tick_labels + + # Per-period bounds keyed on the retained set render fine. + h.event_study_bounds = {-2: {"ci_lb": -0.4, "ci_ub": 0.6}, 0: {"ci_lb": 1.6, "ci_ub": 2.4}} + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + assert plot_honest_event_study(h, show=False) is not None + + # An explicit request for the excluded row fails with a clear + # message, never a KeyError or a fabricated honest interval. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with pytest.raises(ValueError, match="not retained by HonestDiD"): + plot_honest_event_study(h, periods=[-2, 0, 1], show=False) + + def test_honest_plot_reference_anchor(self, cs_universal): + # The container's single reference renders as a hollow + # normalization anchor at 0 with NO standard or honest interval: + # inferred by default, accepted in explicit periods=, on both + # backends. + from diff_diff.visualization import plot_honest_event_study + + surface = cs_universal.aggregate("event_study") + ref = surface.reference_period + h = compute_honest_did(surface, M=0.5) + h.event_study_bounds = { + int(t): {"ci_lb": float(a) - 0.5, "ci_ub": float(a) + 0.5} + for t, a, r in zip(surface.event_time, surface.att, surface.is_reference) + if not r + } + ax = plot_honest_event_study(h, show=False) + labels = [t.get_text() for t in ax.get_xticklabels()] + assert str(ref) in labels + ref_x = float(labels.index(str(ref))) + hollow = [ + line + for line in ax.lines + if len(line.get_xdata()) == 1 and line.get_markerfacecolor() == "white" + ] + assert len(hollow) == 1 + assert float(hollow[0].get_xdata()[0]) == ref_x + # No interval segment (standard or honest) at the reference. + for c in ax.containers: + if hasattr(c, "lines"): + for lc in c.lines[2]: + for s in lc.get_segments(): + if len(s) >= 2: + assert abs(s[0][0] - ref_x) > 1e-9 + # Explicit periods= including the reference are accepted. + all_periods = sorted(int(t) for t in surface.event_time) + assert plot_honest_event_study(h, periods=all_periods, show=False) is not None + plotly = pytest.importorskip("plotly") + assert plotly is not None + fig = plot_honest_event_study(h, show=False, backend="plotly") + assert fig is not None + + def test_honest_plot_off_center_intervals_render(self, cs_universal): + # Honest bounds need not bracket the effects either; the honest + # plotter draws endpoint-based bars for both interval layers. + from diff_diff.visualization import plot_honest_event_study + + surface = cs_universal.aggregate("event_study") + h = compute_honest_did(surface, M=0.5) + h.event_study_bounds = { + int(t): {"ci_lb": float(a) + 0.2, "ci_ub": float(a) + 0.6} + for t, a, r in zip(surface.event_time, surface.att, surface.is_reference) + if not r + } + fig = plot_honest_event_study(h, show=False) + assert fig is not None + + def test_container_alpha_mismatch_warns(self): + # Stored container intervals are at the FIT's level; the plot's + # ``alpha`` cannot re-level bootstrap/t-based intervals from the + # SE, so a mismatch is named, never silently relabeled. + from diff_diff.visualization import plot_event_study + + surface = _tiny_container(alpha=0.10) + with pytest.warns(UserWarning, match="does not apply to an EventStudyResults"): + plot_event_study(surface, show=False) + + def test_container_alpha_match_silent(self): + from diff_diff.visualization import plot_event_study + + surface = _tiny_container() # alpha=0.05 == plot default + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + plot_event_study(surface, show=False) + assert not [w for w in caught if "does not apply" in str(w.message)] + + def test_explicit_normalization_recomputes_at_requested_alpha(self): + # The explicit-reference path discards the stored overrides and + # recomputes pointwise intervals at the requested alpha, so no + # stored-level mismatch remains to warn about. + from diff_diff.visualization import plot_event_study + + surface = _tiny_container(alpha=0.10) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + plot_event_study(surface, reference_period=-1, show=False) + assert not [w for w in caught if "does not apply" in str(w.message)] + + def test_honest_plot_alpha_mismatch_warns(self, cs_universal): + from diff_diff.visualization import plot_honest_event_study + + surface = cs_universal.aggregate("event_study") # alpha=0.05 + h = compute_honest_did(surface, M=0.5, alpha=0.10) + # Bracketing per-period bounds (the established fixture pattern) so + # the honest-CI overlay renders. + h.event_study_bounds = { + int(t): {"ci_lb": float(a) - 0.5, "ci_ub": float(a) + 0.5} + for t, a, r in zip(surface.event_time, surface.att, surface.is_reference) + if not r + } + with pytest.warns(UserWarning, match="container's stored intervals at"): + plot_honest_event_study(h, show=False) + + def test_gapped_cs_container_multi_reference_plot(self): + # End-to-end on a REAL CS universal fit over a gapped grid + # ({1,3,6}, cohorts {3,6}): aggregate("event_study") carries >=2 + # reference rows; automatic plotting hollows every anchor and + # explicit renormalization around a non-anchor fails closed. + from diff_diff.visualization import plot_event_study + + rng = np.random.RandomState(3) + rows = [] + for u in range(40): + g = 3 if u < 15 else (6 if u < 30 else 0) + ufe = rng.randn() * 2 + for t in (1, 3, 6): + post = 1 if (g > 0 and t >= g) else 0 + rows.append( + { + "unit": u, + "period": t, + "outcome": ufe + 0.3 * t + 2.0 * post + rng.randn() * 0.5, + "first_treat": g, + } + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + cs = CallawaySantAnna(n_bootstrap=0, base_period="universal").fit( + pd.DataFrame(rows), + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + ) + surface = cs.aggregate("event_study") + keys = surface.event_time.tolist() + refs = sorted(surface.event_time[surface.is_reference].tolist()) + assert len(refs) >= 2 + + ax = plot_event_study(surface, show=False) + assert self._hollow_marker_xs(ax) == {float(keys.index(r)) for r in refs} + with pytest.raises(ValueError, match="multiple reference rows"): + plot_event_study(surface, reference_period=0, show=False) + + def test_explicit_normalization_discards_stale_bands(self): + # REGISTRY (Event Study Plotting): explicit reference_period= + # normalization recomputes CIs from normalized effects + original + # SEs and NaNs the reference CI - simultaneous bands computed + # around the UN-normalized effects must be discarded, not drawn. + import matplotlib.pyplot as plt + + from diff_diff.visualization import plot_event_study + + surface = _tiny_container( + att=np.array([0.1, 0.05, 2.0, 2.1]), + se=np.array([0.1, 0.1, 0.12, 0.13]), + t_stat=np.array([1.0, 0.5, 15.0, 16.0]), + p_value=np.array([0.3, 0.6, 0.0, 0.0]), + conf_int_lower=np.array([-0.1, -0.15, 1.76, 1.85]), + conf_int_upper=np.array([0.3, 0.25, 2.24, 2.35]), + is_reference=np.array([False, False, False, False]), + n=np.array([10.0, 10.0, 10.0, 10.0]), + cband_lower=np.array([-0.15, -0.2, 1.7, 1.8]), + cband_upper=np.array([0.35, 0.3, 2.3, 2.4]), + ) + ax = plot_event_study(surface, reference_period=-1, show=False) + # Recovered plotted CIs: reference row has NO error bar (NaN CI); + # other rows are recentered pointwise intervals, not the stale + # bands (which bracket the un-normalized effects). Periods map to + # POSITIONAL x (period -1 sits at x=1 for the [-2,-1,0,1] order). + ref_x = float([-2, -1, 0, 1].index(-1)) + containers = [c for c in ax.containers if hasattr(c, "lines")] + assert containers, "expected errorbar containers" + segs = [seg for c in containers for lc in c.lines[2] for seg in lc.get_segments()] + assert segs, "expected drawn interval segments" + for seg in segs: + (x0, y0), (_, y1) = seg + if abs(x0 - ref_x) < 1e-9: + raise AssertionError( + f"reference row at positional x={ref_x} must carry no " + f"error bar, got {(y0, y1)}" + ) + # The drawn intervals are the RECOMPUTED pointwise ones around the + # normalized effects (shifted by the ref effect 0.05), never the + # stale bands. + from scipy import stats as scipy_stats + + z = scipy_stats.norm.ppf(0.975) + expected = set() + for eff, se_v in ((0.1, 0.1), (2.0, 0.12), (2.1, 0.13)): + norm_eff = eff - 0.05 + expected.add((round(norm_eff - z * se_v, 6), round(norm_eff + z * se_v, 6))) + assert self._interval_set(ax) == expected + plt.close("all") + + def test_calendar_container_positional_split(self): + # Calendar labels (str/Timestamp) break numeric p<0 splitting; the + # container branch splits positionally around the reference row. + from diff_diff.visualization._event_study import _extract_plot_data + + surface = _tiny_container( + event_time=np.array(["2018", "2019", "2020", "2021"], dtype=object), + time_scale="calendar", + is_reference=np.array([False, True, False, False]), + ) + out = _extract_plot_data(surface, None, None, None, None) + _, _, _, pre, post, ref, _, _, _, _, _, _ = out + assert ref == "2019" + assert pre == ["2018"] + assert post == ["2020", "2021"] + + def test_stored_pointwise_intervals_preserved(self): + # The container carries the producer's ACTUAL intervals (bootstrap + # percentile / survey-t / Bell-McCaffrey differ from att +/- z*se); + # the plot branch must relay them, not reconstruct normal-style + # intervals from the SE. + from diff_diff.visualization._event_study import _extract_plot_data + + surface = _tiny_container( + # deliberately asymmetric, non +/-1.96*se intervals + conf_int_lower=np.array([-0.4, np.nan, 1.2, 1.4]), + conf_int_upper=np.array([0.2, np.nan, 2.2, 2.5]), + ) + out = _extract_plot_data(surface, None, None, None, None) + _, _, _, _, _, _, _, band_lo, band_hi, pw_lo, pw_hi, _ = out + # No cbands on this surface: the band channel falls back per-row to + # the stored interval; the pointwise channel is the stored interval + # verbatim. Both preserve the reference row's NaN (undefined + # inference is never recomputed from the SE). + for clo, chi in ((band_lo, band_hi), (pw_lo, pw_hi)): + assert clo is not None and chi is not None + assert clo[-2] == -0.4 and chi[-2] == 0.2 + assert clo[0] == 1.2 and chi[0] == 2.2 + assert np.isnan(clo[-1]) and np.isnan(chi[-1]) + + @staticmethod + def _vertical_segments(ax): + segs = [] + for c in ax.containers: + if not hasattr(c, "lines"): + continue + for lc in c.lines[2]: + segs.extend(lc.get_segments()) + return segs + + @classmethod + def _interval_set(cls, ax): + """Drawn vertical error intervals as a set of rounded (lo, hi).""" + out = set() + for s in cls._vertical_segments(ax): + if len(s) < 2: # NaN rows can yield empty segments + continue + ys = sorted((s[0][1], s[1][1])) + if np.isfinite(ys[0]) and np.isfinite(ys[1]): + out.add((round(ys[0], 6), round(ys[1], 6))) + return out + + def test_rendered_intervals_match_stored_ci(self): + # Survey-t / bootstrap-style stored intervals differ from + # att +/- 1.96*se; the RENDERED error bars must be the stored + # endpoints, not normal reconstructions. + from diff_diff.visualization import plot_event_study + + surface = _tiny_container( + conf_int_lower=np.array([-0.4, np.nan, 1.2, 1.4]), + conf_int_upper=np.array([0.2, np.nan, 2.2, 2.5]), + ) + ax = plot_event_study(surface, show=False) + drawn = self._interval_set(ax) + # every drawn interval is a STORED one; the z-reconstruction for + # period 0 (~1.665..2.135) must NOT appear + assert (1.2, 2.2) in drawn and (-0.4, 0.2) in drawn and (1.4, 2.5) in drawn + assert not any(abs(lo - 1.665) < 0.01 for lo, _ in drawn) + + def test_use_cband_false_keeps_stored_intervals(self): + # use_cband=False must select the POINTWISE stored channel, not + # clear all overrides into normal reconstruction (dCDH bootstrap + # percentile intervals are asymmetric). + from diff_diff.visualization import plot_event_study + + surface = _tiny_container( + conf_int_lower=np.array([-0.4, np.nan, 1.2, 1.4]), + conf_int_upper=np.array([0.2, np.nan, 2.2, 2.5]), + cband_lower=np.array([-0.6, np.nan, 1.0, 1.2]), + cband_upper=np.array([0.4, np.nan, 2.4, 2.7]), + ) + ax = plot_event_study(surface, use_cband=False, show=False) + drawn = self._interval_set(ax) + assert (1.2, 2.2) in drawn # stored pointwise, NOT the band + assert (1.0, 2.4) not in drawn + ax2 = plot_event_study(surface, use_cband=True, show=False) + drawn2 = self._interval_set(ax2) + assert (1.0, 2.4) in drawn2 # the band channel + assert (1.2, 2.2) not in drawn2 + + def test_zero_se_row_renders_without_interval(self): + # A non-reference row with se == 0 and NaN stored inference must + # NOT get a zero-width normal interval: its stored NaN bounds are + # preserved, so no error bar is drawn at that period. + from diff_diff.visualization import plot_event_study + + surface = _tiny_container( + se=np.array([0.0, np.nan, 0.12, 0.13]), + t_stat=np.array([np.nan, np.nan, 15.8, 16.2]), + p_value=np.array([np.nan, np.nan, 0.0, 0.0]), + conf_int_lower=np.array([np.nan, np.nan, 1.66, 1.85]), + conf_int_upper=np.array([np.nan, np.nan, 2.14, 2.35]), + ) + ax = plot_event_study(surface, show=False) + drawn = self._interval_set(ax) + # only the two rows with defined stored inference draw intervals; + # the zero-SE row contributes NO zero-width normal interval + assert drawn == {(1.66, 2.14), (1.85, 2.35)} + + def test_plottable_results_membership(self): + from diff_diff.visualization._event_study import PlottableResults + + assert "EventStudyResults" in str(PlottableResults) + + def test_plot_honest_event_study_container_route(self, cs_universal): + # Before this PR the honest plotter's own re-extraction probe fell + # through to TypeError("Cannot extract event study data from + # original_results") for container-route HonestDiD results. + from diff_diff.visualization._event_study import plot_honest_event_study + + surface = cs_universal.aggregate("event_study") + h = compute_honest_did(surface, M=0.5) + # event_study_bounds is populated by the per-period bounds path; + # attach a minimal dict (the established fixture pattern) so the + # plotter's honest-CI overlay has content. + h.event_study_bounds = { + int(t): {"ci_lb": float(a) - 0.5, "ci_ub": float(a) + 0.5} + for t, a, r in zip(surface.event_time, surface.att, surface.is_reference) + if not r + } + fig = plot_honest_event_study(h, show=False) + assert fig is not None + + def test_plot_honest_uses_stored_container_intervals(self): + # The honest plotter's "original" intervals on the container route + # must be the container's STORED intervals (survey-t / bootstrap + # inference), never z-reconstructions; stored NaN (undefined + # inference, e.g. the replicate df_survey=0 sentinel case) is + # preserved rather than replaced with finite bounds. + from diff_diff.visualization._event_study import plot_honest_event_study + + surface = _tiny_container( + conf_int_lower=np.array([np.nan, np.nan, 1.2, 1.4]), + conf_int_upper=np.array([np.nan, np.nan, 2.2, 2.5]), + base_period="universal", + ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + h = compute_honest_did(surface, M=0.5) + h.event_study_bounds = { + int(t): {"ci_lb": float(a) - 0.5, "ci_ub": float(a) + 0.5} + for t, a, r in zip(surface.event_time, surface.att, surface.is_reference) + if not r + } + ax = plot_honest_event_study(h, show=False) + drawn = TestPlotting._interval_set(ax) + # the stored (not z-reconstructed) original interval appears; the + # z-reconstruction for period 0 (~1.665..2.135) does not + assert (1.2, 2.2) in drawn, drawn + assert not any(abs(lo - 1.665) < 0.01 for lo, _ in drawn) diff --git a/tests/test_event_study_surface.py b/tests/test_event_study_surface.py index 8caa47ce8..32cdc4103 100644 --- a/tests/test_event_study_surface.py +++ b/tests/test_event_study_surface.py @@ -160,6 +160,102 @@ def test_to_dict_json_safe_period_labels(): assert d["reference_period"] == "2018" +# =========================================================================== +# Provenance fields (base_period / anticipation / df_survey) +# =========================================================================== + + +def test_provenance_fields_default_none(): + surface = _tiny_surface() + assert surface.base_period is None + assert surface.anticipation is None + assert surface.df_survey is None + + +def test_provenance_round_trips_through_to_dict(): + d = _tiny_surface(base_period="universal", anticipation=1, df_survey=7.0).to_dict() + assert d["base_period"] == "universal" + assert d["anticipation"] == 1 + assert d["df_survey"] == 7.0 + + +def test_builder_threads_cs_provenance(): + # The relative-dict builder reads base_period/anticipation off the + # producer and resolves the scalar df_survey (None here: no survey + # design and no bare-cluster df carrier). + class _FakeCS: + alpha = 0.05 + base_period = "universal" + anticipation = 1 + event_study_effects = { + -2: {"effect": 0.0, "se": np.nan, "n_groups": 0}, + 0: {"effect": 1.0, "se": 0.1, "n_groups": 5}, + } + + surface = build_event_study_surface(_FakeCS()) + assert surface.base_period == "universal" + assert surface.anticipation == 1 + assert surface.df_survey is None + + +def test_empty_surface_threads_provenance(): + # The requested-but-empty early return must carry provenance too - a + # balance_e-emptied aggregation would otherwise read as provenance-free. + class _EmptyCS: + alpha = 0.05 + base_period = "universal" + anticipation = 1 + event_study_effects: dict = {} + + surface = build_event_study_surface(_EmptyCS()) + assert surface.event_time.shape[0] == 0 + assert surface.base_period == "universal" + assert surface.anticipation == 1 + assert surface.df_survey is None + + +def test_df_survey_replicate_undefined_maps_to_zero_sentinel(): + # survey_metadata present, df_survey undefined, replicate design -> + # the 0.0 sentinel (fails closed to NaN critical values downstream). + class _SM: + df_survey = None + replicate_method = "brr" + + class _Fake: + alpha = 0.05 + survey_metadata = _SM() + event_study_effects = {0: {"effect": 1.0, "se": 0.1, "n_groups": 5}} + + surface = build_event_study_surface(_Fake()) + assert surface.df_survey == 0.0 + + +def test_df_survey_prefers_survey_metadata_over_df_inference(): + class _SM: + df_survey = 12 + replicate_method = None + + class _Fake: + alpha = 0.05 + survey_metadata = _SM() + df_inference = 30 + event_study_effects = {0: {"effect": 1.0, "se": 0.1, "n_groups": 5}} + + surface = build_event_study_surface(_Fake()) + assert surface.df_survey == 12.0 + + +def test_df_survey_bare_cluster_falls_back_to_df_inference(): + class _Fake: + alpha = 0.05 + survey_metadata = None + df_inference = 19 + event_study_effects = {0: {"effect": 1.0, "se": 0.1, "n_groups": 5}} + + surface = build_event_study_surface(_Fake()) + assert surface.df_survey == 19.0 + + # =========================================================================== # Producer builders (small analytical fits) # =========================================================================== diff --git a/tests/test_honest_did.py b/tests/test_honest_did.py index 02a16e2ad..b8d71706c 100644 --- a/tests/test_honest_did.py +++ b/tests/test_honest_did.py @@ -315,10 +315,92 @@ def test_extract_from_multiperiod(self, mock_multiperiod_results): pe = mock_multiperiod_results.period_effects[period] assert sigma[i, i] == pytest.approx(pe.se**2, abs=1e-10) + def test_zero_se_period_dropped_from_multiperiod(self, mock_multiperiod_results): + """se == 0 rows drop on the MPD branch (undefined inference). + + safe_inference treats se <= 0 as undefined inference; admitting + such a row would enter Sigma with zero variance and launder NaN + source inference into finite sensitivity bounds. Mirrors the CS, + dCDH and container branches and pretrends. + """ + import dataclasses + + pe0 = mock_multiperiod_results.period_effects[0] + mock_multiperiod_results.period_effects[0] = dataclasses.replace(pe0, se=0.0) + beta_hat, sigma, num_pre, num_post, pre_periods, post_periods, _df = ( + _extract_event_study_params(mock_multiperiod_results) + ) + assert num_pre == 2 and num_post == 4 + assert len(beta_hat) == 6 + assert sigma.shape == (6, 6) + # First retained row is period 1 (period 0 dropped): its variance + # leads the sub-VCV diagonal. + assert sigma[0, 0] == pytest.approx(0.35**2, abs=1e-12) + # Label/index contract: the returned lists are the ESTIMATED + # horizons beta_hat/sigma were built from - the dropped zero-SE + # period and the reference are absent, and lengths match the + # counts. + assert pre_periods == [1, 2] + assert post_periods == [4, 5, 6, 7] + assert len(pre_periods) == num_pre and len(post_periods) == num_post + # End to end: the metadata a user sees relays the same lists. + h = compute_honest_did(mock_multiperiod_results, M=0.5) + assert h.pre_periods_used == [1, 2] + assert h.post_periods_used == [4, 5, 6, 7] + + @pytest.mark.parametrize( + "bad_period", + [1, 2, 4, 5], + ids=["interior-pre", "ref-adjacent-pre", "first-post", "interior-post"], + ) + def test_zero_se_breaking_grid_geometry_fails_closed( + self, mock_multiperiod_results, bad_period + ): + """Interior / reference-adjacent zero-SE drops fail closed. + + The RR constraint builders index retained coefficients + positionally, so dropping an interior or reference-adjacent + horizon would silently treat non-adjacent periods as consecutive + and return wrong bounds. Only leading-pre / trailing-post drops + keep valid geometry (previous test). + """ + import dataclasses + + pe = mock_multiperiod_results.period_effects[bad_period] + mock_multiperiod_results.period_effects[bad_period] = dataclasses.replace(pe, se=0.0) + with pytest.raises(ValueError, match="consecutive estimated horizons"): + _extract_event_study_params(mock_multiperiod_results) + + def test_all_pre_periods_zero_se_fails_closed(self, mock_multiperiod_results): + """Pin: an all-invalid pre block is rejected with a clear message + (never a zero-dimensional restriction build).""" + import dataclasses + + for p in (0, 1, 2): + pe = mock_multiperiod_results.period_effects[p] + mock_multiperiod_results.period_effects[p] = dataclasses.replace(pe, se=0.0) + with pytest.raises(ValueError, match="No pre-period effects"): + _extract_event_study_params(mock_multiperiod_results) + + def test_zero_se_trailing_post_dropped_ok(self, mock_multiperiod_results): + """A trailing post-period drop keeps valid positional geometry.""" + import dataclasses + + pe = mock_multiperiod_results.period_effects[7] + mock_multiperiod_results.period_effects[7] = dataclasses.replace(pe, se=0.0) + beta_hat, sigma, num_pre, num_post, pre_p, post_p, _df = _extract_event_study_params( + mock_multiperiod_results + ) + assert num_pre == 3 and num_post == 3 + assert pre_p == [0, 1, 2] and post_p == [4, 5, 6] + assert len(beta_hat) == 6 and sigma.shape == (6, 6) + def test_extract_unsupported_type_raises(self): - """Test that unsupported types raise TypeError.""" - with pytest.raises(TypeError, match="Unsupported results type"): + """Test that unsupported types raise TypeError, naming the container route.""" + with pytest.raises(TypeError, match="Unsupported results type") as exc_info: _extract_event_study_params("not a results object") + # The expected-types list names the post-fit container route too. + assert "EventStudyResults" in str(exc_info.value) # ============================================================================= @@ -1481,6 +1563,22 @@ def test_dcdh_interior_gap_triggers_trimming_warning(self): # Retained pre should be [-1] only (h=-3 dropped due to gap at -2) assert bounds.pre_periods_used == [-1] + def test_dcdh_zero_se_treated_like_nan(self): + """se == 0 placebo rows drop exactly like NaN-SE rows (undefined + inference; a zero row would enter Sigma with zero variance).""" + import warnings + + results = self._fit_dcdh(n_periods=8, L_max=3) + results.placebo_event_study[-2]["se"] = 0.0 + + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + bounds = compute_honest_did(results) + trim_warns = [x for x in w if "dropping non-consecutive" in str(x.message).lower()] + assert len(trim_warns) >= 1 + # Same outcome as the NaN-SE interior-gap test: only [-1] retained. + assert bounds.pre_periods_used == [-1] + def test_dcdh_missing_boundary_minus1_raises(self): """ValueError when horizon -1 has NaN SE (boundary required).""" import warnings diff --git a/tests/test_pretrends.py b/tests/test_pretrends.py index ce2225766..f7cf74df2 100644 --- a/tests/test_pretrends.py +++ b/tests/test_pretrends.py @@ -863,8 +863,10 @@ def test_unsupported_results_type_raises(self): """Test that unsupported results type raises TypeError.""" pt = PreTrendsPower() - with pytest.raises(TypeError, match="Unsupported results type"): + with pytest.raises(TypeError, match="Unsupported results type") as exc_info: pt.fit("not a results object") + # The expected-types list names the post-fit container route too. + assert "EventStudyResults" in str(exc_info.value) def test_callaway_santanna_universal_base_period(self): """Test that reference period (e=-1) is correctly filtered out with universal base period. diff --git a/tests/test_visualization.py b/tests/test_visualization.py index a6ef132a3..475832e87 100644 --- a/tests/test_visualization.py +++ b/tests/test_visualization.py @@ -280,8 +280,10 @@ def test_error_missing_dataframe_columns(self): def test_error_invalid_results_type(self): """Test error with invalid results type.""" pytest.importorskip("matplotlib") - with pytest.raises(TypeError, match="Cannot extract plot data"): + with pytest.raises(TypeError, match="Cannot extract plot data") as exc_info: plot_event_study("invalid") + # The expected-types list names the unified container too. + assert "EventStudyResults" in str(exc_info.value) def test_plot_with_nan_se_reference_period(self): """Test that reference period with NaN SE is plotted without error bars. @@ -671,7 +673,7 @@ def test_plot_uses_cband_cis_by_default(self, cs_cband_results): assert ax is not None # Verify cband CIs are extracted - _, _, _, _, _, _, _, ci_lower_override, ci_upper_override = _extract_plot_data( + _, _, _, _, _, _, _, ci_lower_override, ci_upper_override, _, _, _ = _extract_plot_data( cs_cband_results, None, None, None, None ) assert ci_lower_override is not None