From dbd7363083529973d94af59e95c7d8138ed38bb1 Mon Sep 17 00:00:00 2001 From: igerber Date: Mon, 3 Aug 2026 15:07:26 -0400 Subject: [PATCH 1/4] feat(v4): StackedDiD post-fit aggregate() + honest/pretrends container admission (2(b) PR-2, M-024) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The event-study surface is always computed at fit (both fit-time gates removed, incl. per-event hc2_bm BM DOFs); fit(aggregate=) deprecates via the NOT_SUPPLIED sentinel (warn-and-still-work; value validation kept). StackedDiDResults adopts AggregationMixin: aggregate('event_study') is a pure builder view with base_period="universal" + singleton reference_event_times provenance; aggregate('simple') relays the stored overall inference bit-exactly (target="att", n=n_treated_units/units, df=inference_df). compute_honest_did/compute_pretrends_power admit Stacked-sourced containers (kappa_pre >= 2; honest needs a non-singular pre-covariance); the seven container guard messages derive the producer from surface.source with producer-conditional remedies; withheld-inference rows admit with a source-scoped warning in both consumers. EventStudyResults.__post_init__ now copies vcov/vcov_index (cross-producer de-aliasing; int index dtype preserved). DiagnosticReport's heterogeneity + parallel_trends (kappa_pre >= 2) checks now run on plain fits with a per-type PT skip remediation; the practitioner balance step moves off the colliding "heterogeneity" step key. Ledger: M-024 -> shimmed (phase 5), M-093 second pre-cut amendment; v4-design §5/§6 amended (view-relay exception + container copy contract); REGISTRY M-024/admission Notes; variance-conventions tail-df pin (309,)*3 + regenerated table; guides/ README/API/choosing/notebook migrated to the post-fit route. --- CHANGELOG.md | 74 +++- README.md | 4 +- TODO.md | 3 +- benchmarks/python/benchmark_stacked_did.py | 6 +- diff_diff/diagnostic_report.py | 20 +- diff_diff/guides/llms-full.txt | 9 +- diff_diff/guides/llms-practitioner.txt | 4 +- diff_diff/guides/llms.txt | 6 +- diff_diff/honest_did.py | 178 +++++--- diff_diff/practitioner.py | 10 +- diff_diff/pretrends.py | 106 +++-- diff_diff/results_base.py | 17 +- diff_diff/stacked_did.py | 256 +++++++----- diff_diff/stacked_did_results.py | 106 ++++- .../diff_diff.StackedDiDResults.rst | 3 + docs/api/stacked_did.rst | 9 +- docs/choosing_estimator.rst | 6 +- docs/methodology/REGISTRY.md | 8 +- docs/methodology/REPORTING.md | 21 +- docs/methodology/variance-conventions.md | 2 +- docs/tutorials/13_stacked_did.ipynb | 8 +- docs/v4-deprecations.yaml | 11 +- docs/v4-design.md | 18 +- tests/test_aggregate_contract.py | 351 +++++++++++++++- tests/test_business_report.py | 47 +++ tests/test_diagnostic_report.py | 66 ++- tests/test_event_study_consumers.py | 390 +++++++++++++++++- tests/test_event_study_surface.py | 4 - tests/test_guides.py | 15 + tests/test_methodology_stacked_did.py | 10 - tests/test_practitioner.py | 23 ++ tests/test_stacked_did.py | 159 ++++--- tests/test_survey_estimator_validation.py | 1 - tests/test_variance_conventions.py | 6 +- 34 files changed, 1631 insertions(+), 326 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 26aa08f49..d80eb7d4f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,70 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- **StackedDiD post-fit `aggregate()` + honest/pretrends container + admission** (v4 program 2(b) PR-2; ledger row [M-024] flips to + shimmed, row [M-093] amended a second time). + - The StackedDiD event-study surface is now ALWAYS computed at fit: + the pooled regression always included the event-time interactions, + so only the extraction and (under `hc2_bm`) the per-event + Bell-McCaffrey contrast DOFs were gated on fit-time `aggregate=` - + both gates are gone. Consequences on plain fits: `summary()` + renders the event-study table, `to_dataframe('event_study')` is + total, `plot_event_study` renders (previously a `TypeError`), + `event_study_effects`/`event_study_vcov`/`event_study_df` are + always populated, and DiagnosticReport's `heterogeneity` (always) + and `parallel_trends` (`kappa_pre >= 2`) checks now run - which + `BusinessReport(auto_diagnostics=True)` surfaces too. The `hc2_bm` + overall-ATT contrast now shares the multi-contrast DOF batch, + whose degeneracy guard is batch-relative: the noise-floor + carve-out on the R `Wald_test` HTZ parity claim widens + quantitatively (fail-closed NaN + a loud warning when tripped - + behavior unchanged, trigger condition wider; REGISTRY M-024 Note). + - `StackedDiD.fit(aggregate=)` is DEPRECATED ([M-024]): supplying it + (any value, `None` included) warns; the parameter is behaviorally + inert (CS-style warn-and-still-work - results are identical to a + plain fit's), and the `group`/`all`/unknown-value `ValueError`s + persist after the warning. The `stacked_did()` convenience wrapper + forwards the shim. `StackedDiDResults.aggregate()` ships as a pure + VIEW: `'event_study'` returns the unified container (with + `base_period='universal'` and singleton `reference_event_times` + provenance - one omitted reference per fit by construction); + `'simple'` relays `overall_att/se/t/p/CI` bit-exactly with + `target='att'` (the post-period average, not the per-event trimmed + aggregate ATT), `n = n_treated_units` (`'units'`; treated-only + scope - the treated and clean-control unit sets overlap across + sub-experiments, so no disjoint total exists) and + `df = inference_df`. Every stored inference mode relays (survey + TSL, replicate refit, `hc2_bm` fail-close); `group`/`calendar` + fail closed, `balance_e` applies to no level, `weights=` rejected. + - `compute_honest_did` and `compute_pretrends_power` now ADMIT + Stacked-sourced containers ([M-093] second pre-cut amendment): + requires `kappa_pre >= 2` (the default grid has no estimated + pre-periods) and, for HonestDiD, a non-singular pre-period + covariance (keep `kappa_pre` small relative to the cluster count). + Analytical Stacked containers carry `df_survey=None` - honest FLCI + critical values are normal-theory there (deliberate; REGISTRY + Note). Rows whose per-row inference the producer withheld (finite + `se`, non-finite `p_value` - the `hc2_bm` BM-DOF fail-close and + replicate-undefined designs) are admitted with a source-scoped + warning in BOTH consumers: bounds/power consume only the point + estimates and covariance, which remain valid (at the replicate + `0.0` sentinel, honest's identified-set bounds stay finite while + its FLCI CI endpoints are NaN). The seven producer-derived guard + messages in both container branches now derive the producer from + `surface.source` with producer-conditional remedies. + - Cross-producer container hardening: `EventStudyResults`'s + `__post_init__` now COPIES `vcov`/`vcov_index` like every other + array field (`np.asarray` aliased the producer's stored matrix on + the post-fit view route, so mutating a container could corrupt the + fitted result; `vcov_index` keeps its native dtype - int labels + never become floats). + - The practitioner "Check sub-experiment balance" step moved off the + shared `step_name='heterogeneity'` key (now + `'sub_experiment_balance'`): with the heterogeneity check running + on every plain fit, the collision silently dropped that unrelated + advice from `next_steps`. Sibling producer-side collisions on + other estimators are recorded in TODO.md. - **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). @@ -20,11 +84,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 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 + HonestDiD/PreTrendsPower is SOURCE-SCOPED (initially CS-only; + widened to Stacked-sourced containers with [M-024] below): 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 diff --git a/README.md b/README.md index 2aa3e054d..c42568f7e 100644 --- a/README.md +++ b/README.md @@ -90,8 +90,8 @@ For rigorous DiD analysis, follow these 8 steps. Skipping diagnostic steps produ 3. **Test parallel trends** - simple 2x2: `check_parallel_trends()`, `equivalence_test_trends()`; staggered: inspect CS event-study pre-period coefficients (generic PT tests are invalid for staggered designs). Insignificant pre-trends do NOT prove PT holds. 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')`/`'event_study'` (post-fit, no refit); SA: `results.event_study_effects` / `to_dataframe(level='cohort')`; subgroup re-estimation. +6. **Sensitivity analysis** - `compute_honest_did(results)` for bounds under PT violations (MultiPeriodDiD, CS, or dCDH natively; a StackedDiD `results.aggregate('event_study')` container also admits - needs `kappa_pre >= 2`), `run_all_placebo_tests()` for 2x2 falsification, specification comparisons for staggered designs. +7. **Heterogeneity** - CS: `results.aggregate('group')`/`'event_study'` (post-fit, no refit); SA: `results.event_study_effects` / `to_dataframe(level='cohort')`; Stacked: `results.aggregate('event_study')`/`'simple'` post-fit views (surface always computed since 3.9); subgroup re-estimation. 8. **Robustness** - compare 2-3 estimators (CS vs SA vs BJS), report with and without covariates (shows whether conditioning drives identification), present pre-trends and sensitivity bounds. Full guide: `diff_diff.get_llm_guide("practitioner")`. diff --git a/TODO.md b/TODO.md index 301e30d29..5d74bb325 100644 --- a/TODO.md +++ b/TODO.md @@ -21,6 +21,7 @@ Related tracking surfaces: | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| +| practitioner `step_name="heterogeneity"` producer-side collisions: three OTHER estimators' advice steps reuse the key with non-heterogeneity labels (`:975` ContinuousDiD dose-response, `:1022` Triple placebo-group, `:1413` LPDiD WAS arrays), so DiagnosticReport's heterogeneity completion silently drops that unrelated advice from `next_steps` via `_filter_steps` - the same latent collision fixed for StackedDiD in M-024 (renamed to `sub_experiment_balance`). Renaming these changes those estimators' report output; audit + rename with per-estimator pins. | `diff_diff/practitioner.py` | 2(b) PR-2 review R9 | Quick | Low | | 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 | @@ -65,7 +66,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | MMM interop PR-B: calibration tutorial notebook (fit DiD/CS -> scope -> `to_pymc_marketing_lift_test` / `to_meridian_roi_prior`) + a `llms-practitioner.txt` Step 8 pointer to the exporters as the MMM hand-off. | `docs/tutorials/`, `diff_diff/guides/llms-practitioner.txt` | mmm-interop | Mid | Low | | Tracking-file contract guard test: reject NEW active deferred-work pointers at `TODO.md` (deferred rows live in `DEFERRED.md`; allowlist for historical/past-tense prose and actionable-row pointers) and assert rows cross-linking a `docs/v4-deprecations.yaml` `M-xxx` id don't restate ledger status. Origin: tracking-split local review R2. | `tests/`, `TODO.md`, `DEFERRED.md` | tracking-split | Quick | Low | | Real-data CI canary for dataset-backed replication tests: `test_methodology_lwdid.py`'s Prop 99 / Walmart goldens skip (visibly) when loaders fall back to synthetic; add a lane or canary asserting `df.attrs["source"] == "lwdid_ssc_ancillary"` in CI so network regressions cannot silently de-gate the replication tests. Follow-on from the loader-fallback repair (#723), which made provenance explicit but deliberately did not add a network-dependent CI lane. | `tests/test_methodology_lwdid.py`, `.github/workflows/` | LWDiD validation suite | Quick | Low | -| df-provenance completion follow-up (PR C descoping): a Wooldridge public per-row df channel (a per-key dict that can also carry the hc2_bm BM contrast dofs the scalar `_df_analytic_fallback` cannot represent — the unified surface's Wooldridge df column stays NaN until then); ImputationDiD lead-horizon `event_study_df` (needs the `_compute_lead_coefficients`/`_aggregate_event_study` return-path change plus bootstrap-clearing semantics for post rows); and the `inference_df` (DiD/MPD/SA/StackedDiD) vs `df_inference` (CS/dCDH) adapter-naming unification. | `diff_diff/wooldridge_results.py`, `diff_diff/imputation.py`, `diff_diff/results_base.py` | #variance-PR-C | Mid | Low | +| df-provenance completion follow-up (PR C descoping): a Wooldridge public per-row df channel (a per-key dict that can also carry the hc2_bm BM contrast dofs the scalar `_df_analytic_fallback` cannot represent — the unified surface's Wooldridge df column stays NaN until then); ImputationDiD lead-horizon `event_study_df` (needs the `_compute_lead_coefficients`/`_aggregate_event_study` return-path change plus bootstrap-clearing semantics for post rows); and the `inference_df` (DiD/MPD/SA/StackedDiD) vs `df_inference` (CS/dCDH) adapter-naming unification. NOTE (M-024): StackedDiD's `df_survey=None` on analytical fits is DELIBERATE - the container resolver's `df_inference` fallback not matching `inference_df` is load-bearing (normal-theory FLCI; see the REGISTRY StackedDiD M-024 Note and the pin in `tests/test_event_study_consumers.py`) - the unification must preserve that outcome or renegotiate it explicitly. | `diff_diff/wooldridge_results.py`, `diff_diff/imputation.py`, `diff_diff/results_base.py` | #variance-PR-C | Mid | Low | | `worktree-rm` safety via a tested argv helper: the prose rewrite (ask-before-remove confirmation gate, detached-HEAD reachability/rescue, tip-identity force-delete guard) was reverted to the main version because editing the prose repeatedly reintroduced shell-injection (last: sourcing a state file built from a git-derived branch name). Restore those guards in a `worktree_rm.py` that takes the name via file ingress, invokes git through argv arrays, and has a metacharacter-branch/path injection regression test — the pattern that worked for `pr_prepare.py`/`premerge_scan.py`. | `.claude/commands/worktree-rm.md`, `.claude/scripts/` | skill-audit | Heavy | Medium | | `premerge_scan.py` should scan the staged blob (`git show :path`) for staged methodology files, not the working-tree copy — a stage-then-revert-working-copy edit currently reads the safe working version and misses the staged violation. Union staged-index findings with unstaged/untracked filesystem findings. | `.claude/scripts/premerge_scan.py` | skill-audit | Mid | Low | | Re-add committed-range methodology scanning to `/push-pr-update` §3b (clean tree, commits ahead) using `premerge_scan.py --range`, with the comparison ref passed as **data** (resolved into a quoted variable in one Bash call, never a raw ``). It was removed to avoid ref interpolation; the helper already implements and tests `--range`. | `.claude/commands/push-pr-update.md` | skill-audit | Quick | Low | diff --git a/benchmarks/python/benchmark_stacked_did.py b/benchmarks/python/benchmark_stacked_did.py index 51de4224d..e070005ef 100644 --- a/benchmarks/python/benchmark_stacked_did.py +++ b/benchmarks/python/benchmark_stacked_did.py @@ -91,10 +91,13 @@ def main(): kappa_pre=kappa_pre, kappa_post=kappa_post, weighting="aggregate", - clean_control="not_yet_treated", + control_group="not_yet_treated", cluster="unit", ) + # The event-study surface is always computed at fit (3.9, row M-024); + # the deprecated aggregate= kwarg would emit a FutureWarning inside + # the timed block and contaminate the measurement. with Timer() as estimation_timer: results = est.fit( df, @@ -102,7 +105,6 @@ def main(): unit="unit", time="time", first_treat="first_treat", - aggregate="event_study", ) estimation_time = estimation_timer.elapsed diff --git a/diff_diff/diagnostic_report.py b/diff_diff/diagnostic_report.py index c33a5d6ab..bb7f02750 100644 --- a/diff_diff/diagnostic_report.py +++ b/diff_diff/diagnostic_report.py @@ -771,6 +771,20 @@ def _instance_skip_reason(self, check: str) -> Optional[str]: "SpilloverDiD(..., event_study=True) to populate the " "per-event-time direct-effect output." ) + # StackedDiD (3.9, row M-024): the event-study surface + # is ALWAYS computed, so ``aggregate='event_study'`` is + # inert (and deprecated - it now warns). Zero estimated + # pre-periods here means the event window has none: + # the default ``kappa_pre=1`` grid minus the reference + # leaves no pre-period columns. Estimator-accurate + # remediation, per the SpilloverDiD precedent above. + if name == "StackedDiDResults": + return ( + "No pre-period event-study coefficients exist on " + "this fit: the event window has no estimated " + "pre-periods. Re-fit with kappa_pre >= 2 so " + "pre-treatment event-study coefficients exist." + ) return ( "No pre-period event-study coefficients are exposed on " "this fit. For staggered estimators, re-fit with " @@ -1774,8 +1788,10 @@ def _infer_cov_source(source_fit: Any) -> str: ``event_study_vcov is None`` (bootstrap or replicate-weight CS / SA / TwoStageDiD fits, plus ImputationDiD / EfficientDiD / etc. which don't yet expose ``event_study_vcov``; StackedDiD - persists its VCV in every inference mode, so it reaches this - fallback only when no event study was requested); + persists its VCV in every inference mode and, since 3.9 + (row M-024), always materializes the event-study surface, so + post-3.9 StackedDiD fits never reach this fallback - only + pre-3.9 pickles can); OR ``MultiPeriodDiDResults`` without ``interaction_indices`` (genuine diag-only path inside ``pretrends.py:_extract_pre_period_params``, no "available but unused" concern, so no downgrade applies). diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 93a748600..bb924895e 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -1032,7 +1032,7 @@ stacked.fit( unit: str, time: str, first_treat: str, - aggregate: str = None, # None, "simple", or "event_study" + aggregate: Any = NOT_SUPPLIED, # DEPRECATED (M-024): warns; surface always computed - use results.aggregate() post-fit population: str = None, # Required when weighting="population" covariates: list[str] = None, # Columns to balance (requires balance="entropy"); values read at t=a-1-anticipation; balanced windows only ) -> StackedDiDResults @@ -1045,10 +1045,11 @@ from diff_diff import StackedDiD, plot_event_study est = StackedDiD(kappa_pre=2, kappa_post=2) results = est.fit(data, outcome='outcome', unit='unit', - time='period', first_treat='first_treat', - aggregate='event_study') + time='period', first_treat='first_treat') results.print_summary() plot_event_study(results) +es = results.aggregate('event_study') # unified container (M-024) +simple = results.aggregate('simple') # one-row overall relay ``` ### EfficientDiD @@ -1760,7 +1761,7 @@ ImputationDiDResults, TwoStageDiDResults, StackedDiDResults, and EfficientDiDRes | `overall_t_stat` | `float` | T-statistic | | `overall_p_value` | `float` | P-value | | `overall_conf_int` | `tuple[float, float]` | CI | -| `event_study_effects` | `dict[int, dict]` | Event study effects (if aggregate includes event_study) | +| `event_study_effects` | `dict[int, dict]` | Event study effects (if aggregate includes event_study; ALWAYS populated for StackedDiD since 3.9 - row M-024) | | `group_effects` | `dict` | Group-level effects (if aggregate includes group) | | `groups` | `list` | Treatment cohorts | | `time_periods` | `list` | All time periods | diff --git a/diff_diff/guides/llms-practitioner.txt b/diff_diff/guides/llms-practitioner.txt index 2650ef568..433e8811e 100644 --- a/diff_diff/guides/llms-practitioner.txt +++ b/diff_diff/guides/llms-practitioner.txt @@ -342,7 +342,9 @@ This step is CRITICAL and most often skipped. Run at least one of: ### HonestDiD (Rambachan & Roth 2023) - recommended Bounds on the treatment effect under violations of parallel trends. Works with MultiPeriodDiD, CallawaySantAnna, and ChaisemartinDHaultfoeuille -(dCDH) results. +(dCDH) results, plus StackedDiD via the post-fit +`results.aggregate('event_study')` container (needs `kappa_pre >= 2` so +estimated pre-periods exist). - For CS: pass the post-fit container - `compute_honest_did( results.aggregate('event_study'))` - no refit needed (the fit-time diff --git a/diff_diff/guides/llms.txt b/diff_diff/guides/llms.txt index 9fc1cfe1d..14301b3a2 100644 --- a/diff_diff/guides/llms.txt +++ b/diff_diff/guides/llms.txt @@ -20,8 +20,8 @@ diagnostic steps produces unreliable results. 3. **Test parallel trends** — simple 2x2: `check_parallel_trends()`, `equivalence_test_trends()`; staggered: inspect CS event-study pre-period coefficients (generic PT tests are invalid for staggered designs). Insignificant pre-trends do NOT prove PT holds. 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; `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. +6. **Sensitivity analysis** — `compute_honest_did(results)` for bounds under PT violations (MultiPeriodDiD, CS, or dCDH natively; a StackedDiD `results.aggregate('event_study')` container also admits - needs `kappa_pre >= 2` so estimated pre-periods exist), `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; `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')`; Stacked: `results.aggregate('event_study')`/`.aggregate('simple')` post-fit views (the surface is ALWAYS computed at fit since 3.9 - row M-024 - and the container admits into `compute_honest_did`/`compute_pretrends_power` with `kappa_pre >= 2`); 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 @@ -88,7 +88,7 @@ The site is organized into 5 sections, each with a landing page: - [Parallel Trends Testing](https://diff-diff.readthedocs.io/en/stable/api/diagnostics.html): Simple and Wasserstein-robust parallel trends tests, equivalence testing (TOST) - [Placebo Tests](https://diff-diff.readthedocs.io/en/stable/api/diagnostics.html): Placebo timing, group, permutation, and leave-one-out diagnostics - [Honest DiD](https://diff-diff.readthedocs.io/en/stable/api/honest_did.html): Rambachan & Roth (2023) sensitivity analysis — robust CI under parallel trends violations, breakdown values -- [Pre-Trends Power Analysis](https://diff-diff.readthedocs.io/en/stable/api/pretrends.html): Roth (2022) Section II.A-B no-individually-significant (NIS) box-probability pretest power + minimum detectable violation; `pretest_form='nis'` (default) implements the paper's primary form, `pretest_form='wald'` retained as paper-supported alternative (Propositions 1+3+4 all apply); linear-violation MDV in Roth's γ units when relative-time labels are threaded through `fit()`; full Σ_22 routing on non-bootstrap CallawaySantAnna and SunAbraham adapters +- [Pre-Trends Power Analysis](https://diff-diff.readthedocs.io/en/stable/api/pretrends.html): Roth (2022) Section II.A-B no-individually-significant (NIS) box-probability pretest power + minimum detectable violation; `pretest_form='nis'` (default) implements the paper's primary form, `pretest_form='wald'` retained as paper-supported alternative (Propositions 1+3+4 all apply); linear-violation MDV in Roth's γ units when relative-time labels are threaded through `fit()`; full Σ_22 routing on non-bootstrap CallawaySantAnna and SunAbraham adapters and on admitted CS-/StackedDiD-sourced `aggregate('event_study')` containers (StackedDiD persists its ES VCV in every inference mode) - [Power Analysis](https://diff-diff.readthedocs.io/en/stable/api/power.html): Analytical and simulation-based power analysis — MDE, sample size, power curves for study design - [MMM Calibration Export](https://diff-diff.readthedocs.io/en/stable/api/mmm.html): Assemble Marketing Mix Model calibration inputs from experiment results (explicit-in / validated-out - the caller passes the already-scoped incremental outcome + SE, the module does NOT rescale a headline ATT). `to_pymc_marketing_lift_test(channel, x, delta_x, delta_y, sigma, dims=, on_wrong_sign=)` builds the PyMC-Marketing/prophetverse lift-test DataFrame with sign/zero/positivity guards. `to_meridian_roi_prior(incremental_outcome, incremental_outcome_se, spend, parameter="roi_m"|"mroi_m", se_widening=)` builds Google Meridian lognormal ROI priors (spend-weighted pooling, lognormal parity with `lognormal_dist_from_mean_std`, channel- and time-scoped `.to_code()` snippet setting `media_prior_type`). Pure numpy/pandas; imports no MMM package; does not introspect result objects. Deriving totals from a fit is deferred to the post-4.0 `results.aggregate()` layer. - Conley spatial HAC SE (`vcov_type="conley"`) on cross-sectional `LinearRegression` / `compute_robust_vcov` PLUS panel `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` (with `conley_lag_cutoff=` for within-unit Bartlett temporal HAC) — Conley (1999) spatial-correlation-aware SEs with haversine/euclidean/callable distance metric and Bartlett/uniform spatial kernel; panel path uses the R `conleyreg`-form block-decomposed sandwich (within-period spatial + within-unit Bartlett serial, same-time excluded); parity vs R `conleyreg` (Düsterhöft 2021) on cross-sectional AND panel `lag_cutoff > 0` fixtures. Combining with explicit `cluster=` applies the combined spatial + cluster product kernel `K_total[i,j] = K_space · 1{c_i = c_j}` (cluster must be constant within each unit across periods on the panel path; validator-enforced). DiD takes `unit=` as a fit-time kwarg when `vcov_type="conley"` (not on `__init__`). Sparse k-d-tree fast path auto-activates for `n > 5_000` with bartlett kernel + haversine/euclidean metric diff --git a/diff_diff/honest_did.py b/diff_diff/honest_did.py index 0f5cae5a4..594adf2af 100644 --- a/diff_diff/honest_did.py +++ b/diff_diff/honest_did.py @@ -583,21 +583,27 @@ def _extract_container_params( """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. + ``CallawaySantAnnaResults.aggregate('event_study')`` or + ``StackedDiDResults.aggregate('event_study')`` (row M-024; Stacked + containers require ``kappa_pre >= 2`` so estimated pre-periods + exist, and a pre-period covariance that is not singular - honest + validates with ``allow_singular=False``, so keep ``kappa_pre`` small + relative to the cluster count). 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": + if surface.source not in ("CallawaySantAnnaResults", "StackedDiDResults"): raise TypeError( "HonestDiD accepts EventStudyResults containers produced by " - "CallawaySantAnnaResults.aggregate('event_study') only " + "CallawaySantAnnaResults.aggregate('event_study') or " + "StackedDiDResults.aggregate('event_study') only " f"(got source={surface.source!r}). For other estimators pass " "the native results object where supported " "(MultiPeriodDiDResults, CallawaySantAnnaResults, or " @@ -605,6 +611,7 @@ def _extract_container_params( "further producers arrives with their own aggregate() " "migrations." ) + _producer = surface.source.replace("Results", "") if surface.time_scale != "relative": raise TypeError( "HonestDiD requires a relative-time event-study container; " @@ -623,11 +630,12 @@ def _extract_container_params( _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 " + f"This {_producer} 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.", + "period cannot be verified (a universal base on a gapped " + "time grid may mix cohort-specific bases). Producer-built " + "containers record it - re-aggregate from the fitted " + "results object.", UserWarning, stacklevel=4, ) @@ -635,15 +643,14 @@ def _extract_container_params( 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." + f"{_producer} container records DISTINCT normalization bases " + f"at event times {sorted(set(_ref_e))}. Coefficients " + "normalized against different bases are not jointly " + "interpretable under Rambachan-Roth's delta_0 = 0 " + "normalization. For CallawaySantAnna this arises from " + "base_period='universal' on a gapped time grid - " + "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 @@ -659,12 +666,21 @@ def _extract_container_params( "use consecutive comparisons (not a common reference period), " "which changes the meaning of the parallel trends restriction. " ) + # Remedy is producer-conditional: only CallawaySantAnna exposes a + # base_period knob (StackedDiD is single-reference by construction + # and always records "universal", so a Stacked-sourced container + # can only land here hand-built/modified). + remedy_clause = ( + "Re-run with CallawaySantAnna(base_period='universal') for " + "methodologically valid HonestDiD bounds." + if surface.source == "CallawaySantAnnaResults" + else "Rebuild the container from the fitted results object " + "(producer-built containers record the true regime)." + ) 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.", + f"HonestDiD sensitivity analysis on {_producer} results " + "requires a universal (common-reference) base for valid " + "interpretation. " + provenance_clause + remedy_clause, UserWarning, stacklevel=4, ) @@ -686,6 +702,32 @@ def _extract_container_params( ) keep = (~surface.is_reference) & np.isfinite(surface.se) & (surface.se > 0) + + # Withheld-inference transparency (row M-024, StackedDiD admission): + # a retained row with finite se but non-finite p_value means the + # producer withheld/never computed its per-row inference (StackedDiD's + # hc2_bm BM-DOF fail-close and replicate-undefined designs both emit + # this shape; a hand-built container can too). The bounds below + # consume only the point estimates and covariance - both valid - so + # the row is ADMITTED, with a warning so the withheld state is never + # silently laundered. Shape-descriptive by design: the container + # cannot prove WHY the producer withheld it. Source-scoped to Stacked + # so no CS-path behavior changes (CS replicate-undefined containers + # keep their pinned silent fail-close via df_survey=0.0). + if surface.source == "StackedDiDResults": + _withheld = keep & ~np.isfinite(surface.p_value) + if bool(np.any(_withheld)): + warnings.warn( + "This container carries non-reference rows whose stored " + "per-row inference is withheld/undefined (finite se, " + "non-finite p_value) at event times " + f"{surface.event_time[_withheld].tolist()}; results are " + "computed from the point estimates and covariance alone, " + "using this consumer's own reference distribution.", + UserWarning, + stacklevel=4, + ) + rel_times = [t for t, k in zip(surface.event_time.tolist(), keep) if k] ref_rows = surface.event_time[surface.is_reference].tolist() @@ -698,11 +740,11 @@ def _extract_container_params( "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." + "reference. For CallawaySantAnna, multiple references arise " + "from 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 @@ -742,14 +784,25 @@ def _extract_container_params( has_gap = True break if has_gap: + # Remedy is producer-conditional: balance_e is CS aggregation + # machinery; StackedDiD's aggregate() has no balance_e level + # (its kappa trimming already balances retained windows), so + # recommending it there would point at a non-existent recovery. + _gap_remedy = ( + "Ensure all event-study periods have valid estimates, " + "or use balance_e to restrict to a balanced subset." + if surface.source == "CallawaySantAnnaResults" + else "Ensure all event-study periods have valid estimates " + "(for StackedDiD, an interior horizon with a non-finite SE " + "indicates a rank-dropped event-time column - inspect the " + "fit's rank_deficient_action warnings)." + ) 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." + "some event-study horizons have non-finite SEs. " + _gap_remedy ) # beta_hat/sigma are subset in EXPLICIT [sorted pre; sorted post] @@ -827,8 +880,11 @@ def _extract_event_study_params( 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``). + ``CallawaySantAnnaResults.aggregate('event_study')`` or + ``StackedDiDResults.aggregate('event_study')`` (CS- and + Stacked-sourced containers only; see + ``_extract_container_params``. Stacked containers need + ``kappa_pre >= 2`` so estimated pre-periods exist). Returns ------- @@ -1393,7 +1449,8 @@ def _largest_consecutive_block(times, boundary_val): f"Unsupported results type: {type(results)}. " "Expected MultiPeriodDiDResults, CallawaySantAnnaResults, " "ChaisemartinDHaultfoeuilleResults, or an EventStudyResults " - "container from CallawaySantAnnaResults.aggregate('event_study')." + "container from CallawaySantAnnaResults.aggregate('event_study') " + "or StackedDiDResults.aggregate('event_study')." ) @@ -2823,13 +2880,20 @@ def fit( 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. + ``CallawaySantAnnaResults.aggregate('event_study')`` or + ``StackedDiDResults.aggregate('event_study')`` (Stacked + containers require ``kappa_pre >= 2`` so estimated + pre-periods exist, and a non-singular pre-period covariance - + keep ``kappa_pre`` small relative to the cluster count). On + the container route the scalar inference df arrives via the + container's ``df_survey`` provenance field (for CS-sourced + containers, bounds and CIs match the native route exactly; + StackedDiD has no native HonestDiD route - the container IS + its route, with normal-theory critical values on analytical + fits), 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. @@ -3201,7 +3265,11 @@ def sensitivity_analysis( Parameters ---------- results : MultiPeriodDiDResults, CallawaySantAnnaResults, ChaisemartinDHaultfoeuilleResults, or EventStudyResults - Results from event study estimation. + Results from event study estimation, or the unified + event-study container from + ``CallawaySantAnnaResults.aggregate('event_study')`` or + ``StackedDiDResults.aggregate('event_study')`` (Stacked + containers require ``kappa_pre >= 2``). M_grid : list of float, optional Grid of M values to evaluate. If None, uses default grid based on method. @@ -3303,7 +3371,11 @@ def breakdown_value( Parameters ---------- results : MultiPeriodDiDResults, CallawaySantAnnaResults, ChaisemartinDHaultfoeuilleResults, or EventStudyResults - Results from event study estimation. + Results from event study estimation, or the unified + event-study container from + ``CallawaySantAnnaResults.aggregate('event_study')`` or + ``StackedDiDResults.aggregate('event_study')`` (Stacked + containers require ``kappa_pre >= 2``). tol : float Tolerance for binary search. @@ -3360,7 +3432,10 @@ def compute_honest_did( Parameters ---------- results : MultiPeriodDiDResults, CallawaySantAnnaResults, ChaisemartinDHaultfoeuilleResults, or EventStudyResults - Results from event study estimation. + Results from event study estimation, or the unified event-study + container from ``CallawaySantAnnaResults.aggregate('event_study')`` + or ``StackedDiDResults.aggregate('event_study')`` (Stacked + containers require ``kappa_pre >= 2``). method : str Type of restriction ("smoothness", "relative_magnitude", "combined"). M : float @@ -3402,7 +3477,10 @@ def sensitivity_plot( Parameters ---------- results : MultiPeriodDiDResults, CallawaySantAnnaResults, ChaisemartinDHaultfoeuilleResults, or EventStudyResults - Results from event study estimation. + Results from event study estimation, or the unified event-study + container from ``CallawaySantAnnaResults.aggregate('event_study')`` + or ``StackedDiDResults.aggregate('event_study')`` (Stacked + containers require ``kappa_pre >= 2``). method : str Type of restriction. M_grid : list of float, optional diff --git a/diff_diff/practitioner.py b/diff_diff/practitioner.py index b0cf41a49..4aba5aa0a 100644 --- a/diff_diff/practitioner.py +++ b/diff_diff/practitioner.py @@ -586,7 +586,15 @@ def _handle_stacked(results: Any): ), code="# Check results.n_sub_experiments and inspect results.stacked_data", priority="medium", - step_name="heterogeneity", + # 3.9 (row M-024): a DISTINCT key, deliberately not + # "heterogeneity". Since the StackedDiD event-study surface is + # always populated, DiagnosticReport's heterogeneity check now + # runs on every plain fit and marks that step completed - which + # used to silently drop this UNRELATED balance advice from + # next_steps (the step_name collision). Like "loo_jackknife", + # this key stays OUT of the STEPS completion vocabulary: no + # diagnostic ever completes it, so the advice always survives. + step_name="sub_experiment_balance", ), _robustness_compare_step("CS, SA, or BJS"), ] diff --git a/diff_diff/pretrends.py b/diff_diff/pretrends.py index 8f39d0af0..e6f1d3d80 100644 --- a/diff_diff/pretrends.py +++ b/diff_diff/pretrends.py @@ -1047,8 +1047,10 @@ def _extract_pre_period_params( 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). + ``CallawaySantAnnaResults.aggregate('event_study')`` or + ``StackedDiDResults.aggregate('event_study')`` (CS- and + Stacked-sourced containers only; Stacked containers require + ``kappa_pre >= 2`` so estimated pre-periods exist). pre_periods : list of int, optional Explicit list of pre-treatment periods. If None, uses results.pre_periods. @@ -1071,8 +1073,11 @@ def _extract_pre_period_params( - ``"full_pre_period_vcov"`` when a full pre-period covariance sub-block was used (MPD with - ``interaction_indices``, or CS/SA with populated - ``event_study_vcov``). + ``interaction_indices``, CS/SA with populated + ``event_study_vcov``, or an admitted CS-/Stacked-sourced + container carrying ``vcov`` - StackedDiD persists its + event-study VCV in every inference mode, so Stacked + containers always take this tier). - ``"diag_fallback"`` when only the per-period standard errors were available (bootstrap / replicate-weight CS or SA fits, MPD without ``interaction_indices``). @@ -1407,7 +1412,8 @@ def _extract_pre_period_params( f"Unsupported results type: {type(results)}. " "Expected MultiPeriodDiDResults, CallawaySantAnnaResults, " "SunAbrahamResults, or an EventStudyResults container from " - "CallawaySantAnnaResults.aggregate('event_study')." + "CallawaySantAnnaResults.aggregate('event_study') or " + "StackedDiDResults.aggregate('event_study')." ) def _extract_container_pre_period_params( @@ -1418,7 +1424,10 @@ def _extract_container_pre_period_params( """Container branch of ``_extract_pre_period_params``. Consumes the unified ``EventStudyResults`` surface produced by - ``CallawaySantAnnaResults.aggregate('event_study')``. Admission is + ``CallawaySantAnnaResults.aggregate('event_study')`` or + ``StackedDiDResults.aggregate('event_study')`` (row M-024; + Stacked containers require ``kappa_pre >= 2`` so estimated + pre-periods exist). 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 @@ -1427,10 +1436,11 @@ def _extract_container_pre_period_params( from pre-trend coefficients, and PreTrendsPower has no native dCDH branch either. """ - if surface.source != "CallawaySantAnnaResults": + if surface.source not in ("CallawaySantAnnaResults", "StackedDiDResults"): raise TypeError( "PreTrendsPower accepts EventStudyResults containers " "produced by CallawaySantAnnaResults.aggregate(" + "'event_study') or StackedDiDResults.aggregate(" "'event_study') only " f"(got source={surface.source!r}). For other estimators " "pass the native results object where supported " @@ -1438,6 +1448,7 @@ def _extract_container_pre_period_params( "SunAbrahamResults); container admission for further " "producers arrives with their own aggregate() migrations." ) + _producer = surface.source.replace("Results", "") if surface.time_scale != "relative": raise TypeError( "PreTrendsPower requires a relative-time event-study " @@ -1462,10 +1473,10 @@ def _extract_container_pre_period_params( _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 " + f"This {_producer} 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 " + "period cannot be verified (a universal base on a gapped " + "time grid may mix cohort-specific bases). Producer-built " "containers record it - re-aggregate from the fitted " "results object.", UserWarning, @@ -1475,14 +1486,16 @@ def _extract_container_pre_period_params( 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." + f"this {_producer} container records DISTINCT " + "normalization bases at event times " + f"{sorted(set(_ref_e))}. The hypothesized violation " + "delta is defined relative to one reference, so " + "power/MDV over mixed-base coefficients target an " + "ill-defined alternative. For CallawaySantAnna this " + "arises from base_period='universal' on a gapped time " + "grid - 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 @@ -1506,12 +1519,19 @@ def _extract_container_pre_period_params( "against a common reference period, so linear power/MDV " "target a different violation shape. " ) + remedy_clause = ( + "Re-run with CallawaySantAnna(base_period='universal') " + "for Roth-faithful linear benchmarks." + if surface.source == "CallawaySantAnnaResults" + else "Rebuild the container from the fitted results " + "object (producer-built containers record the true " + "regime)." + ) warnings.warn( - "PreTrendsPower on a CallawaySantAnna event-study " - "container without base_period='universal'. " + f"PreTrendsPower on a {_producer} event-study " + "container without a universal (common-reference) base. " + provenance_clause - + "Re-run with CallawaySantAnna(base_period='universal') " - "for Roth-faithful linear benchmarks.", + + remedy_clause, UserWarning, stacklevel=4, ) @@ -1529,6 +1549,32 @@ def _extract_container_pre_period_params( & (surface.event_time < _pre_cutoff) ) + # Withheld-inference transparency (row M-024, StackedDiD + # admission; sibling of the HonestDiD container branch): a + # retained row with finite se but non-finite p_value means the + # producer withheld/never computed its per-row inference + # (StackedDiD's hc2_bm BM-DOF fail-close and replicate-undefined + # designs both emit this shape). Power/MDV consume only the + # coefficients and covariance - both valid - so the rows are + # ADMITTED with a warning (no-silent-failures). Deliberately NO + # df_survey conjunct: pretrends never reads df_survey, so its + # output is finite in every df state. Shape-descriptive by + # design; source-scoped so no CS-path behavior changes. + if surface.source == "StackedDiDResults": + _withheld = keep & ~np.isfinite(surface.p_value) + if bool(np.any(_withheld)): + warnings.warn( + "This container carries non-reference rows whose " + "stored per-row inference is withheld/undefined " + "(finite se, non-finite p_value) at event times " + f"{surface.event_time[_withheld].tolist()}; results " + "are computed from the point estimates and " + "covariance alone, using this consumer's own " + "reference distribution.", + UserWarning, + stacklevel=4, + ) + if not keep.any(): raise ValueError("No pre-treatment periods found in event study.") @@ -1891,7 +1937,11 @@ def fit( Parameters ---------- results : MultiPeriodDiDResults, CallawaySantAnnaResults, SunAbrahamResults, or EventStudyResults - Results from an event study estimation. + Results from an event study estimation, or the unified + event-study container from + ``CallawaySantAnnaResults.aggregate('event_study')`` or + ``StackedDiDResults.aggregate('event_study')`` (Stacked + containers require ``kappa_pre >= 2``). M : float, optional Specific violation magnitude to evaluate. If None, evaluates at a default magnitude based on the data. @@ -2146,7 +2196,10 @@ def compute_pretrends_power( Parameters ---------- results : results object - Event study results. + Event study results (MPD / CS / SA), or the unified event-study + container from ``CallawaySantAnnaResults.aggregate('event_study')`` + or ``StackedDiDResults.aggregate('event_study')`` (Stacked + containers require ``kappa_pre >= 2``). M : float, optional Violation magnitude to evaluate. alpha : float, default=0.05 @@ -2210,7 +2263,10 @@ def compute_mdv( Parameters ---------- results : results object - Event study results. + Event study results (MPD / CS / SA), or the unified event-study + container from ``CallawaySantAnnaResults.aggregate('event_study')`` + or ``StackedDiDResults.aggregate('event_study')`` (Stacked + containers require ``kappa_pre >= 2``). alpha : float, default=0.05 Significance level. target_power : float, default=0.80 diff --git a/diff_diff/results_base.py b/diff_diff/results_base.py index 36895777f..92c52d038 100644 --- a/diff_diff/results_base.py +++ b/diff_diff/results_base.py @@ -436,8 +436,17 @@ def __post_init__(self) -> None: "(explicit ordering) or neither." ) if self.vcov is not None: - self.vcov = np.asarray(self.vcov, dtype=float) - self.vcov_index = np.asarray(self.vcov_index) + # COPY, not asarray (matching the _ARRAY_FIELDS convention + # above): np.asarray is a no-op view when the dtype already + # matches, which would alias the PRODUCER's stored matrix - + # e.g. StackedDiDResults.event_study_vcov reaches this + # constructor unmodified via the post-fit aggregate() view + # (row M-024), and a caller mutating the container would + # silently corrupt the fitted result. vcov_index keeps its + # native dtype (int labels must not become floats in + # to_dict()) but is copied for the same isolation. + self.vcov = np.array(self.vcov, dtype=float) + self.vcov_index = np.array(self.vcov_index) k = self.vcov_index.shape[0] if self.vcov.shape != (k, k): raise ValueError( @@ -606,7 +615,9 @@ def summary(self, alpha: Optional[float] = None) -> str: "CallawaySantAnnaResults": "call results.aggregate('event_study')", "ImputationDiDResults": "refit with aggregate='event_study' (or 'all')", "TwoStageDiDResults": "refit with aggregate='event_study' (or 'all')", - "StackedDiDResults": "refit with aggregate='event_study'", + # Absence only possible on pre-3.9 pickles: 3.9+ fits always + # materialize the surface (row M-024). + "StackedDiDResults": "re-fit with diff-diff >= 3.9, which always computes the surface", "StaggeredTripleDiffResults": "refit with aggregate='event_study' (or 'all')", "EfficientDiDResults": "refit with aggregate='event_study' (or 'all')", "ContinuousDiDResults": "refit with aggregate='eventstudy' (or 'all')", diff --git a/diff_diff/stacked_did.py b/diff_diff/stacked_did.py index e7efa79e0..17905e439 100644 --- a/diff_diff/stacked_did.py +++ b/diff_diff/stacked_did.py @@ -174,11 +174,10 @@ class StackedDiD(BaseEstimator): ... time='period', first_treat='first_treat') >>> results.print_summary() - With event study: + The event-study surface is always computed (3.9, row M-024) - view + it post-fit: - >>> results = est.fit(data, outcome='outcome', unit='unit', - ... time='period', first_treat='first_treat', - ... aggregate='event_study') + >>> es = results.aggregate('event_study') >>> from diff_diff import plot_event_study >>> plot_event_study(results) @@ -345,7 +344,7 @@ def fit( unit: str, time: str, first_treat: str, - aggregate: Optional[str] = None, + aggregate: Any = NOT_SUPPLIED, population: Optional[str] = None, survey_design=None, covariates: Optional[List[str]] = None, @@ -367,11 +366,16 @@ def fit( Name of column indicating when unit was first treated. Use 0 or np.inf for never-treated units. aggregate : str, optional - Aggregation mode: None/"simple" (overall ATT only) or - "event_study". Group aggregation is not supported because - the pooled stacked regression cannot produce cohort-specific - effects. Use CallawaySantAnna or ImputationDiD for - cohort-level estimates. + DEPRECATED (3.9, removed in 4.0; ledger row M-024): the + event-study surface is now always computed at fit, so this + parameter is behaviorally inert - passing it (any value, + ``None`` included) emits a ``FutureWarning``. Aggregate as a + post-fit step instead: ``results.aggregate('event_study')`` / + ``results.aggregate('simple')``. Value validation is + unchanged: ``"group"``/``"all"`` raise ``ValueError`` + (the pooled stacked regression cannot produce + cohort-specific effects - use CallawaySantAnna or + ImputationDiD), as do unknown values. population : str, optional Column name for population weights. Required only when weighting="population". @@ -398,6 +402,28 @@ def fit( ValueError If required columns are missing or data validation fails. """ + # ---- fit(aggregate=) shim (row M-024) ---- + # Deprecated in 3.9, removed at 4.0: the event-study surface is now + # ALWAYS computed at fit (see the extraction block below), so the + # param is behaviorally inert and aggregation moves to the post-fit + # results.aggregate(type=...). The sentinel default means the + # warning fires ONLY when the caller supplies the argument (None + # included), never on a plain fit(). Value validation below is + # unchanged, so the deprecated path returns exactly what a plain + # fit does - or raises exactly as it always did. + if aggregate is not NOT_SUPPLIED: + warnings.warn( + "StackedDiD.fit(aggregate=) is deprecated and will be " + "removed in 4.0. The event-study surface is now always " + "computed; aggregate as a post-fit step: " + "results.aggregate('event_study') / " + "results.aggregate('simple').", + FutureWarning, + stacklevel=2, + ) + else: + aggregate = None + # ---- Validate inputs ---- if aggregate in ("group", "all"): raise ValueError( @@ -741,18 +767,24 @@ def fit( # Only build contrasts whose target column is identified; if a # delta_h column itself was dropped, that event-time will get # NaN inference (left to safe_inference's df=None path). - # Per CI codex R1 P3: skip per-event contrast DOFs when the - # event-study surface is not user-visible (aggregate != "event_study"). - # The overall ATT contrast still gets computed below. + # 3.9 (row M-024): the per-event contrasts are ALWAYS built - + # the event-study surface is always materialized so the + # post-fit results.aggregate("event_study") view is total. + # (The pre-M-024 gate skipped them when aggregate != + # "event_study"; that gate's premise died with the post-fit + # surface.) NB: the dof helper's degeneracy guard is + # batch-relative for m > 1 (linalg.py), so the overall-ATT + # contrast now shares a batch with the per-event contrasts; + # the widened noise-floor carve-out is documented in the + # REGISTRY M-024 Note and pinned in test_aggregate_contract. es_keys: List[int] = [] es_cols_full: List[np.ndarray] = [] - if aggregate == "event_study": - for h in event_times: - if h in interaction_indices and _identified[interaction_indices[h]]: - c = np.zeros(k_design) - c[interaction_indices[h]] = 1.0 - es_keys.append(h) - es_cols_full.append(c) + for h in event_times: + if h in interaction_indices and _identified[interaction_indices[h]]: + c = np.zeros(k_design) + c[interaction_indices[h]] = 1.0 + es_keys.append(h) + es_cols_full.append(c) # Overall ATT contrast: average of post-period delta_h columns # (the same 1/K * ones contrast used for overall_se below). Only # construct if ALL post-period delta_h are identified — otherwise @@ -895,96 +927,95 @@ def _refit_stacked(w_r): survey_metadata = compute_survey_metadata(resolved_composed, raw_w_stacked) # ---- Extract event study effects ---- - event_study_effects: Optional[Dict[int, Dict[str, Any]]] = None + # 3.9 (row M-024): the surface is ALWAYS materialized - the + # regression always includes the event-time interactions, and the + # post-fit results.aggregate("event_study") view must be total. + # fit(aggregate=) no longer affects what is computed or stored. + event_study_effects: Dict[int, Dict[str, Any]] = {} es_vcov: Optional[np.ndarray] = None es_vcov_index: Optional[List[int]] = None - es_df_used: Optional[Dict[int, float]] = None - if aggregate == "event_study": - event_study_effects = {} - es_df_used = {} - # Reference period (e = -1 - anticipation) - event_study_effects[ref_period] = { - "effect": 0.0, - "se": 0.0, - "t_stat": np.nan, - "p_value": np.nan, - "conf_int": (np.nan, np.nan), - "n_obs": 0, - } - for h in event_times: - idx = interaction_indices[h] - effect = float(coef[idx]) - se = float(np.sqrt(max(vcov[idx, idx], 0.0))) - _survey_df = ( - max(survey_metadata.df_survey, 1) - if survey_metadata is not None and survey_metadata.df_survey is not None - else (0 if _uses_replicate_sd else None) + es_df_used: Dict[int, float] = {} + # Reference period (e = -1 - anticipation) + event_study_effects[ref_period] = { + "effect": 0.0, + "se": 0.0, + "t_stat": np.nan, + "p_value": np.nan, + "conf_int": (np.nan, np.nan), + "n_obs": 0, + } + for h in event_times: + idx = interaction_indices[h] + effect = float(coef[idx]) + se = float(np.sqrt(max(vcov[idx, idx], 0.0))) + _survey_df = ( + max(survey_metadata.df_survey, 1) + if survey_metadata is not None and survey_metadata.df_survey is not None + else (0 if _uses_replicate_sd else None) + ) + # Override df when replicate replicates were dropped + if _n_valid_rep_sd is not None and resolved_stacked is not None: + if _n_valid_rep_sd < resolved_stacked.n_replicates: + _survey_df = _n_valid_rep_sd - 1 if _n_valid_rep_sd > 1 else 0 + if survey_metadata is not None: + survey_metadata.df_survey = _survey_df if _survey_df > 0 else None + # Use BM contrast DOF for hc2_bm. Fail-closed: when the + # hc2_bm contract is in effect but BM DOF is unavailable (helper + # failed OR noise-floor NaN guard fired), emit all-NaN inference + # rather than fall back to normal-theory CIs/p-values. Mirrors + # the fix in LinearRegression.get_inference() from PR #475 R7 + # (linalg.py:3689-3706). `safe_inference` itself has guarded + # non-finite / <= 0 df since PR #620 (utils.py) with the same + # all-NaN result; this explicit branch predates that guard and + # is kept for explicitness (it also skips the call entirely). + _is_hc2bm_path = self.vcov_type == "hc2_bm" and not _uses_replicate_sd + _bm_df = _bm_contrast_dof_per_event.get(h) + if _is_hc2bm_path and (_bm_df is None or not np.isfinite(_bm_df)): + # BM DOF unavailable on hc2_bm path: NaN-out inference. + t_stat = float("nan") + p_value = float("nan") + conf_int = (float("nan"), float("nan")) + # No safe_inference call happened -> no df provenance. + es_df_used[h] = float("nan") + else: + # BM DOF > survey/replicate df > the knob-resolved + # analytical fallback (residual t by default; the + # replicate 0-sentinel is not None, so it keeps + # precedence over the fallback). + _df_eff = ( + _bm_df + if _bm_df is not None + else (_survey_df if _survey_df is not None else _analytic_fallback_df) ) - # Override df when replicate replicates were dropped - if _n_valid_rep_sd is not None and resolved_stacked is not None: - if _n_valid_rep_sd < resolved_stacked.n_replicates: - _survey_df = _n_valid_rep_sd - 1 if _n_valid_rep_sd > 1 else 0 - if survey_metadata is not None: - survey_metadata.df_survey = _survey_df if _survey_df > 0 else None - # Use BM contrast DOF for hc2_bm. Fail-closed: when the - # hc2_bm contract is in effect but BM DOF is unavailable (helper - # failed OR noise-floor NaN guard fired), emit all-NaN inference - # rather than fall back to normal-theory CIs/p-values. Mirrors - # the fix in LinearRegression.get_inference() from PR #475 R7 - # (linalg.py:3689-3706). `safe_inference` itself has guarded - # non-finite / <= 0 df since PR #620 (utils.py) with the same - # all-NaN result; this explicit branch predates that guard and - # is kept for explicitness (it also skips the call entirely). - _is_hc2bm_path = self.vcov_type == "hc2_bm" and not _uses_replicate_sd - _bm_df = _bm_contrast_dof_per_event.get(h) - if _is_hc2bm_path and (_bm_df is None or not np.isfinite(_bm_df)): - # BM DOF unavailable on hc2_bm path: NaN-out inference. - t_stat = float("nan") - p_value = float("nan") - conf_int = (float("nan"), float("nan")) - # No safe_inference call happened -> no df provenance. - es_df_used[h] = float("nan") - else: - # BM DOF > survey/replicate df > the knob-resolved - # analytical fallback (residual t by default; the - # replicate 0-sentinel is not None, so it keeps - # precedence over the fallback). - _df_eff = ( - _bm_df - if _bm_df is not None - else (_survey_df if _survey_df is not None else _analytic_fallback_df) - ) - t_stat, p_value, conf_int = safe_inference( - effect, se, alpha=self.alpha, df=_df_eff - ) - # Record the df actually handed to safe_inference iff it - # governed a t-reference (finite, > 0); None (normal - # theory) and the df<=0 sentinels record NaN. - es_df_used[h] = ( - float(_df_eff) - if _df_eff is not None and np.isfinite(_df_eff) and _df_eff > 0 - else float("nan") - ) - n_obs_h = int(np.sum((et_vals == h) & (d_vals == 1))) - event_study_effects[h] = { - "effect": effect, - "se": se, - "t_stat": t_stat, - "p_value": p_value, - "conf_int": conf_int, - "n_obs": n_obs_h, - } - - # Persist the event-time sub-block of the pooled-regression VCV - # (the reported ES SEs are exactly its diagonal in every - # inference mode - analytical sandwich, replicate refit, and - # survey TSL all reassign `vcov` before this block). The - # reference period is synthesized, never a regression column, so - # the index is the ESTIMATED event times only. - if event_times: - _delta_cols = [interaction_indices[h] for h in event_times] - es_vcov = vcov[np.ix_(_delta_cols, _delta_cols)] - es_vcov_index = [int(h) for h in event_times] + t_stat, p_value, conf_int = safe_inference(effect, se, alpha=self.alpha, df=_df_eff) + # Record the df actually handed to safe_inference iff it + # governed a t-reference (finite, > 0); None (normal + # theory) and the df<=0 sentinels record NaN. + es_df_used[h] = ( + float(_df_eff) + if _df_eff is not None and np.isfinite(_df_eff) and _df_eff > 0 + else float("nan") + ) + n_obs_h = int(np.sum((et_vals == h) & (d_vals == 1))) + event_study_effects[h] = { + "effect": effect, + "se": se, + "t_stat": t_stat, + "p_value": p_value, + "conf_int": conf_int, + "n_obs": n_obs_h, + } + + # Persist the event-time sub-block of the pooled-regression VCV + # (the reported ES SEs are exactly its diagonal in every + # inference mode - analytical sandwich, replicate refit, and + # survey TSL all reassign `vcov` before this block). The + # reference period is synthesized, never a regression column, so + # the index is the ESTIMATED event times only. + if event_times: + _delta_cols = [interaction_indices[h] for h in event_times] + es_vcov = vcov[np.ix_(_delta_cols, _delta_cols)] + es_vcov_index = [int(h) for h in event_times] # ---- Compute overall ATT ---- # Average of post-treatment delta_h coefficients with delta-method SE @@ -1623,7 +1654,7 @@ def stacked_did( first_treat: str, kappa_pre: int = 1, kappa_post: int = 1, - aggregate: Optional[str] = None, + aggregate: Any = NOT_SUPPLIED, population: Optional[str] = None, survey_design=None, covariates: Optional[List[str]] = None, @@ -1651,7 +1682,10 @@ def stacked_did( kappa_post : int, default=1 Post-treatment event-time periods. aggregate : str, optional - Aggregation mode: None, "simple", or "event_study". + DEPRECATED (3.9, removed in 4.0; row M-024) - the event-study + surface is always computed; passing this warns and changes + nothing. Use ``results.aggregate('event_study')`` / + ``results.aggregate('simple')`` post-fit. population : str, optional Population column for weighting="population". survey_design : SurveyDesign, optional @@ -1674,9 +1708,9 @@ def stacked_did( >>> from diff_diff import stacked_did, generate_staggered_data >>> data = generate_staggered_data(seed=42) >>> results = stacked_did(data, 'outcome', 'unit', 'period', - ... 'first_treat', kappa_pre=2, kappa_post=2, - ... aggregate='event_study') + ... 'first_treat', kappa_pre=2, kappa_post=2) >>> results.print_summary() + >>> es = results.aggregate('event_study') # post-fit view (M-024) """ est = StackedDiD(kappa_pre=kappa_pre, kappa_post=kappa_post, **kwargs) return est.fit( diff --git a/diff_diff/stacked_did_results.py b/diff_diff/stacked_did_results.py index 5e9dc8dd0..89c810063 100644 --- a/diff_diff/stacked_did_results.py +++ b/diff_diff/stacked_did_results.py @@ -6,14 +6,15 @@ """ from dataclasses import dataclass, field -from typing import Any, Dict, List, Optional, Tuple +from typing import Any, ClassVar, Dict, List, Optional, Tuple import numpy as np import pandas as pd from diff_diff._deprecation import deprecated_field_property +from diff_diff.aggregation import AggregationMixin, AggregationResult from diff_diff.results import _format_survey_block, _get_significance_stars -from diff_diff.results_base import BaseResults +from diff_diff.results_base import BaseResults, build_event_study_surface __all__ = [ "StackedDiDResults", @@ -21,7 +22,7 @@ @dataclass -class StackedDiDResults(BaseResults): +class StackedDiDResults(BaseResults, AggregationMixin): """ Results from Stacked DiD estimation (Wing, Freedman & Hollingsworth 2024). @@ -41,6 +42,10 @@ class StackedDiDResults(BaseResults): event_study_effects : dict, optional Dictionary mapping event time h to effect dict with keys: 'effect', 'se', 't_stat', 'p_value', 'conf_int', 'n_obs'. + Always populated on 3.9+ fits (the pooled regression always + includes the event-time interactions and the surface is always + extracted - row M-024); None only on pre-3.9 pickles. View it as + the unified container via ``aggregate('event_study')``. group_effects : dict, optional Dictionary mapping cohort g to effect dict. stacked_data : pd.DataFrame @@ -83,7 +88,9 @@ class StackedDiDResults(BaseResults): (analytical hc1/hc2_bm sandwich, survey replicate refit, and survey TSL all produce the coefficient covariance the SEs are read from). The reference period is synthesized, never a regression column, so - it is absent from the index. None when no event study was requested. + it is absent from the index. Always populated on 3.9+ fits (the + event-study surface is always materialized - row M-024); None + only on pre-3.9 pickles. event_study_vcov_index : list of int, optional Event-time labels ordering ``event_study_vcov``'s rows/columns (the estimated event times, reference excluded). @@ -96,7 +103,8 @@ class StackedDiDResults(BaseResults): fallback otherwise — finite residual df under the 3.9 default, ``G − 1`` under "cluster"), or NaN when the row used normal theory (``df_convention="normal"``), the df was undefined, or hc2_bm - failed closed. None when no event study was requested. + failed closed. Always populated on 3.9+ fits (row M-024); None + only on pre-3.9 pickles. df_convention : str, optional The estimator's ``df_convention`` configuration echoed onto the results ("residual" | "cluster" | "normal"; added 3.9). @@ -177,6 +185,88 @@ def __setstate__(self, state: Dict[str, Any]) -> None: state["control_group"] = state.pop("clean_control") self.__dict__.update(state) + # ------------------------------------------------------------------ + # Container-consumer provenance (rows M-024 / M-093). Class-level + # attrs, deliberately NOT dataclass fields: the generated __init__'s + # positional indexes are public API, and both values are derivable. + # ------------------------------------------------------------------ + #: Every sub-experiment normalizes against the single omitted + #: reference ``e = -1 - anticipation`` - universal-base semantics in + #: the CallawaySantAnna vocabulary. Read by ``_provenance_kwargs`` + #: so honest_did's cannot-verify-universal-base fail-safe stays + #: silent on StackedDiD containers. + base_period: ClassVar[str] = "universal" + + @property + def reference_event_times(self) -> Tuple[int, ...]: + """The singleton common reference event time, ``(-1 - anticipation,)``. + + StackedDiD has exactly one omitted reference shared by every + sub-experiment, so the container consumers' common-reference + guard always sees a single entry (rows M-024 / M-093). + """ + return (-1 - int(self.anticipation),) + + # ------------------------------------------------------------------ + # Post-fit aggregation (row M-024, on the M-122 contract). Both + # levels are pure VIEWS over stored fields - the event-study surface + # is always materialized at fit since 3.9, and "simple" relays the + # stored overall inference bit-exactly. Nothing is recomputed, so + # every inference mode (survey TSL, replicate refit, hc2_bm) relays + # faithfully. + # ------------------------------------------------------------------ + # ClassVar: on a dataclass a bare annotation would turn this routing + # configuration into an ``__init__`` field. + _AGGREGATE_SUPPORTED: ClassVar[Tuple[str, ...]] = ("simple", "event_study") + # StackedDiD has no balance_e machinery on any aggregation level + # (kappa trimming already balances every retained cohort's window). + _AGGREGATE_BALANCE_E_TYPES: ClassVar[Tuple[str, ...]] = () + + 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`` (always + # populated on 3.9+ fits; pre-3.9 pickles raise the absent- + # surface error with a re-fit hint). + return build_event_study_surface(self) + + # level == "simple": one-row view relaying the stored overall + # inference. ``n`` is the TREATED-unit count: StackedDiD's treated + # and control sets OVERLAP (a later-treated unit is treated in its + # own sub-experiment and a clean control in earlier ones), so a + # disjoint total does not exist as a stored scalar and summing the + # two counts would double-count - deliberately narrower in scope + # than CallawaySantAnna's treated+control "units" (REGISTRY + # StackedDiD M-024 Note; cross-container ``n`` comparisons are out + # of contract for this estimator). ``target`` is "att" per the CS + # precedent: ``overall_att`` is the equally-weighted average of + # post-treatment event-study coefficients, NOT the per-event-time + # trimmed aggregate ATT, so weighting-specific target strings + # would misstate the scalar (``describe_target_parameter`` is the + # estimand's prose source of truth). + 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(["att"], 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([float(self.n_treated_units)], dtype=float), + df=np.array( + [float(self.inference_df) if self.inference_df is not None else float("nan")], + dtype=float, + ), + alpha=self.alpha, + n_kind="units", + weight=np.array([1.0], dtype=float), + estimator=type(self).__name__.replace("Results", ""), + ) + # --- Inference-field aliases (balance/external-adapter compatibility) --- @property def att(self) -> float: @@ -462,8 +552,12 @@ def to_dataframe(self, level: str = "event_study") -> pd.DataFrame: """ if level == "event_study": if self.event_study_effects is None: + # Only reachable on pre-3.9 pickles: 3.9+ fits always + # materialize the surface (row M-024). raise ValueError( - "Event study effects not computed. " "Use aggregate='event_study'." + "Event study effects not present on this results object. " + "Re-fit with diff-diff >= 3.9, which always computes the " + "event-study surface." ) rows = [] for h, data in sorted(self.event_study_effects.items()): diff --git a/docs/api/_autosummary/diff_diff.StackedDiDResults.rst b/docs/api/_autosummary/diff_diff.StackedDiDResults.rst index 55e201e46..618b77755 100644 --- a/docs/api/_autosummary/diff_diff.StackedDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.StackedDiDResults.rst @@ -12,6 +12,7 @@ .. autosummary:: ~StackedDiDResults.__init__ + ~StackedDiDResults.aggregate ~StackedDiDResults.print_summary ~StackedDiDResults.summary ~StackedDiDResults.to_dataframe @@ -28,6 +29,7 @@ ~StackedDiDResults.anticipation ~StackedDiDResults.att ~StackedDiDResults.balance + ~StackedDiDResults.base_period ~StackedDiDResults.balance_diagnostics ~StackedDiDResults.cluster_name ~StackedDiDResults.coef_var @@ -49,6 +51,7 @@ ~StackedDiDResults.n_sub_experiments ~StackedDiDResults.n_treated_units ~StackedDiDResults.p_value + ~StackedDiDResults.reference_event_times ~StackedDiDResults.se ~StackedDiDResults.significance_stars ~StackedDiDResults.survey_metadata diff --git a/docs/api/stacked_did.rst b/docs/api/stacked_did.rst index 5c572e46f..e7fa4c2a1 100644 --- a/docs/api/stacked_did.rst +++ b/docs/api/stacked_did.rst @@ -58,6 +58,7 @@ Results container for Stacked DiD estimation. .. autosummary:: + ~StackedDiDResults.aggregate ~StackedDiDResults.summary ~StackedDiDResults.print_summary ~StackedDiDResults.to_dataframe @@ -79,10 +80,14 @@ Basic usage:: est = StackedDiD(kappa_pre=2, kappa_post=2) results = est.fit(data, outcome='outcome', unit='unit', - time='period', first_treat='first_treat', - aggregate='event_study') + time='period', first_treat='first_treat') results.print_summary() + # The event-study surface is always computed (3.9, row M-024); + # view it post-fit as the unified container: + es = results.aggregate('event_study') + simple = results.aggregate('simple') + Accessing the stacked dataset:: # The stacked data is available for custom analysis diff --git a/docs/choosing_estimator.rst b/docs/choosing_estimator.rst index 6d2b0109f..d447cae2e 100644 --- a/docs/choosing_estimator.rst +++ b/docs/choosing_estimator.rst @@ -573,10 +573,12 @@ with clean controls and apply Q-weights to reweight the stacked regression. stk = StackedDiD(kappa_pre=2, kappa_post=3) results = stk.fit(data, outcome='y', unit='unit_id', - time='period', first_treat='first_treat', - aggregate='event_study') + time='period', first_treat='first_treat') results.print_summary() + # The event-study surface is always computed (3.9); view it post-fit: + es = results.aggregate('event_study') + .. note:: The trimmed aggregate ATT may exclude early or late cohorts whose event diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 2ca0d3b70..42bf94021 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2081,6 +2081,8 @@ The pooled estimator is `DID^{CBWSDID}_e = Σ_a (N^D_a/N^D_Ω)(Δ̄^D_{a,e} − - **Note:** Provenance exposure (4.0 program row M-092, no numeric change): the event-time sub-block of the pooled-regression coefficient covariance is now persisted as `StackedDiDResults.event_study_vcov` (+ `event_study_vcov_index` = the estimated event times; the synthesized reference period is not a regression column and is absent). Persisted in every inference mode - analytical hc1/hc2_bm sandwich, survey replicate refit, and survey TSL all reassign the coefficient covariance the reported ES SEs are read from, so the marginal SEs are always this matrix's diagonal. `event_study_df` records, per event time, the df actually passed to `safe_inference` (the per-event Bell-McCaffrey Satterthwaite contrast df under `hc2_bm`, the scalar survey df under survey designs, NaN on the normal-theory hc1 path and on hc2_bm fail-closed rows). DiagnosticReport's parallel-trends check consumes the persisted VCV via its joint-Wald path for hc1/survey fits; `hc2_bm` fits are routed to Bonferroni over the BM-adjusted per-row p-values instead (the chi-square Wald reference would discard the CR2/BM small-sample correction; multi-constraint AHT/HTZ CR2 testing is tracked in DEFERRED.md - see REPORTING.md "hc2_bm parallel-trends policy"). - **Note:** Inference is conditional-on-the-estimated-weights cluster-robust (the existing `hc1`/`hc2_bm` path with `W_{sa}` as the WLS weights) — the paper's default. The paper's weight-re-estimating bootstrap is NOT implemented in v1 (deliberate scope; entropy balancing is smooth so the Abadie–Imbens (2008) nonsmooth-matching bootstrap caveat does not apply). `cluster` is orthogonal to `b_{sa}` (weights conditioned-on); default `unit` matches the paper. - **Note:** v1 scope — only `balance="entropy"` with `weighting="aggregate"`. `balance` + `population`/`sample_share` and `balance` + `survey_design=` raise `NotImplementedError`; matching-based balancing and the repeated `0→1/1→0` episode extension are out of scope. +- **Note (post-fit aggregate() - row M-024):** `fit(aggregate=)` is deprecated (3.9; removed 4.0) and the event-study surface is ALWAYS materialized at fit — the pooled regression always included the event-time interactions, so only the extraction and, under `hc2_bm`, the per-event Bell-McCaffrey contrast DOFs were fit-time-gated; both gates are gone and the param is behaviorally inert (its `group`/`all`/unknown ValueErrors persist after the shim warning). `aggregate('event_study')` returns the unified `EventStudyResults` container (a pure view; no recomputation, so every stored inference mode — survey TSL, replicate refit, hc2_bm fail-close — relays faithfully); `aggregate('simple')` relays `overall_att/se/t/p/CI` bit-exactly with `target="att"` (the scalar is the equally-weighted post-period average, NOT the per-event-time trimmed aggregate ATT `θ^e_κ` — `describe_target_parameter` remains the estimand's prose source of truth), `n = n_treated_units` with `n_kind="units"` (TREATED-only scope: under dynamic clean-control definitions a later-treated unit is treated in its own sub-experiment and a clean control in earlier ones, so the treated and control unit sets OVERLAP and no disjoint total exists — cross-container `n` comparisons vs CallawaySantAnna's treated+control convention are out of contract), and `df = inference_df` (the stored overall inference's exact provenance). Consequences of the always-computed surface: plain-fit `summary()` renders the event-study table, `to_dataframe('event_study')` is total on 3.9+ fits, `plot_event_study` renders plain fits, DiagnosticReport's `heterogeneity` (always) and `parallel_trends` (`kappa_pre >= 2`) checks run on plain fits, and the `hc2_bm` overall-ATT contrast now shares the m>1 batch of `_compute_cr2_bm_contrast_dof`, whose degeneracy guard is batch-relative — the noise-floor carve-out on the `Wald_test(test="HTZ")$df_denom` parity claim in the variance-families bullets above widens quantitatively (unchanged behavior when tripped: fail-closed NaN inference plus linalg's loud noise-floor UserWarning; the absolute floor already made this reachable at m=1). +- **Note (container admission - rows M-024/M-093):** `compute_honest_did` and `compute_pretrends_power` admit Stacked-sourced `aggregate('event_study')` containers. Provenance: `base_period` is the class-level constant `"universal"` (every sub-experiment normalizes against the single omitted reference `e = -1 - anticipation`) and `reference_event_times` is the singleton `(-1 - anticipation,)`. Requirements: `kappa_pre >= 2` (the default `kappa_pre=1` grid minus the reference has ZERO estimated pre-periods — both consumers reject), and for HonestDiD a non-singular pre-period covariance (`allow_singular=False` per Rambachan-Roth's eigenvalues-bounded-away-from-zero assumption — keep `kappa_pre` small relative to the cluster count; PreTrendsPower keeps its documented singular support). The container's scalar `df_survey` is None on analytical fits — the generic resolver's fallback reads `df_inference`, a name StackedDiD does not carry (`inference_df` is the overall POST-AVERAGE contrast df, not a design-level cluster df, and misusing it for FLCI critical values would be worse than the asymptotic default) — so honest FLCI critical values are NORMAL-THEORY on analytical Stacked fits, diverging from bare-`cluster=` CS fits (which thread `G-1`); survey fits thread the survey df and replicate-undefined designs the fail-closed `0.0` sentinel (NaN FLCI CIs; the identified-set `lb`/`ub` bounds are df-independent and stay finite). Withheld-inference rows (finite `se`, non-finite `p_value` — the hc2_bm BM-DOF fail-close and replicate-undefined designs both emit this shape) are ADMITTED with a source-scoped warning in BOTH consumers: bounds/power consume only the point estimates and covariance, which remain valid; the warning is shape-descriptive because containers are publicly constructible. A rank-dropped interior event-time column (NaN-se row) fails closed via honest's consecutive-grid ValueError with a Stacked-specific remedy (no `balance_e` — Stacked's aggregate() has no such level). *Covariate-balancing edge cases:* - Infeasible cohort (treated covariate mean outside the clean-control hull → entropy balancing cannot match the moments): **fail-closed** `ValueError` naming the cohort and worst covariate — NOT silently dropped (dropping a cohort would shift the estimand to an overlap-trimmed ATT, Ustyuzhanin 2026 §3.1). @@ -4648,7 +4650,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. +- **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; producer-built containers always record the field) warns fail-safe that the common reference cannot be verified. Container admission covers CS- and Stacked-sourced containers (row M-024): StackedDiD threads `base_period="universal"` and a singleton `reference_event_times` by construction (one omitted reference per fit), needs `kappa_pre >= 2` for estimated pre-periods, and HonestDiD additionally requires its pre-period covariance non-singular — see the StackedDiD section's container-admission Note for the df_survey and withheld-inference semantics. **Reference implementation(s):** - R: `HonestDiD` package (Rambachan & Roth's official package) @@ -4715,7 +4717,7 @@ 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, 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. +- **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. (Admitted Stacked-sourced containers — row M-024 — are universal-base by construction and never trip this warning from a producer-built container; the guard messages derive the producer from `surface.source`.) 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) @@ -4728,7 +4730,7 @@ Violation types: - Very high power: MDV approaches zero. - Symmetric two-sided pretests under parallel trends: `β̂_post` remains unbiased for `τ_post` (Roth Section II.B paragraph after Prop 1 — `E[β̂_pre | β̂_pre ∈ B] = 0` if B is symmetric and `β_pre = 0`). -- **Note (deviation from paper — diagonal pre-period VCV fallback, bootstrap-only after PR-B):** Roth (2022)'s power and bias objects operate on the full pre-period covariance block Σ_22. After PR-B 2026-05-17, the shipped `compute_pretrends_power` adapter consumes full Σ_22 on the non-bootstrap paths for ALL three result types: +- **Note (deviation from paper — diagonal pre-period VCV fallback, bootstrap-only after PR-B):** Roth (2022)'s power and bias objects operate on the full pre-period covariance block Σ_22. After PR-B 2026-05-17, the shipped `compute_pretrends_power` adapter consumes full Σ_22 on the non-bootstrap paths for ALL three natively supported result types (and, since row M-024, on admitted CS-/StackedDiD-sourced `aggregate('event_study')` containers carrying `vcov` — StackedDiD persists its ES VCV in every inference mode, so Stacked containers always take the full-Σ_22 tier): - `MultiPeriodDiDResults`: full pre-period sub-block from `results.vcov` when `interaction_indices` is populated; diag fallback only when `interaction_indices` is None. - `CallawaySantAnnaResults`: full `event_study_vcov` sub-block on non-bootstrap fits (the matrix is persisted at `staggered_results.py:126-128`). Bootstrap CS fits clear `event_study_vcov` at `staggered.py:2032-2036` to prevent mixing analytical VCV with bootstrap SEs, so they fall through to `diag(ses^2)`. - `SunAbrahamResults`: full `event_study_vcov` sub-block on non-bootstrap fits, constructed in `sun_abraham.py` via `W @ vcov_cohort @ W.T` where W is the cohort-aggregation matrix (PR-B Step 3 SA extension). Bootstrap SA fits and replicate-weight survey fits clear `event_study_vcov` for the same reason as CS. diff --git a/docs/methodology/REPORTING.md b/docs/methodology/REPORTING.md index a8ada36b4..533655e58 100644 --- a/docs/methodology/REPORTING.md +++ b/docs/methodology/REPORTING.md @@ -349,8 +349,11 @@ a library setting. As of the PR-B PreTrendsPower implementation audit (Roth 2022), `compute_pretrends_power()` consumes the full `event_study_vcov` sub-block when it is available — non-bootstrap CS fits - (`staggered_results.py` populates the matrix) and non-bootstrap SA - fits (`sun_abraham.py` builds it via `W @ vcov_cohort @ W.T`). The + (`staggered_results.py` populates the matrix), non-bootstrap SA + fits (`sun_abraham.py` builds it via `W @ vcov_cohort @ W.T`), and + admitted CS-/Stacked-sourced `aggregate('event_study')` containers + carrying `vcov` (row M-024; StackedDiD persists its ES VCV in every + inference mode, so Stacked containers always take this tier). The `PreTrendsPowerResults.covariance_source` field records the actual extraction path (`"full_pre_period_vcov"` vs `"diag_fallback"`), and the `DiagnosticReport.pretrends_power` block surfaces that label @@ -389,12 +392,14 @@ a library setting. diagonal), TwoStageDiD on the analytical paths only (bootstrap and replicate-weight modes clear it). Where the covariance is present, the PT check takes the joint-Wald path (subject to the hc2_bm - policy and rank guard below). Pretrends POWER is unaffected: - `compute_pretrends_power()` supports only MPD / CS / SunAbraham - fits, and Stacked/TwoStage remain outside DR's power applicability - - their new covariance is consumed by the PT check and by the - PRECOMPUTED-power provenance classifier only (a stored power result - on a fit that now exposes a full VCV is labelled + policy and rank guard below). Pretrends POWER: natively + `compute_pretrends_power()` supports MPD / CS / SunAbraham fits; + since row M-024 a Stacked `results.aggregate('event_study')` + CONTAINER also admits (kappa_pre >= 2). Stacked/TwoStage NATIVE + results remain outside DR's power applicability - within + DiagnosticReport their covariance is consumed by the PT check and + by the PRECOMPUTED-power provenance classifier only (a stored power + result on a fit that now exposes a full VCV is labelled `"diag_fallback_available_full_vcov_unused"`). - **Note:** hc2_bm parallel-trends policy (deviation by omission, diff --git a/docs/methodology/variance-conventions.md b/docs/methodology/variance-conventions.md index 0515aac3d..4b1360137 100644 --- a/docs/methodology/variance-conventions.md +++ b/docs/methodology/variance-conventions.md @@ -28,7 +28,7 @@ produced wrong figures three separate times while this inventory was drafted). | `twfe_hc1_cluster_unit_time_post` | 3 | 298 | **legitimate** | K_reference (D2 fixed): 2 visible + rank(post given unit) = 1; matches fixest cluster arm at rel 0 (committed golden) | | `wooldridge_hc1_within` | 15 | 286, 286, 286, 286, 286, 286, 286, 286, 286, 286 | **legitimate** | CR1 k = K_reference (D2 fixed: 9 cells + T = 15, no intercept col -> +1 term; jwdid arms at ratio 1.0); tail df converged in 3.9 (M-127): the default-hc1 arms use t(residual df = n - k_kept - absorbed rank) via the df_convention knob (was silent normal theory); G-1 under 'cluster', z under 'normal' | | `sun_abraham_hc1` | 21 | 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280 | **legitimate** | CR1 k = K_reference (D2 fixed: 15 cells + 6, no intercept col; fixest sunab parity ~5e-15); D4 fixed in 3.9 (M-127): cells AND aggregates share the saturated fit's residual df under the df_convention knob (aggregates previously dropped to normal theory inside the same fit) | -| `stacked_did_hc1` | 6 | 309 | **legitimate** | L1: k_total is clubSandwich CR1S by construction (stacked_did.py pins vcovCR(type='CR1S') at atol=1e-10); tail df converged in 3.9 (M-127): t(pooled residual df n_eff - k_kept) via the df_convention knob (was silent normal theory; stacked rows make the residual df large, so z -> t is a convention alignment) | +| `stacked_did_hc1` | 6 | 309, 309, 309 | **legitimate** | L1: k_total is clubSandwich CR1S by construction (stacked_did.py pins vcovCR(type='CR1S') at atol=1e-10); tail df converged in 3.9 (M-127): t(pooled residual df n_eff - k_kept) via the df_convention knob (was silent normal theory; stacked rows make the residual df large, so z -> t is a convention alignment) | | `lpdid_pre2_post2` | 4, 4, 5, 5, 5, 6 | 59, 59, 59, 59, 59, 59 | **legitimate** | L2: G-1 tail df (Stata/fixest convention) — since 3.9 the df_convention='cluster' DEFAULT on LPDiD (bit-identical; the one surface already at the v4 target) | | `mpd_absorb_hc1_cluster_unit` | 11 | 290, 290, 290, 290, 290, 290 | **legitimate** | K_reference: 6 visible + rank(time given unit) = 5; equals the fixed_effects form's 70 - 59 (MPD absorb/fixed_effects equivalence) | | `mpd_fixed_effects_hc1_cluster_unit` | 11 | 290, 290, 290, 290, 290, 290 | **legitimate** | K_reference: 70 visible (incl. built-in period dummies, MPD's time-FE block) minus the 59 cluster-nested unit dummies = 11 — identical to the absorb form | diff --git a/docs/tutorials/13_stacked_did.ipynb b/docs/tutorials/13_stacked_did.ipynb index 6eef48f91..690daadae 100644 --- a/docs/tutorials/13_stacked_did.ipynb +++ b/docs/tutorials/13_stacked_did.ipynb @@ -94,10 +94,12 @@ "metadata": {}, "outputs": [], "source": [ - "# Fit with event study aggregation\n", + "# The event-study surface is always computed at fit (3.9, ledger row\n", + "# M-024) - no fit-time aggregation flag needed. The unified container\n", + "# view is available post-fit via results_es.aggregate('event_study').\n", "est = StackedDiD(kappa_pre=2, kappa_post=2)\n", "results_es = est.fit(data, outcome='outcome', unit='unit', time='period',\n", - " first_treat='first_treat', aggregate='event_study')\n", + " first_treat='first_treat')\n", "\n", "# Plot event study\n", "if HAS_MATPLOTLIB:\n", @@ -450,4 +452,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} \ No newline at end of file +} diff --git a/docs/v4-deprecations.yaml b/docs/v4-deprecations.yaml index 837808632..20991cbee 100644 --- a/docs/v4-deprecations.yaml +++ b/docs/v4-deprecations.yaml @@ -274,11 +274,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/stacked_did.py] - notes: "StackedDiD's 'simple' value already in-vocabulary." + test_ref: tests/test_aggregate_contract.py + code_refs: [diff_diff/stacked_did.py, diff_diff/stacked_did_results.py, diff_diff/aggregation.py, diff_diff/results_base.py, diff_diff/honest_did.py, diff_diff/pretrends.py, diff_diff/diagnostic_report.py, diff_diff/practitioner.py, diff_diff/guides/llms-practitioner.txt] + notes: "Shimmed in 3.9: fit(aggregate=) warns via the shared NOT_SUPPLIED sentinel (a plain fit() never warns; supplying ANY value, None included, warns - CS-style warn-and-still-work, since the param genuinely worked, unlike dCDH's raise; group/all + unknown values keep their ValueErrors AFTER the warning). The event-study surface is now ALWAYS computed at fit - the pooled regression always included the interactions, and the extraction + per-event Bell-McCaffrey contrast DOFs (hc2_bm) are no longer gated on fit-time aggregate, so the param is behaviorally inert; consequences: plain-fit summary() prints the ES table, to_dataframe('event_study') is total on 3.9+ fits, DiagnosticReport's heterogeneity and (kappa_pre>=2) parallel_trends checks now run on plain fits, and the hc2_bm overall contrast shares the m>1 batch-relative noise floor (fail-closed + warned; REGISTRY Note). Both aggregate() levels are pure VIEWS: 'event_study' -> build_event_study_surface(self); 'simple' relays overall_att/se/t/p/CI bit-exact with target='att' (the scalar is the post-average, NOT the per-event trimmed aggregate ATT), n=n_treated_units/'units' (treated-only scope: control units overlap treated across sub-experiments, so no disjoint total exists), df=inference_df. SUPPORTED SUBSET simple|event_study - group/calendar fail closed via the mixin message; balance_e applies to no level (empty vocabulary); weights= rejected. Every inference mode relays (survey TSL, replicate refit, hc2_bm fail-close). honest/pretrends container admission widened to Stacked-sourced containers (see M-093): base_period='universal' ClassVar + singleton reference_event_times property thread the provenance; analytical fits carry df_survey=None (the resolver reads df_inference, a name Stacked does not use - deliberate: normal-theory FLCI, see the REGISTRY Note and the TODO adapter-naming row's cross-reference)." - id: M-025 kind: param group: aggregate-postfit @@ -995,7 +996,7 @@ rows: 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, 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." + 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. SECOND PRE-CUT AMENDMENT (2(b) PR-2, with M-024): admission widened to STACKED-SOURCED containers in compute_honest_did + compute_pretrends_power (source in {CallawaySantAnnaResults, StackedDiDResults}; kappa_pre >= 2 required for estimated pre-periods; honest additionally needs a non-singular pre-period covariance; withheld-inference rows admitted with a source-scoped warning; the seven producer-derived guard messages now derive the producer from surface.source). Source scoping continues for all other producers. 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 c6bef7cd4..928a09589 100644 --- a/docs/v4-design.md +++ b/docs/v4-design.md @@ -415,7 +415,12 @@ 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. +emitting identical column schemas from every estimator. The container +OWNS its arrays: `__post_init__` copies every array field including the +`vcov`/`vcov_index` pair (amended with [M-024] - `np.asarray` aliased the +producer's stored matrix on the post-fit view route, so mutating a +container could corrupt the fitted result; `vcov_index` keeps its native +dtype, int labels never become floats). **Pickle migration.** Renamed-field classes ship `__setstate__` migration following the existing `SyntheticDiDResults.__setstate__` precedent @@ -457,7 +462,16 @@ canonical 0..n-1 codes, keeping a shared results artifact free of names, emails or administrative IDs. Analytical-vs-bootstrap inference of the aggregated estimand follows the fit's inference method; where bootstrap draws are not retained, `aggregate()` on a bootstrapped fit RAISES rather than silently -returning analytical inference. Estimators whose +returning analytical inference. **View-relay exception (Phase 2b PRs 1-2):** +estimators whose `aggregate()` RELAYS stored fields without recomputation +need no influence-function kit - there is nothing to re-weight. The +retention requirement binds RECOMPUTING estimators (the CallawaySantAnna +class); dCDH [M-026] relays its stored overall/event-study surfaces +verbatim, and StackedDiD [M-024] always materializes its event-study +surface at fit (the pooled regression always includes the interactions) +so both its levels are pure views - which is also why bootstrap-style +fail-closing does not apply to them: every stored inference mode relays +faithfully. Estimators whose estimand is already a single aggregation (SunAbraham's saturated event study, LPDiD's per-horizon design) expose `aggregate()` where meaningful as additive surface; their native params (`only_event`/`only_pooled` etc.) are documented diff --git a/tests/test_aggregate_contract.py b/tests/test_aggregate_contract.py index 249f89770..40b5276e4 100644 --- a/tests/test_aggregate_contract.py +++ b/tests/test_aggregate_contract.py @@ -1106,6 +1106,355 @@ def test_mixin_hooks_are_not_dataclass_fields(self): 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. + # The roster must at least cover the shipped mixin adopters. assert "CallawaySantAnnaResults" in checked assert "ChaisemartinDHaultfoeuilleResults" in checked + assert "StackedDiDResults" in checked + + +# --------------------------------------------------------------------------- # +# StackedDiD (row M-024): fit(aggregate=) shim + the VIEW-based aggregate() +# --------------------------------------------------------------------------- # + +STACKED_KW = dict(outcome="outcome", unit="unit", time="period", first_treat="first_treat") + + +def _stacked_panel(seed=42, n_units=120, n_periods=12): + from diff_diff.prep_dgp import generate_staggered_data + + return generate_staggered_data( + n_units=n_units, n_periods=n_periods, cohort_periods=[4, 6, 8], seed=seed + ) + + +def _fit_stacked(data, *, est_kw=None, **fit_kw): + from diff_diff import StackedDiD + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return StackedDiD(**(est_kw or {"kappa_pre": 2, "kappa_post": 2})).fit( + data, **STACKED_KW, **fit_kw + ) + + +@pytest.fixture(scope="module") +def stacked_panel(): + return _stacked_panel() + + +@pytest.fixture(scope="module") +def stacked_fitted(stacked_panel): + """Plain hc1 fit - the surface is always materialized (M-024).""" + return _fit_stacked(stacked_panel) + + +class TestStackedShim: + def test_plain_fit_does_not_warn(self, stacked_panel): + from diff_diff import StackedDiD + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + StackedDiD(kappa_pre=2, kappa_post=2).fit(stacked_panel, **STACKED_KW) + assert [w for w in caught if issubclass(w.category, FutureWarning)] == [] + + def test_aggregate_kwarg_warns_even_at_none(self, stacked_panel): + from diff_diff import StackedDiD + + with pytest.warns(FutureWarning, match=r"fit\(aggregate=\) is deprecated"): + StackedDiD(kappa_pre=2, kappa_post=2).fit(stacked_panel, aggregate=None, **STACKED_KW) + + def test_deprecated_value_warns_and_still_works(self, stacked_panel, stacked_fitted): + # CS-style warn-and-still-work (the param genuinely worked here, + # unlike dCDH's raise): the deprecated path returns an object whose + # surface equals a plain fit's - the surface is always computed. + from diff_diff import StackedDiD + + with pytest.warns(FutureWarning, match=r"fit\(aggregate=\) is deprecated"): + res = StackedDiD(kappa_pre=2, kappa_post=2).fit( + stacked_panel, aggregate="event_study", **STACKED_KW + ) + assert res.overall_att == stacked_fitted.overall_att + assert res.event_study_effects is not None + assert sorted(res.event_study_effects) == sorted(stacked_fitted.event_study_effects) + np.testing.assert_array_equal(res.event_study_vcov, stacked_fitted.event_study_vcov) + + def test_wrapper_forwarded_aggregate_warns(self, stacked_panel): + # stacked_did() declares aggregate explicitly with its own sentinel + # default and forwards verbatim into fit(). + from diff_diff.stacked_did import stacked_did + + with pytest.warns(FutureWarning, match=r"fit\(aggregate=\) is deprecated"): + stacked_did( + stacked_panel, + "outcome", + "unit", + "period", + "first_treat", + kappa_pre=2, + kappa_post=2, + aggregate="simple", + ) + + def test_plain_wrapper_call_does_not_warn(self, stacked_panel): + from diff_diff.stacked_did import stacked_did + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + res = stacked_did( + stacked_panel, + "outcome", + "unit", + "period", + "first_treat", + kappa_pre=2, + kappa_post=2, + ) + assert [w for w in caught if issubclass(w.category, FutureWarning)] == [] + assert res.event_study_effects is not None + + def test_group_warns_then_raises_educational_error(self, stacked_panel): + from diff_diff import StackedDiD + + with pytest.warns(FutureWarning, match=r"fit\(aggregate=\) is deprecated"): + with pytest.raises(ValueError, match="not supported by StackedDiD"): + StackedDiD(kappa_pre=2, kappa_post=2).fit( + stacked_panel, aggregate="group", **STACKED_KW + ) + + +class TestStackedAggregate: + 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) + 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_event_study_view_matches_builder(self, stacked_fitted): + es = self._assert_surface_matches_builder(stacked_fitted) + # kappa 2/2 grid: {-2, ref -1, 0, 1, 2}; n_kind from n_obs cells. + assert sorted(es.event_time.tolist()) == [-2, -1, 0, 1, 2] + assert es.n_kind == "obs" + assert es.base_period == "universal" + assert es.reference_event_times == (-1,) + + def test_event_study_view_does_not_alias_fit_vcov(self, stacked_fitted): + es = stacked_fitted.aggregate("event_study") + assert not np.shares_memory(es.vcov, stacked_fitted.event_study_vcov) + # int labels survive to_dict (no float coercion of the index) + assert es.to_dict()["vcov_index"] == [-2, 0, 1, 2] + + def test_simple_bit_exact_relay_hc1(self, stacked_fitted): + agg = stacked_fitted.aggregate("simple") + assert isinstance(agg, AggregationResult) + assert agg.level == "simple" + assert list(agg.label) == ["overall"] + # target is "att": overall_att is the equally-weighted post-period + # average, NOT the per-event trimmed aggregate ATT (M-024 Note). + assert list(agg.target) == ["att"] + assert float(agg.att[0]) == stacked_fitted.overall_att + assert float(agg.se[0]) == stacked_fitted.overall_se + assert float(agg.t_stat[0]) == stacked_fitted.overall_t_stat + assert float(agg.p_value[0]) == stacked_fitted.overall_p_value + assert (float(agg.conf_int_lower[0]), float(agg.conf_int_upper[0])) == tuple( + stacked_fitted.overall_conf_int + ) + # Treated-unit count: control units OVERLAP treated across + # sub-experiments, so no disjoint total exists (M-024 Note). + assert float(agg.n[0]) == float(stacked_fitted.n_treated_units) + assert agg.n_kind == "units" + assert float(agg.df[0]) == float(stacked_fitted.inference_df) + assert agg.estimator == "StackedDiD" + + @pytest.mark.parametrize("weighting", ["aggregate", "population", "sample_share"]) + def test_simple_target_att_on_all_weighting_schemes(self, stacked_panel, weighting): + fit_kw = {} + panel = stacked_panel + if weighting == "population": + panel = stacked_panel.copy() + panel["pop"] = 100.0 + (panel["unit"] % 7) + fit_kw["population"] = "pop" + res = _fit_stacked( + panel, + est_kw={"kappa_pre": 2, "kappa_post": 2, "weighting": weighting}, + **fit_kw, + ) + agg = res.aggregate("simple") + assert list(agg.target) == ["att"] + assert agg.n_kind == "units" + + def test_simple_relay_hc2_bm_df_is_overall_bm_dof(self, stacked_panel): + res = _fit_stacked( + stacked_panel, + est_kw={ + "kappa_pre": 2, + "kappa_post": 2, + "vcov_type": "hc2_bm", + "cluster": "unit", + }, + ) + agg = res.aggregate("simple") + assert np.isfinite(agg.df[0]) + assert float(agg.df[0]) == float(res.inference_df) + # per-row BM dfs are present WITHOUT fit-time aggregate (M-024) + assert res.event_study_df is not None + assert all(np.isfinite(v) for v in res.event_study_df.values()) + + def test_simple_relay_bm_failure_nan_inference(self, stacked_panel, monkeypatch): + # The hc2_bm fail-closed state (finite att/se, jointly-NaN t/p/CI, + # inference_df None) must RELAY through the simple view - the df + # comparison is np.isnan, never df == inference_df (nan == None). + import diff_diff.linalg as dl + + def _nan_dof(X, cluster_ids, bread_matrix, contrasts, weights=None): + return np.full(contrasts.shape[1], np.nan) + + monkeypatch.setattr(dl, "_compute_cr2_bm_contrast_dof", _nan_dof) + res = _fit_stacked( + stacked_panel, + est_kw={ + "kappa_pre": 2, + "kappa_post": 2, + "vcov_type": "hc2_bm", + "cluster": "unit", + }, + ) + assert res.inference_df is None + agg = res.aggregate("simple") + assert np.isfinite(agg.att[0]) and np.isfinite(agg.se[0]) + assert np.isnan(agg.t_stat[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_relay_survey_tsl_df(self, stacked_panel): + from diff_diff.survey import SurveyDesign + + panel = stacked_panel.copy() + panel["w"] = 1.0 + 0.1 * (panel["unit"] % 5) + panel["strata"] = panel["unit"] % 4 + panel["psu"] = panel["unit"] + res = _fit_stacked( + panel, + survey_design=SurveyDesign(weights="w", strata="strata", psu="psu"), + ) + agg = res.aggregate("simple") + assert float(agg.att[0]) == res.overall_att + # The stored overall inference used the survey df; the relay + # carries exactly that provenance. + assert float(agg.df[0]) == float(res.inference_df) + + def test_bm_dof_batch_parity_and_overall_reconstruction(self, stacked_panel, monkeypatch): + # Pin (a): per-contrast Satterthwaite dof is column-independent in + # VALUE on a well-conditioned design - the m=1 evaluation of each + # contrast equals its column in the batched call (the batch-relative + # noise floor changes only the DEGENERACY GUARD's scale, documented + # in the REGISTRY M-024 Note). + # Pin (b): the pre-change PLAIN-fit overall inference (m=1 batch) is + # reconstructed in-process from the spy-captured fit-time locals and + # must match the post-change stored overall inference at 1e-14. + import diff_diff.linalg as dl + from diff_diff.utils import safe_inference + + real = dl._compute_cr2_bm_contrast_dof + captured = {} + + def spy(X, cluster_ids, bread_matrix, contrasts, weights=None): + captured.update( + X=X, + cluster_ids=cluster_ids, + bread=bread_matrix, + contrasts=contrasts, + weights=weights, + ) + return real(X, cluster_ids, bread_matrix, contrasts, weights=weights) + + monkeypatch.setattr(dl, "_compute_cr2_bm_contrast_dof", spy) + res = _fit_stacked( + stacked_panel, + est_kw={ + "kappa_pre": 2, + "kappa_post": 2, + "vcov_type": "hc2_bm", + "cluster": "unit", + }, + ) + assert captured, "spy never fired" + contrasts = captured["contrasts"] + batched = real( + captured["X"], + captured["cluster_ids"], + captured["bread"], + contrasts, + weights=captured["weights"], + ) + for j in range(contrasts.shape[1]): + single = real( + captured["X"], + captured["cluster_ids"], + captured["bread"], + contrasts[:, [j]], + weights=captured["weights"], + ) + np.testing.assert_allclose(single[0], batched[j], rtol=1e-14) + # (b) overall contrast is appended LAST at the fit site; its m=1 + # dof + safe_inference reproduce the stored overall inference. + overall_dof_m1 = float( + real( + captured["X"], + captured["cluster_ids"], + captured["bread"], + contrasts[:, [-1]], + weights=captured["weights"], + )[0] + ) + np.testing.assert_allclose(res.inference_df, overall_dof_m1, rtol=1e-14) + t, p, ci = safe_inference( + res.overall_att, res.overall_se, alpha=res.alpha, df=overall_dof_m1 + ) + np.testing.assert_allclose( + [t, p, ci[0], ci[1]], + [ + res.overall_t_stat, + res.overall_p_value, + res.overall_conf_int[0], + res.overall_conf_int[1], + ], + rtol=1e-14, + ) + + def test_legacy_pickle_absent_surface_hint(self, stacked_fitted): + import dataclasses + + legacy = dataclasses.replace( + stacked_fitted, + event_study_effects=None, + event_study_vcov=None, + event_study_vcov_index=None, + event_study_df=None, + ) + with pytest.raises(ValueError, match=r"diff-diff >= 3\.9"): + legacy.aggregate("event_study") + + def test_balance_e_rejected_empty_vocabulary(self, stacked_fitted): + with pytest.raises(ValueError, match="no aggregation type on this estimator"): + stacked_fitted.aggregate("event_study", balance_e=1) + + def test_weights_rejected(self, stacked_fitted): + with pytest.raises(ValueError, match="does not accept a weights selector"): + stacked_fitted.aggregate("simple", weights="cell") + + @pytest.mark.parametrize("bad", ["group", "calendar", "all", "nonsense"]) + def test_unsupported_types_fail_closed(self, stacked_fitted, bad): + with pytest.raises(ValueError, match="Unsupported aggregation type"): + stacked_fitted.aggregate(bad) diff --git a/tests/test_business_report.py b/tests/test_business_report.py index df1a65d0b..1d1a4b816 100644 --- a/tests/test_business_report.py +++ b/tests/test_business_report.py @@ -4958,3 +4958,50 @@ def test_bootstrap_fit_uses_bootstrap_label_on_alpha_override(self): f"allow-list regression: bootstrap fits must not fall through to the " f"analytical label. caveat message: {message!r}" ) + + +class TestStackedAlwaysComputedSurfaceBusinessReport: + """M-024 ripple: BusinessReport(auto_diagnostics=True) lifts DR's + schema, so plain-fit BR output gains the PT/heterogeneity narrative - + and the renamed practitioner step key keeps the sub-experiment + balance advice in next_steps (the step_name collision fix).""" + + @staticmethod + def _panel(seed=42): + rng = np.random.default_rng(seed) + rows = [] + for u in range(60): + g = [4, 6, 0][u % 3] + for t in range(1, 11): + y = 1.0 + 0.1 * t + u * 0.01 + (1.5 if g and t >= g else 0.0) + rng.normal(0, 0.3) + rows.append({"unit": u, "time": t, "outcome": y, "first_treat": g}) + return pd.DataFrame(rows) + + def test_plain_fit_br_runs_checks_and_keeps_balance_step(self): + from diff_diff import StackedDiD + from diff_diff.business_report import BusinessReport + + df = self._panel() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = StackedDiD(kappa_pre=2, kappa_post=2).fit( + df, + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + ) + d = BusinessReport( + res, + data=df, + unit="unit", + time="time", + outcome="outcome", + first_treat="first_treat", + auto_diagnostics=True, + ).to_dict() + schema = d["diagnostics"]["schema"] + assert schema["parallel_trends"]["status"] == "ran" + assert schema["heterogeneity"]["status"] == "ran" + labels = [s["label"] for s in d["next_steps"]] + assert "Check sub-experiment balance" in labels diff --git a/tests/test_diagnostic_report.py b/tests/test_diagnostic_report.py index f65b2ff24..1748c372b 100644 --- a/tests/test_diagnostic_report.py +++ b/tests/test_diagnostic_report.py @@ -3001,7 +3001,6 @@ def test_stacked_pt_check_uses_joint_wald(self): unit="unit", time="time", first_treat="first_treat", - aggregate="event_study", ) rep = DiagnosticReport(res).run_all() pt = rep.to_dict()["parallel_trends"] @@ -3091,7 +3090,6 @@ def _hc2bm_fit(monkeypatch=None): unit="unit", time="time", first_treat="first_treat", - aggregate="event_study", ) def test_stacked_hc2_bm_fail_closed_pt_is_inconclusive(self, monkeypatch): @@ -3222,7 +3220,6 @@ def test_stacked_singular_vcov_downgrades_to_bonferroni(self): unit="unit", time="time", first_treat="first_treat", - aggregate="event_study", ) # Precondition: the retained pre-period block is genuinely singular. V = res.event_study_vcov @@ -3335,3 +3332,66 @@ class _FakeSingular: assert pt["method"] == "bonferroni" stat = pt["test_statistic"] assert stat is None or stat >= 0.0 + + +class TestStackedAlwaysComputedSurfaceCheckFlips: + """M-024: the always-materialized ES surface flips two DR checks on + PLAIN StackedDiD fits (executed-verified during the plan review): + heterogeneity always, parallel_trends at kappa_pre >= 2.""" + + @staticmethod + def _panel(seed=42): + rng = np.random.default_rng(seed) + rows = [] + for u in range(60): + g = [4, 6, 0][u % 3] + for t in range(1, 11): + y = 1.0 + 0.1 * t + u * 0.01 + (1.5 if g and t >= g else 0.0) + rng.normal(0, 0.3) + rows.append({"unit": u, "time": t, "outcome": y, "first_treat": g}) + return pd.DataFrame(rows) + + def _fit(self, kappa_pre): + from diff_diff import StackedDiD + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return StackedDiD(kappa_pre=kappa_pre, kappa_post=2).fit( + self._panel(), + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + ) + + def _report(self, res): + from diff_diff.diagnostic_report import DiagnosticReport + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return DiagnosticReport( + res, + data=self._panel(), + unit="unit", + time="time", + outcome="outcome", + first_treat="first_treat", + ).to_dict() + + def test_kappa_pre_2_runs_heterogeneity_and_parallel_trends(self): + d = self._report(self._fit(kappa_pre=2)) + assert d["parallel_trends"]["status"] == "ran" + assert d["heterogeneity"]["status"] == "ran" + assert d["heterogeneity"]["source"] == "event_study_effects_post" + + def test_kappa_pre_1_default_heterogeneity_only(self): + # The default kappa_pre=1 grid minus the reference has ZERO + # estimated pre-periods, so PT stays skipped - with the M-024 + # per-type remediation, NOT the generic (now-inert for Stacked) + # aggregate='event_study' advice. + d = self._report(self._fit(kappa_pre=1)) + assert d["heterogeneity"]["status"] == "ran" + pt = d["parallel_trends"] + assert pt["status"] == "skipped" + reason = pt.get("reason", "") + " ".join(str(v) for v in pt.values() if isinstance(v, str)) + assert "kappa_pre >= 2" in reason + assert "aggregate='event_study'" not in reason diff --git a/tests/test_event_study_consumers.py b/tests/test_event_study_consumers.py index b626b369c..c593b98d3 100644 --- a/tests/test_event_study_consumers.py +++ b/tests/test_event_study_consumers.py @@ -14,10 +14,11 @@ 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). +- SOURCE-SCOPED ADMISSION: honest/pretrends accept CS- and Stacked-sourced + containers (the latter widened with row M-024; ``kappa_pre >= 2`` + required); dCDH l1 containers, calendar containers, other-producer/e0 + containers and hand-built source=None containers fail closed. The + plotters take no source guard (label-faithful rendering). """ import warnings @@ -960,7 +961,7 @@ def test_multiple_reference_rows_fail_closed_in_honest(self): # consecutive grid instead. msg = str(exc_info.value) assert "native" not in msg - assert "consecutive" in msg and "Re-estimate" in msg + assert "consecutive" in msg and "re-estimate" in msg # --------------------------------------------------------------------------- # @@ -1537,3 +1538,382 @@ def test_plot_honest_uses_stored_container_intervals(self): # 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) + + +# --------------------------------------------------------------------------- # +# StackedDiD-sourced container admission (row M-024, second M-093 pre-cut) +# --------------------------------------------------------------------------- # + +STACKED_FIT_KW = dict(outcome="outcome", unit="unit", time="period", first_treat="first_treat") + + +def _stacked_panel(seed=42, n_units=120, n_periods=12, cohorts=(4, 6, 8)): + from diff_diff.prep_dgp import generate_staggered_data + + return generate_staggered_data( + n_units=n_units, n_periods=n_periods, cohort_periods=list(cohorts), seed=seed + ) + + +def _fit_stacked(data, **est_kw): + from diff_diff import StackedDiD + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return StackedDiD(**est_kw).fit(data, **STACKED_FIT_KW) + + +@pytest.fixture(scope="module") +def stacked_surface(): + """kappa_pre=3 surface: two estimated pre-periods for honest/pretrends.""" + res = _fit_stacked(_stacked_panel(), kappa_pre=3, kappa_post=2) + return res.aggregate("event_study") + + +class TestStackedContainerAdmission: + def test_honest_end_to_end(self, stacked_surface): + with warnings.catch_warnings(): + warnings.simplefilter("error") + h = compute_honest_did(stacked_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) + assert h.ci_lb <= h.lb <= h.ub <= h.ci_ub + + def test_honest_anticipation_fit_behavior(self): + # Behavior pin, not a threading-mechanism test: honest's container + # split runs through the MATERIALIZED reference row (numerically the + # same partition as the anticipation cutoff since ref = -1 - k). + res = _fit_stacked(_stacked_panel(), kappa_pre=3, kappa_post=2, anticipation=1) + surf = res.aggregate("event_study") + assert surf.reference_event_times == (-2,) + with warnings.catch_warnings(): + warnings.simplefilter("error") + h = compute_honest_did(surf, M=0.5) + assert np.isfinite(h.lb) and np.isfinite(h.ub) + + def test_pretrends_end_to_end_full_vcov(self, stacked_surface): + with warnings.catch_warnings(): + warnings.simplefilter("error") + p = compute_pretrends_power(stacked_surface) + assert np.isfinite(p.power) and 0.0 < p.power <= 1.0 + assert np.isfinite(p.mdv) and p.mdv > 0 + # StackedDiD persists its ES VCV in every inference mode, so the + # container always takes the full-covariance tier. + assert p.covariance_source == "full_pre_period_vcov" + + def test_pretrends_anticipation_threading_mechanism(self): + # THIS is where the surface.anticipation channel is exercised: + # pretrends reads it directly for the pre cutoff (t < -k), so with + # anticipation=1 the e=-1 anticipation row must be excluded from + # the pre set while e=-4,-3 stay. + res = _fit_stacked(_stacked_panel(), kappa_pre=3, kappa_post=2, anticipation=1) + surf = res.aggregate("event_study") + assert surf.anticipation == 1 + pt = PreTrendsPower() + effects, ses, vcov, n_pre, rel_times, cov_src = pt._extract_pre_period_params(surf) + # n_pre == 2 proves the e=-1 anticipation row was excluded by the + # t < -k cutoff (with cutoff 0 it would be 3). rel_times are the + # ROTH-ANCHORED offsets t - t_ref (#744 convention): raw pre labels + # {-4, -3} anchored at the ref -2 give [-2, -1]. + assert n_pre == 2 + assert sorted(rel_times.tolist()) == [-2.0, -1.0] + assert cov_src == "full_pre_period_vcov" + + def test_pretrends_last_period_positional_selection(self): + # kappa_pre >= 3 so the eligible pre set has TWO horizons and the + # last_period pattern actually selects one (at kappa_pre=2 the + # weight vector degenerates to [1.] and the test would be vacuous). + res = _fit_stacked(_stacked_panel(), kappa_pre=3, kappa_post=2) + surf = res.aggregate("event_study") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + p_last = compute_pretrends_power(surf, violation_type="last_period") + pt = PreTrendsPower(violation_type="last_period") + effects, ses, vcov, n_pre, rel_times, _ = pt._extract_pre_period_params(surf) + assert n_pre == 2 + w = pt._get_violation_weights(n_pre, rel_times) + # chronologically sorted pre rows [-3, -2]: the LAST pre-period + # (closest to treatment, e=-2) carries the violation mass. + assert w.tolist() == [0.0, 1.0] + assert np.isfinite(p_last.power) + + def test_kappa_pre_1_default_rejected(self): + res = _fit_stacked(_stacked_panel(), kappa_pre=1, kappa_post=2) + surf = res.aggregate("event_study") + with pytest.raises(ValueError, match="No pre-period effects"): + compute_honest_did(surf, M=0.5) + with pytest.raises(ValueError, match="No pre-treatment periods"): + compute_pretrends_power(surf) + + # ---- withheld-inference warning (uniform shape-based predicate) ---- # + + def test_bm_failure_fit_warns_in_both_consumers(self, monkeypatch): + # MANDATORY real-shape fixture: monkeypatched BM-DOF failure on an + # hc2_bm fit produces the joint-NaN rows (t/p/CI all NaN, finite + # se) the estimator actually emits - not a hand-mutated container. + import diff_diff.linalg as dl + + monkeypatch.setattr( + dl, + "_compute_cr2_bm_contrast_dof", + lambda X, cluster_ids, bread_matrix, contrasts, weights=None: np.full( + contrasts.shape[1], np.nan + ), + ) + res = _fit_stacked( + _stacked_panel(), + kappa_pre=3, + kappa_post=2, + vcov_type="hc2_bm", + cluster="unit", + ) + eff = res.event_study_effects + assert any( + np.isnan(v["t_stat"]) + and np.isnan(v["p_value"]) + and np.isfinite(v["se"]) + and np.isnan(v["conf_int"][0]) + and np.isnan(v["conf_int"][1]) + for h, v in eff.items() + if v["n_obs"] > 0 + ) + surf = res.aggregate("event_study") + with pytest.warns(UserWarning, match="withheld/undefined"): + h = compute_honest_did(surf, M=0.5) + assert np.isfinite(h.lb) and np.isfinite(h.ub) + with pytest.warns(UserWarning, match="withheld/undefined"): + p = compute_pretrends_power(surf) + assert np.isfinite(p.power) + + def test_replicate_undefined_state_warns_in_both(self, stacked_surface): + # df_survey=0.0 (replicate-undefined sentinel) + withheld rows: + # BOTH consumers warn (uniform predicate - no df conjunct), and + # honest's identified-set bounds stay FINITE while its FLCI CI + # endpoints are NaN (the executed-verified output split). + import dataclasses + + pv = stacked_surface.p_value.copy() + pv[~stacked_surface.is_reference] = np.nan + ts = stacked_surface.t_stat.copy() + ts[~stacked_surface.is_reference] = np.nan + cl = stacked_surface.conf_int_lower.copy() + cu = stacked_surface.conf_int_upper.copy() + cl[~stacked_surface.is_reference] = np.nan + cu[~stacked_surface.is_reference] = np.nan + surf = dataclasses.replace( + stacked_surface, + p_value=pv, + t_stat=ts, + conf_int_lower=cl, + conf_int_upper=cu, + df_survey=0.0, + ) + with pytest.warns(UserWarning, match="withheld/undefined"): + h = compute_honest_did(surf, M=0.5) + assert np.isfinite(h.lb) and np.isfinite(h.ub) + assert np.isnan(h.ci_lb) and np.isnan(h.ci_ub) + with pytest.warns(UserWarning, match="withheld/undefined"): + p = compute_pretrends_power(surf) + assert np.isfinite(p.power) + + def test_positive_finite_df_survey_still_warns(self, stacked_surface): + # Third state of the "every df_survey state warns" contract - a + # regression guard against re-introducing a df-conditioned + # predicate (superseded round-5 draft). + import dataclasses + + pv = stacked_surface.p_value.copy() + pv[1] = np.nan + surf = dataclasses.replace(stacked_surface, p_value=pv, df_survey=25.0) + with pytest.warns(UserWarning, match="withheld/undefined"): + compute_honest_did(surf, M=0.5) + with pytest.warns(UserWarning, match="withheld/undefined"): + compute_pretrends_power(surf) + + def test_cs_twin_does_not_get_withheld_warning(self, stacked_surface): + import dataclasses + + pv = stacked_surface.p_value.copy() + pv[1] = np.nan + surf = dataclasses.replace(stacked_surface, p_value=pv, source="CallawaySantAnnaResults") + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + compute_honest_did(surf, M=0.5) + compute_pretrends_power(surf) + assert not any("withheld" in str(w.message) for w in caught) + + # ---- source-flip parity (admission is a pure gate) ---- # + + def test_source_flip_parity(self, stacked_surface): + import dataclasses + + cs_surface = dataclasses.replace(stacked_surface, source="CallawaySantAnnaResults") + h_st = compute_honest_did(stacked_surface, M=0.5) + h_cs = compute_honest_did(cs_surface, M=0.5) + np.testing.assert_allclose( + [h_st.lb, h_st.ub, h_st.ci_lb, h_st.ci_ub], + [h_cs.lb, h_cs.ub, h_cs.ci_lb, h_cs.ci_ub], + rtol=1e-12, + ) + pt = PreTrendsPower() + ext_st = pt._extract_pre_period_params(stacked_surface) + ext_cs = pt._extract_pre_period_params(cs_surface) + for a, b in zip(ext_st[:5], ext_cs[:5]): + np.testing.assert_array_equal(np.asarray(a), np.asarray(b)) + assert ext_st[5] == ext_cs[5] + # end-to-end power at the suite's stochastic tolerance (module + # docstring: scipy's Genz MVN CDF is internally randomized) + p_st = compute_pretrends_power(stacked_surface, M=0.1) + p_cs = compute_pretrends_power(cs_surface, M=0.1) + assert abs(p_st.power - p_cs.power) < 1e-3 + + # ---- guard-message reach tests (seven producer-derived sites) ---- # + + def test_honest_missing_provenance_warning_names_producer(self, stacked_surface): + import dataclasses + + surf = dataclasses.replace(stacked_surface, reference_event_times=None) + assert surf.base_period == "universal" + with pytest.warns(UserWarning, match="StackedDiD event-study container carries no"): + compute_honest_did(surf, M=0.5) + + def test_honest_multi_reference_times_error_names_producer(self, stacked_surface): + import dataclasses + + surf = dataclasses.replace(stacked_surface, reference_event_times=(-3, -1)) + with pytest.raises(ValueError, match="StackedDiD container records DISTINCT"): + compute_honest_did(surf, M=0.5) + + def test_honest_varying_base_warning_producer_conditional_remedy(self, stacked_surface): + import dataclasses + + surf = dataclasses.replace( + stacked_surface, base_period="varying", reference_event_times=None + ) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + compute_honest_did(surf, M=0.5) + msgs = " ".join(str(w.message) for w in caught) + assert "Rebuild the container" in msgs + assert "StackedDiD results" in msgs + assert "CallawaySantAnna(base_period='universal')" not in msgs + + def test_honest_multi_reference_rows_error(self, stacked_surface): + import dataclasses + + is_ref = stacked_surface.is_reference.copy() + # mark the earliest pre row as a second reference (att 0, se 0 so + # the zero-count reference convention holds shape-wise) + att = stacked_surface.att.copy() + se = stacked_surface.se.copy() + att[0], se[0] = 0.0, 0.0 + is_ref[0] = True + surf = dataclasses.replace( + stacked_surface, + is_reference=is_ref, + att=att, + se=se, + reference_event_times=(-1,), + ) + with pytest.raises(ValueError, match="multiple reference rows"): + compute_honest_did(surf, M=0.5) + + def test_pretrends_missing_provenance_warning_names_producer(self, stacked_surface): + import dataclasses + + surf = dataclasses.replace(stacked_surface, reference_event_times=None) + with pytest.warns(UserWarning, match="StackedDiD event-study container carries no"): + compute_pretrends_power(surf) + + def test_pretrends_multi_reference_times_error_names_producer(self, stacked_surface): + import dataclasses + + surf = dataclasses.replace(stacked_surface, reference_event_times=(-3, -1)) + with pytest.raises(ValueError, match="StackedDiD container records DISTINCT"): + compute_pretrends_power(surf) + + def test_pretrends_varying_base_warning_producer_conditional_remedy(self, stacked_surface): + import dataclasses + + surf = dataclasses.replace( + stacked_surface, base_period="varying", reference_event_times=None + ) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + compute_pretrends_power(surf) + msgs = " ".join(str(w.message) for w in caught) + assert "Rebuild the container" in msgs + assert "PreTrendsPower on a StackedDiD event-study container" in msgs + + # ---- fail-closed structural paths ---- # + + def test_rank_drop_shape_gap_fail_closed_producer_neutral(self, stacked_surface): + # FIRST pin of honest's retained-grid gap ValueError anywhere in + # the suite: a NaN-se INTERIOR pre row (the rank-drop shape) drops + # from the retained set and the grid gains a hole. The message must + # carry the producer-neutral Stacked remedy, not balance_e (which + # Stacked's aggregate() does not have). + import dataclasses + + se = stacked_surface.se.copy() + # interior pre row (e=-2 on the kappa_pre=3 grid [-3..-2] pre set + # around ref -1): NaN it so pre keeps only -3 next to post 0..2 - + # wait, the gap forms between -3 and the reference split; the + # retained PRE block {-3} is contiguous, so NaN a POST interior + # row (e=1) instead: post block {0, 2} has a gap. + idx = stacked_surface.event_time.tolist().index(1) + se[idx] = np.nan + surf = dataclasses.replace(stacked_surface, se=se) + with pytest.raises(ValueError, match="rank-dropped event-time column"): + compute_honest_did(surf, M=0.5) + + def test_singular_pre_covariance_honest_rejects_pretrends_accepts(self): + # Low-G / high-kappa_pre: the pre-period covariance sub-block is + # singular (vcov rank <= G). Honest validates with + # allow_singular=False (Rambachan-Roth eigenvalues bounded away + # from zero) and REJECTS; pretrends keeps its documented singular + # support and succeeds. + panel = _stacked_panel(n_units=12, n_periods=30, cohorts=(17, 19)) + res = _fit_stacked(panel, kappa_pre=12, kappa_post=2, cluster="unit") + surf = res.aggregate("event_study") + with pytest.raises(ValueError, match="singular"): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + compute_honest_did(surf, M=0.5) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + p = compute_pretrends_power(surf) + assert np.isfinite(p.power) + + # ---- provenance pins ---- # + + def test_provenance_ladder(self, stacked_surface): + from diff_diff.survey import SurveyDesign + + # analytical -> df_survey None (deliberate: the resolver reads + # df_inference, a name Stacked does not use - decision 3; the + # TODO adapter-naming row must preserve this) + assert stacked_surface.base_period == "universal" + assert stacked_surface.reference_event_times == (-1,) + assert stacked_surface.df_survey is None + from diff_diff.results_base import _resolve_scalar_df_survey + + panel = _stacked_panel() + res = _fit_stacked(panel, kappa_pre=3, kappa_post=2) + assert _resolve_scalar_df_survey(res) is None + # survey TSL -> the survey df threads + spanel = panel.copy() + spanel["w"] = 1.0 + 0.1 * (spanel["unit"] % 5) + spanel["strata"] = spanel["unit"] % 4 + spanel["psu"] = spanel["unit"] + from diff_diff import StackedDiD + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + sres = StackedDiD(kappa_pre=3, kappa_post=2).fit( + spanel, + survey_design=SurveyDesign(weights="w", strata="strata", psu="psu"), + **STACKED_FIT_KW, + ) + ssurf = sres.aggregate("event_study") + assert ssurf.df_survey == float(sres.survey_metadata.df_survey) diff --git a/tests/test_event_study_surface.py b/tests/test_event_study_surface.py index 32cdc4103..b8793da06 100644 --- a/tests/test_event_study_surface.py +++ b/tests/test_event_study_surface.py @@ -353,7 +353,6 @@ def surfaces(): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) out["StackedDiD"] = (st, build_event_study_surface(st)) @@ -1366,7 +1365,6 @@ def test_stacked_hc2_bm_per_row_df(): unit="unit", time="time", first_treat="first_treat", - aggregate="event_study", ) with warnings.catch_warnings(): warnings.simplefilter("ignore") @@ -1409,7 +1407,6 @@ def test_stacked_survey_df_broadcast(): unit="unit", time="time", first_treat="first_treat", - aggregate="event_study", survey_design=design, ) assert res.survey_metadata is not None and res.survey_metadata.df_survey is not None @@ -1481,7 +1478,6 @@ def test_stacked_replicate_weight_vcov_diag_matches_ses(): unit="unit", time="time", first_treat="first_treat", - aggregate="event_study", survey_design=design, ) assert res.event_study_vcov is not None diff --git a/tests/test_guides.py b/tests/test_guides.py index 556253c70..bc934a6be 100644 --- a/tests/test_guides.py +++ b/tests/test_guides.py @@ -794,6 +794,21 @@ def test_llms_full_stacked_fit_documents_covariates(self): "StackedDiD.fit() exposes covariates= but the llms-full.txt fit() block " "does not document it." ) + + def test_llms_full_stacked_fit_aggregate_line_documents_shim(self): + # M-024: the documented fit signature must carry the sentinel + # default + deprecation marker (the dCDH precedent's form), not + # the stale `aggregate: str = None` - no other pin covers this + # line, so it could silently go stale. + section = self._stacked_section() + fit_start = section.index("stacked.fit(") + fit_block = section[fit_start : section.index("\n)", fit_start)] + agg_line = next( + line for line in fit_block.splitlines() if line.strip().startswith("aggregate") + ) + assert "NOT_SUPPLIED" in agg_line + assert "DEPRECATED (M-024)" in agg_line + assert "results.aggregate()" in agg_line # balance= must be documented somewhere in the section (constructor param) assert "balance" in section diff --git a/tests/test_methodology_stacked_did.py b/tests/test_methodology_stacked_did.py index c9502462b..acc5f29d6 100644 --- a/tests/test_methodology_stacked_did.py +++ b/tests/test_methodology_stacked_did.py @@ -154,7 +154,6 @@ def _fit(panel, vcov_type, cluster): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) @@ -431,7 +430,6 @@ def test_population_weighting_hc1_matches_clubsandwich_cr1s(self, goldens, popul time="period", first_treat="first_treat", population="pop", - aggregate="event_study", ) py_se = _es_se_vector(res, event_times) np.testing.assert_allclose( @@ -462,7 +460,6 @@ def test_population_weighting_hc2_bm_matches_clubsandwich_cr2(self, goldens, pop time="period", first_treat="first_treat", population="pop", - aggregate="event_study", ) py_se = _es_se_vector(res, event_times) np.testing.assert_allclose( @@ -520,7 +517,6 @@ def test_anticipation1_hc1_matches_clubsandwich_cr1s(self, goldens, panel): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) py_se = _es_se_vector(res, event_times) np.testing.assert_allclose( @@ -550,7 +546,6 @@ def test_anticipation1_hc2_bm_matches_clubsandwich_cr2(self, goldens, panel): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) py_se = _es_se_vector(res, event_times) np.testing.assert_allclose( @@ -609,7 +604,6 @@ def test_sample_share_weighting_hc1_matches_clubsandwich_cr1s(self, goldens, pan unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) py_se = _es_se_vector(res, event_times) np.testing.assert_allclose( @@ -642,7 +636,6 @@ def test_sample_share_weighting_hc2_bm_matches_clubsandwich_cr2(self, goldens, p unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) py_se = _es_se_vector(res, event_times) np.testing.assert_allclose( @@ -777,7 +770,6 @@ def _fit(self, df, balance="entropy"): unit="unit", time="time", first_treat="first_treat", - aggregate="event_study", covariates=["x"] if balance != "none" else None, ) @@ -936,7 +928,6 @@ def _fit(self, df): unit="unit", time="time", first_treat="first_treat", - aggregate="event_study", covariates=["x"], ) @@ -1008,7 +999,6 @@ def test_dynamic_atts_match_r_cbwsdid(self): unit="unit", time="time", first_treat="first_treat", - aggregate="event_study", covariates=["x"], ) for et, r_est, r_se in zip( diff --git a/tests/test_practitioner.py b/tests/test_practitioner.py index 2545b7e30..2faf79c03 100644 --- a/tests/test_practitioner.py +++ b/tests/test_practitioner.py @@ -330,6 +330,29 @@ def test_stacked_results(self, mock_stacked_results): all_text = " ".join(s.get("code", "") + s.get("why", "") for s in output["next_steps"]) assert "not_yet_treated" not in all_text or "control_group" in all_text + def test_stacked_balance_step_uses_distinct_step_name(self, mock_stacked_results): + # M-024: "Check sub-experiment balance" must NOT reuse + # step_name="heterogeneity" - since the ES surface is always + # populated, DiagnosticReport's heterogeneity check runs on every + # plain fit and a shared key silently dropped this unrelated + # advice from next_steps via _filter_steps. + from diff_diff.practitioner import STEPS, _handle_stacked + + steps, _ = _handle_stacked(mock_stacked_results) + balance = [s for s in steps if s["label"] == "Check sub-experiment balance"] + assert len(balance) == 1 + assert balance[0]["_step_name"] == "sub_experiment_balance" + # SURVIVAL: completing heterogeneity (what DiagnosticReport does on + # every surface-populated fit) must no longer drop the balance step. + completed = practitioner_next_steps( + mock_stacked_results, verbose=False, completed_steps=["heterogeneity"] + ) + labels = [s["label"] for s in completed["next_steps"]] + assert "Check sub-experiment balance" in labels + # Like "loo_jackknife", the key deliberately stays OUT of the STEPS + # completion vocabulary: no diagnostic ever completes it. + assert "sub_experiment_balance" not in STEPS + def test_synth_results(self, mock_synth_results): output = practitioner_next_steps(mock_synth_results, verbose=False) assert len(output["next_steps"]) > 0 diff --git a/tests/test_stacked_did.py b/tests/test_stacked_did.py index ff23cd861..e8cfbdf17 100644 --- a/tests/test_stacked_did.py +++ b/tests/test_stacked_did.py @@ -92,7 +92,6 @@ def test_event_study(self, staggered_data): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) assert results.event_study_effects is not None assert -1 in results.event_study_effects # reference period @@ -107,33 +106,35 @@ def test_event_study(self, staggered_data): assert results.event_study_effects[h]["n_obs"] > 0 def test_group_aggregate_raises(self, staggered_data): - """aggregate='group' raises ValueError.""" + """aggregate='group' warns (M-024 shim) then raises ValueError.""" est = StackedDiD(kappa_pre=2, kappa_post=2) - with pytest.raises(ValueError, match="group.*not supported"): - est.fit( - staggered_data, - outcome="outcome", - unit="unit", - time="period", - first_treat="first_treat", - aggregate="group", - ) + with pytest.warns(FutureWarning, match=r"fit\(aggregate=\) is deprecated"): + with pytest.raises(ValueError, match="group.*not supported"): + est.fit( + staggered_data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + aggregate="group", + ) def test_all_aggregate_raises(self, staggered_data): - """aggregate='all' raises ValueError.""" + """aggregate='all' warns (M-024 shim) then raises ValueError.""" est = StackedDiD(kappa_pre=2, kappa_post=2) - with pytest.raises(ValueError, match="all.*not supported"): - est.fit( - staggered_data, - outcome="outcome", - unit="unit", - time="period", - first_treat="first_treat", - aggregate="all", - ) + with pytest.warns(FutureWarning, match=r"fit\(aggregate=\) is deprecated"): + with pytest.raises(ValueError, match="all.*not supported"): + est.fit( + staggered_data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + aggregate="all", + ) def test_simple_att(self, staggered_data): - """aggregate='simple' produces overall ATT only.""" + """Plain fit: overall ATT plus the always-computed surface (M-024).""" est = StackedDiD(kappa_pre=2, kappa_post=2) results = est.fit( staggered_data, @@ -141,10 +142,10 @@ def test_simple_att(self, staggered_data): unit="unit", time="period", first_treat="first_treat", - aggregate="simple", ) assert np.isfinite(results.overall_att) - assert results.event_study_effects is None + # The event-study surface is ALWAYS materialized on 3.9+ fits. + assert results.event_study_effects is not None assert results.group_effects is None def test_known_constant_effect(self, constant_effect_data): @@ -171,7 +172,6 @@ def test_dynamic_effects(self, staggered_data): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) assert results.event_study_effects is not None # Post-treatment effects should generally increase @@ -566,7 +566,6 @@ def test_anticipation_reference_period(self): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) # Reference period is -2 (not -1) @@ -761,7 +760,6 @@ def test_summary(self, staggered_data): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) summary = results.summary() assert "Stacked DiD" in summary @@ -776,7 +774,6 @@ def test_to_dataframe_event_study(self, staggered_data): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) df = results.to_dataframe(level="event_study") assert isinstance(df, pd.DataFrame) @@ -796,8 +793,37 @@ def test_to_dataframe_group_raises(self, staggered_data): with pytest.raises(ValueError, match="Group aggregation is not supported"): results.to_dataframe(level="group") - def test_to_dataframe_no_event_study_raises(self, staggered_data): - """to_dataframe raises when event_study not computed.""" + def test_to_dataframe_event_study_total_on_plain_fit(self, staggered_data): + """Plain fits carry the surface (M-024): to_dataframe is total.""" + est = StackedDiD(kappa_pre=2, kappa_post=2) + results = est.fit( + staggered_data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + ) + df = results.to_dataframe(level="event_study") + assert sorted(df["relative_period"].tolist()) == [-2, -1, 0, 1, 2] + + def test_to_dataframe_legacy_pickle_none_branch(self, staggered_data): + """The effects-None branch survives for pre-3.9 pickles only.""" + import dataclasses + + est = StackedDiD(kappa_pre=2, kappa_post=2) + results = est.fit( + staggered_data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + ) + legacy = dataclasses.replace(results, event_study_effects=None) + with pytest.raises(ValueError, match=r"diff-diff >= 3\.9"): + legacy.to_dataframe(level="event_study") + + def test_summary_prints_event_study_table_on_plain_fit(self, staggered_data): + """summary() renders the ES block on every 3.9+ fit (M-024).""" est = StackedDiD(kappa_pre=2, kappa_post=2) results = est.fit( staggered_data, @@ -806,8 +832,7 @@ def test_to_dataframe_no_event_study_raises(self, staggered_data): time="period", first_treat="first_treat", ) - with pytest.raises(ValueError, match="Event study effects not computed"): - results.to_dataframe(level="event_study") + assert "Event Study (Dynamic) Effects" in results.summary() def test_is_significant(self, staggered_data): """is_significant property works.""" @@ -888,17 +913,18 @@ def test_invalid_cluster(self): StackedDiD(cluster="invalid") def test_invalid_aggregate(self, staggered_data): - """Raises on invalid aggregate parameter.""" + """Warns (M-024 shim) then raises on invalid aggregate parameter.""" est = StackedDiD() - with pytest.raises(ValueError, match="aggregate"): - est.fit( - staggered_data, - outcome="outcome", - unit="unit", - time="period", - first_treat="first_treat", - aggregate="invalid", - ) + with pytest.warns(FutureWarning, match=r"fit\(aggregate=\) is deprecated"): + with pytest.raises(ValueError, match="aggregate"): + est.fit( + staggered_data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + aggregate="invalid", + ) def test_population_required_for_population_weighting(self, staggered_data): """Raises when population col not specified with weighting='population'.""" @@ -1030,7 +1056,6 @@ def test_hc2_bm_finite_and_att_identical_to_hc1(self, staggered_data): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) res_hc1 = StackedDiD(kappa_pre=2, kappa_post=2, vcov_type="hc1").fit( staggered_data, **kwargs @@ -1319,7 +1344,6 @@ def test_hc2_bm_uses_t_distribution_not_normal_on_well_conditioned_design(self, unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) # No fallback warning on a well-conditioned design. fallback_warns = [w for w in caught if "Bell-McCaffrey contrast DOF" in str(w.message)] @@ -1414,7 +1438,6 @@ def _fake_solve_ols(*args, **kwargs): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) # The reduced-design code path should fire without emitting the # fallback warning (the helper succeeds on the kept-column subset). @@ -1480,7 +1503,6 @@ def _fake_contrast_dof(*args, **kwargs): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) # All event_study_effects (except the ref period) should have NaN # inference but finite effect+se. @@ -1528,7 +1550,6 @@ def _fake_raises(*args, **kwargs): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) # Warning fired (informational that DOF was unavailable) warning_msgs = [str(w.message) for w in caught] @@ -1565,7 +1586,6 @@ def test_anticipation_plus_hc2_bm_threads_bm_dof(self, staggered_data): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) fallback_warns = [w for w in caught if "Bell-McCaffrey contrast DOF" in str(w.message)] assert ( @@ -1618,7 +1638,6 @@ def test_population_weighting_plus_hc2_bm_finite_threads_bm_dof(self, staggered_ time="period", first_treat="first_treat", population="pop", - aggregate="event_study", ) assert np.isfinite(res.overall_att) and np.isfinite(res.overall_se) assert res.overall_se > 0 @@ -1644,7 +1663,6 @@ def test_sample_share_weighting_plus_hc2_bm_finite_threads_bm_dof(self, staggere unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) assert np.isfinite(res.overall_att) and np.isfinite(res.overall_se) assert res.overall_se > 0 @@ -1663,7 +1681,6 @@ def test_hc1_vs_hc2_bm_differ_under_anticipation(self, staggered_data): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) res_hc1 = StackedDiD(kappa_pre=2, kappa_post=2, vcov_type="hc1", anticipation=1).fit( staggered_data, **kwargs @@ -1784,7 +1801,6 @@ def _cb_fit(df, **kw): unit="unit", time="time", first_treat="first_treat", - aggregate=kw.get("aggregate", "event_study"), covariates=kw.get("covariates"), survey_design=kw.get("survey_design"), population=kw.get("population"), @@ -1869,13 +1885,15 @@ def test_covariate_scale_invariance(self): b = _cb_fit(df, balance="entropy", covariates=["x2"]) assert abs(a.overall_att - b.overall_att) < 1e-8 - def test_balance_on_default_aggregate_mode(self): - # balancing must work on the overall-ATT (aggregate=None) path, not only - # event_study (the BM-DOF machinery branches on aggregate). + def test_balance_on_plain_fit(self): + # balancing works on a plain fit; since M-024 the surface (and the + # hc2_bm BM-DOF machinery) no longer branches on fit-time + # aggregate, so this is the same code path as every _cb_fit case. df = _balance_panel(seed=8) - res = _cb_fit(df, balance="entropy", covariates=["x"], aggregate=None) + res = _cb_fit(df, balance="entropy", covariates=["x"]) assert np.isfinite(res.overall_att) assert res.balance == "entropy" + assert res.event_study_effects is not None # ---- sklearn surface ---- def test_get_params_includes_balance_and_clone(self): @@ -1894,8 +1912,7 @@ def test_convenience_function_threads_covariates(self): "first_treat", kappa_pre=2, kappa_post=2, - aggregate="event_study", - clean_control="never_treated", + control_group="never_treated", balance="entropy", covariates=["x"], ) @@ -1958,7 +1975,6 @@ def _panel(seed=9, n_units=60): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", ) @staticmethod @@ -2073,3 +2089,28 @@ def test_validation_and_transactional_set_params(self): est.set_params(df_convention="normal", nonexistent_param=1) assert est.get_params() == before assert StackedDiD(df_convention="normal").get_params()["df_convention"] == "normal" + + +class TestPlainFitPlotting: + def test_plot_event_study_renders_on_plain_fit(self, staggered_data): + """M-024: plot_event_study flips from TypeError to rendering on a + plain fit - the duck-typed effects branch now matches because the + surface is always populated.""" + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + from diff_diff import plot_event_study + + est = StackedDiD(kappa_pre=2, kappa_post=2) + results = est.fit( + staggered_data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + ) + ax = plot_event_study(results, show=False) + assert ax is not None + plt.close("all") diff --git a/tests/test_survey_estimator_validation.py b/tests/test_survey_estimator_validation.py index d80bb3c12..ba533345d 100644 --- a/tests/test_survey_estimator_validation.py +++ b/tests/test_survey_estimator_validation.py @@ -217,7 +217,6 @@ def test_s2_stacked_did(self, golden): unit="unit", time="period", first_treat="first_treat", - aggregate="event_study", survey_design=sd, ) diff --git a/tests/test_variance_conventions.py b/tests/test_variance_conventions.py index 2a8f5b7b5..5c9f9c119 100644 --- a/tests/test_variance_conventions.py +++ b/tests/test_variance_conventions.py @@ -274,7 +274,11 @@ def snapshot(self): df, outcome="y", unit="unit", time="time", first_treat="first_treat" ), cr1_k=(6,), - tail_df=(309.0,) * 1, + # 3 safe_inference calls since M-024: the event-study surface is + # ALWAYS materialized (2 estimated event times at the default + # kappa 1/1 grid) alongside the overall ATT, all at the same + # pooled residual df. + tail_df=(309.0,) * 3, status="legitimate", reason=( "L1: k_total is clubSandwich CR1S by construction (stacked_did.py " From a4d6910b490d810a29e724cf7163e24bed14a49d Mon Sep 17 00:00:00 2001 From: igerber Date: Mon, 3 Aug 2026 15:56:23 -0400 Subject: [PATCH 2/4] fix(v4): fail closed on an unobserved stacked reference cell; covariance-wording fix (CI review R1) P1 (verified by execution): on a gapped panel where every retained cohort's calendar period a-1-anticipation is absent, the omitted reference cell is empty, the design rank-drops, and QR pivoting silently re-normalizes against an arbitrary horizon while the surface synthesizes the e=-1-k reference row and the container certifies reference_event_times/base_period="universal" - a fabricated delta_0=0 claim HonestDiD/PreTrendsPower would trust. fit() now raises a fail-closed ValueError when the reference cell lacks treated or control rows stack-wide (both cells are load-bearing: either absence collinears a design block). Partial support - some cohorts' own ref period absent while the pooled omitted category stays populated - still fits, verified by execution: the coefficient-vector normalization is real and the per-cohort raggedness is the documented ragged-window estimand composition, so no legitimate fit changes. Gapped-panel tests pin both sides; REGISTRY/CHANGELOG document the guard. P2: the singularity requirement wording understated the check - honest's container branch validates the FULL retained event-study covariance (pre+post sub-block, Rambachan-Roth Assumption 3), not just the pre-period block; reworded in the docstrings, REGISTRY, ledger notes, and CHANGELOG. --- CHANGELOG.md | 13 +++++-- diff_diff/honest_did.py | 6 ++-- diff_diff/stacked_did.py | 39 +++++++++++++++++++-- docs/methodology/REGISTRY.md | 4 +-- docs/v4-deprecations.yaml | 2 +- tests/test_stacked_did.py | 68 ++++++++++++++++++++++++++++++++++++ 6 files changed, 121 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d80eb7d4f..ae9e9ed78 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,8 +47,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - `compute_honest_did` and `compute_pretrends_power` now ADMIT Stacked-sourced containers ([M-093] second pre-cut amendment): requires `kappa_pre >= 2` (the default grid has no estimated - pre-periods) and, for HonestDiD, a non-singular pre-period - covariance (keep `kappa_pre` small relative to the cluster count). + pre-periods) and, for HonestDiD, a non-singular FULL retained + event-study covariance - the pre+post sub-block, per + Rambachan-Roth Assumption 3 (keep `kappa_pre` small relative to + the cluster count). Analytical Stacked containers carry `df_survey=None` - honest FLCI critical values are normal-theory there (deliberate; REGISTRY Note). Rows whose per-row inference the producer withheld (finite @@ -59,7 +61,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 `0.0` sentinel, honest's identified-set bounds stay finite while its FLCI CI endpoints are NaN). The seven producer-derived guard messages in both container branches now derive the producer from - `surface.source` with producer-conditional remedies. + `surface.source` with producer-conditional remedies. A new + reference-support guard fails `fit()` closed when the omitted + reference cell has no treated or control rows stack-wide (gapped + panel at every cohort's `a - 1 - anticipation`) - previously this + rank-dropped and silently re-normalized while the surface + certified a `delta_0 = 0` reference that never existed. - Cross-producer container hardening: `EventStudyResults`'s `__post_init__` now COPIES `vcov`/`vcov_index` like every other array field (`np.asarray` aliased the producer's stored matrix on diff --git a/diff_diff/honest_did.py b/diff_diff/honest_did.py index 594adf2af..9cc953eaa 100644 --- a/diff_diff/honest_did.py +++ b/diff_diff/honest_did.py @@ -586,7 +586,8 @@ def _extract_container_params( ``CallawaySantAnnaResults.aggregate('event_study')`` or ``StackedDiDResults.aggregate('event_study')`` (row M-024; Stacked containers require ``kappa_pre >= 2`` so estimated pre-periods - exist, and a pre-period covariance that is not singular - honest + exist, and a non-singular FULL retained event-study covariance + (the pre+post sub-block, per Rambachan-Roth's Assumption 3) - honest validates with ``allow_singular=False``, so keep ``kappa_pre`` small relative to the cluster count). Admission is SOURCE-SCOPED: containers from other producers are rejected rather than silently @@ -2883,7 +2884,8 @@ def fit( ``CallawaySantAnnaResults.aggregate('event_study')`` or ``StackedDiDResults.aggregate('event_study')`` (Stacked containers require ``kappa_pre >= 2`` so estimated - pre-periods exist, and a non-singular pre-period covariance - + pre-periods exist, and a non-singular FULL retained event-study + covariance (the pre+post sub-block) - keep ``kappa_pre`` small relative to the cluster count). On the container route the scalar inference df arrives via the container's ``df_survey`` provenance field (for CS-sourced diff --git a/diff_diff/stacked_did.py b/diff_diff/stacked_did.py index 17905e439..3c7a22b57 100644 --- a/diff_diff/stacked_did.py +++ b/diff_diff/stacked_did.py @@ -645,14 +645,47 @@ def fit( # [K+2..2K+1] D_sa * event-time interactions interaction_indices: Dict[int, int] = {} + et_vals = stacked_df["_event_time"].values + d_vals = stacked_df["_D_sa"].values + + # ---- Reference-support guard (row M-024 follow-up) ---- + # The omitted reference event time e = -1 - anticipation is the + # regression's baseline category: every reported delta_h (and the + # synthesized reference row, and the container's + # reference_event_times / base_period="universal" provenance) is + # defined relative to its cells. On a gapped panel where NO + # retained sub-experiment observes that calendar period, the ref + # cell is empty, the design goes rank-deficient, and QR pivoting + # silently re-normalizes against an arbitrary surviving column - + # the surface would then certify a delta_0 = 0 normalization that + # never happened (and HonestDiD / PreTrendsPower would trust it). + # Both cells are required: without TREATED ref rows, D_sa is + # collinear with the delta block (same silent re-normalization); + # without CONTROL ref rows, the baseline time profile is. Fail + # closed - a fabricated reference must never reach consumers. + _ref_mask = et_vals == ref_period + if not np.any(_ref_mask & (d_vals == 1)) or not np.any(_ref_mask & (d_vals == 0)): + _missing_cell = "treated" if not np.any(_ref_mask & (d_vals == 1)) else "control" + raise ValueError( + f"The omitted reference event time e={ref_period} has no " + f"{_missing_cell} observations in the stacked data, so the " + "event-study normalization (delta_0 = 0 at the reference) " + "would be fabricated: the regression's baseline cell is " + "empty and rank handling would silently re-normalize " + "against an arbitrary horizon. This happens on gapped " + "panels where every retained cohort's calendar period " + "a - 1 - anticipation is absent (the IC1 window check " + "only inspects the panel's min/max periods, not interior " + "gaps). Fill the gap, adjust anticipation/kappa_pre so " + "the reference falls on an observed period, or drop the " + "affected cohorts." + ) + # Build design matrix X = np.zeros((n, 2 + 2 * n_event_dummies)) X[:, 0] = 1.0 # intercept X[:, 1] = stacked_df["_D_sa"].values # treatment indicator - et_vals = stacked_df["_event_time"].values - d_vals = stacked_df["_D_sa"].values - for j, h in enumerate(event_times): col_lambda = 2 + j # event-time dummy col_delta = 2 + n_event_dummies + j # interaction diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 42bf94021..24741de62 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2082,7 +2082,7 @@ The pooled estimator is `DID^{CBWSDID}_e = Σ_a (N^D_a/N^D_Ω)(Δ̄^D_{a,e} − - **Note:** Inference is conditional-on-the-estimated-weights cluster-robust (the existing `hc1`/`hc2_bm` path with `W_{sa}` as the WLS weights) — the paper's default. The paper's weight-re-estimating bootstrap is NOT implemented in v1 (deliberate scope; entropy balancing is smooth so the Abadie–Imbens (2008) nonsmooth-matching bootstrap caveat does not apply). `cluster` is orthogonal to `b_{sa}` (weights conditioned-on); default `unit` matches the paper. - **Note:** v1 scope — only `balance="entropy"` with `weighting="aggregate"`. `balance` + `population`/`sample_share` and `balance` + `survey_design=` raise `NotImplementedError`; matching-based balancing and the repeated `0→1/1→0` episode extension are out of scope. - **Note (post-fit aggregate() - row M-024):** `fit(aggregate=)` is deprecated (3.9; removed 4.0) and the event-study surface is ALWAYS materialized at fit — the pooled regression always included the event-time interactions, so only the extraction and, under `hc2_bm`, the per-event Bell-McCaffrey contrast DOFs were fit-time-gated; both gates are gone and the param is behaviorally inert (its `group`/`all`/unknown ValueErrors persist after the shim warning). `aggregate('event_study')` returns the unified `EventStudyResults` container (a pure view; no recomputation, so every stored inference mode — survey TSL, replicate refit, hc2_bm fail-close — relays faithfully); `aggregate('simple')` relays `overall_att/se/t/p/CI` bit-exactly with `target="att"` (the scalar is the equally-weighted post-period average, NOT the per-event-time trimmed aggregate ATT `θ^e_κ` — `describe_target_parameter` remains the estimand's prose source of truth), `n = n_treated_units` with `n_kind="units"` (TREATED-only scope: under dynamic clean-control definitions a later-treated unit is treated in its own sub-experiment and a clean control in earlier ones, so the treated and control unit sets OVERLAP and no disjoint total exists — cross-container `n` comparisons vs CallawaySantAnna's treated+control convention are out of contract), and `df = inference_df` (the stored overall inference's exact provenance). Consequences of the always-computed surface: plain-fit `summary()` renders the event-study table, `to_dataframe('event_study')` is total on 3.9+ fits, `plot_event_study` renders plain fits, DiagnosticReport's `heterogeneity` (always) and `parallel_trends` (`kappa_pre >= 2`) checks run on plain fits, and the `hc2_bm` overall-ATT contrast now shares the m>1 batch of `_compute_cr2_bm_contrast_dof`, whose degeneracy guard is batch-relative — the noise-floor carve-out on the `Wald_test(test="HTZ")$df_denom` parity claim in the variance-families bullets above widens quantitatively (unchanged behavior when tripped: fail-closed NaN inference plus linalg's loud noise-floor UserWarning; the absolute floor already made this reachable at m=1). -- **Note (container admission - rows M-024/M-093):** `compute_honest_did` and `compute_pretrends_power` admit Stacked-sourced `aggregate('event_study')` containers. Provenance: `base_period` is the class-level constant `"universal"` (every sub-experiment normalizes against the single omitted reference `e = -1 - anticipation`) and `reference_event_times` is the singleton `(-1 - anticipation,)`. Requirements: `kappa_pre >= 2` (the default `kappa_pre=1` grid minus the reference has ZERO estimated pre-periods — both consumers reject), and for HonestDiD a non-singular pre-period covariance (`allow_singular=False` per Rambachan-Roth's eigenvalues-bounded-away-from-zero assumption — keep `kappa_pre` small relative to the cluster count; PreTrendsPower keeps its documented singular support). The container's scalar `df_survey` is None on analytical fits — the generic resolver's fallback reads `df_inference`, a name StackedDiD does not carry (`inference_df` is the overall POST-AVERAGE contrast df, not a design-level cluster df, and misusing it for FLCI critical values would be worse than the asymptotic default) — so honest FLCI critical values are NORMAL-THEORY on analytical Stacked fits, diverging from bare-`cluster=` CS fits (which thread `G-1`); survey fits thread the survey df and replicate-undefined designs the fail-closed `0.0` sentinel (NaN FLCI CIs; the identified-set `lb`/`ub` bounds are df-independent and stay finite). Withheld-inference rows (finite `se`, non-finite `p_value` — the hc2_bm BM-DOF fail-close and replicate-undefined designs both emit this shape) are ADMITTED with a source-scoped warning in BOTH consumers: bounds/power consume only the point estimates and covariance, which remain valid; the warning is shape-descriptive because containers are publicly constructible. A rank-dropped interior event-time column (NaN-se row) fails closed via honest's consecutive-grid ValueError with a Stacked-specific remedy (no `balance_e` — Stacked's aggregate() has no such level). +- **Note (container admission - rows M-024/M-093):** `compute_honest_did` and `compute_pretrends_power` admit Stacked-sourced `aggregate('event_study')` containers. Provenance: `base_period` is the class-level constant `"universal"` (every sub-experiment normalizes against the single omitted reference `e = -1 - anticipation`) and `reference_event_times` is the singleton `(-1 - anticipation,)`. Requirements: `kappa_pre >= 2` (the default `kappa_pre=1` grid minus the reference has ZERO estimated pre-periods — both consumers reject), and for HonestDiD a non-singular FULL retained event-study covariance — the validated sub-block spans the retained pre AND post rows (`allow_singular=False` per Rambachan-Roth's eigenvalues-bounded-away-from-zero assumption on the full Σ — keep `kappa_pre` small relative to the cluster count; PreTrendsPower keeps its documented singular support). The container's scalar `df_survey` is None on analytical fits — the generic resolver's fallback reads `df_inference`, a name StackedDiD does not carry (`inference_df` is the overall POST-AVERAGE contrast df, not a design-level cluster df, and misusing it for FLCI critical values would be worse than the asymptotic default) — so honest FLCI critical values are NORMAL-THEORY on analytical Stacked fits, diverging from bare-`cluster=` CS fits (which thread `G-1`); survey fits thread the survey df and replicate-undefined designs the fail-closed `0.0` sentinel (NaN FLCI CIs; the identified-set `lb`/`ub` bounds are df-independent and stay finite). Withheld-inference rows (finite `se`, non-finite `p_value` — the hc2_bm BM-DOF fail-close and replicate-undefined designs both emit this shape) are ADMITTED with a source-scoped warning in BOTH consumers: bounds/power consume only the point estimates and covariance, which remain valid; the warning is shape-descriptive because containers are publicly constructible. A rank-dropped interior event-time column (NaN-se row) fails closed via honest's consecutive-grid ValueError with a Stacked-specific remedy (no `balance_e` — Stacked's aggregate() has no such level). REFERENCE-SUPPORT GUARD: `fit()` raises when the omitted reference cell has no treated or no control rows stack-wide (gapped panel at every retained cohort's `a - 1 - anticipation`; IC1 checks only min/max periods) — without it the design rank-drops and QR pivoting silently re-normalizes against an arbitrary horizon while the surface certifies `delta_0 = 0` at the reference. Partial support (some cohorts' own ref period absent while the POOLED omitted category is populated) still fits: the coefficient-vector normalization is real; the per-cohort baseline raggedness is estimand composition, the documented ragged-window behavior. *Covariate-balancing edge cases:* - Infeasible cohort (treated covariate mean outside the clean-control hull → entropy balancing cannot match the moments): **fail-closed** `ValueError` naming the cohort and worst covariate — NOT silently dropped (dropping a cohort would shift the estimand to an overlap-trimmed ATT, Ustyuzhanin 2026 §3.1). @@ -4650,7 +4650,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; producer-built containers always record the field) warns fail-safe that the common reference cannot be verified. Container admission covers CS- and Stacked-sourced containers (row M-024): StackedDiD threads `base_period="universal"` and a singleton `reference_event_times` by construction (one omitted reference per fit), needs `kappa_pre >= 2` for estimated pre-periods, and HonestDiD additionally requires its pre-period covariance non-singular — see the StackedDiD section's container-admission Note for the df_survey and withheld-inference semantics. +- **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; producer-built containers always record the field) warns fail-safe that the common reference cannot be verified. Container admission covers CS- and Stacked-sourced containers (row M-024): StackedDiD threads `base_period="universal"` and a singleton `reference_event_times` by construction (one omitted reference per fit), needs `kappa_pre >= 2` for estimated pre-periods, and HonestDiD additionally requires the full retained event-study covariance (pre+post sub-block) non-singular — see the StackedDiD section's container-admission Note for the df_survey and withheld-inference semantics. **Reference implementation(s):** - R: `HonestDiD` package (Rambachan & Roth's official package) diff --git a/docs/v4-deprecations.yaml b/docs/v4-deprecations.yaml index 20991cbee..74061748a 100644 --- a/docs/v4-deprecations.yaml +++ b/docs/v4-deprecations.yaml @@ -996,7 +996,7 @@ rows: 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, 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. SECOND PRE-CUT AMENDMENT (2(b) PR-2, with M-024): admission widened to STACKED-SOURCED containers in compute_honest_did + compute_pretrends_power (source in {CallawaySantAnnaResults, StackedDiDResults}; kappa_pre >= 2 required for estimated pre-periods; honest additionally needs a non-singular pre-period covariance; withheld-inference rows admitted with a source-scoped warning; the seven producer-derived guard messages now derive the producer from surface.source). Source scoping continues for all other producers. This row stays planned: its transition remains the 4.0 sentinel retirement + full-producer schema enforcement." + 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. SECOND PRE-CUT AMENDMENT (2(b) PR-2, with M-024): admission widened to STACKED-SOURCED containers in compute_honest_did + compute_pretrends_power (source in {CallawaySantAnnaResults, StackedDiDResults}; kappa_pre >= 2 required for estimated pre-periods; honest additionally needs a non-singular full retained event-study covariance (pre+post sub-block); withheld-inference rows admitted with a source-scoped warning; the seven producer-derived guard messages now derive the producer from surface.source). Source scoping continues for all other producers. 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/tests/test_stacked_did.py b/tests/test_stacked_did.py index e8cfbdf17..9f11361be 100644 --- a/tests/test_stacked_did.py +++ b/tests/test_stacked_did.py @@ -2114,3 +2114,71 @@ def test_plot_event_study_renders_on_plain_fit(self, staggered_data): ax = plot_event_study(results, show=False) assert ax is not None plt.close("all") + + +class TestReferenceSupportGuard: + """M-024 follow-up (CI review): the omitted reference cell must be + populated - a gapped panel where every retained cohort's calendar + period a-1-anticipation is absent would otherwise rank-drop and + silently re-normalize while the surface certifies delta_0 = 0 at the + reference (and the container's reference_event_times provenance + would then mislead HonestDiD/PreTrendsPower).""" + + @staticmethod + def _gapped_panel(periods, cohorts, n_units=90, seed=0): + rng = np.random.default_rng(seed) + rows = [] + for u in range(n_units): + g = cohorts[u % len(cohorts)] + for t in periods: + d = 1 if (g and t >= g) else 0 + rows.append( + { + "unit": u, + "time": t, + "y": u * 0.01 + 0.2 * t + 1.5 * d + rng.normal(0, 0.3), + "first_treat": g, + } + ) + return pd.DataFrame(rows) + + def test_stack_wide_missing_reference_fails_closed(self): + # Calendar period 4 absent; sole cohort a=5 has ref a-1=4 -> the + # reference cell is empty stack-wide. IC1 passes (min/max only), + # so without the guard this rank-drops and re-normalizes against + # an arbitrary horizon while synthesizing a fake e=-1 row. + panel = self._gapped_panel([1, 2, 3, 5, 6, 7], [0, 5]) + est = StackedDiD(kappa_pre=2, kappa_post=2) + with pytest.raises(ValueError, match="would be fabricated"): + est.fit( + panel, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + ) + + def test_partial_reference_support_still_fits(self): + # Cohort 5's own ref period (4) is absent but cohort 7's (6) is + # present, so the POOLED omitted category is populated: the + # coefficient-vector normalization at e=-1 is real (verified by + # execution during the CI-review triage - no rank drop, sane + # surface). Ragged per-cohort baseline support affects estimand + # composition (documented ragged-window behavior), not the + # normalization the container certifies. + panel = self._gapped_panel([1, 2, 3, 5, 6, 7, 8, 9], [0, 5, 7]) + est = StackedDiD(kappa_pre=2, kappa_post=2) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + res = est.fit( + panel, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + ) + assert not [w for w in caught if "Rank-deficient" in str(w.message)] + assert res.event_study_effects is not None + surf = res.aggregate("event_study") + assert surf.reference_event_times == (-1,) + assert np.isfinite(res.overall_att) From 1b847304da11c767d9bbcb132ccfc3f2d03211db Mon Sep 17 00:00:00 2001 From: igerber Date: Mon, 3 Aug 2026 16:06:45 -0400 Subject: [PATCH 3/4] fix(v4): reference-support guard tests positive composed WLS weight, not raw rows (CI review R2) Verified by execution: Equation 3 is Q-weighted WLS, so a reference cell can exist in raw rows yet carry zero effective mass - (a) under aggregate weighting, a sub-experiment's control rows get Q=0 at any event time where that sub-experiment has no treated rows, so a ragged panel whose only control-reference rows sit in such sub-experiments passed the raw-row guard while the effective control baseline was empty (reproduced: rank drop, distorted surface, container still certifying reference_event_times=(-1,)); (b) survey pweights can zero a cell the same way. The guard now runs after weight composition and requires positive total composed weight in BOTH reference cells, with the causes enumerated in the error. Two new tests pin the zero-Q ragged case and the survey-zeroed case; the stack-wide-missing and partial-support scenarios keep their prior behavior (re-verified). REGISTRY/CHANGELOG descriptions updated to the weighted formulation. --- CHANGELOG.md | 11 +++-- diff_diff/stacked_did.py | 78 ++++++++++++++++++++--------------- docs/methodology/REGISTRY.md | 2 +- tests/test_stacked_did.py | 79 ++++++++++++++++++++++++++++++++++++ 4 files changed, 132 insertions(+), 38 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ae9e9ed78..5d29500b5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -62,10 +62,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 its FLCI CI endpoints are NaN). The seven producer-derived guard messages in both container branches now derive the producer from `surface.source` with producer-conditional remedies. A new - reference-support guard fails `fit()` closed when the omitted - reference cell has no treated or control rows stack-wide (gapped - panel at every cohort's `a - 1 - anticipation`) - previously this - rank-dropped and silently re-normalized while the surface + reference-support guard fails `fit()` closed when either omitted- + reference cell (treated / control) carries no positive COMPOSED + WLS weight - raw row presence is not support in a Q-weighted + regression: gapped panels empty the cell outright, aggregate + Q-weights zero a sub-experiment's controls wherever it has no + treated rows, and survey pweights can zero a cell. Previously + these rank-dropped and silently re-normalized while the surface certified a `delta_0 = 0` reference that never existed. - Cross-producer container hardening: `EventStudyResults`'s `__post_init__` now COPIES `vcov`/`vcov_index` like every other diff --git a/diff_diff/stacked_did.py b/diff_diff/stacked_did.py index 3c7a22b57..4d627478d 100644 --- a/diff_diff/stacked_did.py +++ b/diff_diff/stacked_did.py @@ -648,39 +648,6 @@ def fit( et_vals = stacked_df["_event_time"].values d_vals = stacked_df["_D_sa"].values - # ---- Reference-support guard (row M-024 follow-up) ---- - # The omitted reference event time e = -1 - anticipation is the - # regression's baseline category: every reported delta_h (and the - # synthesized reference row, and the container's - # reference_event_times / base_period="universal" provenance) is - # defined relative to its cells. On a gapped panel where NO - # retained sub-experiment observes that calendar period, the ref - # cell is empty, the design goes rank-deficient, and QR pivoting - # silently re-normalizes against an arbitrary surviving column - - # the surface would then certify a delta_0 = 0 normalization that - # never happened (and HonestDiD / PreTrendsPower would trust it). - # Both cells are required: without TREATED ref rows, D_sa is - # collinear with the delta block (same silent re-normalization); - # without CONTROL ref rows, the baseline time profile is. Fail - # closed - a fabricated reference must never reach consumers. - _ref_mask = et_vals == ref_period - if not np.any(_ref_mask & (d_vals == 1)) or not np.any(_ref_mask & (d_vals == 0)): - _missing_cell = "treated" if not np.any(_ref_mask & (d_vals == 1)) else "control" - raise ValueError( - f"The omitted reference event time e={ref_period} has no " - f"{_missing_cell} observations in the stacked data, so the " - "event-study normalization (delta_0 = 0 at the reference) " - "would be fabricated: the regression's baseline cell is " - "empty and rank handling would silently re-normalize " - "against an arbitrary horizon. This happens on gapped " - "panels where every retained cohort's calendar period " - "a - 1 - anticipation is absent (the IC1 window check " - "only inspects the panel's min/max periods, not interior " - "gaps). Fill the gap, adjust anticipation/kappa_pre so " - "the reference falls on an observed period, or drop the " - "affected cohorts." - ) - # Build design matrix X = np.zeros((n, 2 + 2 * n_event_dummies)) X[:, 0] = 1.0 # intercept @@ -714,6 +681,51 @@ def fit( else: composed_weights = Q_weights + # ---- Reference-support guard (row M-024 follow-up, CI R1+R2) ---- + # The omitted reference event time e = -1 - anticipation is the + # regression's baseline category: every reported delta_h (and the + # synthesized reference row, and the container's + # reference_event_times / base_period="universal" provenance) is + # defined relative to its cells. Equation 3 is a Q-WEIGHTED WLS, + # so the check runs on the COMPOSED weights, not raw row + # presence: a cell whose rows all carry zero effective weight + # (aggregate Q-weights zero a sub-experiment's controls at any + # event time where that sub-experiment has no treated rows; + # survey pweights can zero a cell outright) is as empty as a + # missing one. An effectively-empty cell makes the design + # rank-deficient and QR pivoting silently re-normalizes against + # an arbitrary surviving column - the surface would then certify + # a delta_0 = 0 normalization that never happened (and + # HonestDiD / PreTrendsPower would trust it). Both cells are + # load-bearing: without effective TREATED ref weight, D_sa is + # collinear with the delta block; without effective CONTROL ref + # weight, the baseline time profile is. Fail closed - a + # fabricated reference must never reach consumers. + _ref_mask = et_vals == ref_period + _w_arr = np.asarray(composed_weights, dtype=float) + _treated_ref_w = float(np.sum(_w_arr[_ref_mask & (d_vals == 1)])) + _control_ref_w = float(np.sum(_w_arr[_ref_mask & (d_vals == 0)])) + if not (_treated_ref_w > 0) or not (_control_ref_w > 0): + _missing_cell = "treated" if not (_treated_ref_w > 0) else "control" + raise ValueError( + f"The omitted reference event time e={ref_period} has no " + f"{_missing_cell} observations with positive weight in the " + "stacked WLS, so the event-study normalization " + "(delta_0 = 0 at the reference) would be fabricated: the " + "regression's effective baseline cell is empty and rank " + "handling would silently re-normalize against an " + "arbitrary horizon. Causes: a gapped panel where every " + "retained cohort's calendar period a - 1 - anticipation " + "is absent (the IC1 window check only inspects the " + "panel's min/max periods, not interior gaps); ragged " + "panels where the only reference-period control rows sit " + "in sub-experiments without treated reference rows " + "(aggregate Q-weights zero those controls); or survey " + "weights that zero out a reference cell. Fill the gap, " + "adjust anticipation/kappa_pre so the reference falls on " + "an observed period, or drop the affected cohorts." + ) + Y = stacked_df[outcome].values # Cluster IDs diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 24741de62..cea440ba6 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2082,7 +2082,7 @@ The pooled estimator is `DID^{CBWSDID}_e = Σ_a (N^D_a/N^D_Ω)(Δ̄^D_{a,e} − - **Note:** Inference is conditional-on-the-estimated-weights cluster-robust (the existing `hc1`/`hc2_bm` path with `W_{sa}` as the WLS weights) — the paper's default. The paper's weight-re-estimating bootstrap is NOT implemented in v1 (deliberate scope; entropy balancing is smooth so the Abadie–Imbens (2008) nonsmooth-matching bootstrap caveat does not apply). `cluster` is orthogonal to `b_{sa}` (weights conditioned-on); default `unit` matches the paper. - **Note:** v1 scope — only `balance="entropy"` with `weighting="aggregate"`. `balance` + `population`/`sample_share` and `balance` + `survey_design=` raise `NotImplementedError`; matching-based balancing and the repeated `0→1/1→0` episode extension are out of scope. - **Note (post-fit aggregate() - row M-024):** `fit(aggregate=)` is deprecated (3.9; removed 4.0) and the event-study surface is ALWAYS materialized at fit — the pooled regression always included the event-time interactions, so only the extraction and, under `hc2_bm`, the per-event Bell-McCaffrey contrast DOFs were fit-time-gated; both gates are gone and the param is behaviorally inert (its `group`/`all`/unknown ValueErrors persist after the shim warning). `aggregate('event_study')` returns the unified `EventStudyResults` container (a pure view; no recomputation, so every stored inference mode — survey TSL, replicate refit, hc2_bm fail-close — relays faithfully); `aggregate('simple')` relays `overall_att/se/t/p/CI` bit-exactly with `target="att"` (the scalar is the equally-weighted post-period average, NOT the per-event-time trimmed aggregate ATT `θ^e_κ` — `describe_target_parameter` remains the estimand's prose source of truth), `n = n_treated_units` with `n_kind="units"` (TREATED-only scope: under dynamic clean-control definitions a later-treated unit is treated in its own sub-experiment and a clean control in earlier ones, so the treated and control unit sets OVERLAP and no disjoint total exists — cross-container `n` comparisons vs CallawaySantAnna's treated+control convention are out of contract), and `df = inference_df` (the stored overall inference's exact provenance). Consequences of the always-computed surface: plain-fit `summary()` renders the event-study table, `to_dataframe('event_study')` is total on 3.9+ fits, `plot_event_study` renders plain fits, DiagnosticReport's `heterogeneity` (always) and `parallel_trends` (`kappa_pre >= 2`) checks run on plain fits, and the `hc2_bm` overall-ATT contrast now shares the m>1 batch of `_compute_cr2_bm_contrast_dof`, whose degeneracy guard is batch-relative — the noise-floor carve-out on the `Wald_test(test="HTZ")$df_denom` parity claim in the variance-families bullets above widens quantitatively (unchanged behavior when tripped: fail-closed NaN inference plus linalg's loud noise-floor UserWarning; the absolute floor already made this reachable at m=1). -- **Note (container admission - rows M-024/M-093):** `compute_honest_did` and `compute_pretrends_power` admit Stacked-sourced `aggregate('event_study')` containers. Provenance: `base_period` is the class-level constant `"universal"` (every sub-experiment normalizes against the single omitted reference `e = -1 - anticipation`) and `reference_event_times` is the singleton `(-1 - anticipation,)`. Requirements: `kappa_pre >= 2` (the default `kappa_pre=1` grid minus the reference has ZERO estimated pre-periods — both consumers reject), and for HonestDiD a non-singular FULL retained event-study covariance — the validated sub-block spans the retained pre AND post rows (`allow_singular=False` per Rambachan-Roth's eigenvalues-bounded-away-from-zero assumption on the full Σ — keep `kappa_pre` small relative to the cluster count; PreTrendsPower keeps its documented singular support). The container's scalar `df_survey` is None on analytical fits — the generic resolver's fallback reads `df_inference`, a name StackedDiD does not carry (`inference_df` is the overall POST-AVERAGE contrast df, not a design-level cluster df, and misusing it for FLCI critical values would be worse than the asymptotic default) — so honest FLCI critical values are NORMAL-THEORY on analytical Stacked fits, diverging from bare-`cluster=` CS fits (which thread `G-1`); survey fits thread the survey df and replicate-undefined designs the fail-closed `0.0` sentinel (NaN FLCI CIs; the identified-set `lb`/`ub` bounds are df-independent and stay finite). Withheld-inference rows (finite `se`, non-finite `p_value` — the hc2_bm BM-DOF fail-close and replicate-undefined designs both emit this shape) are ADMITTED with a source-scoped warning in BOTH consumers: bounds/power consume only the point estimates and covariance, which remain valid; the warning is shape-descriptive because containers are publicly constructible. A rank-dropped interior event-time column (NaN-se row) fails closed via honest's consecutive-grid ValueError with a Stacked-specific remedy (no `balance_e` — Stacked's aggregate() has no such level). REFERENCE-SUPPORT GUARD: `fit()` raises when the omitted reference cell has no treated or no control rows stack-wide (gapped panel at every retained cohort's `a - 1 - anticipation`; IC1 checks only min/max periods) — without it the design rank-drops and QR pivoting silently re-normalizes against an arbitrary horizon while the surface certifies `delta_0 = 0` at the reference. Partial support (some cohorts' own ref period absent while the POOLED omitted category is populated) still fits: the coefficient-vector normalization is real; the per-cohort baseline raggedness is estimand composition, the documented ragged-window behavior. +- **Note (container admission - rows M-024/M-093):** `compute_honest_did` and `compute_pretrends_power` admit Stacked-sourced `aggregate('event_study')` containers. Provenance: `base_period` is the class-level constant `"universal"` (every sub-experiment normalizes against the single omitted reference `e = -1 - anticipation`) and `reference_event_times` is the singleton `(-1 - anticipation,)`. Requirements: `kappa_pre >= 2` (the default `kappa_pre=1` grid minus the reference has ZERO estimated pre-periods — both consumers reject), and for HonestDiD a non-singular FULL retained event-study covariance — the validated sub-block spans the retained pre AND post rows (`allow_singular=False` per Rambachan-Roth's eigenvalues-bounded-away-from-zero assumption on the full Σ — keep `kappa_pre` small relative to the cluster count; PreTrendsPower keeps its documented singular support). The container's scalar `df_survey` is None on analytical fits — the generic resolver's fallback reads `df_inference`, a name StackedDiD does not carry (`inference_df` is the overall POST-AVERAGE contrast df, not a design-level cluster df, and misusing it for FLCI critical values would be worse than the asymptotic default) — so honest FLCI critical values are NORMAL-THEORY on analytical Stacked fits, diverging from bare-`cluster=` CS fits (which thread `G-1`); survey fits thread the survey df and replicate-undefined designs the fail-closed `0.0` sentinel (NaN FLCI CIs; the identified-set `lb`/`ub` bounds are df-independent and stay finite). Withheld-inference rows (finite `se`, non-finite `p_value` — the hc2_bm BM-DOF fail-close and replicate-undefined designs both emit this shape) are ADMITTED with a source-scoped warning in BOTH consumers: bounds/power consume only the point estimates and covariance, which remain valid; the warning is shape-descriptive because containers are publicly constructible. A rank-dropped interior event-time column (NaN-se row) fails closed via honest's consecutive-grid ValueError with a Stacked-specific remedy (no `balance_e` — Stacked's aggregate() has no such level). REFERENCE-SUPPORT GUARD: `fit()` raises when either omitted-reference cell (treated / control) carries no POSITIVE COMPOSED WLS WEIGHT — Equation 3 is Q-weighted, so raw row presence is not support: a gapped panel can empty the cell outright (IC1 checks only min/max periods), aggregate Q-weights zero a sub-experiment's controls at any event time where that sub-experiment has no treated rows, and survey pweights can zero a cell — without the guard the design rank-drops and QR pivoting silently re-normalizes against an arbitrary horizon while the surface certifies `delta_0 = 0` at the reference. Partial support (some cohorts' own ref period absent while the POOLED omitted category is populated) still fits: the coefficient-vector normalization is real; the per-cohort baseline raggedness is estimand composition, the documented ragged-window behavior. *Covariate-balancing edge cases:* - Infeasible cohort (treated covariate mean outside the clean-control hull → entropy balancing cannot match the moments): **fail-closed** `ValueError` naming the cohort and worst covariate — NOT silently dropped (dropping a cohort would shift the estimand to an overlap-trimmed ATT, Ustyuzhanin 2026 §3.1). diff --git a/tests/test_stacked_did.py b/tests/test_stacked_did.py index 9f11361be..54efde1df 100644 --- a/tests/test_stacked_did.py +++ b/tests/test_stacked_did.py @@ -2182,3 +2182,82 @@ def test_partial_reference_support_still_fits(self): surf = res.aggregate("event_study") assert surf.reference_event_times == (-1,) assert np.isfinite(res.overall_att) + + +class TestReferenceSupportGuardEffectiveWeights: + """CI review R2: Equation 3 is Q-WEIGHTED WLS, so reference support + means positive COMPOSED weight in both baseline cells - raw row + presence is not enough.""" + + def test_zero_q_weight_control_reference_cell_fails_closed(self): + # Raw rows exist in both reference cells, but the only control + # rows at e=-1 sit in sub-experiments whose treated rows are + # absent at that event time - aggregate Q-weights zero them, the + # effective control baseline is empty, and (pre-guard, verified + # by execution) the design rank-dropped and re-normalized while + # the container still certified reference_event_times=(-1,). + rng = np.random.default_rng(3) + rows = [] + for u in range(120): + kind = ["treatA", "ctrlX", "treatB", "ctrlY"][u % 4] + g = {"treatA": 5, "ctrlX": 0, "treatB": 7, "ctrlY": 0}[kind] + for t in range(1, 11): + if kind == "ctrlX" and t == 4: + continue # A's controls gapped at A's ref (cal 4) + if kind == "treatB" and t == 6: + continue # B's treated gapped at B's ref (cal 6) + if kind == "ctrlY" and t == 4: + continue # Y gapped at cal 4 too + d = 1 if (g and t >= g) else 0 + rows.append( + { + "unit": u, + "time": t, + "y": u * 0.01 + 0.2 * t + 1.5 * d + rng.normal(0, 0.3), + "first_treat": g, + } + ) + est = StackedDiD(kappa_pre=2, kappa_post=2) + with pytest.raises(ValueError, match="positive weight.*fabricated"): + est.fit( + pd.DataFrame(rows), + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + ) + + def test_survey_zeroed_reference_cell_fails_closed(self): + from diff_diff.survey import SurveyDesign + + rng = np.random.default_rng(3) + rows = [] + for u in range(90): + g = [0, 5, 7][u % 3] + for t in range(1, 11): + d = 1 if (g and t >= g) else 0 + # controls carry zero pweight at both cohorts' reference + # calendar periods -> the control baseline cell has rows + # but no effective WLS mass + w = 0.0 if (g == 0 and t in (4, 6)) else 1.0 + rows.append( + { + "unit": u, + "time": t, + "y": u * 0.01 + 0.2 * t + 1.5 * d + rng.normal(0, 0.3), + "first_treat": g, + "w": w, + } + ) + est = StackedDiD(kappa_pre=2, kappa_post=2) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + with pytest.raises(ValueError, match="positive weight.*fabricated"): + est.fit( + pd.DataFrame(rows), + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + survey_design=SurveyDesign(weights="w"), + ) From 1c28e7f729b78b8476b19ed12c9f83b3a4a52982 Mon Sep 17 00:00:00 2001 From: igerber Date: Mon, 3 Aug 2026 17:47:52 -0400 Subject: [PATCH 4/4] test(v4): fix CI failures - matplotlib importorskip + deep-tail p-value tolerance Two test-only fixes for the PR CI matrix: - test_plot_event_study_renders_on_plain_fit failed on every leg with ModuleNotFoundError: matplotlib is an optional dependency the CI test env does not install. Guard with pytest.importorskip per the existing convention (test_visualization.py, test_bacon.py, test_event_study_consumers.py). - test_bm_dof_batch_parity_and_overall_reconstruction failed on the Linux x86 legs only: batched-vs-m=1 BLAS kernels differ at ~1 ULP in the Satterthwaite dof (OpenBLAS reassociation), and the overall contrast's deep-tail p-value (~1e-54 at t~29) amplifies that relative noise by ~t^2 (d ln p ~ -t dt), landing at 1.5e-14 - just over the uniform rtol=1e-14. Split the assertion: t/CI stay at rtol=1e-14, the p-value alone loosens to rtol=1e-12 with the amplification documented inline. macOS Accelerate and Windows were bit-stable; no library code changes. --- tests/test_aggregate_contract.py | 9 +++++++-- tests/test_stacked_did.py | 2 +- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/tests/test_aggregate_contract.py b/tests/test_aggregate_contract.py index 40b5276e4..d30720ee8 100644 --- a/tests/test_aggregate_contract.py +++ b/tests/test_aggregate_contract.py @@ -1423,15 +1423,20 @@ def spy(X, cluster_ids, bread_matrix, contrasts, weights=None): res.overall_att, res.overall_se, alpha=res.alpha, df=overall_dof_m1 ) np.testing.assert_allclose( - [t, p, ci[0], ci[1]], + [t, ci[0], ci[1]], [ res.overall_t_stat, - res.overall_p_value, res.overall_conf_int[0], res.overall_conf_int[1], ], rtol=1e-14, ) + # The p-value gets a looser rtol: batched-vs-m=1 BLAS kernels differ + # at ~1 ULP in the dof (platform-dependent - observed on Linux + # x86 OpenBLAS), and the deep tail amplifies that relative noise by + # ~t^2 (d ln p ~ -t dt; t ~ 29 here => ~1e3x), so p carries ~1e-13 + # relative noise while t/CI stay at 1e-14. + np.testing.assert_allclose(p, res.overall_p_value, rtol=1e-12) def test_legacy_pickle_absent_surface_hint(self, stacked_fitted): import dataclasses diff --git a/tests/test_stacked_did.py b/tests/test_stacked_did.py index 54efde1df..25ac9eb8f 100644 --- a/tests/test_stacked_did.py +++ b/tests/test_stacked_did.py @@ -2096,7 +2096,7 @@ def test_plot_event_study_renders_on_plain_fit(self, staggered_data): """M-024: plot_event_study flips from TypeError to rendering on a plain fit - the duck-typed effects branch now matches because the surface is always populated.""" - import matplotlib + matplotlib = pytest.importorskip("matplotlib") matplotlib.use("Agg") import matplotlib.pyplot as plt