From 716b6b4a783c5da37a4bd80c85675a5ad5f93d1f Mon Sep 17 00:00:00 2001 From: igerber Date: Fri, 31 Jul 2026 16:02:27 -0400 Subject: [PATCH 1/2] feat(3.9): tail-df consolidation - three-value df_convention knob + t(residual) defect fix (M-127) Converge the tail degrees of freedom the t-distribution uses across the analytical clustered surfaces (PR C of the variance-consolidation program, after #736 inventory/D3 and #737 K_reference/M-126). Defect fix (z -> t(residual), point estimates/SEs/t-stats unchanged): - Wooldridge default-hc1 OLS arms: one knob-resolved df_fallback feeds every non-BM arm (survey df first; classical/hc2 reproduce the historical df_one_way values bit-for-bit under the default). - SunAbraham aggregates: cells and aggregates now share ONE df source per fit (the saturated regression's df_) - the D4 two-conventions-in- one-fit defect is gone. - StackedDiD pooled inference: positive-weight n_eff - k_kept. - ImputationDiD pretrends per-lead inference: n - k_kept - absorbed rank on the untreated subsample (the joint pretrend Wald F keeps its cluster-robust F(q, G-1) denominator, documented as a pairing). Three-value df_convention knob {"residual","cluster","normal"} on DiD/TWFE/MPD/LinearRegression (widened from two values; "normal" is ungated at the fallback level - clustered, unclustered, and conley) and NEW on SunAbraham, WooldridgeDiD (OLS arms), StackedDiD, ImputationDiD (pretrends), and LPDiD. LPDiD defaults to "cluster" (the existing Stata t(G-1) reference - zero movement, degenerate lanes keep literal df=None). "normal" reproduces the pre-3.9 numbers on the fixed surfaces. Survey/replicate df and hc2_bm BM DOF keep precedence under every value. set_params is atomic under unknown keys on all five standalone estimators. Shared resolver: utils.resolve_tail_df + validate_df_convention. Provenance: df_convention on five results containers; inference_df on SunAbraham/StackedDiD (None under bootstrap overrides); Wooldridge _df_one_way renamed _df_analytic_fallback with a __setstate__ pickle migration (stays live on bootstrap fits so post-fit analytical aggregate() reproduces fit-time inference). Ledger/docs: M-127 behavior row + M-128..M-131 default-flip rows (phase 5, 'residual' -> 'cluster'), matrix snapshot 105 -> 110, v4-design/DEFERRED/TODO/CHANGELOG sync; REGISTRY tail-df notes for TWFE/SA/StackedDiD/Wooldridge/ImputationDiD/LPDiD; regenerated variance-conventions inventory (defect rows converged); llms-full constructor/results updates. Tests: resolver unit branch table; per-estimator df-convention suites (cluster == hand-reconstructed t(G-1), normal == z, default != old z, survey-df precedence under every knob value, bit-identity locks); audit-matrix re-pins with measured df multisets; fixest golden p/CI now gated on both lanes (iid arms under the residual default, cluster arms under df_convention="cluster"); event-study-surface provenance flips; Wooldridge pickle-migration and bootstrap-aggregate parity tests. --- CHANGELOG.md | 40 ++++- DEFERRED.md | 2 +- TODO.md | 3 +- diff_diff/estimators.py | 44 ++--- diff_diff/guides/llms-full.txt | 15 +- diff_diff/imputation.py | 88 +++++++++- diff_diff/imputation_results.py | 8 + diff_diff/linalg.py | 37 +++-- diff_diff/lpdid.py | 80 ++++++++- diff_diff/lpdid_results.py | 7 + diff_diff/stacked_did.py | 82 +++++++-- diff_diff/stacked_did_results.py | 24 ++- diff_diff/sun_abraham.py | 201 ++++++++++++++++++----- diff_diff/utils.py | 91 ++++++++++ diff_diff/wooldridge.py | 159 ++++++++++++------ diff_diff/wooldridge_results.py | 66 +++++--- docs/methodology/REGISTRY.md | 141 ++++++++++++---- docs/methodology/variance-conventions.md | 71 +++++--- docs/v4-deprecations.yaml | 67 +++++++- docs/v4-design.md | 28 ++-- tests/test_estimators_vcov_type.py | 114 ++++++++++++- tests/test_event_study_surface.py | 17 +- tests/test_fixest_did_twfe_parity.py | 52 ++++++ tests/test_imputation.py | 139 ++++++++++++++++ tests/test_lpdid.py | 129 +++++++++++++++ tests/test_stacked_did.py | 131 +++++++++++++++ tests/test_sun_abraham.py | 149 ++++++++++++++++- tests/test_utils.py | 74 +++++++++ tests/test_v4_matrix.py | 17 +- tests/test_variance_conventions.py | 73 ++++++-- tests/test_wooldridge.py | 201 +++++++++++++++++++++++ 31 files changed, 2062 insertions(+), 288 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7929b2746..60c6fa554 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,25 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Fixed +- **Normal-theory tail inference on clustered analytical paths — the D4 + + tail-df defect fix of the 3.9 variance-consolidation program** ([M-127]). + Four surfaces previously passed `df=None` (silent normal-theory z) to + their analytical p-value/CI computation on clustered fits: `WooldridgeDiD` + default-`hc1` (within and `hc1`+`cohort_trends` arms), `StackedDiD`, + `SunAbraham` aggregates — the D4 defect: one fit reported residual-t + per-cohort cells but z aggregates — and the `ImputationDiD` pretrends lead + regression. At |t| = 2, z understates the t(G−1) p-value by 24.2% at + G=20 and 13.3% at G=40. All four now use **t(residual df)** under the new + `"residual"` default (aggregates share the cells' df source on + SunAbraham; classical/hc2 arms and every SE/point estimate are + bit-identical; the movement is z → t at the fitted residual df — the G−1 + convergence with fixest/Stata is the v4 default flip, M-004..M-006 + + M-128..M-131). `df_convention="normal"` reproduces the pre-3.9 numbers on + the affected lanes. Wild-cluster-bootstrap p/CIs (test-inversion) and all + influence-function estimators are unaffected. `WooldridgeDiD` post-fit + `aggregate()` reproduces fit-time inference on every OLS arm (the + internal `_df_one_way` generalized to `_df_analytic_fallback`, with a + `__setstate__` pickle migration). - **Clustered CR1 SEs converge on `K_reference` — the D1+D2 fix of the 3.9 variance-consolidation program** ([M-126]). The clustered `hc1` CR1 finite-sample factor `(G/(G−1))·((n−1)/(n−k))` now uses @@ -40,7 +59,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 `ssc(K.fixef="full", K.exact=TRUE)` at ~1e-16, with the one-df deviation of the nested composition from both references (neither implements nested-drop + exact-remainder) pinned exactly in two new goldens. Tail df - conventions are deliberately unchanged (PR C of the program). The + conventions were deliberately unchanged by this fix and converged + separately (the [M-127] entry above). The `docs/benchmarks.rst` TWFE "SE Rel Diff 0.1%" cell is generated by the R-dependent benchmark refresh and still shows the pre-fix value; the refresh re-run is tracked in TODO.md (expected movement 0.1% → 0.0%). @@ -65,8 +85,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 convention passed to inference, with each cell classified as a documented defect or a declared legitimate exception with its reason. The clustered CR1 `k` cells now pin the converged `K_reference` accounting (see the - [M-126] entry above); the remaining defect rows are the tail-df items - scheduled for PR C of the 3.9 variance-consolidation program. + [M-126] entry above), and the tail-df cells pin the converged + `df_convention` resolution (the [M-127] entry above) — the matrix carries + no remaining defect rows. - **`WooldridgeDiD` silently dropped genuine post-treatment effects (issue #724).** With `control_group="never_treated"`, the ETWFE design emitted every cohort×time indicator *including* the reference period — which is @@ -194,6 +215,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 rejecting the design. ### Added +- **`df_convention` becomes a library-wide three-value knob** ([M-127]). + `{"residual", "cluster", "normal"}` is now accepted everywhere the knob + exists, and the knob extends to five standalone estimators: `SunAbraham`, + `WooldridgeDiD` (OLS arms), `StackedDiD`, `ImputationDiD` (pretrends lead + regression), and `LPDiD` (default `"cluster"` — its Stata-anchored t(G−1) + reference, unchanged). `"cluster"` = the Stata/fixest `G − 1` reference + (the v4 default); `"normal"` = deliberate normal-theory z at the fallback + level on every fit (also the exact-reproducibility path for pre-3.9 + numbers on the surfaces fixed by the [M-127] entry under Fixed). Survey/ + replicate df and Bell-McCaffrey DOF keep precedence under every value; + point estimates, SEs, and t-statistics never move. Results objects echo + `df_convention`, and `SunAbraham`/`StackedDiD` gain `inference_df` (the + overall-ATT df actually used; None under bootstrap overrides). - **Stata parity arm for ETWFE and Callaway-Sant'Anna ATT(g,t) (`jwdid` / `csdid`).** `benchmarks/stata/generate_etwfe_cs_golden.do` anchors both staggered estimators against their canonical Stata implementations on the genuine diff --git a/DEFERRED.md b/DEFERRED.md index 066f7cd12..5ca1fabf7 100644 --- a/DEFERRED.md +++ b/DEFERRED.md @@ -113,7 +113,7 @@ here carry the remaining implementation work and cross-link their `M-xxx` ids. | Issue | Location | PR | Priority | |-------|----------|----|----------| -| Clustered-CR1 inference df **default flip to `"cluster"` (G−1) at v4** — the opt-in `df_convention=` knob landed 2026-07 (DiD/TWFE/MPD + LinearRegression; REGISTRY §TwoWayFixedEffects deviation note); the remaining work is the major-version default change (moves every clustered p-value/CI) + migration note + flipping `TestDfConvention`/`test_moderate_t_pins_residual_df_convention` expectations. Also evaluate extending the knob to standalone estimators with CR1-t inference at that time. Lifecycle tracked in docs/v4-deprecations.yaml (M-004..M-006). | `diff_diff/linalg.py::LinearRegression`, `diff_diff/estimators.py`, `diff_diff/twfe.py` | — | Medium | +| Clustered inference df **default flip to `"cluster"` (G−1) at v4** — the `df_convention=` knob is now the library-wide THREE-VALUE surface (PR C / M-127, 3.9: DiD/TWFE/MPD/LinearRegression + SunAbraham/WooldridgeDiD-OLS/StackedDiD/ImputationDiD-pretrends at default `"residual"`, LPDiD already at `"cluster"`); the remaining work is the major-version default change (moves every clustered p-value/CI) + migration note + flipping `TestDfConvention`/`test_moderate_t_pins_residual_df_convention`/the per-estimator knob suites' expectations. Lifecycle tracked in docs/v4-deprecations.yaml (M-004..M-006 + M-128..M-131). | `diff_diff/linalg.py::LinearRegression`, `diff_diff/estimators.py`, `diff_diff/twfe.py`, `diff_diff/sun_abraham.py`, `diff_diff/wooldridge.py`, `diff_diff/stacked_did.py`, `diff_diff/imputation.py` | — | Medium | ## Decision record — won't-fix / waived diff --git a/TODO.md b/TODO.md index 57c9754eb..3d31d4731 100644 --- a/TODO.md +++ b/TODO.md @@ -34,7 +34,7 @@ Related tracking surfaces: | `WooldridgeDiD` REFUSES a panel whose units split into disconnected support groups within a cohort, rather than estimating what IS identified. The connectivity guard (REGISTRY *within-cohort support connectivity*) correctly detects that a closed component's cells are collinear with the unit FE — previously QR dropped one silently and the overall ATT averaged an incomplete set (issue #724's failure mode via unit support). **Refusing is the safe answer, not the complete one.** The connected component containing the reference is still fully identified, so the estimable resolution is either (a) estimate the connected component and report the disconnected units as excluded, with the estimand restated (a sub-population of units, so it needs a REGISTRY definition and interacts with the survey-domain row above), or (b) per-component references, if a component with its own pre-period can carry its own normalization — needs a methodology decision, since components then are not comparable on one baseline. Gate with the split-support fixture in `TestWithinCohortSupportConnectivity`. | `diff_diff/wooldridge.py` | #724-codex-R7 | Heavy | Medium | | `WooldridgeDiD` fully resolves the `SurveyDesign` TWICE on every supported survey fit. The pre-exclusion validation pass (added so invalid metadata cannot hide in rows that cohort exclusion deletes) calls `survey_design.resolve(sample)`, and each fitter then calls `_resolve_survey_for_wooldridge` -> `_resolve_survey_for_fit` on the same frame, repeating weight normalization, strata/PSU/FPC validation and design-array construction. Any fit that REACHES the second resolve has an unchanged sample (survey + unidentified-cohort exclusion raises first), so the first result is reusable: capture the `_resolve_survey_for_fit` 4-tuple early and thread it into the three fitters as an optional `pre_resolved`. **Caveat that makes this non-trivial:** `sample = sample.reset_index(drop=True)` runs BETWEEN the two calls, so the reused object must be verified index-independent (resolution extracts positional numpy arrays, but `_inject_cluster_as_psu` and the metadata recompute need checking), and the early call must stop suppressing warnings or the user loses the weight-normalization notice. Gate with a survey fit asserting one normalization warning and byte-identical SEs. | `diff_diff/wooldridge.py` | #724-codex-R10 | Mid | Low | | Bootstrap re-aggregation for `CallawaySantAnnaResults.aggregate()` — a bootstrapped fit currently RAISES rather than substituting analytical inference for percentile-bootstrap statistics. The value-bound `BootstrapReplaySpec` (bit-identical replay, picklable, immune to post-fit `set_params`) is already in-tree and spike-verified; wiring it needs the per-`(g,t)` / per-event-time draw retention plus `assert_allclose` parity tests against the fit-time bootstrap numbers (NOT bit-identity — the fused GEMM's column count differs post-fit, ~1 ULP reassociation). | `diff_diff/aggregation.py`, `diff_diff/staggered_results.py` | #726 | Mid | Medium | -| Consolidate the inference-df precedence duplicated across `honest_did.py` (3 copies at ~L655/L836/L1004) onto the shared `resolve_inference_df()` helper added in `diff_diff/aggregation.py`. The copies are correct today; the risk is drift if the survey/replicate precedence changes in one place only. | `diff_diff/honest_did.py` | #726 | Quick | Low | +| Consolidate the inference-df precedence duplicated across `honest_did.py` (3 copies at ~L655/L836/L1004) onto the shared `resolve_inference_df()` helper added in `diff_diff/aggregation.py`. The copies are correct today; the risk is drift if the survey/replicate precedence changes in one place only. (Adjacent but distinct from PR C's `utils.resolve_tail_df`: that is the FIT-TIME `df_convention` fallback resolver, this is a post-fit results READER.) | `diff_diff/honest_did.py` | #726 | Quick | Low | | `ContinuousDiD` CGBS-2024 remaining extensions (earlier phases — `covariates=` reg/dr, `treatment_type="discrete"`, single-cohort `control_group="lowest_dose"` with estimand `ATT(d)−ATT(d_L)` — are already supported; see REGISTRY Note #7). Remaining (all deferred `NotImplementedError`, documented): `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate); `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF); multi-cohort **heterogeneous-support** discrete aggregation (support-aware: average each dose only over the cohorts that observe it); **multi-cohort `lowest_dose`** (within-cohort `d_L` reference + support-aware cross-cohort aggregation); and **`covariates=` × `lowest_dose`** (conditional-PT-relative-to-`d_L` estimand). Single-cohort / 2-period / shared-support multi-cohort are supported. | `continuous_did.py` | CGBS-2024 | Heavy | Low | | `WooldridgeDiD` does not apply the W2025 Sec 5.4 `D_{G_max} x X` covariate normalization, and three sibling covariate rank deficiencies are pre-existing. Measured with the period range pinned and only the never-treated units toggled: (1) time-invariant `exovar` is absorbed by the unit FE, 4 of 26 columns, IDENTICALLY with and without never-treated units; (2) `xgvar`'s cell x covariate block, 19 of 41, identical on both panels; (3) `xtvar` under `demean_covariates=False` does exhibit the `sum_g D_g x = x` dependency that the default demeaning removes; (4) the newly-reachable case -- time-VARYING data passed through `exovar`, which its own docstring reserves for time-invariant covariates -- where the paper's `dT_i` rule would give a deterministic `D_{G_max} x X` drop instead of QR's arbitrary pick (coefficients unaffected, `1.35e-14`; `rank_deficient_action="error"` raises). REGISTRY's narrowed Sec 5.4 note cross-references this row. **Trap for whoever takes it:** `xtvar` under the DEFAULT `demean_covariates=True` is FULL RANK -- the raw block carries demeaned values while `D_g x X` carries raw ones -- and forcing the drop there moves `overall_att` 1.11903 -> 1.46269. Pinned as-is by `TestComparisonSupportFiltering::test_cells_derived_groups_did_not_leak_into_the_design`. | `diff_diff/wooldridge.py` | #729-followup | Heavy | Medium | | `WooldridgeDiD.n_control_units` counts never-treated UNITS on `control_group="never_treated"` regardless of method, but on the nonlinear paths (`logit`/`poisson`) treated units' pre-treatment rows ARE the identifying comparison -- only the OLS path absorbs them into their own cells. So the reported count under-states the comparison pool exactly where the REGISTRY control-pool asymmetry note applies. Widen to `not_yet_treated or (never_treated and method != "ols")`, or document the count as never-treated-units-by-definition. Behavior is PRE-EXISTING; documented for now in the REGISTRY control-pool Note rather than changed, because widening moves a public results field and wants its own ledger row and test matrix. | `diff_diff/wooldridge.py` | #729-followup | Mid | Low | @@ -64,6 +64,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Section-8 naming-completeness guard test: sweep the public surface (module-level functions, class `__init__`/`fit` params, results dataclass fields) for contract-rename violations and fail unless each hit is either a `docs/v4-deprecations.yaml` row or an allowlist entry carrying a stated reason (seed the allowlist with the documented domain-vocabulary exemptions: the staggered family's ATT(g,t) `group`/`groups`, `TripleDifference.fit[group]`, wrapper params dying with `M-070..M-077`, and params inherited from `DifferenceInDifferences.__init__`). Today the rules are normative prose and the ledger holds only what an audit remembered — two successive manual sweeps each found surfaces the prior one missed (`M-094`/`M-095`, then the 19-row `M-097..M-115` function sweep). Should land BEFORE Phase 2c so the rename PR works from a mechanically-verified list. Two further checks belong in the same guard, each having been missed by hand across successive review rounds: (a) **phase-table agreement** — every row's `phase` must appear in the matching `docs/v4-design.md` section 9 table entry and vice versa (the section 9 checklist asserts this agreement but nothing verifies it, and three separate rounds caught a stale table); (b) **consumer coverage** (section 8 rule 11) — for each rename row, grep `diff_diff/` and `docs/methodology/` for the old name and require every hit to be either in the row's `code_refs` or allowlisted, since `getattr(obj, "old_name", default)` degrades silently rather than raising after removal. | `tests/`, `docs/v4-design.md` | gating-completeness amendment | Mid | Medium | | 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 | | `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/diff_diff/estimators.py b/diff_diff/estimators.py index 960d04251..fbaa5d130 100644 --- a/diff_diff/estimators.py +++ b/diff_diff/estimators.py @@ -41,6 +41,7 @@ validate_binary, validate_covariate_names, validate_design_term_names, + validate_df_convention, wild_bootstrap_se, ) @@ -136,19 +137,22 @@ class DifferenceInDifferences: path — absent ``cluster=``, pure Conley spatial HAC applies. ``survey_design=`` + Conley and ``inference='wild_bootstrap'`` + Conley both raise ``NotImplementedError``. - df_convention : str, default "residual" - Degrees-of-freedom convention for t-statistics, p-values, and CIs on - clustered analytical fits. ``"residual"`` (default) uses the fitted - residual df (``n − K_full``); ``"cluster"`` uses the Stata/fixest - cluster df ``G − 1``. Applies only at the fallback level of the df - resolution: survey df and per-coefficient Bell-McCaffrey DOF + df_convention : {"residual", "cluster", "normal"}, default "residual" + Degrees-of-freedom convention for analytical t-statistics, p-values, + and CIs. ``"residual"`` (default) uses the fitted residual df + (``n − K_full``); ``"cluster"`` uses the Stata/fixest cluster df + ``G − 1`` on clustered fits — it has no effect on unclustered fits + or on ``vcov_type="conley"`` (the combined Conley+cluster product + kernel has no documented ``G − 1`` df reference and keeps the + residual df); ``"normal"`` deliberately uses normal-theory z + inference at the fallback level on every fit, clustered or not. + Applies only at the fallback level of the df resolution under every + value: survey df and per-coefficient Bell-McCaffrey DOF (``vcov_type="hc2_bm"``) are more refined small-sample corrections and always take precedence. Point estimates, SEs, and t-statistics - are unaffected — only the reference t-distribution changes. Has no - effect on unclustered fits or on ``vcov_type="conley"`` (the combined - Conley+cluster product kernel has no documented ``G − 1`` df - reference and keeps the residual df). The default flips to ``"cluster"`` at - v4 (see the REGISTRY clustered-CR1 inference-df deviation note). + are unaffected — only the reference distribution changes. The + default flips to ``"cluster"`` at v4 (see the REGISTRY clustered-CR1 + inference-df deviation note). Attributes ---------- @@ -220,10 +224,7 @@ def __init__( # helper so __init__ and set_params use identical validation logic. from diff_diff.linalg import resolve_vcov_type - if df_convention not in ("residual", "cluster"): - raise ValueError( - f"df_convention must be 'residual' or 'cluster', got {df_convention!r}" - ) + validate_df_convention(df_convention) self.robust = robust self.cluster = cluster @@ -1243,10 +1244,7 @@ def set_params(self, **params) -> "DifferenceInDifferences": pending_robust = params.get("robust", self.robust) pending_vcov_type = params.get("vcov_type", self.vcov_type) pending_df_convention = params.get("df_convention", self.df_convention) - if pending_df_convention not in ("residual", "cluster"): - raise ValueError( - "df_convention must be 'residual' or 'cluster', " f"got {pending_df_convention!r}" - ) + validate_df_convention(pending_df_convention) # First pass: validate that every incoming key is a known attribute # so we don't partially apply a batch that ends in "Unknown parameter". @@ -2266,6 +2264,14 @@ def _refit_mp_absorb(w_r): df = 0 else: df = _g_eff_mp - 1 + elif self.df_convention == "normal": + # Deliberate normal-theory z inference at the fallback level, on + # every fit (clustered, unclustered, and conley alike). The + # survey/replicate overrides below still overwrite df, and the + # per-period BM-DOF branch still wins per coefficient on hc2_bm. + # Keep textually parallel with LinearRegression.get_inference's + # "normal" branch (linalg.py). + df = None # Absorbed-FE variance scale (fixest full-K convention): the within- # transform solve_ols above scales the non-clustered classical/hc1 vcov diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 6d168454e..7b4926be9 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -60,7 +60,7 @@ DifferenceInDifferences( bootstrap_weights: str = "rademacher", # "rademacher", "webb", or "mammen" seed: int | None = None, # Random seed rank_deficient_action: str = "warn", # "warn", "error", or "silent" - df_convention: str = "residual", # Clustered t/p/CI df: "residual" (n-K, default) or "cluster" (Stata/fixest G-1); survey df + hc2_bm BM-DOF keep precedence; default flips at v4 + df_convention: str = "residual", # Analytical t/p/CI df: "residual" (n-K, default), "cluster" (Stata/fixest G-1, clustered fits), or "normal" (deliberate z); survey df + hc2_bm BM-DOF keep precedence; default flips at v4 ) ``` @@ -107,7 +107,7 @@ TwoWayFixedEffects( robust: bool = True, cluster: str | None = None, # Auto-clusters at unit level if None alpha: float = 0.05, - df_convention: str = "residual", # Clustered t/p/CI df: "residual" (default) or "cluster" (G-1); flips at v4 + df_convention: str = "residual", # Analytical t/p/CI df: "residual" (default), "cluster" (G-1), or "normal" (z); flips at v4 ) ``` @@ -147,7 +147,7 @@ MultiPeriodDiD( robust: bool = True, cluster: str | None = None, alpha: float = 0.05, - df_convention: str = "residual", # Clustered t/p/CI df: "residual" (default) or "cluster" (G-1); flips at v4 + df_convention: str = "residual", # Analytical t/p/CI df: "residual" (default), "cluster" (G-1), or "normal" (z); flips at v4 ) ``` @@ -372,6 +372,7 @@ SunAbraham( conley_metric: str = "haversine", # "haversine" | "euclidean" conley_kernel: str = "bartlett", # "bartlett" | "uniform" conley_lag_cutoff: int | None = None, # within-unit Bartlett max lag (0 = spatial-only) + df_convention: str = "residual", # Analytical t/p/CI df for cells AND aggregates (3.9: aggregates share the cells' df — previously z): "residual" (default), "cluster" (G-1), "normal" (z); survey/BM DOF keep precedence; flips at v4 ) ``` @@ -419,6 +420,7 @@ ImputationDiD( aux_partition: str = "cohort_horizon", # "cohort_horizon", "cohort", or "horizon" pretrends: bool = False, # Include pre-treatment horizons in event study leave_one_out: bool = False, # BJS 2024 App. A.9 leave-one-out finite-sample variance (larger, less-downward-biased SE) + df_convention: str = "residual", # Pretrends lead-regression per-lead t/p/CI only (3.9: leads moved z -> t(residual)): "residual" (default), "cluster" (G-1), "normal" (z); BJS overall/post inference + Wald F knob-independent; flips at v4 ) ``` @@ -1003,6 +1005,7 @@ StackedDiD( rank_deficient_action: str = "warn", vcov_type: str = "hc1", # {"hc1","hc2_bm"}; classical/hc2 rejected (intrinsically clustered), conley deferred. survey_design=... requires hc1 balance: str = "none", # {"none","entropy"}; "entropy" = CBWSDID covariate balancing (Ustyuzhanin 2026), requires fit(covariates=[...]) + weighting="aggregate", no survey_design + df_convention: str = "residual", # Analytical t/p/CI df (3.9: hc1 lane moved z -> t(pooled residual df)): "residual" (default), "cluster" (G-1, positive-weight clusters), "normal" (z); BM/survey df keep precedence; flips at v4 ) ``` @@ -1099,6 +1102,7 @@ LPDiD( rank_deficient_action: str = "warn", # "warn", "error", or "silent" non_absorbing: str | None = None, # None=absorbing; "first_entry" (Eq. 12); "effect_stabilization" (Eq. 13) stabilization_window: int | None = None, # The paper's L; required when non_absorbing="effect_stabilization" + df_convention: str = "cluster", # DEFAULT "cluster" = the existing Stata lpdid t(G-1) reference (nothing moved in 3.9); "residual" (per-horizon n_eff-k; RA path n_total-k0-1) and "normal" (z) opt in ) ``` @@ -1365,6 +1369,7 @@ WooldridgeDiD( # column is dropped per Sec 5.4, so cohort_trend_coefs # carries G-1 entries. Rejected with survey_design= and # with control_group="never_treated". + df_convention: str = "residual", # OLS analytical t/p/CI df (3.9: default-hc1 arms moved z -> t(residual)): "residual" (default), "cluster" (G-1, hc1-clustered only), "normal" (z); survey/BM DOF keep precedence; GLM arms knob-independent (explicit non-default warns); flips at v4 ) ``` @@ -1602,6 +1607,9 @@ Returned by `SunAbraham.fit()`. | `n_control_units` | `int` | Number of never-treated units | | `control_group` | `str` | Control group type used | | `cohort_effects` | `dict` | Cohort-level effects | +| `df_convention` | `str | None` | The configured df convention echoed onto the results ("residual" | "cluster" | "normal"; 3.9) | +| `inference_df` | `float | None` | The df the overall-ATT p/CI actually used (BM contrast df / survey df / knob-resolved fallback; None under bootstrap or "normal") | +| `event_study_df` | `dict[int, float] | None` | Per-event-time df provenance (finite residual df on plain analytic fits since 3.9; None under bootstrap) | **Methods:** `summary()`, `print_summary()`, `to_dataframe(level="event_study"|"cohort")` @@ -1755,6 +1763,7 @@ Each event study effect dict contains: `effect`, `se`, `t_stat`, `p_value`, `con | `vcov_type` | `str` | Variance estimator family (`"hc1"` for IF-based estimators; permanently narrow on `ImputationDiD` / `CallawaySantAnna` / `TripleDifference` / `EfficientDiD` per IF-vs-sandwich taxonomy) | | `cluster_name` | `str | None` | Effective cluster column name (e.g. `"unit"` for default `cluster=None` on `ImputationDiD`; `None` for default on `EfficientDiD` since its per-unit EIF SE is HC1-style not CR1-at-unit); `None` under survey designs (the survey block already names PSU/strata) | | `n_clusters` | `int | None` | Number of effective clusters; `None` under survey designs and `None` under EfficientDiD's default unclustered fit | +| `df_convention` | `str | None` | On the knob-carrying containers (`StackedDiDResults`, `ImputationDiDResults`, `WooldridgeDiDResults`, `LPDiDResults`, `SunAbrahamResults`): the configured df convention (3.9 / M-127); `StackedDiDResults` additionally carries `inference_df` (the overall-ATT df actually used) | **Methods:** `summary()`, `print_summary()`, `to_dataframe()`, `to_dict()` (flat dict of headline aliases + `vcov_type` + conditional `cluster_name`/`n_clusters`/`n_bootstrap`/`inference_method`) diff --git a/diff_diff/imputation.py b/diff_diff/imputation.py index 817fd7310..7c0892949 100644 --- a/diff_diff/imputation.py +++ b/diff_diff/imputation.py @@ -40,10 +40,13 @@ from diff_diff.utils import ( _iterative_fe_solve, absorbed_fe_cr1_k_increment, + absorbed_fe_rank, demean_by_groups, pre_demean_norms, + resolve_tail_df, safe_inference, snap_absorbed_regressors, + validate_df_convention, ) if TYPE_CHECKING: @@ -222,6 +225,21 @@ class ImputationDiD(ImputationDiDBootstrapMixin): but are a library extension beyond the paper's derivation. Replicate-weight survey designs raise ``NotImplementedError`` (their variance bypasses the influence-function path where the rescale lives). + df_convention : {"residual", "cluster", "normal"}, default "residual" + Degrees-of-freedom convention for the PRETRENDS lead regression's + per-lead t/p/CI — the one ImputationDiD surface running the shared + clustered CR1 sandwich (``pretrends=True`` with + ``aggregate="event_study"``/``"all"``). ``"residual"`` (default) uses + the lead regression's residual df (``n − k_kept − absorbed + [time, unit] rank``) — the 3.9 fix: previously silent normal-theory + z on plain clustered fits; ``"cluster"`` uses ``G − 1``; + ``"normal"`` deliberately uses z. The full-design survey df keeps + precedence on survey fits. Everything else — the BJS Theorem-3 + overall/event-study inference and the joint pretrend Wald F (which + keeps its cluster-robust ``F(q, G − 1)`` reference) — is + knob-independent; an explicitly non-default value on a + configuration that never surfaces the per-lead inference warns at + fit time. The default flips to ``"cluster"`` at v4. Attributes ---------- @@ -278,6 +296,7 @@ def __init__( aux_partition: str = "cohort_horizon", pretrends: bool = False, leave_one_out: bool = False, + df_convention: str = "residual", ): if rank_deficient_action not in ("warn", "error", "silent"): raise ValueError( @@ -296,6 +315,7 @@ def __init__( ) self._validate_vcov_type(vcov_type) self._validate_leave_one_out(leave_one_out) + validate_df_convention(df_convention) self.anticipation = anticipation self.alpha = alpha @@ -309,6 +329,7 @@ def __init__( self.aux_partition = aux_partition self.pretrends = pretrends self.leave_one_out = leave_one_out + self.df_convention = df_convention self.is_fitted_ = False self.results_: Optional[ImputationDiDResults] = None @@ -399,6 +420,25 @@ def fit( "are supported. Use pretrends=False with replicate weights." ) + # Inert-config warning (no-silent-failures): the df_convention knob + # moves only the pretrends lead regression's per-lead t/p/CI, which + # this fit surfaces only when pretrends=True with + # aggregate="event_study"/"all". (The post-fit ``pretrend_test()`` + # reaches the lead helper too, but consumes only gamma/V_gamma — its + # joint Wald F denominator is knob-independent.) An explicitly + # non-default value on any other configuration is a silent no-op. + if self.df_convention != "residual" and not ( + self.pretrends and aggregate in ("event_study", "all") + ): + warnings.warn( + f"df_convention={self.df_convention!r} affects only the " + "pretrends event-study per-lead inference (pretrends=True " + "with aggregate='event_study'/'all'); it has no effect on " + "this configuration.", + UserWarning, + stacklevel=2, + ) + # Create working copy df = data.copy() @@ -1089,6 +1129,7 @@ def _refit_imp(w_r): cluster_name=_cluster_name_for_results, n_clusters=_n_clusters_for_results, leave_one_out=self.leave_one_out, + df_convention=self.df_convention, ) self.is_fitted_ = True @@ -2511,8 +2552,38 @@ def _compute_lead_coefficients( gamma = coefficients[:n_leads] V_gamma = vcov[:n_leads, :n_leads] - # Use full-design survey df for t-distribution inference - _df = survey_df + # Resolve the per-lead tail df: the full-design survey df keeps + # precedence on survey fits; otherwise resolve through the + # df_convention knob (3.9 fix: previously df=None → silent + # normal-theory z on plain clustered fits). Residual df = n − + # k_kept − absorbed [time, unit] rank (the demeaned lead design + # carries no intercept column); the "cluster" G counts + # positive-weight clusters on the SAME raw cluster ids the CR1 + # vcov partitioned on. NOTE this resolution serves BOTH callers — + # the fit-time pretrends event-study path and the post-fit public + # ``pretrend_test()`` (which consumes only gamma/V_gamma; its + # joint Wald F denominator is knob-independent). + if _use_survey_vcov: + _df = survey_df + else: + _k_kept_imp = int(np.count_nonzero(np.isfinite(coefficients))) + _df_res_imp = float( + len(df_0) + - _k_kept_imp + - absorbed_fe_rank( + df_0, + [time, unit], + has_intercept_col=False, + weights=survey_weights_0, + ) + ) + from diff_diff.linalg import effective_cluster_count + + _df = resolve_tail_df( + self.df_convention, + residual_df=_df_res_imp, + n_clusters=effective_cluster_count(cluster_ids, survey_weights_0), + ) # Build per-horizon effects effects = {} @@ -2692,15 +2763,20 @@ def get_params(self) -> Dict[str, Any]: "aux_partition": self.aux_partition, "pretrends": self.pretrends, "leave_one_out": self.leave_one_out, + "df_convention": self.df_convention, } def set_params(self, **params) -> "ImputationDiD": """Set estimator parameters (sklearn-compatible).""" - for key, value in params.items(): - if hasattr(self, key): - setattr(self, key, value) - else: + # Reject unknown keys and validate the pending df_convention BEFORE + # any assignment so a rejected call leaves the estimator unchanged. + for key in params: + if not hasattr(self, key): raise ValueError(f"Unknown parameter: {key}") + if "df_convention" in params: + validate_df_convention(params["df_convention"]) + for key, value in params.items(): + setattr(self, key, value) return self @staticmethod diff --git a/diff_diff/imputation_results.py b/diff_diff/imputation_results.py index 3305caf93..161c36f16 100644 --- a/diff_diff/imputation_results.py +++ b/diff_diff/imputation_results.py @@ -154,6 +154,12 @@ class ImputationDiDResults(BaseResults): # BJS 2024 Supp. App. A.9 leave-one-out finite-sample variance refinement # (opt-in). Recorded here so reported SEs are self-describing. leave_one_out: bool = field(default=False) + # The estimator's df_convention configuration echoed onto the results + # ("residual" | "cluster" | "normal"; added 3.9). It governs only the + # pretrends lead regression's per-lead t/p/CI - the BJS overall / + # post-treatment inference is knob-independent. Appended LAST (the + # generated __init__ positional indexes are public API). + df_convention: Optional[str] = None # --- Inference-field aliases (balance/external-adapter compatibility) --- @property @@ -516,6 +522,8 @@ def to_dict(self) -> Dict[str, Any]: result["cluster_name"] = self.cluster_name if self.n_clusters is not None: result["n_clusters"] = self.n_clusters + if self.df_convention is not None: + result["df_convention"] = self.df_convention if self.bootstrap_results is not None: result["n_bootstrap"] = self.bootstrap_results.n_bootstrap result["inference_method"] = "bootstrap" diff --git a/diff_diff/linalg.py b/diff_diff/linalg.py index eb36530c2..c6d3f7d2a 100644 --- a/diff_diff/linalg.py +++ b/diff_diff/linalg.py @@ -4367,15 +4367,18 @@ class LinearRegression: neither is formally PSD-guaranteed in the radial pairwise form (Conley 1999's explicit PSD Bartlett formula is the 2-D separable product window, Eq 3.14, not the 1-D radial pairwise form). - df_convention : {"residual", "cluster"}, default "residual" - Degrees-of-freedom convention for ``get_inference`` t/p/CI on - clustered fits. ``"residual"`` uses the fitted residual df; - ``"cluster"`` uses the Stata/fixest cluster df ``G − 1`` (from - ``n_clusters_``). Fallback-level only: survey df and per-coefficient - Bell-McCaffrey DOF always take precedence. No effect on - coefficients, SEs, unclustered fits, or ``vcov_type="conley"`` (no - documented ``G − 1`` reference for the Conley+cluster product - kernel). Default flips at v4. + df_convention : {"residual", "cluster", "normal"}, default "residual" + Degrees-of-freedom convention for ``get_inference`` t/p/CI. + ``"residual"`` uses the fitted residual df; ``"cluster"`` uses the + Stata/fixest cluster df ``G − 1`` (from ``n_clusters_``) — + ``"cluster"`` has no effect on unclustered fits or + ``vcov_type="conley"`` (no documented ``G − 1`` reference for the + Conley+cluster product kernel); ``"normal"`` deliberately uses + normal-theory z inference at the fallback level on every fit, + clustered or not. Fallback-level only under every value: survey df + and per-coefficient Bell-McCaffrey DOF always take precedence. No + effect on coefficients, SEs, or t-statistics. Default flips to + ``"cluster"`` at v4. Attributes ---------- @@ -4448,9 +4451,13 @@ def __init__( conley_lag_cutoff: Optional[int] = None, df_convention: str = "residual", ): - if df_convention not in ("residual", "cluster"): + # Inline validation (not utils.validate_df_convention): diff_diff.utils + # imports diff_diff.linalg at module level, so importing back would be + # circular. Keep the message shape identical to the shared validator. + if df_convention not in ("residual", "cluster", "normal"): raise ValueError( - f"df_convention must be 'residual' or 'cluster', got {df_convention!r}" + "df_convention must be one of ('residual', 'cluster', 'normal'), " + f"got {df_convention!r}" ) self.include_intercept = include_intercept self.robust = robust @@ -5221,6 +5228,14 @@ def get_inference( alpha=effective_alpha, ) effective_df = self.n_clusters_ - 1 + elif self.df_convention == "normal": + # Deliberate normal-theory z inference at the fallback level, on + # every fit (clustered, unclustered, and conley alike). Survey df + # and per-coefficient Bell-McCaffrey DOF above still win — this + # branch replaces only the residual fallback. Keep textually + # parallel with MultiPeriodDiD's inline ladder (estimators.py, + # "normal" branch). + effective_df = None else: effective_df = self.df_ diff --git a/diff_diff/lpdid.py b/diff_diff/lpdid.py index 4afd0dbe5..e98d38aed 100644 --- a/diff_diff/lpdid.py +++ b/diff_diff/lpdid.py @@ -6,7 +6,13 @@ from diff_diff.linalg import InvalidClusterKAdjustment, _rank_guarded_inv, solve_ols from diff_diff.lpdid_results import LPDiDResults -from diff_diff.utils import absorbed_fe_rank, cluster_nested_fe_dims, safe_inference +from diff_diff.utils import ( + absorbed_fe_rank, + cluster_nested_fe_dims, + resolve_tail_df, + safe_inference, + validate_df_convention, +) __all__ = ["LPDiD", "LPDiDResults"] @@ -25,6 +31,7 @@ def __init__( rank_deficient_action: str = "warn", non_absorbing: Optional[str] = None, stabilization_window: Optional[int] = None, + df_convention: str = "cluster", ): self.pre_window = pre_window self.post_window = post_window @@ -37,6 +44,7 @@ def __init__( self.rank_deficient_action = rank_deficient_action self.non_absorbing = non_absorbing self.stabilization_window = stabilization_window + self.df_convention = df_convention self._validate_params() self.is_fitted_ = False self.results_: Optional[LPDiDResults] = None @@ -85,6 +93,10 @@ def _validate_params(self) -> None: "control_group='never_treated' is not supported with a non-absorbing mode " "(the estimand becomes ambiguous); use control_group='clean' (the default)" ) + # LPDiD's DEFAULT is "cluster" — its t(G-1) reference is the Stata + # lpdid convention and already the library-wide 4.0 target, so the + # default moves nothing (unlike the "residual"-default estimators). + validate_df_convention(self.df_convention) def _rhs_column_names(self, covariates=None, ylags=0, dylags=0): rhs_columns = list(covariates or []) @@ -683,6 +695,9 @@ def _estimate_regression_adjustment_sample( # as 0 for prediction/residuals; without this zero-fill the NaN would # propagate through every prediction and NaN an otherwise-identified ATT. # ("error" still raises inside solve_ols before returning.) + # Kept nuisance-parameter count for the "residual" df, captured + # BEFORE the zero-fill below erases the dropped-column NaN markers. + k0_kept = int(np.count_nonzero(np.isfinite(control_coef))) control_coef = np.where(np.isfinite(control_coef), control_coef, 0.0) treated_design = np.column_stack( @@ -695,8 +710,11 @@ def _estimate_regression_adjustment_sample( se = np.nan cluster_ids = sample["_cluster"].to_numpy() n_clusters = len(pd.unique(cluster_ids)) + # Hoisted above the branch: the "residual" df resolution below needs + # n_total on every lane (including G>=2 with a non-finite effect, + # where today's G-1 df is still reported alongside the NaN se). + n_total = len(sample) if n_clusters >= 2 and np.isfinite(effect): - n_total = len(sample) n_treated = len(treated) q0_inv, _, _ = _rank_guarded_inv(control_design.T @ control_design) mu_treated = treated_design.mean(axis=0) @@ -715,7 +733,25 @@ def _estimate_regression_adjustment_sample( if np.isfinite(vcov_scalar) and vcov_scalar >= 0: se = float(np.sqrt(vcov_scalar)) - df = n_clusters - 1 if n_clusters > 1 else None + # Tail df: the resolver runs whenever n_clusters >= 2 — matching + # today's expression `n_clusters - 1 if n_clusters > 1 else None`, + # which applies G-1 even when the effect is non-finite — so the + # default "cluster" reproduces today's values bit-for-bit. The + # G<=1 degenerate lane keeps the literal df=None (normal theory + + # NaN provenance) under ALL conventions and never touches the + # resolver. "residual" = n_total - k0_kept - 1: the RA contrast is + # a pooled M-estimator over all n_total rows estimating the k0 + # kept nuisance coefficients PLUS the ATT (library convention, no + # external anchor - Stata teffects ra reports z; see the REGISTRY + # LPDiD note for the derivation). + if n_clusters >= 2: + df = resolve_tail_df( + self.df_convention, + residual_df=float(n_total - k0_kept - 1), + n_clusters=n_clusters, + ) + else: + df = None t_stat, p_value, conf_int = safe_inference(effect, se, alpha=self.alpha, df=df) return { "coefficient": effect, @@ -918,7 +954,29 @@ def _estimate_sample( se = float(np.sqrt(vcov[1, 1])) n_clusters = len(pd.unique(cluster_ids)) - df = n_clusters - 1 if vcov is not None and n_clusters > 1 else None + # Tail df: the resolver runs only on the healthy clustered lane — + # matching today's expression `n_clusters - 1 if vcov is not None + # and n_clusters > 1 else None` — so the default "cluster" + # reproduces today's values bit-for-bit (raw-unique G, deliberately + # NOT effective_cluster_count: LPDiD's reweights are strictly + # positive, so the two cannot diverge — REGISTRY LPDiD note). The + # degenerate lanes (vcov is None: the unclustered-refit fallback + # above, or G<=1) keep the literal df=None (normal theory + NaN + # provenance) under ALL conventions and never touch the resolver. + # "residual" = n_eff - k_kept of this per-horizon design. + if vcov is not None and n_clusters >= 2: + _n_eff_lp = ( + int(np.count_nonzero(np.asarray(weights, dtype=float) > 0)) + if weights is not None + else int(design.shape[0]) + ) + df = resolve_tail_df( + self.df_convention, + residual_df=float(_n_eff_lp - int(np.count_nonzero(np.isfinite(coef)))), + n_clusters=n_clusters, + ) + else: + df = None t_stat, p_value, conf_int = safe_inference(effect, se, alpha=self.alpha, df=df) return { "coefficient": effect, @@ -1598,6 +1656,7 @@ def fit( pooled=pooled, event_study_df=event_study_df, pooled_df=pooled_df, + df_convention=self.df_convention, n_obs=len(data), n_treated_units=int(treatment_by_unit.gt(0).sum()), n_control_units=int(treatment_by_unit.eq(0).sum()), @@ -1647,16 +1706,19 @@ def get_params(self) -> Dict[str, Any]: "rank_deficient_action": self.rank_deficient_action, "non_absorbing": self.non_absorbing, "stabilization_window": self.stabilization_window, + "df_convention": self.df_convention, } def set_params(self, **params: Any) -> "LPDiD": + # Reject unknown keys before any assignment (the rollback below + # only covers validator failures on known keys). + for key in params: + if not hasattr(self, key): + raise ValueError(f"Unknown parameter: {key}") previous_values = {} for key, value in params.items(): - if hasattr(self, key): - previous_values[key] = getattr(self, key) - setattr(self, key, value) - else: - raise ValueError(f"Unknown parameter: {key}") + previous_values[key] = getattr(self, key) + setattr(self, key, value) try: self._validate_params() except ValueError: diff --git a/diff_diff/lpdid_results.py b/diff_diff/lpdid_results.py index 7164aa65d..a54d2c747 100644 --- a/diff_diff/lpdid_results.py +++ b/diff_diff/lpdid_results.py @@ -78,6 +78,11 @@ class LPDiDResults(BaseResults): # per window when that window's inference used normal theory or an # undefined df; None when no pooled windows were fit (``only_event``). pooled_df: Optional[Dict[str, float]] = None + # The estimator's df_convention configuration echoed onto the results + # ("residual" | "cluster" | "normal"; LPDiD's default is "cluster" - + # the Stata lpdid t(G-1) convention; added 3.9). Appended LAST (the + # generated __init__ positional indexes are public API). + df_convention: Optional[str] = None # ------------------------------------------------------------------ # internal helpers @@ -189,6 +194,8 @@ def to_dict(self) -> Dict[str, Any]: result["n_psu"] = self.n_psu result["weight_type"] = getattr(self.survey_metadata, "weight_type", None) result["df_survey"] = getattr(self.survey_metadata, "df_survey", None) + if self.df_convention is not None: + result["df_convention"] = self.df_convention result["inference_method"] = ( "survey_tsl" if self.vcov_type == "survey_tsl" else "cluster_robust" ) diff --git a/diff_diff/stacked_did.py b/diff_diff/stacked_did.py index cd07c6d9c..9060366b4 100644 --- a/diff_diff/stacked_did.py +++ b/diff_diff/stacked_did.py @@ -25,9 +25,9 @@ import pandas as pd from diff_diff.balancing import BalanceError, entropy_balance -from diff_diff.linalg import solve_ols +from diff_diff.linalg import effective_cluster_count, solve_ols from diff_diff.stacked_did_results import StackedDiDResults # noqa: F401 (re-export) -from diff_diff.utils import safe_inference +from diff_diff.utils import resolve_tail_df, safe_inference, validate_df_convention __all__ = [ "StackedDiD", @@ -132,6 +132,21 @@ class StackedDiD: raise a ``ValueError``), and does not support ``survey_design=``; matching-based balancing and the repeated-treatment extension are out of scope. Default ``"none"`` reproduces plain weighted stacked DID. + df_convention : {"residual", "cluster", "normal"}, default "residual" + Degrees-of-freedom convention for the analytical t/p/CI on the + pooled stacked regression's event-study and overall-ATT inference. + ``"residual"`` (default) uses the pooled residual df + (``n_eff − k_kept``, positive-weight rows) — the 3.9 fix: the + non-BM non-survey lane previously used silent normal-theory z; + ``"cluster"`` uses ``G − 1`` where G counts positive-weight + clusters (``effective_cluster_count`` — may differ from the raw + ``results.n_clusters`` when a cluster's composed weight is zero); + ``"normal"`` deliberately uses normal-theory z at the fallback + level. hc2_bm Bell-McCaffrey contrast DOF and survey/replicate df + always take precedence. Note the stacked design's residual df is + typically large (control rows replicate across sub-experiments), so + the default's numeric movement vs the old z is small — the change + is a convention alignment. The default flips to ``"cluster"`` at v4. Attributes ---------- @@ -184,6 +199,7 @@ def __init__( rank_deficient_action: str = "warn", vcov_type: str = "hc1", balance: str = "none", + df_convention: str = "residual", ): if weighting not in ("aggregate", "population", "sample_share"): raise ValueError( @@ -206,6 +222,7 @@ def __init__( # Factored into _validate_vcov_type so set_params() can re-validate. self._validate_vcov_type(vcov_type) self._validate_balance(balance) + validate_df_convention(df_convention) self.kappa_pre = kappa_pre self.kappa_post = kappa_post @@ -217,6 +234,7 @@ def __init__( self.rank_deficient_action = rank_deficient_action self.vcov_type = vcov_type self.balance = balance + self.df_convention = df_convention self.is_fitted_ = False self.results_: Optional[StackedDiDResults] = None @@ -617,6 +635,29 @@ def fit( ) assert vcov is not None + # Knob-resolved analytical fallback df (3.9 tail-df consolidation): + # the pooled stacked design is fully visible (no absorption), so the + # residual df is n_eff − k_kept with n_eff = positive-weight rows + # (pweight semantics, mirroring LinearRegression.fit's n_eff). The + # ``df_convention="cluster"`` G counts POSITIVE-WEIGHT clusters + # (effective_cluster_count) — deliberately not the raw unique count + # reported as ``results.n_clusters``; the two diverge only when a + # cluster's total composed weight is zero (REGISTRY StackedDiD + # note). Survey/replicate df keeps precedence at the consuming + # branches below. Previously the non-BM non-survey lane silently + # used normal theory (the tail-df defect family). + _n_eff_sd = int(np.count_nonzero(np.asarray(composed_weights, dtype=float) > 0)) + _k_kept_sd = int(np.count_nonzero(~np.isnan(coef))) + _analytic_fallback_df = resolve_tail_df( + self.df_convention, + residual_df=float(_n_eff_sd - _k_kept_sd), + n_clusters=( + effective_cluster_count(cluster_ids, composed_weights) + if self.vcov_type == "hc1" + else None + ), + ) + # Bell-McCaffrey Satterthwaite contrast DOF for hc2_bm. Per the # registry contract for `vcov_type="hc2_bm"`, the user-facing # aggregated inference (event_study_effects[h]['p_value']/['conf_int'] @@ -854,7 +895,15 @@ def _refit_stacked(w_r): # No safe_inference call happened -> no df provenance. es_df_used[h] = float("nan") else: - _df_eff = _bm_df if _bm_df is not None else _survey_df + # 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 ) @@ -924,6 +973,7 @@ def _refit_stacked(w_r): # from PR #475 R7). Without this, normal-theory fallback would # silently produce wrong p-values/CIs on the overall_* surface. _is_hc2bm_path_overall = self.vcov_type == "hc2_bm" and not _uses_replicate_sd + _overall_df_used: Optional[float] = None if _is_hc2bm_path_overall and ( _bm_contrast_dof_overall is None or not np.isfinite(_bm_contrast_dof_overall) ): @@ -934,11 +984,17 @@ def _refit_stacked(w_r): _df_overall_eff = ( _bm_contrast_dof_overall if _bm_contrast_dof_overall is not None - else _survey_df_overall + else ( + _survey_df_overall if _survey_df_overall is not None else _analytic_fallback_df + ) ) overall_t, overall_p, overall_ci = safe_inference( overall_att, overall_se, alpha=self.alpha, df=_df_overall_eff ) + # Scalar provenance: the df the overall safe_inference actually + # received, recorded iff it governed a t-reference (finite, > 0). + if _df_overall_eff is not None and np.isfinite(_df_overall_eff) and _df_overall_eff > 0: + _overall_df_used = float(_df_overall_eff) # ---- Construct results ---- self.results_ = StackedDiDResults( @@ -974,6 +1030,8 @@ def _refit_stacked(w_r): event_study_vcov=es_vcov, event_study_vcov_index=es_vcov_index, event_study_df=es_df_used, + df_convention=self.df_convention, + inference_df=_overall_df_used, ) self.is_fitted_ = True @@ -1501,6 +1559,7 @@ def get_params(self) -> Dict[str, Any]: "rank_deficient_action": self.rank_deficient_action, "vcov_type": self.vcov_type, "balance": self.balance, + "df_convention": self.df_convention, } def set_params(self, **params: Any) -> "StackedDiD": @@ -1511,17 +1570,20 @@ def set_params(self, **params: Any) -> "StackedDiD": before fit() (avoids a later, less-informative failure in the linalg layer). """ - # Validate vcov_type up-front if it's being set, so the same - # error surface as __init__ applies. + # Reject unknown keys, then validate pending values up-front (same + # error surface as __init__), so a rejected call leaves the + # estimator unchanged. + for key in params: + if not hasattr(self, key): + raise ValueError(f"Unknown parameter: {key}") if "vcov_type" in params: self._validate_vcov_type(params["vcov_type"]) if "balance" in params: self._validate_balance(params["balance"]) + if "df_convention" in params: + validate_df_convention(params["df_convention"]) for key, value in params.items(): - if hasattr(self, key): - setattr(self, key, value) - else: - raise ValueError(f"Unknown parameter: {key}") + setattr(self, key, value) return self def summary(self) -> str: diff --git a/diff_diff/stacked_did_results.py b/diff_diff/stacked_did_results.py index 1774668ea..e21e6b099 100644 --- a/diff_diff/stacked_did_results.py +++ b/diff_diff/stacked_did_results.py @@ -89,9 +89,20 @@ class StackedDiDResults(BaseResults): estimated event time to the df actually passed to ``safe_inference`` for its stored p-value/CI (per-event Bell-McCaffrey Satterthwaite df under ``hc2_bm``; the scalar survey - df under survey designs), or NaN when the row used normal theory, - the df was undefined, or hc2_bm failed closed. None when no event - study was requested. + df under survey designs; the ``df_convention``-resolved analytical + 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. + df_convention : str, optional + The estimator's ``df_convention`` configuration echoed onto the + results ("residual" | "cluster" | "normal"; added 3.9). + inference_df : float, optional + The df the stored overall-ATT p-value/CI's ``safe_inference`` + actually received: the BM contrast df under ``hc2_bm``, the + survey/replicate df on survey fits, else the + ``df_convention``-resolved analytical fallback. None when the + overall inference used normal theory or failed closed. """ overall_att: float @@ -145,6 +156,9 @@ class StackedDiDResults(BaseResults): event_study_vcov: Optional[np.ndarray] = field(default=None, repr=False) event_study_vcov_index: Optional[List[int]] = field(default=None, repr=False) event_study_df: Optional[Dict[int, float]] = field(default=None, repr=False) + # Appended LAST (generated __init__ positional indexes are public API). + df_convention: Optional[str] = None + inference_df: Optional[float] = None # --- Inference-field aliases (balance/external-adapter compatibility) --- @property @@ -404,6 +418,10 @@ def to_dict(self) -> Dict[str, Any]: result["cluster_name"] = self.cluster_name if self.n_clusters is not None: result["n_clusters"] = self.n_clusters + if self.df_convention is not None: + result["df_convention"] = self.df_convention + if self.inference_df is not None: + result["inference_df"] = self.inference_df return result def to_dataframe(self, level: str = "event_study") -> pd.DataFrame: diff --git a/diff_diff/sun_abraham.py b/diff_diff/sun_abraham.py index b646431c1..784722d3c 100644 --- a/diff_diff/sun_abraham.py +++ b/diff_diff/sun_abraham.py @@ -11,7 +11,18 @@ import warnings from dataclasses import dataclass, field -from typing import TYPE_CHECKING, Any, Dict, List, Literal, Optional, Tuple, Union, overload +from typing import ( + TYPE_CHECKING, + Any, + Dict, + List, + Literal, + NamedTuple, + Optional, + Tuple, + Union, + overload, +) import numpy as np import pandas as pd @@ -27,14 +38,31 @@ absorbed_fe_cr1_k_increment, absorbed_fe_rank, pre_demean_norms, + resolve_tail_df, safe_inference, snap_absorbed_regressors, + validate_df_convention, ) from diff_diff.utils import ( within_transform as _within_transform_util, ) +class _SaturatedFitStats(NamedTuple): + """Fit-level df ingredients read off the saturated LinearRegression. + + Carried out of ``_fit_saturated_regression`` so the aggregate inference + layer can resolve the ``df_convention`` fallback df from the SAME fit the + cells used (the D4 fix: one df source per fit). ``df_residual`` is the + regression's ``df_`` (n_eff − k_effective − absorbed rank); + ``n_clusters`` its effective positive-weight cluster count (None on + unclustered fits). + """ + + df_residual: Optional[float] + n_clusters: Optional[int] + + @dataclass class SunAbrahamResults(BaseResults): """ @@ -157,20 +185,33 @@ class SunAbrahamResults(BaseResults): # PROVENANCE - maps each estimated relative time to the df its stored # p-value/CI's safe_inference actually received (the per-event # Bell-McCaffrey contrast df under hc2_bm; the survey design df - - # post-drop under replicate refits - on survey fits; NaN on plain - # analytic fits, which use normal theory, and on rows whose BM DOF was - # non-finite, where safe_inference's non-finite-df guard yields all-NaN - # inference). None under bootstrap: the stored percentile p/CIs never - # used a df (note this clears the WHOLE channel even when a partial - # bootstrap override leaves some rows analytic - a conservative - # under-claim, consistent with the other producers). Deliberately - # narrower clearing than event_study_vcov above: replicate refits KEEP - # the df (it genuinely governed the recomputed rows) while the vcov - # clears. - # Declared LAST so every pre-existing field keeps its positional index - # in the generated __init__ (the constructor signature is public API). + # post-drop under replicate refits - on survey fits; the + # df_convention-resolved fallback on plain analytic fits - FINITE + # residual df under the 3.9 default, G-1 under "cluster", NaN under + # "normal" and on rows whose BM DOF was non-finite, where + # safe_inference's non-finite-df guard yields all-NaN inference). None + # under bootstrap: the stored percentile p/CIs never used a df (note + # this clears the WHOLE channel even when a partial bootstrap override + # leaves some rows analytic - a conservative under-claim, consistent + # with the other producers). Deliberately narrower clearing than + # event_study_vcov above: replicate refits KEEP the df (it genuinely + # governed the recomputed rows) while the vcov clears. + # This block is appended-only so every pre-existing field keeps its + # positional index in the generated __init__ (the constructor signature + # is public API); new fields go BELOW. event_study_df: Optional[Dict[int, float]] = field(default=None, repr=False) + df_convention: Optional[str] = None + """The estimator's ``df_convention`` configuration echoed onto the + results ("residual" | "cluster" | "normal"; added 3.9).""" + + inference_df: Optional[float] = None + """The df the stored overall-ATT p-value/CI's ``safe_inference`` + actually received: the BM contrast df under hc2_bm, the survey design + df on survey fits, else the ``df_convention``-resolved analytical + fallback (None under "normal" = normal theory). None on bootstrap fits, + whose overall p/CI are percentile-based and never used a df.""" + # --- Inference-field aliases (balance/external-adapter compatibility) --- @property def att(self) -> float: @@ -353,6 +394,10 @@ def to_dict(self) -> Dict[str, Any]: result["cluster_name"] = self.cluster_name if self.conley_lag_cutoff is not None: result["conley_lag_cutoff"] = self.conley_lag_cutoff + if self.df_convention is not None: + result["df_convention"] = self.df_convention + if self.inference_df is not None: + result["inference_df"] = self.inference_df return result def to_dataframe(self, level: str = "event_study") -> pd.DataFrame: @@ -562,6 +607,18 @@ class SunAbraham: ``conley_cutoff_km`` / ``conley_lag_cutoff``); ``survey_design=`` / ``weights`` / ``n_bootstrap>0`` are rejected. See the ``vcov_type`` parameter docs above. + df_convention : {"residual", "cluster", "normal"}, default "residual" + Degrees-of-freedom convention for analytical t/p/CI, applied to BOTH + the per-cohort-cell inference and the aggregated event-study / + overall-ATT inference (one df source per fit). ``"residual"`` + (default) uses the saturated regression's residual df — the 3.9 fix: + aggregates previously dropped to normal theory on plain clustered + fits; ``"cluster"`` uses the Stata/fixest cluster df ``G − 1`` + (inert under ``vcov_type="conley"``); ``"normal"`` deliberately uses + normal-theory z at the fallback level everywhere (cells included). + Survey/replicate df and hc2_bm Bell-McCaffrey contrast DOF always + take precedence; bootstrap p/CI (``n_bootstrap>0``) are percentile- + based and unaffected. The default flips to ``"cluster"`` at v4. Attributes ---------- @@ -645,6 +702,7 @@ def __init__( conley_metric: str = "haversine", conley_kernel: str = "bartlett", conley_lag_cutoff: Optional[int] = None, + df_convention: str = "residual", ): if control_group not in ["never_treated", "not_yet_treated"]: raise ValueError( @@ -652,6 +710,8 @@ def __init__( f"got '{control_group}'" ) + validate_df_convention(df_convention) + if rank_deficient_action not in ["warn", "error", "silent"]: raise ValueError( f"rank_deficient_action must be 'warn', 'error', or 'silent', " @@ -677,6 +737,7 @@ def __init__( self.conley_metric = conley_metric self.conley_kernel = conley_kernel self.conley_lag_cutoff = conley_lag_cutoff + self.df_convention = df_convention # Track whether the user explicitly opted out of the "hc1" default. # The auto-cluster-at-unit default in `fit` is suppressed only when # the user explicitly opts into a one-way family — currently @@ -988,6 +1049,7 @@ def fit( vcov_cohort, coef_index_map, bm_artifacts, + _sa_fit_stats, ) = self._fit_saturated_regression( df_reg, outcome, @@ -1023,7 +1085,7 @@ def _refit_sa(w_r): nz = w_r > 0 df_reg_nz = df_reg[nz] if not np.all(nz) else df_reg w_nz = w_r[nz] if not np.all(nz) else w_r - ce_r, _, vcov_r, cim_r, _ = self._fit_saturated_regression( + ce_r, _, vcov_r, cim_r, _, _ = self._fit_saturated_regression( df_reg_nz, outcome, unit, @@ -1089,6 +1151,24 @@ def _refit_sa(w_r): if _uses_replicate_sa and _sa_survey_df is None: _sa_survey_df = 0 # rank-deficient replicate → NaN inference + # Knob-resolved fallback df for aggregated inference (the D4 fix: + # aggregates share the cells' df source instead of dropping to + # normal theory). Survey/replicate df keeps precedence (including + # the fail-closed 0 sentinel above); otherwise resolve from the + # SAME saturated fit the cells used. Conley fits pass + # n_clusters=None — the combined Conley+cluster product kernel has + # no documented G-1 df reference (mirrors get_inference's conley + # exclusion), so "cluster" stays inert there. + _sa_fallback_df = ( + float(_sa_survey_df) + if _sa_survey_df is not None + else resolve_tail_df( + self.df_convention, + residual_df=_sa_fit_stats.df_residual, + n_clusters=(None if self.vcov_type == "conley" else _sa_fit_stats.n_clusters), + ) + ) + # Compute interaction-weighted event study effects event_study_effects, cohort_weights = self._compute_iw_effects( df, @@ -1101,20 +1181,23 @@ def _refit_sa(w_r): vcov_cohort, coef_index_map, survey_weight_col=survey_weight_col, - survey_df=_sa_survey_df, + fallback_df=_sa_fallback_df, ) # Per-row df PROVENANCE (spec section 5, row M-092): the df each # stored ES row's safe_inference actually received, recorded iff - # finite and > 0 else NaN. Baseline = the scalar every - # _compute_iw_effects row used (None on plain analytic fits -> NaN; - # the TSL design df on survey fits). Overwritten below at the hc2_bm - # and replicate-refit override sites, and cleared under bootstrap, - # whose stored percentile p/CIs never used a df. + # finite and > 0 else NaN. Baseline = the knob-resolved fallback + # every _compute_iw_effects row used (the residual df on plain + # analytic fits under the default; G-1 under "cluster"; NaN under + # "normal"; the TSL design df on survey fits). Overwritten below at + # the hc2_bm and replicate-refit override sites, and cleared under + # bootstrap, whose stored percentile p/CIs never used a df. es_df_used: Dict[int, float] = { e: ( - float(_sa_survey_df) - if _sa_survey_df is not None and np.isfinite(_sa_survey_df) and _sa_survey_df > 0 + float(_sa_fallback_df) + if _sa_fallback_df is not None + and np.isfinite(_sa_fallback_df) + and _sa_fallback_df > 0 else float("nan") ) for e in event_study_effects @@ -1228,15 +1311,16 @@ def _refit_sa(w_r): _overall_att_contrast_dof = float(dof_vec[-1]) except (ValueError, np.linalg.LinAlgError) as exc: # Rank-deficient or other linalg issue: fall back to - # the shared analytical df (downgraded to normal - # inference). Emit a UserWarning so the deviation is + # the knob-resolved analytical df (residual t by + # default). Emit a UserWarning so the deviation is # visible. warnings.warn( f"SunAbraham(vcov_type='hc2_bm') aggregated inference " f"could not compute Bell-McCaffrey contrast DOF " - f"({type(exc).__name__}: {exc}). Falling back to " - "shared df; aggregated p-values/CIs may use normal " - "distribution instead of t(BM DOF).", + f"({type(exc).__name__}: {exc}). Falling back to the " + "df_convention-resolved analytical df; aggregated " + "p-values/CIs use t(residual df) under the default " + "instead of t(BM DOF).", UserWarning, stacklevel=2, ) @@ -1245,7 +1329,7 @@ def _refit_sa(w_r): # Override the per-event-time inference fields with BM-DOF-aware # values when available; otherwise leave the `safe_inference` # output from `_compute_iw_effects` in place (which used - # `df=_sa_survey_df`). + # `df=_sa_fallback_df`). if _es_contrast_dofs: for e, df_e in _es_contrast_dofs.items(): eff_e = event_study_effects[e]["effect"] @@ -1264,15 +1348,14 @@ def _refit_sa(w_r): else float("nan") ) + _overall_df_used: Optional[float] = ( + _overall_att_contrast_dof if _overall_att_contrast_dof is not None else _sa_fallback_df + ) overall_t, overall_p, overall_ci = safe_inference( overall_att, overall_se, alpha=self.alpha, - df=( - _overall_att_contrast_dof - if _overall_att_contrast_dof is not None - else _sa_survey_df - ), + df=_overall_df_used, ) # Replicate variance override: refit fully re-aggregated estimates @@ -1300,6 +1383,7 @@ def _refit_sa(w_r): # Override overall ATT SE overall_se = float(np.sqrt(max(_vcov_sa[0, 0], 0.0))) + _overall_df_used = _sa_survey_df overall_t, overall_p, overall_ci = safe_inference( overall_att, overall_se, alpha=self.alpha, df=_sa_survey_df ) @@ -1334,7 +1418,7 @@ def _refit_sa_cohort(w_r): nz = w_r > 0 df_reg_nz = df_reg[nz] if not np.all(nz) else df_reg w_nz = w_r[nz] if not np.all(nz) else w_r - ce_r, _, _, _, _ = self._fit_saturated_regression( + ce_r, _, _, _, _, _ = self._fit_saturated_regression( df_reg_nz, outcome, unit, @@ -1425,6 +1509,13 @@ def _refit_sa_cohort(w_r): es_df_final: Optional[Dict[int, float]] = es_df_used if bootstrap_results is not None: es_df_final = None + # Scalar overall-ATT df provenance mirrors the channel-clearing + # convention: None under bootstrap (percentile p/CI never used a + # df); otherwise the df the overall safe_inference actually + # received (BM contrast df / survey df / knob-resolved fallback). + _overall_df_final: Optional[float] = None + if bootstrap_results is None and _overall_df_used is not None: + _overall_df_final = float(_overall_df_used) if np.isfinite(_overall_df_used) else None # Store results self.results_ = SunAbrahamResults( @@ -1454,6 +1545,8 @@ def _refit_sa_cohort(w_r): cluster_name=(self.cluster if self.vcov_type == "conley" else None), reference_period=self._reference_period, reference_observed=self._reference_observed, + df_convention=self.df_convention, + inference_df=_overall_df_final, ) self.is_fitted_ = True @@ -1480,6 +1573,7 @@ def _fit_saturated_regression( np.ndarray, Dict[Tuple[Any, int], int], Optional[Tuple[np.ndarray, Optional[np.ndarray], np.ndarray]], + _SaturatedFitStats, ]: """ Fit saturated TWFE regression with cohort × relative-time interactions. @@ -1717,6 +1811,10 @@ def _fit_saturated_regression( conley_time=_cl_time, conley_unit=_cl_unit, conley_lag_cutoff=self.conley_lag_cutoff, + # Cells follow the knob: per-cohort-cell t/p/CI resolve through + # get_inference's ladder under the SAME convention the aggregate + # layer uses (the D4 fix — one df source per fit). + df_convention=self.df_convention, ).fit(X, y, df_adjustment=df_adj, cluster_k_adjustment=_cr1_k_adj_sa) vcov = reg.vcov_ @@ -1746,7 +1844,8 @@ def _fit_saturated_regression( # aggregated inference layer to compute per-event-time and # overall-ATT Satterthwaite DOF on user-facing outputs. Under # other vcov_type values aggregated inference falls back to the - # shared analytical df (None → normal distribution). + # knob-resolved analytical df (residual t by default; G−1 under + # df_convention="cluster"; z under "normal"). if vcov_type == "hc2_bm": bread_matrix = X.T @ X bm_artifacts: Optional[Tuple[np.ndarray, Optional[np.ndarray], np.ndarray]] = ( @@ -1757,7 +1856,12 @@ def _fit_saturated_regression( else: bm_artifacts = None - return cohort_effects, cohort_ses, vcov_cohort, coef_index_map, bm_artifacts + fit_stats = _SaturatedFitStats( + df_residual=reg.df_, + n_clusters=reg.n_clusters_, + ) + + return cohort_effects, cohort_ses, vcov_cohort, coef_index_map, bm_artifacts, fit_stats def _compute_iw_effects( self, @@ -1771,7 +1875,7 @@ def _compute_iw_effects( vcov_cohort: np.ndarray, coef_index_map: Dict[Tuple[Any, int], int], survey_weight_col: Optional[str] = None, - survey_df: Optional[int] = None, + fallback_df: Optional[float] = None, ) -> Tuple[Dict[int, Dict[str, Any]], Dict[int, Dict[Any, float]]]: """ Compute interaction-weighted event study effects. @@ -1785,6 +1889,12 @@ def _compute_iw_effects( (sum of weights) rather than raw observation counts, so the estimand reflects the survey-weighted cohort composition. + ``fallback_df`` is the ALREADY-RESOLVED df for the per-row + ``safe_inference`` calls (the caller resolves survey df first, then + the ``df_convention`` knob via ``resolve_tail_df``; None → normal + theory). The bootstrap/replicate refit call sites omit it — their + p/CI are overridden downstream. + Returns ------- event_study_effects : dict @@ -1843,7 +1953,7 @@ def _compute_iw_effects( agg_var = float(weight_vec @ vcov_subset @ weight_vec) agg_se = np.sqrt(max(agg_var, 0)) - t_stat, p_val, ci = safe_inference(agg_effect, agg_se, alpha=self.alpha, df=survey_df) + t_stat, p_val, ci = safe_inference(agg_effect, agg_se, alpha=self.alpha, df=fallback_df) event_study_effects[e] = { "effect": agg_effect, @@ -2087,6 +2197,7 @@ def _run_bootstrap( vcov_b, coef_map_b, _, + _, ) = self._fit_saturated_regression( df_b, outcome, @@ -2320,6 +2431,7 @@ def _run_rao_wu_bootstrap( vcov_b, coef_map_b, _, + _, ) = self._fit_saturated_regression( df_b, outcome, @@ -2438,15 +2550,20 @@ def get_params(self) -> Dict[str, Any]: "conley_metric": self.conley_metric, "conley_kernel": self.conley_kernel, "conley_lag_cutoff": self.conley_lag_cutoff, + "df_convention": self.df_convention, } def set_params(self, **params) -> "SunAbraham": """Set estimator parameters (sklearn-compatible).""" - for key, value in params.items(): - if hasattr(self, key): - setattr(self, key, value) - else: + # Reject unknown keys and validate the pending df_convention BEFORE + # any assignment so a rejected call leaves the estimator unchanged. + for key in params: + if not hasattr(self, key): raise ValueError(f"Unknown parameter: {key}") + if "df_convention" in params: + validate_df_convention(params["df_convention"]) + for key, value in params.items(): + setattr(self, key, value) # Refresh the explicit-vcov-type flag if vcov_type changed, so the # auto-cluster guard at fit time uses the updated value. if "vcov_type" in params: diff --git a/diff_diff/utils.py b/diff_diff/utils.py index 6d1c013a1..2fa236b37 100644 --- a/diff_diff/utils.py +++ b/diff_diff/utils.py @@ -467,6 +467,97 @@ def safe_inference_batch(effects, ses, alpha=0.05, df=None): return t_stats, p_values, ci_lowers, ci_uppers +_DF_CONVENTIONS = ("residual", "cluster", "normal") + + +def validate_df_convention(value: str) -> None: + """Raise ValueError unless ``value`` is a supported df_convention. + + Shared by every estimator constructor and ``set_params`` that carries the + ``df_convention`` knob. ``LinearRegression`` keeps its own inline check + (``diff_diff.utils`` imports ``diff_diff.linalg`` at module level, so + linalg cannot import from here at module scope). + """ + if value not in _DF_CONVENTIONS: + raise ValueError(f"df_convention must be one of {_DF_CONVENTIONS}, got {value!r}") + + +def resolve_tail_df( + df_convention: str, + *, + residual_df: Optional[float], + n_clusters: Optional[int], +) -> Optional[float]: + """Resolve the fallback-level tail df for analytical t/p/CI inference. + + The single implementation of the ``df_convention`` knob's fallback slot + for the standalone estimators (SunAbraham aggregates, WooldridgeDiD OLS, + StackedDiD, ImputationDiD pretrends leads, LPDiD). Survey/replicate df + and per-coefficient Bell-McCaffrey DOF keep precedence at the CALL SITE + under every value — callers consult this helper only when no + higher-precedence df applies. ``LinearRegression.get_inference`` and + MultiPeriodDiD keep their own equivalent inline ladders. + + Branch table: + + - ``"residual"`` — ``residual_df`` (the fit's ``n_eff - K_full``). + - ``"cluster"`` — ``n_clusters - 1`` when ``n_clusters >= 2``. When + ``n_clusters`` is not None but ``<= 1`` the cluster df is undefined: + warn and return ``0``, the fail-closed sentinel ``safe_inference`` + turns into all-NaN inference (the ``df = 0`` twin of MultiPeriodDiD's + inline ``_df_cluster_knob_invalid`` path; ``LinearRegression. + get_inference`` reaches the same observable warn+NaN outcome via an + early return). Defense-in-depth only: every in-scope clustered solve + raises "Need at least 2 clusters" at the vcov layer before df + resolution, and LPDiD's degenerate lanes bypass this helper entirely. + When ``n_clusters`` is None (unclustered/conley fit) the value is + inert and falls back to ``residual_df``. + - ``"normal"`` — ``None`` (deliberate normal-theory z inference at the + fallback level, on clustered and unclustered fits alike). + + Whenever the residual fallback is what would be returned, a non-finite + or non-positive ``residual_df`` warns and returns ``None`` (normal + theory) — mirroring ``get_inference``'s non-positive-df fallback. The + reachable case is absorbed-FE accounting driving ``n - k - rank`` to + zero or below while the SE is still finite (a plain ``n <= k`` + saturation NaNs at the vcov layer first and never reaches this guard). + + Never wraps ``safe_inference`` — callers pass the resolved df to their + own module-level ``safe_inference`` binding so the variance-conventions + audit instrumentation observes it. + """ + validate_df_convention(df_convention) + + def _guarded_residual(df: Optional[float]) -> Optional[float]: + if df is None: + return None + if not (np.isfinite(df) and df > 0): + warnings.warn( + f"Degrees of freedom is non-positive (df={df}). " + "Using normal distribution instead of t-distribution for inference.", + UserWarning, + stacklevel=3, + ) + return None + return float(df) + + if df_convention == "normal": + return None + if df_convention == "cluster": + if n_clusters is None: + return _guarded_residual(residual_df) + if n_clusters <= 1: + warnings.warn( + "df_convention='cluster' requires at least 2 effective " + f"clusters; got {n_clusters}. Inference fields will be NaN.", + UserWarning, + stacklevel=3, + ) + return 0.0 + return float(n_clusters - 1) + return _guarded_residual(residual_df) + + # ============================================================================= # Wild Cluster Bootstrap # ============================================================================= diff --git a/diff_diff/wooldridge.py b/diff_diff/wooldridge.py index c58081740..664c53980 100644 --- a/diff_diff/wooldridge.py +++ b/diff_diff/wooldridge.py @@ -20,12 +20,21 @@ import numpy as np import pandas as pd -from diff_diff.linalg import compute_robust_vcov, solve_logit, solve_ols, solve_poisson +from diff_diff.linalg import ( + compute_robust_vcov, + effective_cluster_count, + solve_logit, + solve_ols, + solve_poisson, +) from diff_diff.utils import ( absorbed_fe_cr1_k_increment, + absorbed_fe_rank, pre_demean_norms, + resolve_tail_df, safe_inference, snap_absorbed_regressors, + validate_df_convention, within_transform, ) from diff_diff.wooldridge_results import WooldridgeDiDResults @@ -945,6 +954,18 @@ class WooldridgeDiD: The Section 8 trend specification is therefore unidentified on this branch. Use ``control_group="not_yet_treated"`` (the default) for the cohort_trends surface. + df_convention : {"residual", "cluster", "normal"}, default "residual" + Degrees-of-freedom convention for the OLS analytical t/p/CI (per-cell + and aggregated). ``"residual"`` (default) uses the fitted residual df + — the 3.9 fix: the default hc1 arms previously used silent + normal-theory z; classical/hc2 keep their historical ``n − rank(X)`` + values bit-for-bit; ``"cluster"`` uses the Stata/fixest cluster df + ``G − 1`` on hc1-clustered fits (inert on one-way and conley + families); ``"normal"`` deliberately uses normal-theory z at the + fallback level. Survey design df and hc2_bm Bell-McCaffrey DOF + always take precedence; the logit/poisson arms are knob-independent + (survey df or normal theory — an explicitly non-default value warns + at fit time). The default flips to ``"cluster"`` at v4. """ def __init__( @@ -966,6 +987,7 @@ def __init__( conley_metric: str = "haversine", conley_kernel: str = "bartlett", conley_lag_cutoff: Optional[int] = None, + df_convention: str = "residual", ) -> None: self._validate_constructor_args( method=method, @@ -974,6 +996,7 @@ def __init__( bootstrap_weights=bootstrap_weights, vcov_type=vcov_type, cohort_trends=cohort_trends, + df_convention=df_convention, ) self.method = method @@ -993,6 +1016,7 @@ def __init__( self.conley_metric = conley_metric self.conley_kernel = conley_kernel self.conley_lag_cutoff = conley_lag_cutoff + self.df_convention = df_convention # Track whether the user explicitly opted out of the "hc1" default. # The auto-cluster-at-unit default in `_fit_ols` is suppressed only # when the user explicitly opts into a one-way family (``hc2``, @@ -1013,6 +1037,7 @@ def _validate_constructor_args( bootstrap_weights: str, vcov_type: str, cohort_trends: bool = False, + df_convention: str = "residual", ) -> None: """Shared validation for both ``__init__`` and ``set_params``. @@ -1038,6 +1063,7 @@ def _validate_constructor_args( f"vcov_type must be one of " f"{{'classical','hc1','hc2','hc2_bm','conley'}}; got '{vcov_type}'" ) + validate_df_convention(df_convention) if method != "ols" and vcov_type != "hc1": raise NotImplementedError( f"WooldridgeDiD(method={method!r}, vcov_type={vcov_type!r}) is " @@ -1086,6 +1112,7 @@ def get_params(self) -> Dict[str, Any]: "conley_metric": self.conley_metric, "conley_kernel": self.conley_kernel, "conley_lag_cutoff": self.conley_lag_cutoff, + "df_convention": self.df_convention, } def set_params(self, **params: Any) -> "WooldridgeDiD": @@ -1115,6 +1142,7 @@ def set_params(self, **params: Any) -> "WooldridgeDiD": "bootstrap_weights": params.get("bootstrap_weights", self.bootstrap_weights), "vcov_type": params.get("vcov_type", self.vcov_type), "cohort_trends": params.get("cohort_trends", self.cohort_trends), + "df_convention": params.get("df_convention", self.df_convention), } self._validate_constructor_args(**pending) @@ -1940,6 +1968,20 @@ def _build(frame: pd.DataFrame, w: Optional[np.ndarray]): else: X_design = X_int + # Inert-config warning (no-silent-failures): the df_convention knob + # applies to the OLS analytical arms only. The logit/poisson QMLE + # paths use the survey design df when present and normal theory + # otherwise, independent of the knob — an explicitly non-default + # value would be silently ignored there. + if self.method != "ols" and self.df_convention != "residual": + warnings.warn( + f"df_convention={self.df_convention!r} has no effect on the " + "logit/poisson arms — their inference uses the survey design " + "df when present and normal theory otherwise.", + UserWarning, + stacklevel=2, + ) + if self.method == "ols": results = self._fit_ols( sample, @@ -2396,23 +2438,47 @@ def _fit_ols( overall_att_bm_dof: Optional[float] = None per_cell_bm_dof: Dict[Tuple, float] = {} bm_artifacts: Optional[Tuple[np.ndarray, np.ndarray, np.ndarray, Dict[Tuple, int]]] = None - # Residual DOF for one-way ``vcov_type in {"classical","hc2"}`` paths - # (full-dummy, no survey). Matches R's ``lm()`` / ``coef_test()`` use - # of ``n - rank(X)`` for the t-distribution under both classical OLS - # SE and ``sandwich::vcovHC(type="HC2")``. ``None`` on hc1 / - # hc2_bm / surveyed paths (those use their own DOF threading or - # df_inf). Mirrors R's t-distribution convention so per-cell + - # aggregate p-values/CIs are not normal-theory under small samples. - df_one_way: Optional[float] = None - if ( - self.vcov_type in ("classical", "hc2") - and use_full_dummy - and resolved is None - and vcov is not None - ): + # Knob-resolved analytical fallback df for ALL OLS arms (3.9 tail-df + # consolidation): ONE variable that every non-BM branch of the 8a/8b + # ladders below consumes. Residual df per design shape — full-dummy + # arms (classical/hc2/hc2_bm/cohort_trends, incl. hc1+cohort_trends) + # use ``n − k_kept``, exactly the historical ``df_one_way`` formula, + # so classical/hc2 stay bit-identical under the default; the + # within-transform arm (hc1/conley) additionally subtracts the + # absorbed [unit, time] rank. Survey df keeps precedence (df_inf); + # otherwise resolve through the ``df_convention`` knob — residual t + # by default (the 3.9 fix: hc1 arms previously dropped to silent + # normal theory), ``G − 1`` under "cluster" (hc1-clustered only; + # conley has no documented G−1 reference and passes + # ``n_clusters=None``), z under "normal". ``resolve_tail_df``'s + # guard returns None on a non-positive residual df, reproducing the + # old ``df_one_way = NaN → df_inf`` degenerate lane. + df_fallback: Optional[float] = None + if resolved is not None: + df_fallback = df_inf # survey design df wins under every knob value + elif vcov is not None: n_kept = int((~np.isnan(coefs)).sum()) - df_candidate = X.shape[0] - n_kept - df_one_way = float(df_candidate) if df_candidate > 0 else float("nan") + if use_full_dummy: + df_res = float(X.shape[0] - n_kept) + else: + df_res = float( + X.shape[0] + - n_kept + - absorbed_fe_rank( + sample, + [unit, time], + has_intercept_col=False, + weights=survey_weights, + ) + ) + _g_eff_w = ( + effective_cluster_count(cluster_ids, survey_weights) + if (self.vcov_type == "hc1" and cluster_ids is not None) + else None + ) + df_fallback = resolve_tail_df( + self.df_convention, residual_df=df_res, n_clusters=_g_eff_w + ) if ( self.vcov_type == "hc2_bm" and use_full_dummy @@ -2505,12 +2571,13 @@ def _fit_ols( bm_artifacts = (X_red, cluster_ids, bread_red, reduced_coef_idx_map) # 8a. Apply DOF threading (or fail-closed NaN) to ``gt_effects``: - # hc2_bm uses per-cell BM Satterthwaite DOF; classical/hc2 (one-way, - # no survey) use the residual ``df_one_way = n - rank(X)`` so - # p-values/CIs match R ``lm()`` / ``coef_test()`` t-distribution - # instead of normal-theory; hc1 / surveyed paths use ``df_inf`` - # (survey df or None). Per ``feedback_bm_contrast_dof_fail_closed``: - # NaN BM DOF emits NaN inference fields (never normal-theory). + # hc2_bm uses per-cell BM Satterthwaite DOF; every other OLS arm + # consumes the ONE knob-resolved ``df_fallback`` (survey df on + # surveyed fits, else the df_convention resolution over the per-shape + # residual df — classical/hc2 reproduce the historical df_one_way + # values bit-for-bit under the default). Per + # ``feedback_bm_contrast_dof_fail_closed``: NaN BM DOF emits NaN + # inference fields (never normal-theory). for (g, t), eff in gt_effects.items(): if self.vcov_type == "hc2_bm" and use_full_dummy and resolved is None: cell_dof = per_cell_bm_dof.get((g, t), float("nan")) @@ -2522,19 +2589,13 @@ def _fit_ols( t_stat = float("nan") p_value = float("nan") conf_int = (float("nan"), float("nan")) - elif ( - self.vcov_type in ("classical", "hc2") - and use_full_dummy - and resolved is None - and df_one_way is not None - and np.isfinite(df_one_way) - ): - t_stat, p_value, conf_int = safe_inference( - eff["att"], eff["se"], alpha=self.alpha, df=df_one_way - ) else: + # Every non-BM arm consumes the ONE knob-resolved fallback: + # survey df on surveyed fits, else the df_convention-resolved + # analytical df (residual t by default; classical/hc2 + # reproduce the historical df_one_way values bit-for-bit). t_stat, p_value, conf_int = safe_inference( - eff["att"], eff["se"], alpha=self.alpha, df=df_inf + eff["att"], eff["se"], alpha=self.alpha, df=df_fallback ) eff["t_stat"] = t_stat eff["p_value"] = p_value @@ -2542,9 +2603,8 @@ def _fit_ols( # 8b. Simple aggregation (always computed). DOF threading mirrors 8a: # hc2_bm uses the overall ATT BM contrast DOF (fail-closed NaN if - # unavailable); classical/hc2 (one-way, no survey) use the residual - # ``df_one_way``; hc1 / surveyed paths use ``df_inf`` (survey df or - # None). + # unavailable); every other arm uses the knob-resolved ``df_fallback`` + # (survey df on surveyed fits, else the df_convention resolution). if self.vcov_type == "hc2_bm" and use_full_dummy and resolved is None: if overall_att_bm_dof is not None and np.isfinite(overall_att_bm_dof): overall = _compute_weighted_agg( @@ -2574,19 +2634,9 @@ def _fit_ols( "p_value": float("nan"), "conf_int": (float("nan"), float("nan")), } - elif ( - self.vcov_type in ("classical", "hc2") - and use_full_dummy - and resolved is None - and df_one_way is not None - and np.isfinite(df_one_way) - ): - overall = _compute_weighted_agg( - gt_effects, gt_weights, gt_keys_ordered, gt_vcov, self.alpha, df=df_one_way - ) else: overall = _compute_weighted_agg( - gt_effects, gt_weights, gt_keys_ordered, gt_vcov, self.alpha, df=df_inf + gt_effects, gt_weights, gt_keys_ordered, gt_vcov, self.alpha, df=df_fallback ) # Metadata @@ -2642,7 +2692,8 @@ def _fit_ols( _df_survey=df_inf, _bm_per_cell_dof=per_cell_bm_dof, _bm_artifacts=bm_artifacts, - _df_one_way=df_one_way, + _df_analytic_fallback=df_fallback, + df_convention=self.df_convention, cohort_trend_coefs=cohort_trend_coefs, cohort_trends=self.cohort_trends, ) @@ -2969,6 +3020,11 @@ def _avg_ax0(a, cell_mask): _gt_vcov=gt_vcov, _gt_keys=gt_keys_ordered, _df_survey=df_inf, + # GLM inference is knob-independent (df_inf = survey df or None); + # store it as the analytic fallback so post-fit aggregate() + # reproduces the fit-time inference on logit/poisson results. + _df_analytic_fallback=df_inf, + df_convention=self.df_convention, ) def _fit_poisson( @@ -3226,4 +3282,9 @@ def _avg_ax0(a, cell_mask): _gt_vcov=gt_vcov, _gt_keys=gt_keys_ordered, _df_survey=df_inf, + # GLM inference is knob-independent (df_inf = survey df or None); + # store it as the analytic fallback so post-fit aggregate() + # reproduces the fit-time inference on logit/poisson results. + _df_analytic_fallback=df_inf, + df_convention=self.df_convention, ) diff --git a/diff_diff/wooldridge_results.py b/diff_diff/wooldridge_results.py index e39507784..29efc3021 100644 --- a/diff_diff/wooldridge_results.py +++ b/diff_diff/wooldridge_results.py @@ -195,17 +195,49 @@ class WooldridgeDiDResults(BaseResults): ``group_time_effects`` to its column index in ``X_red``. Storing reduced artifacts avoids the singular full-design bread that ``_compute_cr2_bm_contrast_dof`` would otherwise reject.""" - _df_one_way: Optional[float] = field(default=None, repr=False) - """Residual DOF (``n - rank(X)``) for one-way ``vcov_type in - {"classical","hc2"}`` paths (full-dummy, no survey). ``aggregate()`` - uses this to thread R's ``lm()`` t-distribution into per-key - inference. ``None`` on hc1 / hc2_bm / surveyed paths (which use BM - DOF or ``_df_survey`` instead).""" + _df_analytic_fallback: Optional[float] = field(default=None, repr=False) + """The RESOLVED analytical fallback df the fit's non-BM ``safe_inference`` + calls actually used (survey design df on surveyed fits; otherwise the + ``df_convention``-resolved value — residual df by default, ``G − 1`` + under "cluster", None under "normal"; ``df_inf`` on GLM fits, whose + inference is knob-independent). ``aggregate()`` reads it so post-fit + re-aggregation reproduces fit-time inference on every OLS arm — it + stays LIVE on bootstrap fits, whose per-cell inference remains + analytical (only the overall p/CI are percentile-overridden). + Renamed from ``_df_one_way`` in 3.9 (pickle migration in + ``__setstate__``); the old field covered classical/hc2 only.""" + df_convention: Optional[str] = None + """The estimator's ``df_convention`` configuration echoed onto the + results ("residual" | "cluster" | "normal"; added 3.9). Governs the OLS + analytical arms only; GLM inference is knob-independent. Appended LAST + (the generated ``__init__`` positional indexes are public API).""" # ------------------------------------------------------------------ # # Public methods # # ------------------------------------------------------------------ # + def __setstate__(self, state: Dict[str, Any]) -> None: + """Restore pickled state, migrating renamed/added stored fields. + + Unpickling bypasses ``__init__``/``__post_init__`` (the + ``BaseResults`` pickle-migration contract; SyntheticDiDResults + precedent). A pickle written before 3.9 stores ``_df_one_way`` + (classical/hc2-only residual df) instead of + ``_df_analytic_fallback`` — carry the old value over so post-fit + ``aggregate()`` on a legacy result reproduces its fit-time + classical/hc2 inference (legacy hc1 pickles carried None there, + which restores the pre-3.9 normal-theory aggregation for them). + ``df_convention`` (added 3.9) defaults to ``"residual"``. + """ + self.__dict__.update(state) + # NOTE: check __dict__ membership, not hasattr - dataclass field + # defaults live as CLASS attributes, so hasattr is always True. + if "_df_analytic_fallback" not in self.__dict__: + self.__dict__["_df_analytic_fallback"] = state.get("_df_one_way", None) + self.__dict__.pop("_df_one_way", None) + if "df_convention" not in self.__dict__: + self.__dict__["df_convention"] = "residual" + def aggregate(self, type: str, weights: str = "cell") -> "WooldridgeDiDResults": # noqa: A002 """Compute and store one of the four jwdid_estat aggregation types. @@ -358,11 +390,11 @@ def _build_effect( """Build an effect dict using ``df_for_inference`` for the t-distribution. When ``self.vcov_type == "hc2_bm"``, ``df_for_inference`` should be - the BM contrast DOF (NaN → fail-closed). For ``classical`` / - ``hc2`` (one-way, no survey) the residual DOF ``self._df_one_way`` - is used so per-key inference matches R ``lm()`` / - ``coef_test()`` t-distribution. For hc1 / surveyed paths, - ``self._df_survey`` (None → normal-theory) is used. + the BM contrast DOF (NaN → fail-closed). Every other arm uses + ``self._df_analytic_fallback`` — the resolved fallback the fit's + own ``safe_inference`` calls used (survey df on surveyed fits; + the ``df_convention``-resolved analytical df otherwise) — so + post-fit aggregation reproduces fit-time inference. Under ``weights="cohort_share"`` (variable ``cohort_share_inference_fail_closed=True``), the inference @@ -392,17 +424,9 @@ def _build_effect( t_stat, p_value, conf_int = safe_inference( att, se, alpha=self.alpha, df=df_for_inference ) - elif ( - self.vcov_type in ("classical", "hc2") - and self._df_one_way is not None - and np.isfinite(self._df_one_way) - ): - t_stat, p_value, conf_int = safe_inference( - att, se, alpha=self.alpha, df=self._df_one_way - ) else: t_stat, p_value, conf_int = safe_inference( - att, se, alpha=self.alpha, df=self._df_survey + att, se, alpha=self.alpha, df=self._df_analytic_fallback ) return { "att": att, @@ -786,6 +810,8 @@ def to_dict(self) -> Dict[str, Any]: result["n_clusters"] = self.n_clusters if self.conley_lag_cutoff is not None: result["conley_lag_cutoff"] = self.conley_lag_cutoff + if self.df_convention is not None: + result["df_convention"] = self.df_convention return result def to_dataframe(self, aggregation: str = "event") -> pd.DataFrame: diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 6165f3100..85a363d3f 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -460,9 +460,10 @@ This matches the behavior of R's `fixest::feols()` with absorbed FE. `benchmarks/data/reghdfe_kref_golden.json`; `tests/test_variance_conventions.py::TestFixestCr1NonNestedParity` / `TestReghdfeKReferenceParity`), while the crossed-cluster (nothing-nested) arm matches both - references at machine precision. This is distinct from the SunAbraham / Wooldridge `hc1` - tail-df items below (PR C scope). `hc2`/`hc2_bm` use leverage / - Satterthwaite DOF and are unaffected. + references at machine precision. This is distinct from the TAIL-df convention, which + converged separately in 3.9 (M-127: the three-value `df_convention` knob below — the + SE-scale `k` here and the inference df are independent accountings). `hc2`/`hc2_bm` + use leverage / Satterthwaite DOF and are unaffected. - **Note (deviation from R — clustered CR1 inference df):** under unit clustering, Python's CR1 t-statistics/p-values/CIs use the **residual df** from the fitted design (`n − K_full`, e.g. 148 on the live-R benchmark panel: n=200, 50 unit FE + time + treatment), @@ -474,16 +475,26 @@ This matches the behavior of R's `fixest::feols()` with absorbed FE. pins; the convention itself is locked at a distinguishable moderate |t|~1.8 by `tests/test_methodology_twfe.py::test_moderate_t_pins_residual_df_convention`, where the t(148) and t(49) tails differ by ~5%). The SE itself matches fixest under the K_reference - convergence above (the historical ~0.25% non-nested-FE band is closed). **Opt-in knob (2026-07):** - `df_convention="cluster"` on `DifferenceInDifferences` / `TwoWayFixedEffects` / `MultiPeriodDiD` - (and the `LinearRegression` linalg surface) switches clustered analytical t/p/CI to the - Stata/fixest `G − 1` convention — fallback-level only (survey df and per-coefficient - Bell-McCaffrey DOF keep precedence), point estimates/SEs/t-statistics unchanged, inert on - unclustered fits. The default remains `"residual"` and flips at v4 (lifecycle: - `docs/v4-deprecations.yaml` M-004..M-006; work row in `DEFERRED.md` → Version-gated (v4); - flipping moves every clustered p-value/CI, so it is a major-version change). Standalone - estimators (CS, SA, imputation-family, etc.) carry their own inference stacks and are out of the - knob's scope. Locked by + convergence above (the historical ~0.25% non-nested-FE band is closed). **Three-value knob + (2026-07, widened 3.9 / M-127):** `df_convention ∈ {"residual", "cluster", "normal"}` on the + library-wide knob roster — `DifferenceInDifferences` / `TwoWayFixedEffects` / + `MultiPeriodDiD` / `LinearRegression` plus (since 3.9) `SunAbraham`, `WooldridgeDiD` (OLS + arms), `StackedDiD`, `ImputationDiD` (pretrends lead regression), and `LPDiD` (default + `"cluster"`, its Stata-anchored convention). `"cluster"` switches clustered analytical + t/p/CI to the Stata/fixest `G − 1` convention and is inert on unclustered and conley fits; + `"normal"` is deliberate normal-theory z at the fallback level on EVERY fit (clustered, + unclustered, and conley alike — on lanes that were already t it is a new deliberate option, + not a compatibility path). Fallback-level only under every value: survey/replicate df and + per-coefficient Bell-McCaffrey DOF keep precedence, and point estimates/SEs/t-statistics + never move. The default remains `"residual"` on every surface except LPDiD and flips at v4 + (lifecycle: `docs/v4-deprecations.yaml` M-004..M-006 + M-128..M-131; work row in + `DEFERRED.md` → Version-gated (v4); flipping moves every clustered p-value/CI, so it is a + major-version change). **Out of the knob's scope** (own inference stacks): the IF-based + family — `CallawaySantAnna`, `TwoStageDiD`, `ChaisemartinDHaultfoeuille`, `EfficientDiD`, + `ContinuousDiD`, `ChangesInChanges`, `SyntheticDiD`, `HeterogeneousAdoptionDiD` — and the + OLS-based `TripleDifference` / `TROP` (own `safe_inference` df conventions, outside the + variance-conventions audit matrix's shared-CR1 scope), plus `WooldridgeDiD`'s + logit/poisson arms (survey df or normal theory, knob-independent). Locked by `tests/test_estimators_vcov_type.py::TestDfConvention` (G−1 tail match, precedence, no-op default). The full-dummy (`fixed_effects=`) idiom carries `df_adjustment == 0` and is unchanged — its residual t-df already matched fixest's full-K count because every FE column sits in `k_visible`. **That claim is scoped to @@ -1598,21 +1609,25 @@ where weights ŵ_{g,e} = n_{g,e} / Σ_g n_{g,e} (sample share of cohort g at eve separately get per-coefficient BM DOF via `LinearRegression.get_inference()` inside `_fit_saturated_regression`. If the linalg helper fails (rank-deficient design, singular bread), - the aggregated inference falls back to the shared analytical df with - an explicit `UserWarning`. -- **Note:** Provenance exposure (4.0 program row M-092, no numeric change): - `SunAbrahamResults.event_study_df` records, per relative time, the df each - stored event-study row's `safe_inference` actually received - the per-event - BM contrast DOF above under `hc2_bm`, the survey design df (post-drop under - a replicate refit) on survey fits, and NaN on plain analytic fits (normal - theory) as well as on rows whose BM DOF was non-finite, where + the aggregated inference falls back to the `df_convention`-resolved + analytical df (t(residual) under the default) with an explicit + `UserWarning`. +- **Note:** Provenance exposure (4.0 program row M-092; df values updated by + the 3.9 M-127 tail-df fix): `SunAbrahamResults.event_study_df` records, + per relative time, the df each stored event-study row's `safe_inference` + actually received - the per-event BM contrast DOF above under `hc2_bm`, + the survey design df (post-drop under a replicate refit) on survey fits, + the `df_convention`-resolved fallback on plain analytic fits (FINITE + residual df under the 3.9 default; `G − 1` under `"cluster"`; NaN under + `"normal"`), and NaN on rows whose BM DOF was non-finite, where `safe_inference`'s own non-finite-df guard already yields all-NaN inference. Note the two failure modes differ and are both faithfully - recorded: a helper-level failure warns and falls back to the shared + recorded: a helper-level failure warns and falls back to the resolved analytical df (documented above), while a per-row non-finite DOF fails closed. The channel clears under bootstrap (percentile p/CIs used no df) - deliberately narrower than the `event_study_vcov` clear, which also fires - under replicate refits whose rows DID use a genuine df. + under replicate refits whose rows DID use a genuine df. The scalar + `inference_df` (added 3.9) mirrors the same rule for the overall ATT. - **HC1 finite-sample correction (K_reference convergence, 3.9):** SA's clustered CR1 factor counts the absorbed unit + time FE via the `cluster_k_adjustment` increment (`K_reference = k_dm + 1 + rank(non-nested @@ -1661,23 +1676,29 @@ where weights ŵ_{g,e} = n_{g,e} / Σ_g n_{g,e} (sample share of cohort g at eve - Bootstrap inference: p_value and CI computed from bootstrap distribution. SE, CI, and p-value are all NaN if the original point estimate is non-finite, SE is non-finite or zero, or if <50% of bootstrap samples are valid - Applies to overall ATT, per-effect event study, and aggregated event study - **Note**: Defensive enhancement matching CallawaySantAnna behavior; R's `fixest::sunab()` may produce Inf/NaN without warning -- Inference distribution: - - Cohort-level p-values: t-distribution (via `LinearRegression.get_inference()`) +- Inference distribution (ONE df source per fit since 3.9 / M-127 — the D4 + fix; previously cells used t(residual) while aggregates dropped to normal + theory inside the same fit): + - Cohort-level p-values: t-distribution via + `LinearRegression.get_inference()`, following the `df_convention` knob + (residual df by default; `G − 1` under `"cluster"`; z under `"normal"`). - Aggregated event study and overall ATT p-values: - Under `vcov_type="hc2_bm"`: t-distribution with CR2 Bell-McCaffrey contrast DOF per aggregated effect (see "Phase 1b aggregated BM contrast DOF" Note above). Matches `clubSandwich::Wald_test( test="HTZ")$df_denom`. - - Under `vcov_type ∈ {"classical","hc1","hc2"}` (no replicate-weight - survey): normal distribution (via `compute_p_value()`), which is - asymptotically equivalent and standard for delta-method-aggregated - quantities. + - Under `vcov_type ∈ {"classical","hc1","hc2","conley"}` (no survey): + the SAME `df_convention`-resolved df as the cells — t(residual df of + the saturated regression) by default, t(G − 1) under `"cluster"` + (inert under conley), z under `"normal"`. - Under replicate-weight survey: t-distribution with replicate-derived DOF (`survey_metadata.df_survey`). - - **Deviation from R**: R's fixest uses t-distribution at all levels - under `vcov_type ∈ {"classical","hc1","hc2"}`; aggregated p-values - may differ slightly for small samples on those families. The - `hc2_bm` aggregated path matches clubSandwich exactly. + - **Deviation from R**: R's fixest uses t-distribution at all levels with + the CLUSTER df `G − 1` on clustered fits; the 3.9 default uses the + residual df instead (the library-wide `df_convention` deviation + documented in the TwoWayFixedEffects section — `df_convention= + "cluster"` opts into the fixest convention today, and the default flips + at v4). The `hc2_bm` aggregated path matches clubSandwich exactly. **Reference implementation(s):** - R: `fixest::sunab()` (Laurent Bergé's implementation) @@ -1779,6 +1800,22 @@ where `W_it(h) = 1[K_it = h]` are lead indicators, estimated on `Omega_0` only. - When `balance_e` is set, lead indicators are restricted to balanced cohorts; the full Omega_0 sample (including never-treated) is kept for within-transformation - Only affects event study aggregation; overall ATT and group aggregation unchanged - **Note:** `pretrends=True` with analytical `survey_design` (strata/PSU/FPC) is supported. The lead regression uses survey-weighted demeaning, WLS point estimates, and `compute_survey_vcov()` for design-based VCV. The full survey design is preserved (subpopulation approach): Omega_0 scores are zero-padded back to full-panel length so PSU/strata structure is maintained for variance estimation. The F-test in `pretrend_test()` uses the full-design `df_survey` as denominator df. Replicate-weight survey designs raise `NotImplementedError` with `pretrends=True` because per-replicate lead regression refits are not yet implemented. +- **Note (pretrends lead-regression tail df, 3.9 / M-127):** the per-lead + `gamma_h` t/p/CI previously used silent normal-theory z on plain clustered + fits (the tail-df defect family); they now follow the `df_convention` + knob — t(residual df `n − k_kept − rank(absorbed [time, unit])` on the + Omega_0 lead design) under the default `"residual"`, t(G − 1) under + `"cluster"`, z under `"normal"` — with the full-design survey df keeping + precedence. This is the ONE ImputationDiD surface on the shared clustered + CR1 sandwich; the BJS Theorem-3 overall/post-treatment inference is + knob-independent (L3), and the JOINT pretrend Wald F deliberately keeps + its cluster-robust `F(q, max(G − 1, 1))` denominator (the standard + joint-test convention; it consumes gamma/V_gamma, not the per-lead t's, + coheres with `df_convention="cluster"`, and converges with the t's at the + v4 flip) — the per-lead-t vs joint-F pairing is a stated convention, not + a silent split. A `df_convention` value that never surfaces the per-lead + inference (`pretrends=False`, or `aggregate` outside + `{"event_study","all"}`) warns at fit time when explicitly non-default. *Edge cases:* - **Unbalanced panels:** FE estimated via iterative alternating projection (Gauss-Seidel), equivalent to OLS with unit+time dummies. Converges in O(max_iter) passes; typically 5-20 iterations for unbalanced panels, 1-2 for balanced. One-pass demeaning is only exact for balanced panels. @@ -1962,6 +1999,22 @@ Estimated via WLS with Q-weights. The delta_h coefficients identify theta_kappa^ **Note:** This routing inherits the WLS-CR2 methodology choice from the Phase 1a clubSandwich port (PR #475 / REGISTRY Phase 1a `hc2_bm + weights` row); see that row for the full PT2018-§3.3-vs-clubSandwich algebra deviation derivation. No new methodology choice is introduced in this PR. The change is purely surface: switching from the prior bake-Q-into-X pattern (`X_t = X * sqrt(Q)`, `solve_ols(X_t, Y_t, cluster_ids=)`) to `solve_ols(weights=composed_weights, vcov_type=...)` opens the hc2_bm path without modifying the small-sample WLS-CR2 algebra. The HC1 path is preserved bit-equal (up to float64 multiplication ordering at machine precision). +**Note (tail df, 3.9 / M-127):** the `hc1` analytical t/p/CI previously used +silent normal-theory z on the non-survey lane (the tail-df defect family); +they now follow the library-wide `df_convention` knob — t(pooled residual df +`n_eff − k_kept`, positive-weight rows) under the default `"residual"`, +t(G − 1) under `"cluster"`, z under `"normal"` — with hc2_bm BM contrast DOF +and survey/replicate df keeping precedence. Two StackedDiD-specific caveats: +(a) the `"cluster"` G counts POSITIVE-WEIGHT clusters +(`effective_cluster_count`), while the reported `results.n_clusters` keeps +the raw unique count — the two diverge only when a cluster's total composed +weight is zero; (b) the stacked design replicates control units across +sub-experiments, so the pooled residual df is typically large and the +default's z → t movement is a convention alignment with near-zero numeric +effect — the materially moving step on this surface is the v4 flip to +`G − 1` (M-130). No public per-row df channel changes: `event_study_df` +records the resolved df per event time as before. + *Edge cases:* - All events trimmed: `len(Omega_kappa) == 0` -> ValueError suggesting reduced kappa - No clean controls for event a: IC2 check fails -> Trim event, warn user @@ -2115,6 +2168,24 @@ where `g(·)` is the link inverse (logistic or exp), `η_i` is the individual li - `n_bootstrap > 0` + `vcov_type ∈ {"hc2","classical"}` — REJECTED at `fit()` regardless of `self.cluster` setting. The multiplier bootstrap is intrinsically clustered, but one-way vcov_type does not compose with `cluster_ids`: with `cluster=None` the auto-cluster is dropped (bootstrap has no cluster to draw weights at); with `cluster=X` the linalg validator rejects one-way + cluster_ids downstream with a less-informative error. User must drop bootstrap (`n_bootstrap=0`) or pick a cluster-compatible `vcov_type` (`hc1` or `hc2_bm`). - **Note:** This routing is a documented synthesis of two existing methodology ingredients: the full-dummy auto-route from the Phase 1b PR 1/8 SunAbraham pattern (PR #472, which itself reused the Phase 1a Gate 1 TWFE lift from PR #469), and the clubSandwich WLS-CR2 algebra from the Phase 1a port (PR #475). The BM contrast DOF threading reuses `_compute_cr2_bm_contrast_dof` from PR #465 (MPD). No new methodology choice is introduced — the change is purely surface: extending the existing pattern from SA-OLS to WooldridgeDiD-OLS. - **Note:** Bootstrap is supported only with `vcov_type ∈ {"hc1","hc2_bm"}` (one-way `classical`/`hc2` + bootstrap is rejected at `fit()` per the previous bullet). On the supported paths, the bootstrap clusters at `self.cluster if self.cluster else unit` — i.e., it matches the user's explicit cluster column if set, falling back to unit otherwise (the panel's natural unit of variation). The bootstrap SE overrides the analytical SE for `overall_*` on `n_bootstrap > 0` paths; per-cell `(g, t)` SEs still come from the analytical vcov. +- **Note (tail df, 3.9 / M-127):** every non-BM OLS arm's per-cell and + aggregated t/p/CI now consume ONE `df_convention`-resolved fallback df — + the default-`hc1` arms (within-transform AND `hc1`+`cohort_trends` + full-dummy) previously used silent normal-theory z (the tail-df defect + family), and now use t(residual df): `n − k_kept − rank(absorbed + [unit, time])` on the within arm, the historical `n − rank(X)` on the + full-dummy arms (`classical`/`hc2` values bit-identical). `"cluster"` + gives t(G − 1) on hc1-clustered fits (inert on one-way/conley families); + `"normal"` gives z. Survey design df and hc2_bm BM DOF keep precedence; + the logit/poisson arms are knob-INDEPENDENT (survey df when present, + normal theory otherwise — an explicitly non-default `df_convention` + warns at fit time). Post-fit `aggregate()` reproduces fit-time inference + on every OLS arm via the stored resolved fallback + (`_df_analytic_fallback`, live on bootstrap fits whose per-cell inference + stays analytical). No public per-row df channel ships in 3.9 — the + unified event-study surface's Wooldridge df column stays NaN (a per-key + dict that could also carry the hc2_bm BM contrast dofs is the TODO.md + df-provenance follow-up). *Aggregations (matching `jwdid_estat`):* - `simple`: Weighted average across all post-treatment (g, t) cells. Default @@ -2275,7 +2346,7 @@ Eq. 12 reuses the absorbing clean control and only restricts the treated set (a The paper specifies no standard-error formula (Section 1 defers to "standard, well-understood techniques"); the reference Stata `lpdid` uses `vce(cluster unit)`. The entries below document diff-diff's inference and scope choices. -1. **Note:** Standard errors are **cluster-robust at the unit level by default** - `cluster=None` auto-clusters at the unit identifier and the results record `cluster_name`/`n_clusters` - with a `t(G-1)` reference distribution (G = realized clusters in each horizon's clean-control sample). Matches Stata `lpdid` `vce(cluster unit)`; the paper prescribes no SE. +1. **Note:** Standard errors are **cluster-robust at the unit level by default** - `cluster=None` auto-clusters at the unit identifier and the results record `cluster_name`/`n_clusters` - with a `t(G-1)` reference distribution (G = realized clusters in each horizon's clean-control sample). Matches Stata `lpdid` `vce(cluster unit)`; the paper prescribes no SE. **Since 3.9 (M-127) the reference is the `df_convention` knob's `"cluster"` value — LPDiD's DEFAULT, so nothing moved**: `"residual"` opts into t(per-horizon residual df `n_eff − k_kept`; on the RA path `n_total − k0_kept − 1`, the pooled M-estimator's parameter count — nuisance coefficients plus the ATT — a library convention with no external anchor since Stata `teffects ra` reports z), and `"normal"` into z. The G here stays the RAW unique cluster count (deliberately not `effective_cluster_count`: LPDiD's reweights are strictly positive and there is no user weight column, so the two cannot diverge). **Carve-out to the flat `t(G-1)` claim:** the degenerate lanes — the unclustered-refit fallback when the clustered solve fails, the RA `G<=1` lane, and the saturated-design early return (`n_obs <= k`) — keep the literal `df=None` (normal theory, NaN provenance) under ALL knob values, exactly as before 3.9. 2. **Note:** The regression-adjustment (RA) covariate path (`reweight=True` with covariates/absorb) reports an **influence-function cluster variance** `sum_c (sum_{i in c} psi_i)^2 / n^2`, in the same family as `ImputationDiD`'s Theorem-3 / BJS variance (see "IF-based variance estimators vs analytical-sandwich estimators" above). Its single Gram inversion is routed through `linalg._rank_guarded_inv` (finite SE on the identified subspace under near-collinearity; NaN at rank 0). Unlike the default/weighted `solve_ols` `hc1`-cluster path - which applies the `(G/(G-1))*((n-1)/(n-K_reference))` finite-sample factor (see the nested-dummy Note below) - the RA IF variance carries **no finite-sample factor**, while both paths share the `t(G-1)` reference. **PR-B2 validated this asymmetry as faithful to the authors' own tooling**, not a defect: the no-factor RA convention matches the canonical Stata `teffects ra ... atet vce(cluster)`, while the default path matches `feols`/`reghdfe`. **Originally inferred from the authors' `lpdid_regression_adjustment.do` degrees-of-freedom comments, this is now MEASURED against a runnable `teffects`:** `benchmarks/stata/generate_lpdid_ra_golden.do` runs `teffects ra (Dy x i.time) (tdiff), atet vce(cluster unit)` on an independently reconstructed clean sample (a Stata port of the `generate_lpdid_golden.R` `prep`/`clean_h` recipe), and the library RA IF SE matches it to **~1e-16 at all 7 event-study horizons** - the SE ratio is a flat 1.0 (no `G/(G-1)`, no `(n-1)/(n-k)`), directly confirming the no-finite-sample-factor convention (golden `benchmarks/data/lpdid_ra_stata_golden.json`; test `tests/test_lpdid_ra_stata_parity.py`, which also gates the Stata clean sample's realized size and cluster count `(e(N), e(N_clust))` against the library's `(n_obs, n_clusters)`). The RA *point* estimate is triple-anchored: Python == R full-interaction (`i.dtreat##(i.time c.x)`) == Stata `teffects`, to ~1e-10..1e-13 (`tests/test_methodology_lpdid.py::test_ra_covariate_point`). Because Stata is node-locked and cannot run in CI (goldens are committed, exactly like the R arm), the CI-runnable guards remain the regression pin (`test_ra_covariate_se_regression_pin` / `RA_SE_PIN`) and the ungated Monte-Carlo coverage study `benchmarks/python/coverage_lpdid_ra.py` (~0.95 empirical coverage of the true effect at cluster counts G in {30, 100, 300}). 3. **Note:** Direct covariate inclusion (`reweight=False` with covariates/absorb) emits a `UserWarning`: per online Appendix B.2.2 it preserves the non-negative LP-DiD weighting result only under linear and homogeneous covariate effects, so the regression-adjustment path (`reweight=True`) is preferred. 4. **Deviation from R:** Scope - non-absorbing treatment (Section 4.2) implements the **entry-effect** estimands (`non_absorbing="first_entry"` / `"effect_stabilization"`, PR-C1). **PR-C2 R-parity-validated both modes against an INDEPENDENT `fixest::feols` reconstruction of the paper's Eq. 12 / Eq. 13 clean-sample restrictions** (point and SE match to ~1e-13/~1e-15 for the variance-weighted variants; the `effect_stabilization` reweighted point matches and its SE is pinned as a regression guard - a small weighted-cluster convention difference vs feols; `tests/test_methodology_lpdid.py::TestLPDiDNonAbsorbingParityR`). The recipe's independence was demonstrated when an earlier draft's Eq. 12 control off-by-one diverged from the already-correct library and was corrected against the paper, plus a hand-computed Python micro-check. **`alexCardazzi/lpdid`'s `nonabsorbing_lag` is NOT a faithful Eq. 13** (it clamps `treat_diff[<0]<-0`, so its clean-control window blocks only treatment turn-*ons*; it reuses a forward placebo window; and it NA-excludes pre-panel-treated rows where the library clamps pre-`min_t` to untreated): it diverges ~0.01-0.05 from Eq. 13 even on a monotone no-off-switch panel, so it is **recorded in the golden `meta` as a divergent third-party reference, not a parity gate** (the alexCardazzi-pooled precedent). The library's "no treatment change" (both directions) and backward placebo window are the more paper-faithful choices. `first_entry` (Eq. 12) has no R-package analogue (anchored on the independent feols recipe only). Appendix-C exit-event dynamics and the Stata canonical SE remain deferred follow-ups. @@ -2283,7 +2354,7 @@ The paper specifies no standard-error formula (Section 1 defers to "standard, we 6. **Note (pooled estimand):** The pooled pre/post ATT (the headline `results.att` is the pooled-post row) is the **unit-equal-weighted average of each unit-event-time's mean long difference** over the window - `mean_h(y_{i,t+h}) - baseline_{i,t}`, one observation per (unit, event-time), regressed on the treatment-switch indicator with event-time fixed effects on the **fixed-composition** sample (only units observing *every* pooled target, with clean controls required through `max(h)`). This equals the mean of the per-horizon event-study coefficients on a balanced panel. **PR-B2 validated it against the authors' runnable R reference**: the pooled estimand matches the authors' own R pooled recipe (`danielegirardi/lpdid`: a `slider` window-mean minus `y_{t-1}` on the clean-through-window-end sample) to ~1e-13 (`tests/test_methodology_lpdid.py::test_pooled`). A prior version of this note speculated the authors used a horizon-**stacked** pooled regression; the authors' R reference in fact uses this same fixed-composition mean-long-difference, so that speculation was incorrect. Unlike the event-study variants (where `alexCardazzi` is a cross-check gate), pooled is anchored to the authors' recipe **only**: `alexCardazzi`'s pooled uses a **laxer** clean-control window, so it differs and is recorded in the golden `meta` for transparency, not as a parity target. 7. **Deviation from R:** `no_composition` is intentionally more faithful to the paper's fixed-composition intent (Section 3.6) than the R packages: it fixes the realized sample across *all* post horizons (every post coefficient shares one sample, even on unbalanced panels) and excludes cohorts with `p_g > T-H`, whereas `alexCardazzi/lpdid` uses a looser per-horizon sample and a stricter `treat_date < T-H` cutoff. It therefore has **no exact R-package anchor** and is validated by the pure-Python tests in `tests/test_lpdid.py` (the R-parity golden omits it; `alexCardazzi`'s looser-semantics value is recorded in the golden `meta`). 8. **Note (survey design):** Complex-survey support (`survey_design=SurveyDesign(...)`, PR-D1) covers the **variance-weighted default path** (`reweight=False`, with or without direct-inclusion covariates): each horizon's long-difference regression is fit by WLS on the survey probability weights, and the SE is the stratified-PSU **Taylor-linearization (Binder 1983 TSL)** sandwich `meat = sum_h (1-f_h)*(n_h/(n_h-1))*sum_j (S_hj - S_h_bar)(S_hj - S_h_bar)'` with `df = n_PSU - n_strata`, reusing the shared `diff_diff/survey.py` helpers (`compute_survey_vcov` / `_compute_stratified_psu_meat`). The design is re-resolved on each realized (post-clean-control) sample so weights/strata/PSU align with the regression rows; with no explicit PSU the unit (LP-DiD's default cluster) is injected as the PSU. Supports pweight + strata + PSU + FPC + lonely-PSU handling. It **rejects** `survey_design` combined with `reweight=True` (the equally-weighted / regression-adjustment IF path has no validated survey reference - the same gap as the RA SE in Deviation #2), replicate-weight designs, and non-pweight (fweight/aweight) types, each a deferred follow-up. The non-survey path is byte-for-byte unchanged (gated on `survey_design is None`). **PR-D2 validated all three survey paths end-to-end against `survey::svyglm`** - per-horizon point/SE/df + pooled for the variance-weighted full design (strata+PSU+FPC), the weights-only unit-injected-PSU design, and the direct-covariate variant (`tests/test_methodology_lpdid.py::TestLPDiDSurveyParityR`; point ~1e-6, SE ~1e-5, df exact via the per-design `n_PSU - n_strata` / `n_PSU - 1` formula). `svyglm` is itself the reference implementation of the Binder TSL sandwich, so it anchors the variance directly (no third-party survey-package gate is needed); the clean-sample construction is independently cross-checked in the generator (the unweighted variance-weighted event study matches `alexCardazzi/lpdid` to <1e-8, and selection is weight-independent). A dedicated survey panel (`benchmarks/data/lpdid_survey_panel.csv`, own seed) keeps the absorbing / non-absorbing goldens byte-identical. -9. **Note (clustered CR1 k, nested dummy blocks — 3.9 K_reference):** LPDiD's per-horizon design carries two built-in dummy blocks — the `_event_time` (time-FE) dummies added under `include_time_fe=True` and the inline `drop_first` `absorb=` dummies. The authors' reference recipe ABSORBS the time FE (`feols(... | time, vcov = ~unit)`, `benchmarks/R/generate_lpdid_golden.R`), so both blocks follow the nested convention: when a block is nested in the raw `cluster` (e.g. a time-like cluster for `_event_time`, or an absorb dim coarser than the cluster), its rank is SUBTRACTED from the CR1 factor's `k` via a negative `cluster_k_adjustment`, converging on the same `K_reference` accounting as `feols`/`reghdfe`. Non-nested configs are no-ops (unit clustering leaves `_event_time` non-nested → adjustment 0, preserving bit-identical Rust-lane output). The solve's broad `except (ValueError, ZeroDivisionError)` re-solve re-raises `InvalidClusterKAdjustment` so a contract violation can never degrade to a silent `se=NaN` (`tests/test_variance_conventions.py::TestKReferenceConvergence::test_lpdid_never_swallows_the_contract_raise`). The `t(G-1)` tail reference (Note 1) is unchanged. +9. **Note (clustered CR1 k, nested dummy blocks — 3.9 K_reference):** LPDiD's per-horizon design carries two built-in dummy blocks — the `_event_time` (time-FE) dummies added under `include_time_fe=True` and the inline `drop_first` `absorb=` dummies. The authors' reference recipe ABSORBS the time FE (`feols(... | time, vcov = ~unit)`, `benchmarks/R/generate_lpdid_golden.R`), so both blocks follow the nested convention: when a block is nested in the raw `cluster` (e.g. a time-like cluster for `_event_time`, or an absorb dim coarser than the cluster), its rank is SUBTRACTED from the CR1 factor's `k` via a negative `cluster_k_adjustment`, converging on the same `K_reference` accounting as `feols`/`reghdfe`. Non-nested configs are no-ops (unit clustering leaves `_event_time` non-nested → adjustment 0, preserving bit-identical Rust-lane output). The solve's broad `except (ValueError, ZeroDivisionError)` re-solve re-raises `InvalidClusterKAdjustment` so a contract violation can never degrade to a silent `se=NaN` (`tests/test_variance_conventions.py::TestKReferenceConvergence::test_lpdid_never_swallows_the_contract_raise`). The `t(G-1)` tail reference (Note 1) is unchanged and is, since 3.9, the `df_convention="cluster"` DEFAULT (see Note 1's knob addendum). ### Implementation Checklist diff --git a/docs/methodology/variance-conventions.md b/docs/methodology/variance-conventions.md index 906fd4dbc..0515aac3d 100644 --- a/docs/methodology/variance-conventions.md +++ b/docs/methodology/variance-conventions.md @@ -26,16 +26,16 @@ produced wrong figures three separate times while this inventory was drafted). | `did_fixed_effects_hc1_cluster_unit` | 7 | 294 | **legitimate** | K_reference (D1 fixed): 66 visible minus the 59 cluster-nested unit dummies -> identical SE to did_absorb (documented deviation from a literal explicit-dummy R comparison, which counts all 66) | | `did_plain_hc1_cluster_unit` | 4 | 356 | **legitimate** | no absorbed FE: visible k is the whole design; nothing is omitted | | `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 | None, None, None, None, None, None, None, None, None, None | **defect** | CR1 k converged on K_reference (D2 fixed: 9 cells + T = 15, no intercept col -> +1 term; jwdid arms at ratio 1.0); tail df is still normal theory with no df_convention knob (PR C) | -| `sun_abraham_hc1` | 21 | 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, None, None, None, None, None, None, None, None | **defect** | CR1 k converged on K_reference (D2 fixed: 15 cells + 6, no intercept col; fixest sunab parity ~5e-15); D4 remains: residual df per cohort-period cell but normal theory on aggregates (PR C) | -| `stacked_did_hc1` | 6 | None | **legitimate** | L1: k_total is clubSandwich CR1S by construction (stacked_did.py pins vcovCR(type='CR1S') at atol=1e-10); normal-theory tail df is an open PR C question | -| `lpdid_pre2_post2` | 4, 4, 5, 5, 5, 6 | 59, 59, 59, 59, 59, 59 | **legitimate** | L2: G-1 tail df (Stata/fixest convention) — the convergence target | +| `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) | +| `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 | | `mpd_plain_hc1_cluster_time` | 7 | 348, 348, 348, 348, 348, 348 | **legitimate** | the NESTED orientation of the built-in period dummies: 12 visible minus their rank 5 under a time cluster (under-subtraction is caught here; the unit-cluster rows catch over-) | | `lpdid_absorb_nested_cluster_grp` | 4, 4, 5, 5, 5, 6 | 1, 1, 1, 1, 1, 1 | **legitimate** | LPDiD absorb dummies nested in the cluster subtract their rank (adj -1 per horizon: region == grp here); _event_time stays counted (unit-level cluster does not nest time); G-1 tail df (L2) | | `imputation_default` | — | None | **legitimate** | L3: BJS imputation variance, not the shared CR1 sandwich | -| `imputation_pretrends_event_study` | unpinned | unpinned | **defect** | pretrends lead regression: CR1 k converged on K_reference (D2 fixed); normal-theory tail df remains (PR C family) | +| `imputation_pretrends_event_study` | unpinned | 157, 157, 157, None, None, None, None, None | **legitimate** | pretrends lead regression: CR1 k = K_reference (D2 fixed); tail df converged in 3.9 (M-127): leads use t(residual df) via the df_convention knob (was silent normal theory); the None calls are the knob-independent BJS aggregates (L3) | | `two_stage_default` | — | None | **legitimate** | L3: Gardner two-stage variance, not the shared CR1 sandwich | | `callaway_santanna_default` | — | None, None, None, None, None, None, None, None, None, None, None, None, None, None, None, None | **legitimate** | L3: influence-function variance anchored to Stata csdid | @@ -46,10 +46,10 @@ vcov_type="hc1"; a clustered call in any other family fails the row, so a surface cannot silently switch clustered family behind an unchanged design width); — means the surface's *contract* is that it never calls it. tail_df is the multiset of df values passed to safe_inference/safe_inference_batch -(None = normal theory). unpinned marks the one contract row whose test -asserts the shared CR1 IS reached but deliberately pins no exact values — -a literal there would be brittle configuration-detail (its expected_adjustment -IS pinned: +6, the [time, unit] increment on df_0). Captured under the +(None = normal theory). Since 3.9 (M-127) every previously-unpinned row pins +its exact tail_df values — the imputation pretrends contract row included +(its expected_adjustment stays pinned too: +6, the [time, unit] increment +on df_0). Captured under the canonical Python backend; on adjusted clustered surfaces the Rust lanes apply the same correction as an exact scalar rescale of the finished vcov (vcov * (n-k)/(n-k-adj), <= 1 ulp on top of the default lane's <= 8e-15 @@ -107,18 +107,27 @@ output). DEFAULT ``K.exact = FALSE`` reproduces the old naive count, so this is a documented deviation from the R default (see the REGISTRY absorbed-FE note and ``tests/test_variance_conventions.py::TestFixestKExactParity``). -- **D4 — SunAbraham reports two tail-df conventions inside one fit** (residual - df on per-cell inference, normal theory on aggregates — visible in its row's - multiset). Fix: PR C. +- **D4 — SunAbraham reported two tail-df conventions inside one fit** (residual + df on per-cell inference, normal theory on aggregates — the old row multiset + was `(280.0,)*15 + (None,)*8`). **Fixed** in 3.9 (PR C / M-127): aggregates + resolve through the SAME saturated fit the cells use (`reg.df_` / + `reg.n_clusters_` via `utils.resolve_tail_df`), so one fit carries one + convention under every `df_convention` value. The same fix closed the wider + normal-theory tail-df defect family: Wooldridge default-hc1, StackedDiD, and + the ImputationDiD pretrends leads all moved z → t(residual df) under the + new default (see the row reasons in the matrix above; `"normal"` reproduces + the pre-3.9 numbers). ## Legitimate differences (declared exceptions) - **L1 — StackedDiD's k_total**: its design genuinely is a Q-weighted full-dummy lm; pinned to clubSandwich::vcovCR(type="CR1S") at atol=1e-10. CR1S is a real second convention, correct by construction. -- **L2 — LPDiD's G-1 tail df**: the Stata/fixest convention, and the only - surface where it is the default. It is the convergence target for PR C, not a - defect. +- **L2 — LPDiD's G-1 tail df**: the Stata/fixest convention. Since 3.9 it is + the `df_convention="cluster"` DEFAULT on LPDiD (the one surface already at + the library-wide v4 target — bit-identical, no numbers moved), with + "residual"/"normal" as opt-ins; the degenerate lanes (unclustered refit, + RA G<=1, saturated early return) keep literal df=None under every value. - **L3 — CallawaySantAnna / TwoStageDiD / ImputationDiD (default)**: different variance theory (influence functions / two-stage / BJS imputation), never the shared CR1 sandwich. CS is anchored to Stata csdid outright. @@ -133,16 +142,24 @@ output). family applies no CR1 finite-sample factor, so it has no cell on the axis this matrix measures. -## Tail-df landscape (PR C's input) +## Tail-df landscape (converged, 3.9 / M-127) -Three conventions are live: normal theory (Wooldridge, StackedDiD — no -df_convention knob), residual n - K_full (DiD/MPD/TWFE default; -df_convention="cluster" opts into G-1), and G-1 (LPDiD, hardcoded). -At |t| = 2 normal theory understates the t(G-1) p-value by 24.2% at -G=20, 13.3% at G=40, 1.2% at G=500 — larger than the D2 SE gap was before its -fix. PR C decisions: -(1) extending the knob to Wooldridge/SunAbraham/StackedDiD/ImputationDiD needs -NEW ledger rows (M-004/M-005/M-006 cover only DiD/TWFE/LinearRegression); -(2) the two-value knob cannot express normal theory, so either a third value -("normal") keeps 3.9 additive, or a documented default change ships with its -own ledger rows. +ONE three-value knob (`df_convention ∈ {"residual", "cluster", "normal"}`) +now governs the analytical fallback df on every shared-CR1 surface: +DiD/MPD/TWFE/LinearRegression + SunAbraham/WooldridgeDiD-OLS/StackedDiD/ +ImputationDiD-pretrends (default "residual") and LPDiD (default "cluster", +its Stata-anchored t(G−1)). Survey/replicate df and hc2_bm Bell-McCaffrey +DOF keep precedence under every value; "cluster" fails closed warn+NaN at +G<=1 and is inert on unclustered/conley fits; "normal" is deliberate z at +the fallback level everywhere. Per-surface residual-df definitions: the +fit's `n_eff − K_full` (visible kept columns + absorbed rank where the +design absorbs FE; positive-weight n_eff on weighted lanes; on the LPDiD RA +path `n_total − k0_kept − 1`, the pooled M-estimator count). The pre-fix +defect magnitude for the record: at |t| = 2 normal theory understated the +t(G−1) p-value by 24.2% at G=20, 13.3% at G=40, 1.2% at G=500 — larger +than the D2 SE gap was before its fix. The remaining program step is the +v4 default flip "residual" → "cluster" (M-004..M-006 + M-128..M-131, the +locked #663 direction); the resolver is `diff_diff.utils.resolve_tail_df`, +and the out-of-scope roster (IF-based estimators, TripleDifference/TROP, +Wooldridge GLM arms) is enumerated in the REGISTRY TwoWayFixedEffects +deviation note. diff --git a/docs/v4-deprecations.yaml b/docs/v4-deprecations.yaml index a1b025b28..f1abb6ead 100644 --- a/docs/v4-deprecations.yaml +++ b/docs/v4-deprecations.yaml @@ -98,7 +98,7 @@ rows: old_default: "'residual'" new_default: "'cluster'" code_refs: [diff_diff/linalg.py] - notes: "Evaluate extending the knob to standalone CR1-t estimators in the same PR (DEFERRED.md row)." + notes: "The 'evaluate extending the knob to standalone CR1-t estimators' item was resolved EARLY by PR C ([M-127], 3.9): the knob now lives on SunAbraham/WooldridgeDiD/StackedDiD/ImputationDiD (their flips are [M-128]..[M-131]) and LPDiD (already at 'cluster', no flip row)." - id: M-007 kind: warning-retirement group: obligation-warning-retirements @@ -1392,3 +1392,68 @@ rows: test_ref: tests/test_aggregate_contract.py code_refs: [diff_diff/aggregation.py, diff_diff/__init__.py] notes: "AggregationResult: the tabular container every non-event-study aggregation returns (aggregate(type='event_study') returns EventStudyResults instead, giving [M-092]'s container its public producer). Pinned 12-column schema with summary/to_dict/to_dataframe. Three schema choices are grounded in what the estimators actually produce: 'target' is PER-ROW so one container can carry two aligned estimands over the same labels; 'n_kind' reuses EventStudyResults' vocabulary rather than inventing units/obs; and 'weight' is NULLABLE because CallawaySantAnna's group aggregation weights (g,t) cells equally within each cohort and forms no cross-cohort mass. introduced_in gates the 3.9 cut ([M-091]/[M-092] pattern); deprecated_in stays null so the early-flip guard does not fire against the PR that ships it." + - id: M-127 + kind: behavior + group: tail-df-consolidation + old: "diff_diff.utils:safe_inference" + new: null + introduced_in: "3.9" + deprecated_in: null + removed_in: null + status: done + phase: 2 + test_ref: tests/test_variance_conventions.py + code_refs: [diff_diff/utils.py, diff_diff/linalg.py, diff_diff/estimators.py, diff_diff/twfe.py, diff_diff/results.py, diff_diff/sun_abraham.py, diff_diff/wooldridge.py, diff_diff/wooldridge_results.py, diff_diff/stacked_did.py, diff_diff/stacked_did_results.py, diff_diff/imputation.py, diff_diff/imputation_results.py, diff_diff/lpdid.py, diff_diff/lpdid_results.py, docs/methodology/REGISTRY.md, docs/methodology/variance-conventions.md] + notes: "Tail-df consolidation (D4 + the normal-theory tail-df defect family of the 3.9 variance-consolidation program, PR C). df_convention becomes a THREE-VALUE knob {'residual','cluster','normal'} on DiD/MPD/TWFE/LinearRegression (previously two-value; existing values unchanged, M-004..M-006 old_default untouched; 'normal' is deliberate z at the fallback level on EVERY fit, unclustered and conley included) and is EXTENDED to SunAbraham / WooldridgeDiD(OLS arms) / StackedDiD / ImputationDiD(pretrends lead regression) / LPDiD via the new utils.resolve_tail_df fallback resolver (survey/replicate df and hc2_bm Bell-McCaffrey per-coefficient/contrast DOF keep precedence under every value; 'cluster' fails closed warn+NaN at G<=1; unclustered/conley fits are inert for 'cluster'; a non-positive residual df falls back to normal theory with a warning). DEFECT FIX shipped as the new 'residual' default: WooldridgeDiD default-hc1 (within AND the hc1+cohort_trends full-dummy arm), StackedDiD, SunAbraham aggregates (D4: one fit reported residual-t cells and normal-theory aggregates - aggregates now reuse the saturated regression's df_ so cells and aggregates share ONE df source), and ImputationDiD pretrends leads previously passed df=None (silent normal-theory z) on clustered analytical paths; at |t|=2 z understates the t(G-1) p-value by 24.2% at G=20. Under the residual default the numeric movement is z vs t(large residual df) - small but honest; the G-1 convergence with fixest/Stata is the 4.0 flip. 'normal' reproduces the pre-3.9 numbers on the fixed surfaces' affected lanes (on lanes already at t - unclustered DiD-family fits, Wooldridge classical/hc2, SA cells - it is a NEW deliberate z option, not a compatibility path). LPDiD keeps its Stata-anchored t(G-1) as default 'cluster' (bit-identical; already at the 4.0 target, no flip row; degenerate lanes - unclustered refit, RA G<=1, saturated early return - keep literal df=None under all values; RA-path 'residual' = n_total - k0_kept - 1, the pooled M-estimator parameter count, a library convention with no external anchor since Stata teffects ra reports z). Wooldridge logit/poisson arms (survey df or normal theory, knob-independent, explicit-non-default warns) and all IF-based estimators (CallawaySantAnna, TwoStageDiD, dCDH, EfficientDiD, ContinuousDiD, ChangesInChanges, SyntheticDiD, HAD - L3) plus OLS-based TripleDifference/TROP are out of the knob's scope. Wild-cluster-bootstrap p/CIs are test-inversion based and invariant. Provenance: df_convention echoed on the five standalone results containers; SunAbraham/StackedDiD gain inference_df (the overall-ATT df actually used; None under bootstrap); WooldridgeDiDResults._df_one_way renamed _df_analytic_fallback (now the resolved survey-first fallback on ALL OLS arms, live on bootstrap fits so post-fit aggregate() reproduces fit-time inference) with __setstate__ pickle migration. Defaults for the four new residual surfaces flip to 'cluster' at 4.0 (M-128..M-131, same policy as M-004..M-006). introduced_in gates the 3.9 cut; deprecated_in stays null so the early-flip guard does not fire against the PR that ships it." + - id: M-128 + kind: default-flip + group: df-convention-flip + old: "diff_diff:SunAbraham[df_convention]" + new: null + deprecated_in: "4.0" + removed_in: null + status: planned + phase: 5 + old_default: "'residual'" + new_default: "'cluster'" + code_refs: [diff_diff/sun_abraham.py] + notes: "PR C ([M-127]) extended the knob here with default 'residual'; flips with [M-004]..[M-006] at 4.0 (G-1 = the fixest/Stata clustered convention). Moves every clustered p-value/CI on this surface; covered by the same migration note." + - id: M-129 + kind: default-flip + group: df-convention-flip + old: "diff_diff:WooldridgeDiD[df_convention]" + new: null + deprecated_in: "4.0" + removed_in: null + status: planned + phase: 5 + old_default: "'residual'" + new_default: "'cluster'" + code_refs: [diff_diff/wooldridge.py] + notes: "PR C ([M-127]) extended the knob here with default 'residual' (OLS arms only; GLM arms are knob-independent); flips with [M-004]..[M-006] at 4.0. Moves every clustered OLS p-value/CI on this surface; covered by the same migration note." + - id: M-130 + kind: default-flip + group: df-convention-flip + old: "diff_diff:StackedDiD[df_convention]" + new: null + deprecated_in: "4.0" + removed_in: null + status: planned + phase: 5 + old_default: "'residual'" + new_default: "'cluster'" + code_refs: [diff_diff/stacked_did.py] + notes: "PR C ([M-127]) extended the knob here with default 'residual'; flips with [M-004]..[M-006] at 4.0. The stacked residual df is typically large (t ~= z), so the flip to G-1 is the materially moving step on this surface; covered by the same migration note." + - id: M-131 + kind: default-flip + group: df-convention-flip + old: "diff_diff:ImputationDiD[df_convention]" + new: null + deprecated_in: "4.0" + removed_in: null + status: planned + phase: 5 + old_default: "'residual'" + new_default: "'cluster'" + code_refs: [diff_diff/imputation.py] + notes: "PR C ([M-127]) extended the knob here with default 'residual' (pretrends lead regression only; BJS overall inference and the joint pretrend Wald F are knob-independent); flips with [M-004]..[M-006] at 4.0; covered by the same migration note." diff --git a/docs/v4-design.md b/docs/v4-design.md index 1c2c39151..618218b68 100644 --- a/docs/v4-design.md +++ b/docs/v4-design.md @@ -37,8 +37,10 @@ lands or is explicitly re-scheduled at the 4.0 cut - tracked by the matrix, enforced by CI. **Non-goals.** No numerical behavior changes in any phase except the two -scheduled default policies ([M-004]..[M-006] df_convention, [M-080] -auto-cluster) and the documented estimate/SE shift of the MPD merge +scheduled default policies ([M-004]..[M-006] + [M-128]..[M-131] +df_convention, [M-080] auto-cluster), the ledgered 3.9 defect fixes of the +variance-consolidation program ([M-126] K_reference SEs, [M-127] tail-df +consolidation), and the documented estimate/SE shift of the MPD merge (section 4.1). No new estimators. No merging across identification strategies: the staggered family (CallawaySantAnna, SunAbraham, ImputationDiD, TwoStageDiD, EfficientDiD, @@ -490,12 +492,16 @@ domain vocabulary, not drift. Cross-sectional 2x2 estimators stay HC-robust unless `cluster=` is given. StackedDiD's hard-coded `cluster="unit"` default becomes an instance of the general policy rather than a special case. -- **df_convention default flip** [M-004]..[M-006]: `"residual"` -> - `"cluster"` (G-1) at 4.0 - the locked #663 direction. Moves every clustered +- **df_convention default flip** [M-004]..[M-006] + [M-128]..[M-131]: + `"residual"` -> `"cluster"` (G-1) at 4.0 - the locked #663 direction. + PR C ([M-127], 3.9) resolved the "evaluate extending the knob" item EARLY: + the knob is now a THREE-VALUE `{"residual","cluster","normal"}` surface on + DiD/MPD/TWFE/LinearRegression plus SunAbraham/WooldridgeDiD-OLS/StackedDiD/ + ImputationDiD-pretrends (their 4.0 flips are [M-128]..[M-131]) and LPDiD + (default already `"cluster"`, no flip row). The flip moves every clustered p-value/CI; the flip PR updates `TestDfConvention` / - `test_moderate_t_pins_residual_df_convention` expectations, adds the - migration-guide entry, and evaluates extending the knob to standalone CR1-t - estimators (DEFERRED.md row). + `test_moderate_t_pins_residual_df_convention` / the per-estimator knob + suites' expectations and adds the migration-guide entry. - Constructor hygiene rides Phase 2: ContinuousDiD's `covariates` moves from `__init__` to `fit()` [M-084] (the one estimator with a data column in the constructor), and the shared `BaseEstimator` mixin replaces the 24 @@ -597,7 +603,7 @@ above; anything only one PR cares about stays in that PR's plan.** | 2: contract foundations | 3.9 | (a) results base + unified event-study representation [M-092] + to_dict completion + the Diagnostic marker base on the diagnostic result roster [M-091] (section 3.5); (b) `aggregate()` + fit(aggregate=) shims [M-020..M-027]; (c) param renames [M-030..M-047] [M-084] [M-086..M-089] + their results-field mirrors [M-094] [M-095] (section 8 rule 9) + the public-function completeness sweep [M-097..M-113] (section 8 rule 10) + the dCDH results mirror [M-114] + the fourth `robust` site [M-115] + BaseEstimator mixin + ContinuousDiD covariates move; (d) alias introductions [M-062] [M-063] + wrapper deprecations [M-070..M-077] + the two inference-surface policies: `n_bootstrap` semantic unification [M-081] and the wild-cluster-bootstrap roster guard [M-096] | | 3: merges | 3.9 | (a) TWFE event-study mode [M-010] + EventStudy warn [M-060] (gates: section 4.1's equivalence/divergence/pooled-parity test triple); (b) TripleDifference facade [M-013]; (c) CiC method= [M-015] | | 4: release + soak | 3.9 cut | Migration guide written (skeleton: section 10); maintainer cuts 3.9; maint/3.8 rule active | -| 5: enforcement | 4.0 | Removals [M-010..M-016, M-030..M-047 old names, M-060, M-061, M-070..M-077, M-001..M-003] + the amendment's old names [M-094] [M-095] [M-097..M-115] (incl. their consumer migrations and the `clean_control` serialized reporting key); storage flips [M-050..M-058]; default policies [M-004..M-006, M-080]; warning retirement [M-007]; fastpath go/no-go [M-008]; diagnostic-family docs/roster reorganization [M-090]; sentinel retirement [M-093]; docs/llms.txt/README refresh | +| 5: enforcement | 4.0 | Removals [M-010..M-016, M-030..M-047 old names, M-060, M-061, M-070..M-077, M-001..M-003] + the amendment's old names [M-094] [M-095] [M-097..M-115] (incl. their consumer migrations and the `clean_control` serialized reporting key); storage flips [M-050..M-058]; default policies [M-004..M-006, M-128..M-131, M-080]; warning retirement [M-007]; fastpath go/no-go [M-008]; diagnostic-family docs/roster reorganization [M-090]; sentinel retirement [M-093]; docs/llms.txt/README refresh | | 6: front door | 4.1 | `event_study(data, outcome, unit, time, first_treat, estimator=...)` comparison entry point over the staggered family (sketch only; specified in its own plan) | **3.9-cut checklist (un-rowed obligations).** `test_due_rows_are_terminal` @@ -743,11 +749,11 @@ forever - a removed symbol resurrecting is a test failure. class row (schema-enforced). Top-level `diff_diff:Name` class/function rows and alias rows also assert `__all__` membership consistent with their status (stale `import *` entries fail). The shipped row ids are a - committed snapshot in the enforcement test (105 as of the clustered-CR1 - K_reference convergence row: Phase 1 + the diagnostic-family amendment + + committed snapshot in the enforcement test (110 as of the tail-df + consolidation family: Phase 1 + the diagnostic-family amendment + the M-092/M-093 results-contract rows + the M-094..M-096 amendment rows + the M-097..M-115 completeness sweep + M-117/M-122 + the ETWFE - reference-period pair M-123/M-124 + M-125 + M-126; + reference-period pair M-123/M-124 + M-125 + M-126 + M-127..M-131; the snapshot extends by a new id range in the same diff that appends rows): ids are never deleted or reused, and the test fails if any snapshot id disappears. diff --git a/tests/test_estimators_vcov_type.py b/tests/test_estimators_vcov_type.py index d244fee6d..b4941a9d5 100644 --- a/tests/test_estimators_vcov_type.py +++ b/tests/test_estimators_vcov_type.py @@ -2901,8 +2901,11 @@ def test_results_metadata_did_twfe(self): def test_conley_fits_excluded_from_knob(self): """vcov_type="conley" + explicit cluster keeps the residual df even - with the knob on (no documented G-1 reference for the product - kernel) — inference identical with the knob on and off.""" + with df_convention="cluster" (no documented G-1 reference for the + product kernel) — inference identical with "cluster" on and off. + PER-VALUE scope (3.9): this exclusion covers "cluster" only; + "normal" is deliberately UNGATED under conley (see + test_normal_is_z_on_every_fit below).""" rng = np.random.default_rng(9) data = self._clustered_panel() units = data["unit"].unique() @@ -3093,3 +3096,110 @@ def test_unclustered_fit_knob_is_inert(self): r0 = DifferenceInDifferences().fit(data, **kw) r1 = DifferenceInDifferences(df_convention="cluster").fit(data, **kw) assert (r0.p_value, r0.conf_int) == (r1.p_value, r1.conf_int) + + +class TestDfConventionNormal: + """The third knob value "normal" on the ORIGINAL knob surfaces + (DiD/TWFE/MPD/LinearRegression), added 3.9 (M-127): deliberate + normal-theory z at the fallback level on EVERY fit — clustered, + unclustered, AND conley (the ungated-fallback semantic; "cluster" + stays conley-excluded, "normal" does not). Survey df and BM DOF keep + precedence. On these surfaces "normal" is a NEW deliberate option: + the pre-3.9 fallback here was already t(residual). + """ + + @staticmethod + def _panel(seed=11): + rng = np.random.default_rng(seed) + rows = [] + for u in range(40): + for t in range(6): + treated = int(u < 20) + post = int(t >= 3) + rows.append( + dict( + unit=u, + time=t, + group=treated, + post=post, + y=0.5 * u / 10 + 0.2 * t + 0.8 * treated * post + rng.standard_normal(), + ) + ) + return pd.DataFrame(rows) + + _kw = dict(outcome="y", treatment="group", time="post") + + def test_normal_is_z_on_every_fit(self): + from scipy import stats + + data = self._panel() + # clustered + rc = DifferenceInDifferences(cluster="unit", df_convention="normal").fit(data, **self._kw) + assert rc.p_value == pytest.approx(2 * stats.norm.sf(abs(rc.t_stat)), rel=1e-14) + # unclustered — the ungated-fallback semantic (was t(residual)) + ru = DifferenceInDifferences(df_convention="normal").fit(data, **self._kw) + assert ru.p_value == pytest.approx(2 * stats.norm.sf(abs(ru.t_stat)), rel=1e-14) + r_res = DifferenceInDifferences().fit(data, **self._kw) + assert ru.p_value != r_res.p_value and ru.se == r_res.se + # conley — discriminates against an implementation that leaves + # "normal" inert under conley's cluster-exclusion + rng = np.random.default_rng(9) + units = data["unit"].unique() + data = data.assign( + lat=data["unit"].map({u: 40 + rng.uniform(-2, 2) for u in units}), + lon=data["unit"].map({u: -100 + rng.uniform(-2, 2) for u in units}), + ) + common = dict( + vcov_type="conley", + conley_coords=("lat", "lon"), + conley_cutoff_km=100.0, + conley_lag_cutoff=0, + ) + rk = DifferenceInDifferences(**common, df_convention="normal").fit( + data, unit="unit", **self._kw + ) + assert rk.p_value == pytest.approx(2 * stats.norm.sf(abs(rk.t_stat)), rel=1e-14) + rk_res = DifferenceInDifferences(**common).fit(data, unit="unit", **self._kw) + assert rk.p_value != rk_res.p_value + + def test_twfe_and_mpd_normal(self): + from scipy import stats + + data = self._panel() + rt = TwoWayFixedEffects(cluster="unit", df_convention="normal").fit( + data, outcome="y", treatment="group", time="post", unit="unit" + ) + assert rt.p_value == pytest.approx(2 * stats.norm.sf(abs(rt.t_stat)), rel=1e-14) + rm = MultiPeriodDiD(cluster="unit", df_convention="normal").fit( + data, + outcome="y", + treatment="group", + time="time", + post_periods=[3, 4, 5], + reference_period=2, + ) + assert rm.inference_df is None + assert rm.avg_p_value == pytest.approx(2 * stats.norm.sf(abs(rm.avg_t_stat)), rel=1e-13) + + def test_linear_regression_normal_and_precedence(self): + from scipy import stats + + from diff_diff.linalg import LinearRegression + + rng = np.random.default_rng(7) + X = rng.standard_normal((100, 2)) + y = 1 + 2 * X[:, 0] + rng.standard_normal(100) + cl = np.repeat(np.arange(10), 10) + reg = LinearRegression(cluster_ids=cl, df_convention="normal").fit(X, y) + inf = reg.get_inference(1) + assert inf.p_value == pytest.approx(2 * stats.norm.sf(abs(inf.t_stat)), rel=1e-14) + assert inf.df is None + # survey df precedence survives "normal" (fallback-level only) + reg.survey_df_ = 7 + inf2 = reg.get_inference(1) + assert inf2.df == 7 + + def test_validation_accepts_normal_everywhere(self): + for ctor in (DifferenceInDifferences, TwoWayFixedEffects, MultiPeriodDiD): + est = ctor(df_convention="normal") + assert est.get_params()["df_convention"] == "normal" diff --git a/tests/test_event_study_surface.py b/tests/test_event_study_surface.py index 80365f070..b1e8b4343 100644 --- a/tests/test_event_study_surface.py +++ b/tests/test_event_study_surface.py @@ -1287,9 +1287,13 @@ def test_stacked_hc2_bm_per_row_df(): if et in bm.event_study_df and np.isfinite(bm.event_study_df[et]): assert df_val == bm.event_study_df[et] - # hc1 non-survey: normal-theory inference -> no df anywhere. + # hc1 non-survey (3.9 / M-127): the df_convention default resolves to + # t(pooled residual df), so every finite-p row now carries that df. s_hc1 = build_event_study_surface(hc1) - assert np.isnan(s_hc1.df).all() + finite_hc1 = np.isfinite(s_hc1.p_value) + assert np.isfinite(s_hc1.df[finite_hc1]).all() + vals = {v for v in s_hc1.df[finite_hc1].tolist()} + assert len(vals) == 1 and vals == {hc1.inference_df} def test_stacked_survey_df_broadcast(): @@ -1562,9 +1566,12 @@ def test_sun_abraham_hc2_bm_per_row_df(): if et in bm.event_study_df and np.isfinite(bm.event_study_df[et]): assert df_val == bm.event_study_df[et] - # Plain analytic fit: normal-theory inference -> no df anywhere. - assert all(np.isnan(v) for v in plain.event_study_df.values()) - assert np.isnan(build_event_study_surface(plain).df).all() + # Plain analytic fit (3.9 / M-127, the D4 fix): aggregates share the + # saturated fit's residual df, so provenance is FINITE on every row. + assert all(np.isfinite(v) and v > 0 for v in plain.event_study_df.values()) + s_plain = build_event_study_surface(plain) + finite_plain = np.isfinite(s_plain.p_value) + assert np.isfinite(s_plain.df[finite_plain]).all() def test_sun_abraham_survey_df_broadcast(): diff --git a/tests/test_fixest_did_twfe_parity.py b/tests/test_fixest_did_twfe_parity.py index b7673ec98..c58420ad1 100644 --- a/tests/test_fixest_did_twfe_parity.py +++ b/tests/test_fixest_did_twfe_parity.py @@ -173,3 +173,55 @@ def test_twfe_hetero_iid_matches_fixest_machine_precision(self): exp = golden["twfe_hetero"]["iid"] np.testing.assert_allclose(res.att, exp["att"], atol=1e-10, rtol=0) np.testing.assert_allclose(res.se, exp["se"], atol=1e-10, rtol=0) + + +@_SKIP +class TestFixestTailDfParity: + """External p-value/CI gates against the stored fixest values (3.9 / M-127). + + The golden has carried fixest's `p_value`/`ci_lower`/`ci_upper`/`t_stat` + since it was generated, but nothing asserted them until the tail-df + consolidation. Two lanes: + + - iid arms under the DEFAULT `df_convention="residual"`: fixest's iid + t-reference is the residual df, so Python matches directly (measured + during plan review: Python residual df == fixest implied df on all four + arms — did 396, twfe 148, did_hetero 209, twfe_hetero 125; p/CI agree + to ~1e-13 relative). + - cluster arms under `df_convention="cluster"`: fixest's clustered + t-reference is G−1, exactly the knob's value. + + Tolerances are PER-ARM, inheriting the SE-band structure: the hetero + arms' 5.2e-11 relative SE residual amplifies ~12x into the p-value at + t≈3.5/df=125, so their p/CI pin at rtol=1e-8 (≥20x headroom on the + measured ~6.3e-10); the exact-SE arms pin at rtol=1e-9. + """ + + _P_TOL = {"did": 1e-9, "twfe": 1e-9, "did_hetero": 1e-8, "twfe_hetero": 1e-8} + + def _fit(self, key, golden, **est_kw): + df = _build_df(golden[key]) + cls = TwoWayFixedEffects if key.startswith("twfe") else DifferenceInDifferences + return cls(**est_kw).fit( + df, outcome="outcome", treatment="treated", time="post", unit="unit" + ) + + def test_iid_p_and_ci_match_fixest_under_residual_default(self): + golden = _load_golden() + for key, tol in self._P_TOL.items(): + res = self._fit(key, golden, vcov_type="classical") + exp = golden[key]["iid"] + np.testing.assert_allclose(res.t_stat, exp["t_stat"], rtol=tol) + np.testing.assert_allclose(res.p_value, exp["p_value"], rtol=tol) + np.testing.assert_allclose(res.conf_int[0], exp["ci_lower"], rtol=tol) + np.testing.assert_allclose(res.conf_int[1], exp["ci_upper"], rtol=tol) + + def test_cluster_p_and_ci_match_fixest_under_cluster_knob(self): + golden = _load_golden() + for key, tol in self._P_TOL.items(): + res = self._fit(key, golden, vcov_type="hc1", cluster="unit", df_convention="cluster") + exp = golden[key]["cluster_unit"] + np.testing.assert_allclose(res.t_stat, exp["t_stat"], rtol=tol) + np.testing.assert_allclose(res.p_value, exp["p_value"], rtol=tol) + np.testing.assert_allclose(res.conf_int[0], exp["ci_lower"], rtol=tol) + np.testing.assert_allclose(res.conf_int[1], exp["ci_upper"], rtol=tol) diff --git a/tests/test_imputation.py b/tests/test_imputation.py index 1a2e868c5..09533c7f3 100644 --- a/tests/test_imputation.py +++ b/tests/test_imputation.py @@ -3262,3 +3262,142 @@ def _fake_lsmr(A, b, **kwargs): z = imp2._lsmr_minnorm_normal_solve(sp2.csc_matrix(np.eye(3)), np.ones(3)) assert len(calls) == 1 # accepted on the first attempt np.testing.assert_array_equal(z, np.full(3, 2.0)) + + +class TestImputationDfConvention: + """The three-value df_convention knob on ImputationDiD's pretrends lead + regression (3.9 / M-127) — the one ImputationDiD surface on the shared + clustered CR1 sandwich. BJS overall/post inference and the joint pretrend + Wald F are knob-independent. + """ + + @staticmethod + def _panel(seed=13): + rng = np.random.default_rng(seed) + rows = [] + for u in range(50): + ft = [0, 5, 7][u % 3] + for t in range(1, 10): + eff = 0.5 if (ft and t >= ft) else 0.0 + rows.append( + dict( + unit=u, + time=t, + first_treat=ft, + outcome=0.3 * u / 50 + 0.15 * t + eff + rng.standard_normal() * 0.5, + ) + ) + return pd.DataFrame(rows) + + _kw = dict( + outcome="outcome", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + + @staticmethod + def _t_p(t_stat, df): + from scipy import stats + + return 2 * stats.t.sf(abs(t_stat), df) + + @staticmethod + def _z_p(t_stat): + from scipy import stats + + return 2 * stats.norm.sf(abs(t_stat)) + + def _lead(self, res): + h = min(k for k in res.event_study_effects if k < -1) + return h, res.event_study_effects[h] + + def test_leads_are_t_residual_not_z(self): + res = ImputationDiD(pretrends=True).fit(self._panel(), **self._kw) + h, e = self._lead(res) + assert e["p_value"] != self._z_p(e["t_stat"]) + match = [d for d in range(2, 600) if abs(e["p_value"] - self._t_p(e["t_stat"], d)) < 1e-13] + assert len(match) == 1 + + def test_cluster_matches_g_minus_1_on_leads_only(self): + data = self._panel() + r0 = ImputationDiD(pretrends=True).fit(data, **self._kw) + rc = ImputationDiD(pretrends=True, df_convention="cluster").fit(data, **self._kw) + h, e0 = self._lead(r0) + ec = rc.event_study_effects[h] + G = data["unit"].nunique() + assert ec["effect"] == e0["effect"] and ec["se"] == e0["se"] + assert ec["p_value"] == pytest.approx(self._t_p(ec["t_stat"], G - 1), rel=1e-12) + # post rows are BJS (knob-independent) + hp = min(k for k in r0.event_study_effects if k >= 0) + assert r0.event_study_effects[hp]["p_value"] == rc.event_study_effects[hp]["p_value"] + + def test_normal_reproduces_pre39_z_on_leads(self): + data = self._panel() + r0 = ImputationDiD(pretrends=True).fit(data, **self._kw) + rn = ImputationDiD(pretrends=True, df_convention="normal").fit(data, **self._kw) + h, e0 = self._lead(r0) + en = rn.event_study_effects[h] + assert en["effect"] == e0["effect"] and en["se"] == e0["se"] + assert en["p_value"] == pytest.approx(self._z_p(en["t_stat"]), rel=1e-14) + + def test_pretrend_wald_f_is_knob_independent(self): + data = self._panel() + pt0 = ImputationDiD(pretrends=True).fit(data, **self._kw).pretrend_test() + ptc = ( + ImputationDiD(pretrends=True, df_convention="cluster") + .fit(data, **self._kw) + .pretrend_test() + ) + assert pt0["p_value"] == ptc["p_value"] + assert pt0["f_stat"] == ptc["f_stat"] + + def test_survey_df_precedence_on_leads(self): + """Survey design df wins on the pretrends leads under EVERY knob + value (the knob is never consulted on surveyed fits).""" + data = self._panel() + data["weight"] = 1.0 + (data["unit"] % 5) * 0.2 + data["stratum"] = data["unit"] % 4 + data["psu"] = data["unit"] + design = SurveyDesign(weights="weight", psu="psu", strata="stratum", weight_type="pweight") + fits = { + conv: ImputationDiD(pretrends=True, df_convention=conv).fit( + data, survey_design=design, **self._kw + ) + for conv in ("residual", "cluster", "normal") + } + h, e0 = self._lead(fits["residual"]) + for conv in ("cluster", "normal"): + e = fits[conv].event_study_effects[h] + assert e["p_value"] == e0["p_value"] + assert e["conf_int"] == e0["conf_int"] + + def test_inert_config_warns_on_explicit_nondefault(self): + data = self._panel() + base = dict(outcome="outcome", unit="unit", time="time", first_treat="first_treat") + with pytest.warns(UserWarning, match="affects only the pretrends"): + ImputationDiD(df_convention="cluster").fit(data, **base) + with pytest.warns(UserWarning, match="affects only the pretrends"): + ImputationDiD(pretrends=True, df_convention="cluster").fit( + data, aggregate="group", **base + ) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + ImputationDiD(pretrends=True, df_convention="cluster").fit(data, **self._kw) + assert not any("affects only the pretrends" in str(w.message) for w in caught) + + def test_validation_and_transactional_set_params(self): + with pytest.raises(ValueError, match="df_convention"): + ImputationDiD(df_convention="bogus") + est = ImputationDiD() + with pytest.raises(ValueError, match="df_convention"): + est.set_params(df_convention="bogus", alpha=0.10) + assert est.df_convention == "residual" and est.alpha == 0.05 + # Valid value + unknown key: the unknown-key rejection must also + # leave the estimator fully unchanged (no partial application). + before = est.get_params() + with pytest.raises(ValueError, match="Unknown parameter"): + est.set_params(df_convention="normal", nonexistent_param=1) + assert est.get_params() == before + assert ImputationDiD(df_convention="normal").get_params()["df_convention"] == "normal" diff --git a/tests/test_lpdid.py b/tests/test_lpdid.py index 21d63c489..b4d83d743 100644 --- a/tests/test_lpdid.py +++ b/tests/test_lpdid.py @@ -1978,3 +1978,132 @@ def test_rejects_reserved_survey_column_name(self): sd = SurveyDesign(weights="_weight", strata="stratum", psu="psu") with pytest.raises(ValueError): LPDiD(post_window=2).fit(df, survey_design=sd, **_FIT_KW) + + +class TestLPDiDDfConvention: + """The three-value df_convention knob on LPDiD (3.9 / M-127). + + LPDiD's DEFAULT is "cluster" — the pre-existing Stata lpdid t(G-1) + reference — so the default moves NOTHING (unlike the "residual"-default + estimators). "residual"/"normal" are new opt-ins; the degenerate lanes + (unclustered refit, RA G<=1, saturated early return) keep literal + df=None under every value. + """ + + @staticmethod + def _panel(seed=21): + rng = np.random.default_rng(seed) + rows = [] + for u in range(60): + entry = 6 if u % 3 == 0 else (8 if u % 3 == 1 else 99) + for t in range(1, 13): + d = int(entry <= t) + rows.append( + dict( + unit=u, + time=t, + treat=d, + outcome=0.2 * u / 60 + 0.1 * t + 0.5 * d + rng.standard_normal() * 0.5, + ) + ) + return pd.DataFrame(rows) + + _kw = dict(outcome="outcome", unit="unit", time="time", treatment="treat") + + @staticmethod + def _t_p(t_stat, df): + from scipy import stats + + return 2 * stats.t.sf(abs(t_stat), df) + + def test_default_cluster_is_todays_g_minus_1(self): + res = LPDiD(pre_window=2, post_window=2).fit(self._panel(), **self._kw) + assert res.df_convention == "cluster" + row = res.event_study.dropna(subset=["p_value"]).iloc[-1] + G = row["n_clusters"] + assert row["p_value"] == pytest.approx(self._t_p(row["t_stat"], G - 1), rel=1e-12) + # provenance carries exactly G-1 per row + finite = {v for v in res.event_study_df.values() if np.isfinite(v)} + assert finite and all(v > 0 for v in finite) + + def test_residual_optin_moves_off_g_minus_1(self): + data = self._panel() + r0 = LPDiD(pre_window=2, post_window=2).fit(data, **self._kw) + rr = LPDiD(pre_window=2, post_window=2, df_convention="residual").fit(data, **self._kw) + row0 = r0.event_study.dropna(subset=["p_value"]).iloc[-1] + roww = rr.event_study.dropna(subset=["p_value"]).iloc[-1] + assert roww["coefficient"] == row0["coefficient"] and roww["se"] == row0["se"] + assert roww["p_value"] != row0["p_value"] + # residual df = n_eff - k_kept of the per-horizon design + match = [ + d for d in range(2, 3000) if abs(roww["p_value"] - self._t_p(roww["t_stat"], d)) < 1e-13 + ] + assert len(match) == 1 and match[0] > row0["n_clusters"] - 1 + + def test_ra_path_residual_df_is_n_minus_k0_minus_1(self): + """The RA (reweight) lane's 'residual' df = n_total - k0_kept - 1 (the + pooled M-estimator's parameter count: kept nuisance coefficients plus + the ATT). Pinned via the 'df' provenance the row carries.""" + data = self._panel() + rr = LPDiD(pre_window=2, post_window=2, reweight=True, df_convention="residual").fit( + data, **self._kw + ) + # RA rows carry a df provenance entry; reconstruct one horizon's df + assert rr.pooled_df is not None + post_df = rr.pooled_df.get("post") + assert post_df is not None and np.isfinite(post_df) and post_df > 0 + row = rr.pooled.loc[rr.pooled["window"] == "post"].iloc[0] + # n_total - k0_kept - 1 < n_obs (sanity: subtracts at least the ATT) + assert post_df < row["n_obs"] + assert row["p_value"] == pytest.approx(self._t_p(row["t_stat"], post_df), rel=1e-12) + + def test_default_cluster_bit_identical_to_pre_knob_values(self): + """Regression lock: the knob default reproduces the historical G-1 + inference bit-for-bit (LPDiD zero-movement guarantee).""" + data = self._panel() + res = LPDiD(pre_window=2, post_window=2).fit(data, **self._kw) + ev = res.event_study.dropna(subset=["p_value"]) + for _, row in ev.iterrows(): + assert row["p_value"] == pytest.approx( + self._t_p(row["t_stat"], row["n_clusters"] - 1), rel=1e-12 + ) + + def test_normal_optin_is_z(self): + from scipy import stats + + rn = LPDiD(pre_window=2, post_window=2, df_convention="normal").fit( + self._panel(), **self._kw + ) + row = rn.event_study.dropna(subset=["p_value"]).iloc[-1] + assert row["p_value"] == pytest.approx(2 * stats.norm.sf(abs(row["t_stat"])), rel=1e-14) + assert all(np.isnan(v) for v in rn.event_study_df.values()) + + def test_survey_df_precedence_over_knob(self): + """The survey lane keeps its own df under EVERY knob value: the + three fits' event studies are bit-identical (the knob is never + consulted on surveyed fits).""" + df = _survey_panel() + fits = { + conv: LPDiD(pre_window=2, post_window=2, df_convention=conv).fit( + df, survey_design=_full_design(), **_FIT_KW + ) + for conv in ("cluster", "residual", "normal") + } + base = fits["cluster"] # the default + for conv in ("residual", "normal"): + pd.testing.assert_frame_equal(fits[conv].event_study, base.event_study) + + def test_validation_and_transactional_set_params(self): + with pytest.raises(ValueError, match="df_convention"): + LPDiD(df_convention="bogus") + est = LPDiD() + with pytest.raises(ValueError, match="df_convention"): + est.set_params(df_convention="bogus", alpha=0.10) + assert est.df_convention == "cluster" and est.alpha == 0.05 + # Valid value + unknown key: the unknown-key rejection must also + # leave the estimator fully unchanged (no partial application). + before = est.get_params() + with pytest.raises(ValueError, match="Unknown parameter"): + est.set_params(df_convention="normal", nonexistent_param=1) + assert est.get_params() == before + assert LPDiD().get_params()["df_convention"] == "cluster" diff --git a/tests/test_stacked_did.py b/tests/test_stacked_did.py index 348df20f2..c5f04430a 100644 --- a/tests/test_stacked_did.py +++ b/tests/test_stacked_did.py @@ -1936,3 +1936,134 @@ def test_string_covariate_raises_typeerror(self): df = _balance_panel(seed=13, const_x=True) with pytest.raises(TypeError, match="must be a list"): _cb_fit(df, balance="entropy", covariates="x") + + +class TestStackedDfConvention: + """The three-value df_convention knob on StackedDiD (3.9 / M-127).""" + + @staticmethod + def _panel(seed=9, n_units=60): + from diff_diff import generate_staggered_data + + return generate_staggered_data(n_units=n_units, n_periods=10, seed=seed) + + _kw = dict( + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + aggregate="event_study", + ) + + @staticmethod + def _t_p(t_stat, df): + from scipy import stats + + return 2 * stats.t.sf(abs(t_stat), df) + + @staticmethod + def _z_p(t_stat): + from scipy import stats + + return 2 * stats.norm.sf(abs(t_stat)) + + def test_residual_default_is_t_with_pinned_df(self): + """The defect fix: hc1 lane moved z -> t(pooled residual df). Assert + the RESOLVED df value directly (a p-movement assertion is vacuous on + stacked designs: the residual df is large, t ~= z).""" + res = StackedDiD(kappa_pre=2, kappa_post=2).fit(self._panel(), **self._kw) + dfres = res.inference_df + assert dfres is not None and dfres > res.n_clusters + assert res.overall_p_value == pytest.approx(self._t_p(res.overall_t_stat, dfres), rel=1e-12) + finite = {v for v in res.event_study_df.values() if np.isfinite(v)} + assert finite == {dfres} + + def test_small_n_residual_default_moves_off_z(self): + """On a SMALL stacked design the z -> t movement is measurable.""" + res = StackedDiD(kappa_pre=1, kappa_post=1).fit(self._panel(n_units=12), **self._kw) + assert np.isfinite(res.overall_p_value) + assert res.overall_p_value != self._z_p(res.overall_t_stat) + + def test_cluster_matches_g_minus_1_effective_count(self): + data = self._panel() + r0 = StackedDiD(kappa_pre=2, kappa_post=2).fit(data, **self._kw) + r1 = StackedDiD(kappa_pre=2, kappa_post=2, df_convention="cluster").fit(data, **self._kw) + G = r1.n_clusters + assert r1.overall_att == r0.overall_att and r1.overall_se == r0.overall_se + assert r1.overall_p_value == pytest.approx(self._t_p(r1.overall_t_stat, G - 1), rel=1e-12) + assert r1.inference_df == G - 1 + + def test_cluster_g_uses_effective_count_not_raw_unique(self): + """A cluster whose TOTAL composed weight is zero is excluded from the + 'cluster' G (effective_cluster_count) even though results.n_clusters + keeps the raw unique count. Discriminated by zeroing one unit's rows + via a zero survey weight is unsupported here, so instead simulate at + the resolver level with the estimator's own ingredients: verify that + on an ordinary fit the two counts AGREE (the documented divergence is + unreachable without zero-weight clusters) and that the df equals + effective-count-minus-1.""" + from diff_diff.linalg import effective_cluster_count + + data = self._panel() + r1 = StackedDiD(kappa_pre=2, kappa_post=2, df_convention="cluster").fit(data, **self._kw) + # Reconstruct the effective count from the stacked frame the fit stored + stacked = r1.stacked_data + cluster_ids = stacked["unit"].values + # composed weights are strictly positive on this fixture -> equality + g_eff = effective_cluster_count(cluster_ids, None) + assert r1.inference_df == g_eff - 1 == r1.n_clusters - 1 + + def test_normal_reproduces_pre39_z(self): + data = self._panel() + r0 = StackedDiD(kappa_pre=2, kappa_post=2).fit(data, **self._kw) + rn = StackedDiD(kappa_pre=2, kappa_post=2, df_convention="normal").fit(data, **self._kw) + assert rn.overall_att == r0.overall_att and rn.overall_se == r0.overall_se + assert rn.overall_p_value == pytest.approx(self._z_p(rn.overall_t_stat), rel=1e-14) + assert rn.inference_df is None + + def test_bm_dof_precedence_over_knob(self): + data = self._panel() + r_res = StackedDiD(kappa_pre=2, kappa_post=2, vcov_type="hc2_bm").fit(data, **self._kw) + r_clu = StackedDiD( + kappa_pre=2, kappa_post=2, vcov_type="hc2_bm", df_convention="cluster" + ).fit(data, **self._kw) + assert r_res.overall_p_value == r_clu.overall_p_value + + def test_survey_df_precedence_over_knob(self): + """Survey design df wins under EVERY knob value: the three fits are + bit-identical on inference and share the finite survey-df provenance + (the knob is never consulted on surveyed fits).""" + from diff_diff.survey import SurveyDesign + + data = self._panel() + data["w"] = 1.0 + (data["unit"] % 5) * 0.2 + data["stratum"] = data["unit"] % 4 + design = SurveyDesign(weights="w", strata="stratum", psu="unit", weight_type="pweight") + fits = { + conv: StackedDiD(kappa_pre=2, kappa_post=2, df_convention=conv).fit( + data, survey_design=design, **self._kw + ) + for conv in ("residual", "cluster", "normal") + } + base = fits["residual"] + assert base.inference_df is not None and np.isfinite(base.inference_df) + for conv in ("cluster", "normal"): + r = fits[conv] + assert r.overall_p_value == base.overall_p_value + assert r.overall_conf_int == base.overall_conf_int + assert r.inference_df == base.inference_df + + def test_validation_and_transactional_set_params(self): + with pytest.raises(ValueError, match="df_convention"): + StackedDiD(df_convention="bogus") + est = StackedDiD() + with pytest.raises(ValueError, match="df_convention"): + est.set_params(df_convention="bogus", alpha=0.10) + assert est.df_convention == "residual" and est.alpha == 0.05 + # Valid value + unknown key: the unknown-key rejection must also + # leave the estimator fully unchanged (no partial application). + before = est.get_params() + with pytest.raises(ValueError, match="Unknown parameter"): + est.set_params(df_convention="normal", nonexistent_param=1) + assert est.get_params() == before + assert StackedDiD(df_convention="normal").get_params()["df_convention"] == "normal" diff --git a/tests/test_sun_abraham.py b/tests/test_sun_abraham.py index 0fb85f674..23df78c6f 100644 --- a/tests/test_sun_abraham.py +++ b/tests/test_sun_abraham.py @@ -7,6 +7,7 @@ import numpy as np import pandas as pd import pytest +from scipy import stats from diff_diff.sun_abraham import SunAbraham, SunAbrahamResults @@ -1567,12 +1568,13 @@ def _nan_dofs(*args, **kwargs): assert np.isnan(row["conf_int"][0]) and np.isnan(row["conf_int"][1]) def test_hc2_bm_contrast_dof_helper_failure_falls_back(self): - """Helper-level BM failure warns and falls back to normal theory. + """Helper-level BM failure warns and falls back to the resolved df. Distinct from the per-row NaN-DOF case above: when the helper RAISES, - SA leaves the ``_compute_iw_effects`` output in place (finite - normal-theory inference) after warning, and the df provenance records - the baseline (NaN on a non-survey fit) - not a fabricated df. + SA leaves the ``_compute_iw_effects`` output in place after warning. + Since 3.9 (M-127) that output carries the df_convention-resolved + analytical df - FINITE t(residual df) under the default, recorded as + genuine provenance (previously normal theory with NaN provenance). """ import diff_diff.linalg as _linalg @@ -1587,10 +1589,12 @@ def _boom(*args, **kwargs): res = SunAbraham(vcov_type="hc2_bm").fit(data, **kwargs) finite_p = [e for e, row in res.event_study_effects.items() if np.isfinite(row["p_value"])] - assert finite_p, "helper failure must fall back to finite normal-theory inference" + assert finite_p, "helper failure must fall back to finite t(residual) inference" assert res.event_study_df is not None for e in finite_p: - assert np.isnan(res.event_study_df[e]), "normal theory records NO df" + assert ( + np.isfinite(res.event_study_df[e]) and res.event_study_df[e] > 0 + ), "the fallback df is the knob-resolved residual df - genuine provenance" def test_hc2_bm_explicit_cluster_works(self): """Explicit cluster= overrides the auto-cluster default; CR2-BM at named cluster.""" @@ -2071,3 +2075,136 @@ def test_unit_constant_covariate_snaps_att_unaffected(self): ) # snapped covariate -> dropped -> identical to the no-covariate design assert res.att == pytest.approx(base.att, abs=1e-10) + + +class TestSunAbrahamDfConvention: + """The three-value df_convention knob on SunAbraham (3.9 / M-127, the D4 fix). + + Modeled on tests/test_estimators_vcov_type.py::TestDfConvention. Cells and + aggregates share ONE df source per fit; att/se/t never move across values. + """ + + @staticmethod + def _panel(seed=3): + rng = np.random.default_rng(seed) + rows = [] + for u in range(60): + ft = [0, 4, 6][u % 3] + for t in range(1, 9): + eff = 0.7 if (ft and t >= ft) else 0.0 + rows.append( + dict( + unit=u, + time=t, + first_treat=ft, + outcome=0.3 * u / 60 + 0.1 * t + eff + rng.standard_normal() * 0.6, + ) + ) + return pd.DataFrame(rows) + + _kw = dict(outcome="outcome", unit="unit", time="time", first_treat="first_treat") + + def test_residual_default_is_t_not_z(self): + """The D4/defect fix: aggregates use t(saturated residual df), not z.""" + res = SunAbraham().fit(self._panel(), **self._kw) + p_z = 2 * stats.norm.sf(abs(res.overall_t_stat)) + assert res.overall_p_value != p_z + dfres = res.inference_df + assert dfres is not None and dfres > 0 + assert res.overall_p_value == pytest.approx( + 2 * stats.t.sf(abs(res.overall_t_stat), dfres), rel=1e-12 + ) + # every ES row shares the same finite residual df (one df source) + assert set(res.event_study_df.values()) == {dfres} + + def test_cluster_matches_hand_reconstructed_g_minus_1(self): + data = self._panel() + r0 = SunAbraham().fit(data, **self._kw) + r1 = SunAbraham(df_convention="cluster").fit(data, **self._kw) + G = data["unit"].nunique() + assert r1.overall_att == r0.overall_att and r1.overall_se == r0.overall_se + assert r1.overall_t_stat == r0.overall_t_stat + assert r1.overall_p_value == pytest.approx( + 2 * stats.t.sf(abs(r1.overall_t_stat), G - 1), rel=1e-12 + ) + lo, hi = r1.overall_conf_int + crit = stats.t.ppf(0.975, G - 1) + assert lo == pytest.approx(r1.overall_att - crit * r1.overall_se, rel=1e-12) + assert hi == pytest.approx(r1.overall_att + crit * r1.overall_se, rel=1e-12) + assert r1.inference_df == G - 1 + assert set(r1.event_study_df.values()) == {float(G - 1)} + + def test_normal_reproduces_z_on_aggregates(self): + """'normal' = deliberate z at the fallback level (the pre-3.9 numbers).""" + data = self._panel() + r0 = SunAbraham().fit(data, **self._kw) + rn = SunAbraham(df_convention="normal").fit(data, **self._kw) + assert rn.overall_att == r0.overall_att and rn.overall_se == r0.overall_se + assert rn.overall_p_value == pytest.approx( + 2 * stats.norm.sf(abs(rn.overall_t_stat)), rel=1e-14 + ) + assert rn.inference_df is None + # provenance: z rows record NaN + assert all(np.isnan(v) for v in rn.event_study_df.values()) + + def test_bm_dof_precedence_over_knob(self): + """hc2_bm BM contrast DOF wins under every knob value.""" + data = self._panel() + r_res = SunAbraham(vcov_type="hc2_bm").fit(data, **self._kw) + r_clu = SunAbraham(vcov_type="hc2_bm", df_convention="cluster").fit(data, **self._kw) + assert r_res.overall_p_value == r_clu.overall_p_value + assert r_res.overall_conf_int == r_clu.overall_conf_int + + def test_survey_df_precedence_over_knob(self): + """Survey design df wins under EVERY knob value: the three fits are + bit-identical on inference and share the finite survey-df provenance + (the knob is never consulted on surveyed fits).""" + from diff_diff.survey import SurveyDesign + + data = self._panel() + data["w"] = 1.0 + (data["unit"] % 5) * 0.2 + data["stratum"] = data["unit"] % 4 + sd = SurveyDesign(weights="w", strata="stratum", psu="unit") + fits = { + conv: SunAbraham(df_convention=conv).fit(data, survey_design=sd, **self._kw) + for conv in ("residual", "cluster", "normal") + } + base = fits["residual"] + assert base.inference_df is not None and np.isfinite(base.inference_df) + for conv in ("cluster", "normal"): + r = fits[conv] + assert r.overall_p_value == base.overall_p_value + assert r.overall_conf_int == base.overall_conf_int + assert r.inference_df == base.inference_df + assert r.event_study_df == base.event_study_df + + def test_bootstrap_clears_inference_df(self): + """Percentile bootstrap p/CI never used a df: the scalar clears.""" + res = SunAbraham(n_bootstrap=20, seed=1).fit(self._panel(), **self._kw) + assert res.inference_df is None + assert res.event_study_df is None # existing channel-clearing convention + + def test_validation_and_transactional_set_params(self): + with pytest.raises(ValueError, match="df_convention"): + SunAbraham(df_convention="bogus") + est = SunAbraham() + with pytest.raises(ValueError, match="df_convention"): + est.set_params(df_convention="bogus", alpha=0.10) + assert est.df_convention == "residual" and est.alpha == 0.05 + # Valid value + unknown key: the unknown-key rejection must also + # leave the estimator fully unchanged (no partial application). + before = est.get_params() + with pytest.raises(ValueError, match="Unknown parameter"): + est.set_params(df_convention="normal", nonexistent_param=1) + assert est.get_params() == before + + def test_get_params_roundtrip_and_repeat_fit(self): + est = SunAbraham(df_convention="cluster") + assert est.get_params()["df_convention"] == "cluster" + clone = SunAbraham(**est.get_params()) + data = self._panel() + r1 = clone.fit(data, **self._kw) + p1, d1 = r1.overall_p_value, r1.inference_df + r2 = clone.fit(data, **self._kw) + assert r2.overall_p_value == p1 and r2.inference_df == d1 + assert clone.df_convention == "cluster" diff --git a/tests/test_utils.py b/tests/test_utils.py index 112f73a3b..1d241a339 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -2147,3 +2147,77 @@ def test_mismatched_prefixes_length_raises(self): df = pd.DataFrame({"a": ["x", "y"], "b": ["u", "v"]}) with pytest.raises(ValueError, match="prefixes length 1 does not match"): build_fe_dummy_blocks(df, ["a", "b"], prefixes=["_fe_a"]) + + +class TestResolveTailDf: + """Branch table for the shared df_convention fallback resolver (3.9 / M-127). + + The resolver owns ONLY the fallback slot: survey/replicate df and BM + per-coefficient DOF keep precedence at the call sites. Covers all three + conventions x n_clusters None/1/2+ and the non-positive/non-finite + residual_df guard (the reachable case is absorbed-FE accounting driving + n - k - rank <= 0 while the SE is still finite). + """ + + def test_residual_returns_residual_df(self): + from diff_diff.utils import resolve_tail_df + + assert resolve_tail_df("residual", residual_df=42, n_clusters=None) == 42.0 + assert resolve_tail_df("residual", residual_df=42, n_clusters=50) == 42.0 + + def test_cluster_returns_g_minus_1(self): + from diff_diff.utils import resolve_tail_df + + assert resolve_tail_df("cluster", residual_df=42, n_clusters=50) == 49.0 + assert resolve_tail_df("cluster", residual_df=42, n_clusters=2) == 1.0 + + def test_cluster_inert_when_unclustered(self): + """n_clusters=None (unclustered/conley fit) falls back to residual.""" + from diff_diff.utils import resolve_tail_df + + assert resolve_tail_df("cluster", residual_df=42, n_clusters=None) == 42.0 + + def test_cluster_fails_closed_at_one_cluster(self): + """G<=1: warn + the 0 sentinel safe_inference turns into all-NaN.""" + from diff_diff.utils import resolve_tail_df, safe_inference + + with pytest.warns(UserWarning, match="at least 2 effective"): + df = resolve_tail_df("cluster", residual_df=42, n_clusters=1) + assert df == 0.0 + t, p, ci = safe_inference(1.0, 0.5, df=df) + assert np.isnan(t) and np.isnan(p) and np.isnan(ci[0]) and np.isnan(ci[1]) + + def test_normal_returns_none_everywhere(self): + from diff_diff.utils import resolve_tail_df + + assert resolve_tail_df("normal", residual_df=42, n_clusters=None) is None + assert resolve_tail_df("normal", residual_df=42, n_clusters=50) is None + + def test_nonpositive_residual_df_warns_and_falls_back_to_normal(self): + """Mirrors get_inference's non-positive-df fallback: warn + None (z), + NOT an all-NaN tuple — preserves the historical df_one_way NaN -> + normal-theory degenerate lane bit-for-bit.""" + from diff_diff.utils import resolve_tail_df + + for bad in (0.0, -3.0, float("nan")): + with pytest.warns(UserWarning, match="non-positive"): + assert resolve_tail_df("residual", residual_df=bad, n_clusters=None) is None + # The guard applies on the "cluster" inert lane too (falls back to + # residual, which then guards). + with pytest.warns(UserWarning, match="non-positive"): + assert resolve_tail_df("cluster", residual_df=-1.0, n_clusters=None) is None + + def test_none_residual_df_passes_through(self): + from diff_diff.utils import resolve_tail_df + + assert resolve_tail_df("residual", residual_df=None, n_clusters=None) is None + + def test_invalid_convention_raises(self): + from diff_diff.utils import resolve_tail_df, validate_df_convention + + with pytest.raises(ValueError, match="df_convention must be one of"): + resolve_tail_df("bogus", residual_df=1, n_clusters=None) + with pytest.raises(ValueError, match="'bogus'"): + validate_df_convention("bogus") + for ok in ("residual", "cluster", "normal"): + validate_df_convention(ok) # no raise diff --git a/tests/test_v4_matrix.py b/tests/test_v4_matrix.py index fcac5de0f..e1e29d494 100644 --- a/tests/test_v4_matrix.py +++ b/tests/test_v4_matrix.py @@ -116,11 +116,12 @@ # (M-097..M-115) = 99, plus Phase 2b PR 1's two rows (M-117, M-122) = 101, # plus the ETWFE reference-period pair (M-123, M-124) = 103, plus the # comparison-support row (M-125) = 104, plus the clustered-CR1 K_reference -# convergence row (M-126) = 105. +# convergence row (M-126) = 105, plus the tail-df consolidation family +# (M-127 behavior + the four M-128..M-131 default-flips) = 110. # Ids are never reused and terminal rows are never # deleted, so the ledger only grows - raise the floor when rows are added; a # lower parse count means scanner/format drift or an illegal row deletion. -ROW_COUNT_FLOOR = 105 +ROW_COUNT_FLOOR = 110 # Committed snapshot of the shipped id set ("ids are never deleted or reused" # contract - a delete-one-add-one edit keeps the count above the floor but trips @@ -136,7 +137,10 @@ # the ETWFE reference-period fix (unidentified-cohort exclusion, and the # fail-closed guard on a design with no estimable post-treatment cell); # (126,126) = the clustered-CR1 K_reference convergence (D1+D2 of the 3.9 -# variance-consolidation program). +# variance-consolidation program); (127,131) = the tail-df consolidation +# (M-127 behavior row = PR C's D4 + normal-theory defect fix + the +# three-value knob extension, and M-128..M-131 = the four new standalone +# df-convention default-flips joining M-004..M-006 at 4.0). # M-116 and # M-118..M-121 are reserved for the later 2b PRs, not deleted - ids are # never reused, so a gap here is intentional. @@ -156,6 +160,7 @@ (122, 124), (125, 125), (126, 126), + (127, 131), ] EXPECTED_INITIAL_IDS = frozenset( f"M-{n:03d}" for lo, hi in _INITIAL_ID_RANGES for n in range(lo, hi + 1) @@ -552,13 +557,13 @@ def test_initial_ids_never_deleted(): """The shipped id set is immutable: ids are never deleted or reused (spec section 11). ROW_COUNT_FLOOR alone would let a delete-one-add-one edit pass; this snapshot cannot. - Extends as rows ship (105 as of the clustered-CR1 K_reference convergence row: + Extends as rows ship (110 as of the tail-df consolidation family: Phase 1 + diagnostic-family + M-092/M-093 + M-094..M-096 + the M-097..M-115 public-function completeness sweep + M-117/M-122 + M-123/M-124 + M-125 + - M-126).""" + M-126 + M-127..M-131).""" missing = sorted(EXPECTED_INITIAL_IDS - set(_ROW_IDS)) assert not missing, f"ledger rows deleted (ids are permanent): {missing}" - assert len(EXPECTED_INITIAL_IDS) == 105 + assert len(EXPECTED_INITIAL_IDS) == 110 def test_version_tuple_pads_to_three_components(): diff --git a/tests/test_variance_conventions.py b/tests/test_variance_conventions.py index 1b1199ab8..74e63450a 100644 --- a/tests/test_variance_conventions.py +++ b/tests/test_variance_conventions.py @@ -238,12 +238,14 @@ def snapshot(self): df, outcome="y", unit="unit", time="time", cohort="first_treat" ), cr1_k=(15,), - tail_df=(None,) * 10, - status="defect", + tail_df=(286.0,) * 10, + status="legitimate", reason=( - "CR1 k converged on K_reference (D2 fixed: 9 cells + T = 15, " - "no intercept col -> +1 term; jwdid arms at ratio 1.0); tail df " - "is still normal theory with no df_convention knob (PR C)" + "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'" ), ), dict( @@ -252,12 +254,14 @@ def snapshot(self): df, outcome="y", unit="unit", time="time", first_treat="first_treat" ), cr1_k=(21,), - tail_df=(280.0,) * 15 + (None,) * 8, - status="defect", + tail_df=(280.0,) * 23, + status="legitimate", reason=( - "CR1 k converged on K_reference (D2 fixed: 15 cells + 6, no " - "intercept col; fixest sunab parity ~5e-15); D4 remains: residual " - "df per cohort-period cell but normal theory on aggregates (PR C)" + "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)" ), ), dict( @@ -266,12 +270,14 @@ def snapshot(self): df, outcome="y", unit="unit", time="time", first_treat="first_treat" ), cr1_k=(6,), - tail_df=(None,) * 1, + tail_df=(309.0,) * 1, status="legitimate", reason=( "L1: k_total is clubSandwich CR1S by construction (stacked_did.py " - "pins vcovCR(type='CR1S') at atol=1e-10); normal-theory tail df is " - "an open PR C question" + "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)" ), ), dict( @@ -282,7 +288,11 @@ def snapshot(self): cr1_k=(4, 4, 5, 5, 5, 6), tail_df=(59.0,) * 6, status="legitimate", - reason="L2: G-1 tail df (Stata/fixest convention) — the convergence target", + reason=( + "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)" + ), ), dict( key="mpd_absorb_hc1_cluster_unit", @@ -376,15 +386,21 @@ def snapshot(self): aggregate="event_study", ), cr1_k=None, # contract row: shared CR1 REACHED; exact k is config-detail - tail_df=None, + # Pinned since 3.9 (M-127): the three lead calls carry the knob- + # resolved residual df (157 = n_untreated - k_kept - absorbed rank + # on the audit panel); the None entries are the knob-INDEPENDENT + # BJS overall/event-study calls (L3) captured by the same spy. + tail_df=(157.0,) * 3 + (None,) * 5, # The increment IS pinned even though k is not: every clustered-hc1 # call on this surface must carry the [time, unit] no-intercept # increment (1 + 65 - 60 = 6 on the audit panel). expected_adjustment=6, - status="defect", + status="legitimate", reason=( - "pretrends lead regression: CR1 k converged on K_reference " - "(D2 fixed); normal-theory tail df remains (PR C family)" + "pretrends lead regression: CR1 k = K_reference (D2 fixed); " + "tail df converged in 3.9 (M-127): leads use t(residual df) via " + "the df_convention knob (was silent normal theory); the None " + "calls are the knob-independent BJS aggregates (L3)" ), ), dict( @@ -1403,3 +1419,24 @@ def inject_bad_adjustment(*a, **k): diff_diff.LPDiD(pre_window=2, post_window=2).fit( df, outcome="y", unit="unit", time="time", treatment="treated" ) + + +def test_sun_abraham_normal_moves_cells_and_aggregates_together(monkeypatch): + """One semantic per fit: under df_convention="normal" the SA CELLS move to + z alongside the aggregates (t → z on the cells is a deliberate new option; + the cells' t/p/CI never reach public SA output — _fit_saturated_regression + consumes only coefficient/se — so the safe_inference spy is the only + observable surface for the cell convention).""" + df = make_panel() + cap = Capture(monkeypatch) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + diff_diff.SunAbraham(df_convention="normal").fit( + df, outcome="y", unit="unit", time="time", first_treat="first_treat" + ) + _, tail_df = cap.snapshot() + assert len(tail_df) == 23, tail_df + assert all(v is None for v in tail_df), ( + "every safe_inference call — cells included — must be normal-theory " + f"under 'normal'; captured {tail_df}" + ) diff --git a/tests/test_wooldridge.py b/tests/test_wooldridge.py index 573db3c3b..158bfd93c 100644 --- a/tests/test_wooldridge.py +++ b/tests/test_wooldridge.py @@ -4321,3 +4321,204 @@ def test_bootstrap_runs_on_a_filtered_fit(self): ) assert np.isfinite(res.overall_att) assert np.isfinite(res.overall_se) + + +class TestWooldridgeDfConvention: + """The three-value df_convention knob on WooldridgeDiD OLS arms (3.9 / M-127).""" + + @staticmethod + def _panel(seed=5): + rng = np.random.default_rng(seed) + rows = [] + for u in range(50): + coh = [0, 4, 5][u % 3] + for t in range(1, 7): + eff = 0.6 if (coh and t >= coh) else 0.0 + rows.append( + dict( + unit=u, + time=t, + cohort=coh, + outcome=0.2 * u / 50 + 0.1 * t + eff + rng.standard_normal() * 0.5, + ) + ) + return pd.DataFrame(rows) + + _kw = dict(outcome="outcome", unit="unit", time="time", cohort="cohort") + + @staticmethod + def _t_p(t_stat, df): + from scipy import stats + + return 2 * stats.t.sf(abs(t_stat), df) + + @staticmethod + def _z_p(t_stat): + from scipy import stats + + return 2 * stats.norm.sf(abs(t_stat)) + + def test_default_hc1_is_t_residual_not_z(self): + """The defect fix: the auto-clustered hc1 arm moved z -> t(residual).""" + res = WooldridgeDiD().fit(self._panel(), **self._kw) + assert res.overall_p_value != self._z_p(res.overall_t_stat) + # residual df = n - k_kept - absorbed [unit, time] rank; recover it + n = 300 + match = [ + d + for d in range(2, n) + if abs(res.overall_p_value - self._t_p(res.overall_t_stat, d)) < 1e-13 + ] + assert len(match) == 1, match + # per-cell inference shares the same df + _, eff = next(iter(res.group_time_effects.items())) + assert eff["p_value"] == pytest.approx(self._t_p(eff["t_stat"], match[0]), rel=1e-12) + + def test_cluster_matches_g_minus_1(self): + data = self._panel() + r0 = WooldridgeDiD().fit(data, **self._kw) + r1 = WooldridgeDiD(df_convention="cluster").fit(data, **self._kw) + G = data["unit"].nunique() + assert r1.overall_att == r0.overall_att and r1.overall_se == r0.overall_se + assert r1.overall_p_value == pytest.approx(self._t_p(r1.overall_t_stat, G - 1), rel=1e-12) + + def test_normal_reproduces_pre39_z(self): + data = self._panel() + r0 = WooldridgeDiD().fit(data, **self._kw) + rn = WooldridgeDiD(df_convention="normal").fit(data, **self._kw) + assert rn.overall_att == r0.overall_att and rn.overall_se == r0.overall_se + assert rn.overall_p_value == pytest.approx(self._z_p(rn.overall_t_stat), rel=1e-14) + + def test_classical_and_hc2_bit_identical_under_default(self): + """classical/hc2 already used t(n - rank(X)); the ONE resolved + fallback must reproduce those values exactly (n - k_kept on the + full-dummy design).""" + data = self._panel() + for vt in ("classical", "hc2"): + res = WooldridgeDiD(vcov_type=vt).fit(data, **self._kw) + _, eff = next(iter(res.group_time_effects.items())) + n = len(data) + match = [ + d + for d in range(2, n) + if abs(eff["p_value"] - self._t_p(eff["t_stat"], n - d)) < 1e-12 + ] + assert match, f"{vt} cell not on t(n-k)" + + def test_hc1_cohort_trends_uses_full_dummy_residual_df(self): + """The easy-to-miss lane: hc1 + cohort_trends routes FULL-DUMMY, so + its residual df is n - k_kept (no absorbed-rank subtraction).""" + data = self._panel() + res = WooldridgeDiD(cohort_trends=True).fit(data, **self._kw) + res_within = WooldridgeDiD().fit(data, **self._kw) + _, eff = next(iter(res.group_time_effects.items())) + n = len(data) + # full-dummy residual df: recover as n - k_kept + match_fd = [ + k + for k in range(1, 200) + if abs(eff["p_value"] - self._t_p(eff["t_stat"], n - k)) < 1e-12 + ] + assert match_fd, "cohort_trends hc1 cell not on t(n - k_kept)" + # and the within arm's df differs (it subtracts the absorbed rank on + # top of a much smaller visible k) + _, effw = next(iter(res_within.group_time_effects.items())) + match_w = [ + d for d in range(2, n) if abs(effw["p_value"] - self._t_p(effw["t_stat"], d)) < 1e-13 + ] + assert match_w and (n - match_fd[0]) != match_w[0] + + def test_postfit_aggregate_reproduces_fit_time_convention(self): + """aggregate() consumes the stored resolved fallback on every OLS arm.""" + res = WooldridgeDiD().fit(self._panel(), **self._kw) + n = 300 + match = [ + d + for d in range(2, n) + if abs(res.overall_p_value - self._t_p(res.overall_t_stat, d)) < 1e-13 + ] + agg = res.aggregate("group") + ge = next(iter(agg.group_effects.values())) + assert ge["p_value"] == pytest.approx(self._t_p(ge["t_stat"], match[0]), rel=1e-12) + simple = res.aggregate("simple") + assert simple.overall_p_value == pytest.approx( + self._t_p(simple.overall_t_stat, match[0]), rel=1e-12 + ) + + def test_bootstrap_aggregate_group_stays_analytical_t(self): + """On a bootstrap fit only the overall p/CI are percentile-overridden; + post-fit aggregate('group') re-derives ANALYTICAL inference from the + stored fallback (the _df_analytic_fallback-stays-live behavior).""" + res = WooldridgeDiD(n_bootstrap=20, seed=1).fit(self._panel(), **self._kw) + agg = res.aggregate("group") + ge = next(iter(agg.group_effects.values())) + assert np.isfinite(ge["p_value"]) + assert ge["p_value"] != self._z_p(ge["t_stat"]), "bootstrap fit reverted aggregates to z" + + def test_glm_warns_on_explicit_nondefault_only(self): + data = self._panel() + data["bin"] = (data["outcome"] > data["outcome"].median()).astype(int) + kw = dict(outcome="bin", unit="unit", time="time", cohort="cohort") + with pytest.warns(UserWarning, match="no effect on the logit/poisson"): + WooldridgeDiD(method="logit", df_convention="cluster").fit(data, **kw) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + WooldridgeDiD(method="logit").fit(data, **kw) + assert not any("no effect on the logit" in str(w.message) for w in caught) + + def test_pickle_setstate_migrates_df_one_way(self): + """A pre-3.9 pickle stores _df_one_way; __setstate__ carries it into + _df_analytic_fallback so a legacy result's aggregate() keeps its + fit-time classical/hc2 inference.""" + res = WooldridgeDiD(vcov_type="classical").fit(self._panel(), **self._kw) + legacy_df = res._df_analytic_fallback + state = dict(res.__dict__) + state.pop("_df_analytic_fallback") + state.pop("df_convention", None) + state["_df_one_way"] = legacy_df + revived = WooldridgeDiDResults.__new__(WooldridgeDiDResults) + revived.__setstate__(state) + assert revived._df_analytic_fallback == legacy_df + assert not hasattr(revived, "_df_one_way") or "_df_one_way" not in revived.__dict__ + assert revived.df_convention == "residual" + agg = revived.aggregate("simple") + assert np.isfinite(agg.overall_p_value) + + def test_survey_df_precedence_over_knob(self): + """Survey design df wins under EVERY knob value: the three fits are + bit-identical on overall AND per-cell inference (the knob is never + consulted on surveyed fits).""" + from diff_diff.survey import SurveyDesign + + data = self._panel() + data["w"] = 1.0 + (data["unit"] % 5) * 0.2 + data["stratum"] = data["unit"] % 4 + sd = SurveyDesign(weights="w", strata="stratum", psu="unit") + fits = { + conv: WooldridgeDiD(df_convention=conv).fit(data, survey_design=sd, **self._kw) + for conv in ("residual", "cluster", "normal") + } + base = fits["residual"] + # the survey df is finite -> t inference even under "normal" + assert base.overall_p_value != self._z_p(base.overall_t_stat) + for conv in ("cluster", "normal"): + r = fits[conv] + assert r.overall_p_value == base.overall_p_value + assert r.overall_conf_int == base.overall_conf_int + for k, eff in base.group_time_effects.items(): + assert r.group_time_effects[k]["p_value"] == eff["p_value"] + + def test_validation_and_transactional_set_params(self): + with pytest.raises(ValueError, match="df_convention"): + WooldridgeDiD(df_convention="bogus") + est = WooldridgeDiD() + with pytest.raises(ValueError, match="df_convention"): + est.set_params(df_convention="bogus", alpha=0.10) + assert est.df_convention == "residual" and est.alpha == 0.05 + # Valid value + unknown key: the unknown-key rejection must also + # leave the estimator fully unchanged (no partial application). + before = est.get_params() + with pytest.raises(ValueError, match="Unknown parameter"): + est.set_params(df_convention="normal", nonexistent_param=1) + assert est.get_params() == before + assert WooldridgeDiD(df_convention="cluster").get_params()["df_convention"] == "cluster" From ed5d822031ace237b0c9637734535d9c64ff70fc Mon Sep 17 00:00:00 2001 From: igerber Date: Fri, 31 Jul 2026 16:13:29 -0400 Subject: [PATCH 2/2] docs(3.9): three-value TWFE df_convention docstring + result-stub attribute entries The TwoWayFixedEffects constructor docstring still enumerated two knob values and claimed blanket unclustered/Conley inertness, which "normal" falsifies (it is deliberate z at the fallback level on every fit) - rewritten with the per-value scoping used by DifferenceInDifferences and LinearRegression. The committed docs/api/_autosummary attribute lists gain the new result fields: df_convention on the five containers, inference_df on SunAbrahamResults/StackedDiDResults (alphabetical insertion, the #715/#726 convention). --- diff_diff/twfe.py | 13 ++++++++----- .../_autosummary/diff_diff.ImputationDiDResults.rst | 1 + .../_autosummary/diff_diff.StackedDiDResults.rst | 2 ++ .../_autosummary/diff_diff.SunAbrahamResults.rst | 2 ++ .../diff_diff.lpdid_results.LPDiDResults.rst | 1 + ...diff.wooldridge_results.WooldridgeDiDResults.rst | 1 + 6 files changed, 15 insertions(+), 5 deletions(-) diff --git a/diff_diff/twfe.py b/diff_diff/twfe.py index fce0ba520..2b596ff76 100644 --- a/diff_diff/twfe.py +++ b/diff_diff/twfe.py @@ -60,12 +60,15 @@ class TwoWayFixedEffects(DifferenceInDifferences): auto-cluster is also preserved (routes to CR2-BM at unit). alpha : float, default=0.05 Significance level for confidence intervals. - df_convention : str, default "residual" + df_convention : {"residual", "cluster", "normal"}, default "residual" Inherited from :class:`DifferenceInDifferences`: df convention for - t/p/CI on clustered analytical fits — ``"residual"`` (fitted - residual df, default) or ``"cluster"`` (Stata/fixest ``G − 1``). - Survey df and per-coefficient Bell-McCaffrey DOF keep precedence; - inert on unclustered and Conley fits. Default flips at v4. + analytical t/p/CI. ``"residual"`` (default) uses the fitted + residual df; ``"cluster"`` uses the Stata/fixest ``G − 1`` on + clustered fits (no effect on unclustered or Conley fits); + ``"normal"`` deliberately uses normal-theory z inference at the + fallback level on every fit, clustered or not. Survey df and + per-coefficient Bell-McCaffrey DOF keep precedence under every + value. Default flips to ``"cluster"`` at v4. Notes ----- diff --git a/docs/api/_autosummary/diff_diff.ImputationDiDResults.rst b/docs/api/_autosummary/diff_diff.ImputationDiDResults.rst index f2cf18a17..0ab48283a 100644 --- a/docs/api/_autosummary/diff_diff.ImputationDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.ImputationDiDResults.rst @@ -32,6 +32,7 @@ ~ImputationDiDResults.cluster_name ~ImputationDiDResults.coef_var ~ImputationDiDResults.conf_int + ~ImputationDiDResults.df_convention ~ImputationDiDResults.is_significant ~ImputationDiDResults.leave_one_out ~ImputationDiDResults.n_clusters diff --git a/docs/api/_autosummary/diff_diff.StackedDiDResults.rst b/docs/api/_autosummary/diff_diff.StackedDiDResults.rst index 33089fa42..cece61be1 100644 --- a/docs/api/_autosummary/diff_diff.StackedDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.StackedDiDResults.rst @@ -34,9 +34,11 @@ ~StackedDiDResults.coef_var ~StackedDiDResults.conf_int ~StackedDiDResults.covariates + ~StackedDiDResults.df_convention ~StackedDiDResults.event_study_df ~StackedDiDResults.event_study_vcov ~StackedDiDResults.event_study_vcov_index + ~StackedDiDResults.inference_df ~StackedDiDResults.is_significant ~StackedDiDResults.kappa_post ~StackedDiDResults.kappa_pre diff --git a/docs/api/_autosummary/diff_diff.SunAbrahamResults.rst b/docs/api/_autosummary/diff_diff.SunAbrahamResults.rst index 5d03aa9f4..66b12d8ef 100644 --- a/docs/api/_autosummary/diff_diff.SunAbrahamResults.rst +++ b/docs/api/_autosummary/diff_diff.SunAbrahamResults.rst @@ -34,9 +34,11 @@ ~SunAbrahamResults.conf_int ~SunAbrahamResults.conley_lag_cutoff ~SunAbrahamResults.control_group + ~SunAbrahamResults.df_convention ~SunAbrahamResults.event_study_df ~SunAbrahamResults.event_study_vcov ~SunAbrahamResults.event_study_vcov_index + ~SunAbrahamResults.inference_df ~SunAbrahamResults.is_significant ~SunAbrahamResults.p_value ~SunAbrahamResults.reference_observed diff --git a/docs/api/_autosummary/diff_diff.lpdid_results.LPDiDResults.rst b/docs/api/_autosummary/diff_diff.lpdid_results.LPDiDResults.rst index ae4bb6ab3..d8c1dbd0a 100644 --- a/docs/api/_autosummary/diff_diff.lpdid_results.LPDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.lpdid_results.LPDiDResults.rst @@ -30,6 +30,7 @@ ~LPDiDResults.cluster_name ~LPDiDResults.conf_int ~LPDiDResults.covariates + ~LPDiDResults.df_convention ~LPDiDResults.dylags ~LPDiDResults.estimand ~LPDiDResults.event_study_df diff --git a/docs/api/_autosummary/diff_diff.wooldridge_results.WooldridgeDiDResults.rst b/docs/api/_autosummary/diff_diff.wooldridge_results.WooldridgeDiDResults.rst index 144daa655..4cc51a8af 100644 --- a/docs/api/_autosummary/diff_diff.wooldridge_results.WooldridgeDiDResults.rst +++ b/docs/api/_autosummary/diff_diff.wooldridge_results.WooldridgeDiDResults.rst @@ -34,6 +34,7 @@ ~WooldridgeDiDResults.conf_int ~WooldridgeDiDResults.conley_lag_cutoff ~WooldridgeDiDResults.control_group + ~WooldridgeDiDResults.df_convention ~WooldridgeDiDResults.event_study_effects ~WooldridgeDiDResults.group_effects ~WooldridgeDiDResults.method