From 993b2d0a3e19db61a0ed48c465dbef6769d6b928 Mon Sep 17 00:00:00 2001 From: igerber Date: Thu, 30 Jul 2026 08:14:24 -0400 Subject: [PATCH] feat(3.9): variance-conventions inventory + component-aware absorbed-FE rank (D3) The measured map behind the 3.9 variance-consolidation program, plus the first of its fixes. Inventory (docs/methodology/variance-conventions.md, repo-internal): - Measured matrix of every surface reaching the shared clustered-CR1 denominator (linalg._compute_robust_vcov_numpy, vcov_type="hc1" only - a clustered call in any other family fails its row) and every tail-df convention passed to safe_inference/safe_inference_batch. The table is generated from expected literals in tests/test_variance_conventions.py; the parametrized tests assert those literals against live instrumentation (fast subset by default, full sweep under -m slow) and a byte-equality test keeps the committed doc in sync. - Classifies the defects scheduled for the consolidation PRs - D1 (absorb= vs fixed_effects=: 10.35% SE split on the same model, ratio exactly sqrt((360-2)/(360-66))), D2 (clustered CR1 never counts absorbed FE not nested in the cluster; anti-conservative, up to 5.51% at G=60), D4 (SunAbraham reports residual df per cohort-period cell but normal theory on aggregates) - and the legitimate exceptions L1-L4 (StackedDiD CR1S, LPDiD G-1, influence-function paths, hc2/hc2_bm + survey TSL) with reasons and external anchors. D3 fix (this PR's behavior change): - New diff_diff.utils.absorbed_fe_rank: two-way absorbed-FE df from the connected components of the bipartite level graph (sum(levels) - C, minus 1 when the visible design carries an intercept column), exact per Abowd-Creecy-Kramarz. Levels and edges use positive-weight rows only, per the REGISTRY zero-weight-padding guarantee; NaN group keys raise an actionable ValueError. N>=3 keeps sum(levels-1) with the over-count limitation documented and tracked in TODO.md. - Consumers: DiD(absorb=) and MultiPeriodDiD(absorb=) (intercept form, computed on the pre-transform frame), within-transform TwoWayFixedEffects (intercept form), SunAbraham (no-intercept form - its saturated design has coef_offset=0). Connected independent panels are bit-identical to the old count; disconnected and hierarchical designs (absorb=["state","state_year"]: true rank 29, old count 34) get corrected residual df and non-clustered classical/HC1 SE scaling, with the fail-closed NaN boundary moving consistently in both directions. - External anchor: matches fixest::ssc(K.exact = TRUE) at <=1e-12 on a committed hierarchical R golden (benchmarks/data/fixest_kexact_golden.json, fixest 0.14.2, generator script included); fixest's DEFAULT K.exact=FALSE reproduces the old approximate count, recorded as a labeled deviation from the R default in REGISTRY.md. - demean_by_groups/demean_by_group's second return value is now documented as the raw level count, NOT a valid df adjustment, pointing to absorbed_fe_rank (no callers consume it as df any more; returning the rank from the demeaner instead was rejected because within_transform delegates there and would compute the graph twice per TWFE/SA fit). - REGISTRY.md: the TwoWayFixedEffects absorbed-FE df note is now component-aware, with cross-references from the DiD, MultiPeriodDiD, and SunAbraham sections; both sides of the D1 split and the Wooldridge hc1 mechanism status are documented. doc-deps.yaml registers the inventory doc against its 11 source modules; performance-plan.md records the helper's measured 1.9 ms / 186k-row cost with the codes-reuse optimization tracked as its own TODO row. --- CHANGELOG.md | 23 + DEFERRED.md | 4 +- TODO.md | 4 +- benchmarks/R/generate_fixest_kexact_golden.R | 59 ++ benchmarks/data/fixest_kexact_golden.json | 24 + diff_diff/estimators.py | 27 +- diff_diff/sun_abraham.py | 12 +- diff_diff/twfe.py | 12 +- diff_diff/utils.py | 122 ++- docs/conf.py | 1 + docs/doc-deps.yaml | 33 + docs/methodology/REGISTRY.md | 73 +- docs/methodology/variance-conventions.md | 118 +++ docs/performance-plan.md | 17 + tests/test_utils.py | 24 +- tests/test_variance_conventions.py | 812 +++++++++++++++++++ 16 files changed, 1348 insertions(+), 17 deletions(-) create mode 100644 benchmarks/R/generate_fixest_kexact_golden.R create mode 100644 benchmarks/data/fixest_kexact_golden.json create mode 100644 docs/methodology/variance-conventions.md create mode 100644 tests/test_variance_conventions.py diff --git a/CHANGELOG.md b/CHANGELOG.md index d7d678f59..8344f0c23 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,29 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Fixed +- **Absorbed-FE degrees of freedom over-counted on disconnected and hierarchical + panels.** The df adjustment for absorbed fixed effects (`TwoWayFixedEffects`, + `SunAbraham`, `DifferenceInDifferences(absorb=)`, `MultiPeriodDiD(absorb=)`) + used `sum_d (levels_d − 1)`, which assumes the FE dimensions are mutually + independent and the unit×time incidence graph is connected. On disconnected + panels, and on nested specs like `absorb=["state", "state_year"]` (measured: + 6 states × 5 years → true absorbed rank 29, old count 34), it over-stated the + absorbed rank, making the reported residual df too small and the non-clustered + classical/hc1 SEs too large. The new `diff_diff.utils.absorbed_fe_rank` + computes the two-way rank from the connected components of the bipartite level + graph; connected panels are bit-identical (`rtol=0, atol=1e-14` across all ten + estimator surfaces). Levels and connectivity are now evaluated over + positive-weight rows only, restoring the documented zero-weight-padding + inference-invariance on weighted fits with inert rows. +- **New CI-enforced variance-convention audit matrix** + (`tests/test_variance_conventions.py` + + `docs/methodology/variance-conventions.md`): pins, per estimator surface, the + visible `k` reaching the shared clustered CR1 denominator and the tail-df + convention passed to inference, with each cell classified as a documented + defect (scheduled for the 3.9 variance-consolidation program) or a declared + legitimate exception with its reason. This is the map for the follow-up PRs + that converge the clustered CR1 `k` on the reghdfe/fixest nested convention + and the tail df on `t(G−1)`. - **`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 diff --git a/DEFERRED.md b/DEFERRED.md index 19827426c..10d99e00c 100644 --- a/DEFERRED.md +++ b/DEFERRED.md @@ -50,8 +50,8 @@ exists but parity can't be verified without a local toolchain. | Extend `WooldridgeDiD` `method ∈ {logit, poisson}` with `vcov_type ∈ {classical, hc2, hc2_bm}`: composing HC2 leverage + Bell-McCaffrey DOF with the QMLE pseudo-residual sandwich needs derivation + R parity vs `clubSandwich::vcovCR(glm, type="CR2")`. Rejected at `__init__`. | `wooldridge.py` | follow-up | Medium | | Multi-constraint CR2 parallel-trends test (AHT/HTZ) for `hc2_bm` fits: DiagnosticReport's PT check routes `vcov_type="hc2_bm"` sources to Bonferroni over the BM-adjusted per-row p-values because the generic chi-square joint Wald would discard the CR2 small-sample correction (see REPORTING.md "hc2_bm parallel-trends policy"). The proper joint test is the AHT/HTZ Wald with a Satterthwaite-style denominator df over the pre-period contrast block; needs derivation for the stacked/pooled WLS-CR2 layout + parity vs `clubSandwich::Wald_test(..., test="HTZ")`. | `diagnostic_report.py`, `linalg.py` | vcov/df round-trip PR | Low | | `PreTrendsPower` CS/SA `anticipation=1` R-parity fixture: R `pretrends` has no anticipation parameter, so the Python `_extract_pre_period_params` anticipation filter isn't R-parity-locked. Build a synthetic CS/SA result with `anticipation=1` and assert γ_p matches R's `slope_for_power()`. (Mechanism already covered by MC + full-VCV tests.) | `tests/test_methodology_pretrends.py`, `generate_pretrends_golden.R` | PR-C | Low | -| Harmonize SunAbraham's HC1 within-transform finite-sample correction with `fixest::sunab()` — SA applies `n/(n-k_dm)`, fixest applies `n/(n-k_total)` (counts absorbed FE); ~1-2% SE difference, documented as a "Deviation from R" and pinned at `atol=5e-3`. Either thread `df_adjustment` or keep as an intentional, R-verified difference. | `sun_abraham.py`, `linalg.py` | follow-up | Low | -| Absorbed-FE **clustered** CR1 with *non-nested* FE: for `absorb=[FE1,FE2], cluster=FE1` (e.g. `absorb=["unit","time"], cluster="unit"`), `fixest` counts the non-nested FE (time) in the CR1 `(n-1)/(n-k)` finite-sample denominator, but the clustered path uses only `k_visible`. D4 harmonized the *non-clustered* classical/hc1 full-K scale (`_absorbed_fe_vcov_scale`) and left the clustered path unchanged — correct for FE nested in the cluster, a small deviation for non-nested FE (documented in REGISTRY within-transform note). Thread a non-nested `df_adjustment` into the clustered CR1 factor; verify vs `fixest::feols(..., cluster=)`. | `linalg.py`, `estimators.py` | SE-audit D4 | Low | +| Harmonize SunAbraham's HC1 within-transform finite-sample correction with `fixest::sunab()` — SA applies `n/(n-k_dm)`, fixest applies `n/(n-k_total)` (counts absorbed FE); ~1-2% SE difference, documented as a "Deviation from R" and pinned at `atol=5e-3`. Either thread `df_adjustment` or keep as an intentional, R-verified difference. Mechanism + measured inventory now in `docs/methodology/variance-conventions.md` (defect D2 family; 3.9 consolidation PR B). | `sun_abraham.py`, `linalg.py` | follow-up | Low | +| Absorbed-FE **clustered** CR1 with *non-nested* FE: for `absorb=[FE1,FE2], cluster=FE1` (e.g. `absorb=["unit","time"], cluster="unit"`), `fixest` counts the non-nested FE (time) in the CR1 `(n-1)/(n-k)` finite-sample denominator, but the clustered path uses only `k_visible`. D4 harmonized the *non-clustered* classical/hc1 full-K scale (`_absorbed_fe_vcov_scale`) and left the clustered path unchanged — correct for FE nested in the cluster, a small deviation for non-nested FE (documented in REGISTRY within-transform note). Thread a non-nested `df_adjustment` into the clustered CR1 factor; verify vs `fixest::feols(..., cluster=)`. Closed form now derived and externally verified (Stata ~1e-15 / fixest ~1e-12) — see `docs/methodology/variance-conventions.md` defect D2; scheduled as PR B of the 3.9 consolidation program. | `linalg.py`, `estimators.py` | SE-audit D4 | Low | | Rust multiplier-bootstrap weight RNG (`generate_bootstrap_weights_batch`) seeds `Xoshiro256PlusPlus::seed_from_u64(seed+i)` per row; audit Python callers (`sdid.py`, `efficient_did_bootstrap.py`, `bootstrap_utils.py`) for parity-test gaps and, where a numpy-canonical equivalent exists, pre-generate in Python and pass through PyO3 (same fix shape as TROP RNG parity #354). | `rust/src/bootstrap.rs`, `bootstrap_utils.py` | follow-up | Medium | | `SyntheticDiD` bootstrap cross-language parity anchor vs R `synthdid::vcov(method="bootstrap")` or Julia `Synthdid.jl` (refit-native). Same-library validation is in place; Julia is the cleanest target. Tolerance ~1e-6 (BLAS+RNG paths preclude 1e-10). | `benchmarks/R/`, `benchmarks/julia/`, `tests/` | follow-up | Low | | CS R helpers hard-code `xformla = ~1`; no covariate-adjusted R benchmark for the IRLS path. | `tests/test_methodology_callaway.py` | #202 | Low | diff --git a/TODO.md b/TODO.md index a472e2401..77421c89f 100644 --- a/TODO.md +++ b/TODO.md @@ -21,9 +21,10 @@ Related tracking surfaces: | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| +| `absorbed_fe_rank` N>=3 general rank: the helper keeps `sum(levels-1)` for 3+ absorbed dims, exact for independent connected dims but an over-count for duplicated/nested triples (measured `a(5),b(4),c==b(4)`: true 7 vs formula 10) and for disconnected N-way graphs. Two-way is component-exact. Deriving general N-way FE rank is a hypergraph problem; do it with a reference (fixest's `fixef.rm`/reghdfe df_a) rather than a guess. See `docs/methodology/variance-conventions.md` D3. | `diff_diff/utils.py` | #variance-inventory | Mid | Low | | `SyntheticControl` conformal (CWZ 2021) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies — the residual-permutation shortcut is only valid for time-permutation-invariant proxies (SC/Lasso/DiD); an AR proxy needs innovation permutation. | `diff_diff/conformal.py`, `diff_diff/synthetic_control_results.py` | CWZ-2021 | Heavy | Low | | Make the post-fit `results.aggregate("event_study")` container consumable downstream. `EventStudyResults` is rejected by all THREE consumers that read a CS event study — `compute_honest_did` (`honest_did.py`, dispatches on `CallawaySantAnnaResults` and raises `TypeError`), `compute_pretrends_power` (`pretrends.py`, same), and `plot_event_study` (`visualization`, same) — so `fit(aggregate="event_study")` is still the only route for them and their error messages say so explicitly. Needs an `EventStudyResults` branch in each extraction path (consuming `event_time` / `is_reference` / `vcov` / `vcov_index` / per-row `df`) PLUS `base_period` and `anticipation` provenance, which the unified container does not carry and HonestDiD needs for its universal-base-period warning and pre-period classification. Gate with end-to-end tests: `compute_honest_did(res.aggregate("event_study"))` at `base_period="universal"`, and `compute_pretrends_power(...)` at `anticipation=1`. | `diff_diff/honest_did.py`, `diff_diff/pretrends.py`, `diff_diff/results_base.py` | #726 | Mid | Medium | -| Derive and fix the `WooldridgeDiD` `hc1` SE gap vs Stata `jwdid`. Measured: every SE is uniformly SMALLER than `jwdid`'s, by 1.0280 at G=20 / 1.0132 at G=40 / 1.00264 at G=191 / 1.0010 at G=500 (ATT(g,t) points match exactly). The gap tracks `sqrt(G/(G-1))` but sits consistently above it, and `solve_ols` already applies the full CR1 `(G/(G-1))*((n-1)/(n-k))`, so a missing cluster factor is ruled out -- the likely source is the within-transform `k` accounting vs `hdfe`/`reghdfe`'s. Derive the exact factor FIRST; a `sqrt(G/(G-1))` patch would match only approximately and would force a loose tolerance on the parity test, defeating its purpose. Move the REGISTRY note and the pinned ratio in `tests/test_etwfe_cs_stata_parity.py` together with the fix. **Required artifact:** a committed subsample LADDER in the Stata golden (rosters = first N units per `first_treat` by sorted `countyreal`; rungs spanning G≈20..500, storing Stata `G`/`n`/`df_a`/`rank`/`df_r` and per-cell `att`/`se`), with parameterized ratio assertions. Only G=500 (full panel) and G=191 (all-eventually-treated arm) are pinned today, so the few-cluster behavior — where the gap is materially largest (~2.8%) — is ungated, and the ladder is also the instrument for comparing Stata's `df_a`/`rank` against the library's within-transform `k`, which is the leading hypothesis for the factor. | `diff_diff/wooldridge.py`, `diff_diff/linalg.py` | #723-followup | Mid | Medium | +| Derive and fix the `WooldridgeDiD` `hc1` SE gap vs Stata `jwdid`. **MECHANISM NOW DERIVED** (see `docs/methodology/variance-conventions.md`, defect D2): the clustered CR1 denominator uses `k_visible` and never counts absorbed FE not nested in the cluster; the closed form `K_reference = explicit cols + (1 if no intercept col) + rank(non-nested FE | nested)` reproduces jwdid/reghdfe to ~1e-15 on three arms and fixest to ~1e-12 on two, retrodicting the G=20/G=40 rungs. Remaining work (PR B of the 3.9 consolidation program) is the threading + Rust port + Stata ladder, not derivation. Measured: every SE is uniformly SMALLER than `jwdid`'s, by 1.0280 at G=20 / 1.0132 at G=40 / 1.00264 at G=191 / 1.0010 at G=500 (ATT(g,t) points match exactly). The gap tracks `sqrt(G/(G-1))` but sits consistently above it, and `solve_ols` already applies the full CR1 `(G/(G-1))*((n-1)/(n-k))`, so a missing cluster factor is ruled out -- the likely source is the within-transform `k` accounting vs `hdfe`/`reghdfe`'s. Derive the exact factor FIRST; a `sqrt(G/(G-1))` patch would match only approximately and would force a loose tolerance on the parity test, defeating its purpose. Move the REGISTRY note and the pinned ratio in `tests/test_etwfe_cs_stata_parity.py` together with the fix. **Required artifact:** a committed subsample LADDER in the Stata golden (rosters = first N units per `first_treat` by sorted `countyreal`; rungs spanning G≈20..500, storing Stata `G`/`n`/`df_a`/`rank`/`df_r` and per-cell `att`/`se`), with parameterized ratio assertions. Only G=500 (full panel) and G=191 (all-eventually-treated arm) are pinned today, so the few-cluster behavior — where the gap is materially largest (~2.8%) — is ungated, and the ladder is also the instrument for comparing Stata's `df_a`/`rank` against the library's within-transform `k`, which is the leading hypothesis for the factor. | `diff_diff/wooldridge.py`, `diff_diff/linalg.py` | #723-followup | Mid | Medium | | `SunAbraham`: a cohort not observed at its own reference relative period (`e = -1 - anticipation`) makes that cohort's block collinear, so QR drops an unnamed column (`dropping 1 of 12 columns (column 9)`) and `overall_att` comes back **NaN**. Found by auditing the sibling estimator while fixing the ETWFE analogue (#724); PRE-EXISTING, not introduced there. Lower severity than #724 — that returned a silently WRONG finite number, this returns NaN with a rank warning — but the event-study surface still looks complete, so a user may not notice the loss. SA already omits its reference explicitly and tracks `_reference_observed`, so the fix is per-cohort support for that flag rather than the ETWFE-style redesign. | `diff_diff/sun_abraham.py` | #724-audit | Mid | Low | | Define `N_g` (W2025 Eqs. 7.4/7.6) for UNBALANCED panels where comparison-support filtering removes every observation of some units in an estimated cohort, then replace the fail-closed guard with the defined behavior. `_n_g_per_cohort` is read off the final sample, so those units vanish from the cohort-share weights; measured on a cohort supplied with 100 units of which 90 appear only at a dropped period, `aggregate(weights="cohort_share")` moves 1.8078 -> 3.8157. The paper assumes a balanced panel and does not say whether `N_g` counts the supplied cohort or the surviving units, and the two disagree materially, so `aggregate` currently raises naming the cohorts and counts ([M-125]); `weights="cell"` is unaffected and balanced panels never trip it. Settle the estimand (likely: count the supplied cohort, since ATT(g,t) is a cohort-level quantity, but that weights units with no retained observation) and gate with a test computing Eq. 7.4 by hand on unequal cohort sizes. | `diff_diff/wooldridge_results.py`, `diff_diff/wooldridge.py` | #729-followup | Mid | Medium | | `WooldridgeDiD` + `survey_design=` does not support DOMAIN ESTIMATION, so BOTH row-deleting paths are currently REFUSED (`NotImplementedError`, all three methods) rather than performed: unidentified-cohort exclusion ([M-123]) and comparison-support period filtering ([M-125]). One fix unblocks both. Implementing it properly means zero-padding the excluded rows' weights while retaining strata/PSU/FPC, per REGISTRY *Subpopulation Analysis (Phase 6)* / Lumley (2004) 3.4, so TSL variance and `df_survey = n_PSU - n_strata` use the full design (naive deletion measured 22 -> 14 on a two-stratum panel). `SurveyDesign.subpopulation()` already implements the contract and SpilloverDiD Wave E.3 is the in-repo precedent; the blocker is that the weighted within-transform rejects zero-weight units, shared machinery behind 7 estimators. Landing it would turn both refusals back into supported fits. Gate with a `SurveyDesign.subpopulation()` parity test on ATT, TSL SE and survey df where the excluded cohort exhausts a PSU. | `diff_diff/wooldridge.py`, `diff_diff/utils.py` | #724-codex-R4/R5 | Heavy | Medium | @@ -48,6 +49,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| +| Reuse the demeaner's factorized codes in `absorbed_fe_rank` instead of re-factorizing: the helper adds ~1.9 ms per absorbed fit at 186k rows (7.7% of the fastest Rust-served TWFE fit; see `docs/performance-plan.md` "Component-aware absorbed-FE rank"), and both the helper and `demean_by_groups` factorize the same group columns. Threading the codes through the four call sites halves the factorize work; the `connected_components` call itself is ~1.1 ms. Deliberately not done in the correctness PR. | `diff_diff/utils.py` | #variance-inventory | Quick | Low | | `EfficientDiD` conditional path: the largest remaining O(n) stage is the sieve/nuisance construction outside the tiled pass (~9s at 10k). (The `_ridge_solve_weights` Python-prep shave landed 2026-07-07 — the `omega_stack[rest]` fancy-index copy and tail scatter are skipped when no row is zero-masked, byte-identical outputs; the `zero_mask` abs scan itself remains, needed for correctness.) | `efficient_did_covariates.py` | CS-scaling | Mid | Low | | `_rq_fit` LP assembly is dense (`A_eq = [X, I, -I]` with dense identity blocks, rebuilt per cell fit): a `scipy.sparse` construction would cut memory and likely HiGHS time for large cells / bootstrap-heavy covariate CiC/QDiD fits. CAVEAT before doing it: a different matrix representation can change HiGHS's vertex selection at degenerate/tied QR optima - end-to-end covariate goldens are tie-selection-gated (fine), but the `qr_cases` tight coefficient matches may shift to the equal-loss branch; re-run the parity suite and re-calibrate if needed. | `diff_diff/changes_in_changes.py::_rq_fit` | covariates PR | Quick | Low | | Evaluate flipping `DIFF_DIFF_SOLVE_OLS_FASTPATH` default-ON after an opt-in soak (the 2026-07 certified normal-equations Cholesky fast path, both backends). A flip needs: golden/parity-suite recapture at the tol-bounded posture (fitted ~1e-8 abs / SE ~1e-6 rel — the default today is byte-pinned in several benchmark conventions), certification-rate telemetry across real workloads (any decline is silent-correct but forfeits the speedup), and the staged default-flip protocol used for `df_convention` (v4-class change). Lifecycle tracked in docs/v4-deprecations.yaml (M-008). | `diff_diff/linalg.py::_resolve_solve_ols_fastpath`, `rust/src/linalg.rs::solve_ols_chol` | CS-scaling | Mid | Low | diff --git a/benchmarks/R/generate_fixest_kexact_golden.R b/benchmarks/R/generate_fixest_kexact_golden.R new file mode 100644 index 000000000..1bd89905d --- /dev/null +++ b/benchmarks/R/generate_fixest_kexact_golden.R @@ -0,0 +1,59 @@ +# Golden: fixest exact-vs-default FE counting on a hierarchical two-way design. +# +# absorb = [state, state_year] with state_year nested in state splits the +# bipartite level graph into one component per state (C = 6), so the absorbed +# dummy-space rank is 29 beyond the intercept -- not the naive +# sum(levels - 1) = 34. fixest's default ssc(K.exact = FALSE) uses the naive +# count (df.K = 36 here); ssc(K.exact = TRUE) computes the exact rank +# (df.K = 31). diff-diff's component-aware absorbed_fe_rank matches the EXACT +# side at machine precision (a documented deviation from the R *default*). +# +# Regenerate: Rscript benchmarks/R/generate_fixest_kexact_golden.R +suppressMessages(library(fixest)) +suppressMessages(library(jsonlite)) + +set.seed(7) +d <- expand.grid(s = 0:5, y = 0:4, r = 1:4) +d$state <- d$s +d$state_year <- d$s * 100 + d$y +d$x <- rnorm(nrow(d)) +d$out <- rnorm(nrow(d)) + 0.5 * d$x + 0.3 * d$s + +m_default <- feols(out ~ x | state + state_year, data = d, vcov = "iid") +m_exact <- feols(out ~ x | state + state_year, data = d, vcov = "iid", + ssc = ssc(K.exact = TRUE)) + +golden <- list( + meta = list( + generator = "benchmarks/R/generate_fixest_kexact_golden.R", + r_version = paste(R.version$major, R.version$minor, sep = "."), + fixest_version = as.character(packageVersion("fixest")), + description = paste( + "Hierarchical two-way FE (state_year nested in state, C=6):", + "fixest default ssc(K.exact=FALSE) vs exact FE-rank counting.", + "diff-diff absorbed_fe_rank matches the K.exact=TRUE side." + ) + ), + data = list( + state = d$state, + state_year = d$state_year, + x = d$x, + out = d$out + ), + n_obs = nrow(d), + coef = unname(coef(m_default)[["x"]]), + iid_default = list( + se = unname(se(m_default)[["x"]]), + df_k = degrees_freedom(m_default, "k") + ), + iid_k_exact = list( + se = unname(se(m_exact)[["x"]]), + df_k = degrees_freedom(m_exact, "k") + ) +) + +path <- "benchmarks/data/fixest_kexact_golden.json" +write_json(golden, path, digits = NA, auto_unbox = TRUE, pretty = TRUE) +cat("wrote", path, "\n") +cat(sprintf("coef=%.15f default_se=%.15f exact_se=%.15f\n", + golden$coef, golden$iid_default$se, golden$iid_k_exact$se)) diff --git a/benchmarks/data/fixest_kexact_golden.json b/benchmarks/data/fixest_kexact_golden.json new file mode 100644 index 000000000..5fdf223ab --- /dev/null +++ b/benchmarks/data/fixest_kexact_golden.json @@ -0,0 +1,24 @@ +{ + "meta": { + "generator": "benchmarks/R/generate_fixest_kexact_golden.R", + "r_version": "4.5.2", + "fixest_version": "0.14.2", + "description": "Hierarchical two-way FE (state_year nested in state, C=6): fixest default ssc(K.exact=FALSE) vs exact FE-rank counting. diff-diff absorbed_fe_rank matches the K.exact=TRUE side." + }, + "data": { + "state": [0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5], + "state_year": [0, 100, 200, 300, 400, 500, 1, 101, 201, 301, 401, 501, 2, 102, 202, 302, 402, 502, 3, 103, 203, 303, 403, 503, 4, 104, 204, 304, 404, 504, 0, 100, 200, 300, 400, 500, 1, 101, 201, 301, 401, 501, 2, 102, 202, 302, 402, 502, 3, 103, 203, 303, 403, 503, 4, 104, 204, 304, 404, 504, 0, 100, 200, 300, 400, 500, 1, 101, 201, 301, 401, 501, 2, 102, 202, 302, 402, 502, 3, 103, 203, 303, 403, 503, 4, 104, 204, 304, 404, 504, 0, 100, 200, 300, 400, 500, 1, 101, 201, 301, 401, 501, 2, 102, 202, 302, 402, 502, 3, 103, 203, 303, 403, 503, 4, 104, 204, 304, 404, 504], + "x": [2.28724716134052, -1.19677168222235, -0.694292510435459, -0.412292951136803, -0.970673341119483, -0.947279945228107, 0.748139340290551, -0.116955225887152, 0.152657626282234, 2.18997810732938, 0.356986230329022, 2.71675178313072, 2.28145192598956, 0.324020540138516, 1.89606706680993, 0.467680511321698, -0.893800723085444, -0.307328299537195, -0.00482242226757041, 0.988164149499945, 0.839750359624071, 0.7053418309055, 1.30596472081169, -1.38799621659285, 1.27291686425524, 0.184192771235767, 0.752279895740033, 0.591745052462727, -0.983052595771021, -0.276063955112006, -0.870851022568591, 0.718710553084245, 0.110652877769336, -0.0784667679717042, -0.420490459341998, -0.562125876285266, 0.997513444755305, -1.10513005881326, -0.142287830774585, 0.314994904887913, 1.21855053450735, -0.699317078685514, -0.285432751528726, -1.31155267260939, -0.391012431449258, -0.401526613094972, 1.35051758092295, 0.591190027089221, 0.100525455628569, 0.931071995520097, -0.262742348566532, -0.00766810471266396, 0.367153006545634, 1.70716254513761, 0.723740262528379, 0.481036048707917, -1.56786824422525, 0.318250283480828, 0.16599145067735, -0.899907629628172, 0.0763714738605741, 0.159155278262728, 0.543674184709699, 0.704807352613047, 0.31896914255711, 1.10924978897106, 0.769154194657112, 1.15347367477894, 1.26068350268094, 0.700623506572325, 0.432627160845955, -0.922601718256921, -0.615584206630919, -0.866659688251375, -1.63951708718114, -1.32583924384341, -0.88903672763797, -0.557602330302113, -0.0624023088383481, 2.42269297715943, 0.342585350147299, 0.00424823625213038, 0.0292198420016504, -0.393423429121302, -0.792704562788928, -0.311701865295172, -0.34606859178086, -0.304607588245324, -1.78589348744452, 0.58727467185797, 1.63579443444659, -0.645423473634397, 0.61899216878734, 0.236393598401322, 0.846500898751643, -0.573645738849693, 1.11799320399617, -1.54000113193302, -0.438123899300085, -0.150672970896448, 0.519058365272953, 0.587539704901891, -0.0793330610552102, -1.17436101486938, 0.308722117638371, -1.60387854268349, 0.991289625052712, 1.02322044470037, 0.840145438883483, 0.12007860795572, -0.426255055445707, 0.458926243504027, 0.645047947693915, 0.611530549290199, -0.889211293868794, 1.54389234869222, -1.24176360488504, 1.10344734039128, 0.982772356675575, 0.304327174033201], + "out": [-0.411094635525043, 1.27150493528914, 0.94130343806635, 0.516249865480014, 1.44386459536935, 2.55961095883309, 0.880648118150659, 0.274850063319732, -0.791221901996461, 3.01414676959371, 0.78515378194107, 3.67078578152616, 2.00689379330002, 0.830318573310043, 2.68285263202372, 0.376740975592928, 1.19832020735595, 2.26226879408419, 0.267347531284372, 1.80162252363394, -0.441558085110783, 0.377885304713433, 2.0182112987635, 1.01320891845771, 1.11882032292138, 0.320438010809833, 0.0872115224938079, 1.68181609644557, 1.04897045935062, 0.26404382845742, -0.791119954062928, 1.75665568726567, -0.251365531107804, 0.653309991876617, 1.66861620437602, 0.421149320632615, -1.09278221559987, 0.927783095111009, 1.75142537445086, 1.04658826588411, 2.15479773300139, 1.24453061870585, -0.135923007890613, 0.386846541794514, 1.44668032043699, 0.379899555222139, 2.19790218668496, 2.48064799237163, 0.370505502159387, -1.1494020262643, -1.87133407368194, 1.37903024122261, 2.56210655065419, 1.06080423906594, 0.97764321650855, 0.796560105240529, 0.543904880277871, 2.38323606027677, 1.43015421551732, 0.451506792431283, 2.23071566370947, 2.70669043470961, 2.00670669212431, 1.06151912011815, 1.83432913527221, 1.50914995922734, 1.44022564219937, 1.30575866988233, -0.589614662815086, 0.558335664246801, -0.515514667252586, 2.24837115351643, -0.586792151622676, -1.22911307106722, -0.321944586053909, 0.5941033362345, -0.126820929569639, 1.16455104089066, -0.612220105840004, 2.15543467764524, 2.79463672831966, 1.76461661568623, 1.18970042952239, 1.90392323442221, 0.820128453653892, -1.03238248236585, -0.182374330210397, 1.1349587790873, -1.09208688820508, 3.0261291866446, 0.833404717872278, -1.64367078878929, 0.244031420089011, 0.443356289627561, 0.72135747390784, 2.70477650094062, 0.421717284804278, -0.3617188678887, -0.654277052965218, 0.379919595830044, 1.26347991694893, 0.524405360701437, 0.913854654123115, 0.185366146391167, 0.195582957428515, 1.34012528839471, 1.6803863463402, 1.2198218331505, 0.0184539768943683, -1.53666198555819, 1.35884740670205, 0.615508265844485, 1.53708564653686, 1.54482970443691, 1.07775274507472, -0.40142357143868, -0.0374398164962227, 1.47626427907725, 1.6901291894783, 1.21823088039179] + }, + "n_obs": 120, + "coef": 0.405853325711128, + "iid_default": { + "se": 0.115064780324081, + "df_k": 36 + }, + "iid_k_exact": { + "se": 0.111785906348027, + "df_k": 31 + } +} diff --git a/diff_diff/estimators.py b/diff_diff/estimators.py index 1f4043ecd..272d0d8f2 100644 --- a/diff_diff/estimators.py +++ b/diff_diff/estimators.py @@ -29,6 +29,7 @@ from diff_diff.results import DiDResults, MultiPeriodDiDResults, PeriodEffect from diff_diff.utils import ( WildBootstrapResults, + absorbed_fe_rank, build_fe_dummy_blocks, demean_by_groups, fe_dummy_names, @@ -492,10 +493,20 @@ def fit( vars_to_demean = [outcome, treatment, time, "_treat_time"] + (covariates or []) _absorb_regressors = vars_to_demean[1:] # everything except outcome _pre_norms = pre_demean_norms(working_data, _absorb_regressors, weights=survey_weights) + # Absorbed df MUST be measured before the in-place demean below + # overwrites the group columns with demeaned floats. Equals + # demean_by_groups' historical `sum_d (n_d - 1)` on a connected panel; + # smaller when the incidence graph splits (disconnected/hierarchical). + _absorbed_df = absorbed_fe_rank( + working_data, + list(absorb), + has_intercept_col=True, + weights=survey_weights, + ) # Method of alternating projections: for N > 1 absorbed dimensions a # single sequential sweep is only exact on balanced (orthogonal-FE) # panels; demean_by_groups iterates to the exact (W)LS-FWL residual. - working_data, n_fe = demean_by_groups( + working_data, _ = demean_by_groups( # count superseded by absorbed_fe_rank above working_data, vars_to_demean, list(absorb), @@ -515,7 +526,7 @@ def fit( display_names={"_treat_time": f"{treatment}:{time}"}, weights=survey_weights, ) - n_absorbed_effects += n_fe + n_absorbed_effects += _absorbed_df absorbed_vars = list(absorb) # Extract variables (may be demeaned if absorb was used) @@ -1751,9 +1762,17 @@ def fit( # type: ignore[override] ) _absorb_regressors = vars_to_demean[1:] # everything except outcome _pre_norms = pre_demean_norms(working_data, _absorb_regressors, weights=survey_weights) + # Absorbed df MUST be measured before the in-place demean below + # overwrites the group columns (see the DiD path for the rationale). + _absorbed_df = absorbed_fe_rank( + working_data, + list(absorb), + has_intercept_col=True, + weights=survey_weights, + ) # Method of alternating projections (exact for unbalanced panels; a # single sequential sweep is exact only on balanced orthogonal-FE panels). - working_data, n_fe = demean_by_groups( + working_data, _ = demean_by_groups( # count superseded by absorbed_fe_rank above working_data, vars_to_demean, list(absorb), @@ -1778,7 +1797,7 @@ def fit( # type: ignore[override] }, weights=survey_weights, ) - n_absorbed_effects += n_fe + n_absorbed_effects += _absorbed_df # Extract outcome and treatment (may be demeaned if absorb was used) y = working_data[outcome].values.astype(float) diff --git a/diff_diff/sun_abraham.py b/diff_diff/sun_abraham.py index 9826fa4a1..a199647a4 100644 --- a/diff_diff/sun_abraham.py +++ b/diff_diff/sun_abraham.py @@ -24,6 +24,7 @@ from diff_diff.results import _format_survey_block, _get_significance_stars from diff_diff.results_base import BaseResults from diff_diff.utils import ( + absorbed_fe_rank, pre_demean_norms, safe_inference, snap_absorbed_regressors, @@ -1653,7 +1654,16 @@ def _fit_saturated_regression( cluster_ids = df_demeaned[cluster_var].values else: cluster_ids = None - df_adj = n_units_fe + n_times_fe - 1 + # Absorbed df from the PRE-transform frame. This design carries NO + # intercept column (coef_offset = 0), so it takes the raw FE rank — + # equal to the historical `n_units_fe + n_times_fe - 1` on a + # connected panel, smaller when the incidence graph splits. + df_adj = absorbed_fe_rank( + df, + [unit, time], + has_intercept_col=False, + weights=survey_weights, + ) # Interactions occupy columns 0..n_interactions-1 (no intercept) coef_offset = 0 diff --git a/diff_diff/twfe.py b/diff_diff/twfe.py index fc83dc6f3..84f42024b 100644 --- a/diff_diff/twfe.py +++ b/diff_diff/twfe.py @@ -16,6 +16,7 @@ from diff_diff.linalg import LinearRegression from diff_diff.results import DiDResults from diff_diff.utils import ( + absorbed_fe_rank, build_fe_dummy_blocks, fe_dummy_names, pre_demean_norms, @@ -414,7 +415,16 @@ def fit( # type: ignore[override] for cov in covariates: X_list.append(data_demeaned[f"{cov}_demeaned"].values) X = np.column_stack([np.ones(len(y))] + X_list) - df_adjustment = n_units + n_times - 2 + # Absorbed df from the PRE-transform frame. Equals the historical + # `n_units + n_times - 2` on a connected panel; smaller when the + # unit x time incidence graph splits into components (disconnected + # or hierarchical FE), where the old count over-stated rank. + df_adjustment = absorbed_fe_rank( + data, + [unit, time], + has_intercept_col=True, + weights=survey_weights, + ) # Within-transform path: preserve the historical # `{"ATT": att}` user-facing `result.coefficients` contract. # Broadening this dict here would silently change the diff --git a/diff_diff/utils.py b/diff_diff/utils.py index a250873fb..46ee594a8 100644 --- a/diff_diff/utils.py +++ b/diff_diff/utils.py @@ -2712,7 +2712,9 @@ def demean_by_group( data : pd.DataFrame DataFrame with demeaned variables. n_effects : int - Number of absorbed fixed effects (nunique - 1). + Raw level count (``nunique - 1``). **NOT a valid degrees-of-freedom + adjustment** on weighted fits — it counts levels carried only by + zero-weight rows; for df adjustments use :func:`absorbed_fe_rank`. Examples -------- @@ -2918,6 +2920,115 @@ def _demean_map_rust( return demeaned, iters_all +def absorbed_fe_rank( + data: pd.DataFrame, + group_vars: List[str], + *, + has_intercept_col: bool, + weights: Optional[np.ndarray] = None, +) -> int: + """Degrees of freedom absorbed by a set of fixed-effect dimensions. + + Returns the rank of the absorbed FE dummy space, minus one when the caller's + design matrix already carries its own intercept column (so the shared + constant is not counted twice). + + The rank INCLUDING an implicit intercept is:: + + N == 2: sum(levels) - C C = connected components of the bipartite + level graph (Abowd-Creecy-Kramarz) + else: sum(levels) - N + 1 + + For a single connected component the two-dimensional form collapses to the + historical ``sum_d (n_d - 1) + 1``, so every existing caller is unchanged on + a connected panel. It differs — correctly — in two supported cases the old + count over-stated: + + - **Disconnected panels** (``C > 1``): the FE are not jointly identified + across components, so the dummy space is rank ``sum(levels) - C``. + - **Nested / hierarchical dimensions** (e.g. ``["state", "state_year"]``): + each parent level forms its own component, so ``C = n_parents`` and the + absorbed rank is that of the finer dimension alone. + + ``N >= 3`` keeps the ``sum(levels) - N + 1`` form, which is exact for + mutually independent, connected dimensions and over-states rank otherwise + (a duplicated third dimension is the simplest counterexample). Computing the + general N-way rank is a hypergraph problem and is deliberately not attempted + here; see the TODO backlog row. + + Parameters + ---------- + data : DataFrame + The frame the design is built from. Must be pre-transform: after an + in-place demean the group columns may hold demeaned floats. + group_vars : list of str + Absorbed FE columns. + has_intercept_col : bool + True when the caller's ``X`` already contains an intercept column. + weights : ndarray, optional + Observation weights. Levels and connectivity are evaluated over + POSITIVE-weight rows only, per the REGISTRY guarantee that zero-weight + padding is inference-invariant on the generic HC1/classical paths. + + Returns + ------- + int + Absorbed degrees of freedom for the caller's df adjustment. + """ + if not group_vars: + return 0 + + frame = data + if weights is not None: + w = np.asarray(weights, dtype=np.float64) + if w.shape[0] != len(data): + raise ValueError(f"weights length ({w.shape[0]}) must match data rows ({len(data)})") + positive = w > 0 + if not positive.any(): + return 0 + if not positive.all(): + frame = data.loc[positive] + + codes_list = [] + for g in group_vars: + codes = pd.factorize(frame[g].values, sort=False)[0] + if codes.size and codes.min() < 0: + # pd.factorize assigns NaN keys code -1; a negative index would + # otherwise surface as a cryptic sparse-graph error downstream. + raise ValueError( + f"Absorbed fixed-effect column '{g}' contains NaN group keys; " + "drop or impute those rows before fitting." + ) + codes_list.append(codes) + levels = [int(c.max()) + 1 if c.size else 0 for c in codes_list] + total_levels = sum(levels) + if total_levels == 0: + return 0 + + n_dims = len(group_vars) + if n_dims == 2: + from scipy.sparse import coo_matrix + from scipy.sparse.csgraph import connected_components + + left, right = codes_list + offset = levels[0] + adjacency = coo_matrix( + ( + np.ones(left.shape[0], dtype=np.int8), + (left, right + offset), + ), + shape=(total_levels, total_levels), + ) + # Weak connectivity of the directed one-way bipartite graph equals + # undirected connectivity, and skips materializing A + A.T (~2.4x). + n_components = int(connected_components(adjacency, directed=True, connection="weak")[0]) + rank = total_levels - n_components + else: + rank = total_levels - n_dims + 1 + + return rank - 1 if has_intercept_col else rank + + def demean_by_groups( data: pd.DataFrame, variables: List[str], @@ -2976,8 +3087,13 @@ def demean_by_groups( data : pd.DataFrame DataFrame with demeaned variables. n_effects : int - Number of absorbed fixed effects, ``sum_d (nunique_d - 1)`` over - ``group_vars`` (the standard DOF-accounting convention). + Raw level count ``sum_d (nunique_d - 1)`` over ``group_vars``. + **NOT a valid degrees-of-freedom adjustment**: it over-counts the + absorbed dummy-space rank on disconnected or nested/hierarchical + dimensions and counts levels carried only by zero-weight rows. No + library caller consumes it as df any more — for df adjustments use + :func:`absorbed_fe_rank`, which computes the component-aware rank + over positive-weight rows. Raises ------ diff --git a/docs/conf.py b/docs/conf.py index 5633878fc..0342ff032 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -61,6 +61,7 @@ "v4-design.md", "methodology/continuous-did.md", "methodology/survey-theory.md", + "methodology/variance-conventions.md", # Internal paper-review notes (methodology validation artifacts). "methodology/papers/*", "tutorials/README.md", diff --git a/docs/doc-deps.yaml b/docs/doc-deps.yaml index f5f53335f..fb0d5daa3 100644 --- a/docs/doc-deps.yaml +++ b/docs/doc-deps.yaml @@ -88,6 +88,9 @@ sources: diff_diff/estimators.py: drift_risk: medium docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "DifferenceInDifferences, TwoWayFixedEffects, MultiPeriodDiD" type: methodology @@ -122,6 +125,9 @@ sources: diff_diff/twfe.py: drift_risk: low docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "TwoWayFixedEffects" type: methodology @@ -133,6 +139,9 @@ sources: diff_diff/staggered.py: drift_risk: high docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "CallawaySantAnna" type: methodology @@ -198,6 +207,9 @@ sources: diff_diff/sun_abraham.py: drift_risk: medium docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "SunAbraham" type: methodology @@ -227,6 +239,9 @@ sources: diff_diff/imputation.py: drift_risk: medium docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "ImputationDiD" type: methodology @@ -257,6 +272,9 @@ sources: diff_diff/two_stage.py: drift_risk: medium docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "TwoStageDiD" type: methodology @@ -604,6 +622,9 @@ sources: diff_diff/stacked_did.py: drift_risk: medium docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "StackedDiD" type: methodology @@ -630,6 +651,9 @@ sources: diff_diff/wooldridge.py: drift_risk: medium docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "WooldridgeDiD" type: methodology @@ -656,6 +680,9 @@ sources: diff_diff/lpdid.py: drift_risk: medium docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "LPDiD" type: methodology @@ -990,6 +1017,9 @@ sources: diff_diff/linalg.py: drift_risk: medium docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "Variance Estimation, Cluster-Robust SE" type: methodology @@ -1013,6 +1043,9 @@ sources: diff_diff/utils.py: drift_risk: medium docs: + - path: docs/methodology/variance-conventions.md + type: methodology + note: "Clustered-variance conventions inventory (CR1 k + tail-df map, D/L classification); table generated by tests/test_variance_conventions.py" - path: docs/methodology/REGISTRY.md section: "Inference, safe_inference NaN gating" type: methodology diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 408b5afa4..6676c60e1 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -73,6 +73,11 @@ where τ is the ATT. - Default: HC1 heteroskedasticity-robust - Optional: Cluster-robust (specify `cluster` parameter) - Optional: Wild cluster bootstrap for small number of clusters +- With `absorb=`, the absorbed-FE degrees-of-freedom adjustment uses the + component-aware rank (`diff_diff.utils.absorbed_fe_rank`) — see the + TwoWayFixedEffects section's absorbed-FE degrees-of-freedom note for the + shared convention and its fixest `K.exact` anchor. Consumed by the + classical/HC1 variance rescale and the reported residual df. *Edge cases:* - Empty cells (e.g., no treated-pre observations) cause rank deficiency, handled per `rank_deficient_action` setting @@ -279,7 +284,11 @@ where V is the VCV sub-matrix for post-treatment δ_e coefficients. variance. Regression: `tests/test_methodology_wls_cr2.py::TestLinearRegressionFENanGuardEndToEnd`. - Optional: Wild cluster bootstrap (complex for multi-coefficient testing; requires joint bootstrap distribution) -- Degrees of freedom adjusted for absorbed fixed effects +- Degrees of freedom adjusted for absorbed fixed effects: component-aware rank + via `diff_diff.utils.absorbed_fe_rank` (own `_absorbed_fe_vcov_scale` gate at + its fit site — MultiPeriodDiD is a second implementation, not an alias of + DifferenceInDifferences) — see the TwoWayFixedEffects section's absorbed-FE + degrees-of-freedom note for the shared convention *Edge cases:* - Reference period: omitted from design matrix; coefficient is zero by construction. @@ -364,12 +373,46 @@ This matches the behavior of R's `fixest::feols()` with absorbed FE. *Standard errors:* - Default: Cluster-robust at unit level (accounts for serial correlation) -- Degrees of freedom adjusted for absorbed fixed effects: `df_adjustment = n_units + n_times - 2` +- Degrees of freedom adjusted for absorbed fixed effects: component-aware rank via + `diff_diff.utils.absorbed_fe_rank` — equal to the historical + `n_units + n_times - 2` on a connected panel, and `sum(levels) - C - 1` in + general, where `C` is the number of connected components of the unit x time + incidence graph (Abowd-Creecy-Kramarz). +- **Note:** the pre-2026-07 count `sum_d (levels_d - 1)` over-stated absorbed rank + in two supported cases: disconnected panels (the FE are not jointly identified + across components, true rank `sum(levels) - C`) and nested/hierarchical + dimensions (e.g. `absorb=["state", "state_year"]`, where each parent level is + its own component: measured 6 states x 5 years gives true rank 29 vs the old + count 34). On those inputs the reported residual df was too small and the + non-clustered classical/hc1 SEs too large; connected panels are bit-identical. + Levels and connectivity are evaluated over POSITIVE-weight rows only, per the + zero-weight-padding inference-invariance guarantee (see "Weight Type Effects on + Inference"); the old all-rows `nunique()` count violated that guarantee on + weighted fits with inert rows. For 3+ absorbed dimensions the historical + `sum(levels - 1)` form is retained (exact for independent connected dims; + the general N-way rank is tracked in TODO.md). Full measured inventory: + `docs/methodology/variance-conventions.md`. - **Note (absorbed-FE variance scale = fixest full-K):** for the *non-clustered* `classical` and `hc1` (hetero) variance families, the finite-sample scale (`sse/(n-k)` / `n/(n-k)`) now counts the absorbed FE in `k` -- i.e. `K_full = k_visible + df_adjustment` -- matching `fixest::feols(vcov="iid"/"hetero")` and the reported t-`df` (`linalg._absorbed_fe_vcov_scale`, a single scalar rescale of the `k_visible` vcov, fail-closed when `n - K_full <= 0`). + **Note (deviation from R default — irregular FE designs):** with the 2026-07 + component-aware `absorbed_fe_rank`, `df_adjustment` on disconnected or + nested/hierarchical two-way designs is the EXACT dummy-space rank + (`sum(levels) - C`), which matches `fixest::feols(..., ssc = ssc(K.exact = TRUE))` + at machine precision but deviates from the fixest DEFAULT + (`ssc(K.exact = FALSE)`), whose approximate `sum(levels - 1)` count reproduces the + library's own pre-2026-07 over-count on exactly those designs. Measured on the + committed hierarchical golden (`benchmarks/data/fixest_kexact_golden.json`, + 6 states x 5 nested state-years, C = 6): fixest default `df.K = 36` / iid SE + 0.115064780324081 vs `K.exact = TRUE` `df.K = 31` / iid SE 0.111785906348027; + the library matches the exact side to <= 1e-12 + (`tests/test_variance_conventions.py::TestFixestKExactParity`). Exact rank is + chosen deliberately: `K.exact` exists in fixest precisely because its default is + an approximation, and the approximate count misstates the residual df the + reported t-distribution uses. Connected, independent designs are identical under + both conventions. Previously the within-transform SE used `k_visible`, sitting ~6.5% below fixest even though the t-`df` already used `K_full` (an internal inconsistency). Applies to `TwoWayFixedEffects(vcov_type="classical")`, `DifferenceInDifferences(absorb=..., vcov_type in {classical,hc1})`, @@ -405,7 +448,23 @@ This matches the behavior of R's `fixest::feols()` with absorbed FE. estimators (CS, SA, imputation-family, etc.) carry their own inference stacks and are out of the knob's scope. 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 (it already matched fixest). + `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 + the residual df and the non-clustered variance families; it does NOT extend to the + clustered CR1 finite-sample factor** — see the deviation note below. +- **Deviation from R (clustered CR1, full-dummy `fixed_effects=` path):** with + `fixed_effects=[unit, time]` and `cluster=unit`, the CR1 factor's `k` counts every + FE dummy including the cluster-NESTED unit FE (measured `k = 66` on the audit + panel), which fixest's default `ssc(fixef.K = "nested")` drops. The full-dummy + clustered SE is therefore CONSERVATIVE relative to fixest — the opposite direction + from the `absorb=` path's anti-conservative `k_visible` (measured `k = 2` on the + same model). The two documented-equivalent idioms consequently return clustered + SEs differing by 10.35% on an identical model with an identical ATT + (`sqrt((360−2)/(360−66))`, pinned by + `tests/test_variance_conventions.py::test_d1_divergence_is_pinned`), and NEITHER + matches the reghdfe/fixest nested convention. Full measured map: + `docs/methodology/variance-conventions.md` (defect D1); convergence of both paths + on `K_reference` is scheduled as PR B of the 3.9 variance-consolidation program. *Edge cases:* - Singleton units/periods are automatically dropped @@ -1520,6 +1579,12 @@ where weights ŵ_{g,e} = n_{g,e} / Σ_g n_{g,e} (sample share of cohort g at eve aggregation step is otherwise identical. Tracked as a follow-up (harmonizing the correction or documenting it as an intentional difference). +- The within-transform residual df subtracts the absorbed unit + time FE via + the component-aware rank (`diff_diff.utils.absorbed_fe_rank`, no-intercept + form — the saturated design carries no intercept column) — see the + TwoWayFixedEffects section's absorbed-FE degrees-of-freedom note for the + shared convention. On disconnected panels the absorbed count falls by C-1, + raising the per-cell residual df used for cohort-period t-inference. - Survey designs (`survey_design=`) + `vcov_type ∈ {"classical","hc2", "hc2_bm"}` are rejected at fit-time: the survey-design Taylor Series Linearization (or replicate-weight refit) variance overrides the @@ -1998,7 +2063,7 @@ where `g(·)` is the link inverse (logistic or exp), `η_i` is the individual li - **Note:** QMLE sandwich uses `weight_type="aweight"` which applies `(G/(G-1)) * ((n-1)/(n-k))` small-sample adjustment. Stata `jwdid` uses `G/(G-1)` only. The `(n-1)/(n-k)` term is conservative (inflates SEs slightly). For typical ETWFE panels where n >> k, the difference is negligible. *Variance families (`vcov_type`, OLS path only):* -- `hc1` (default) — CR1 Liang-Zeger cluster-robust on the within-transformed design. Bit-equal to prior behavior (FWL preserves the score). The natural R anchor is `fixest::feols(y ~ | unit + time, cluster=~unit)` or Stata `jwdid` (both within-transform). **Deviation from Stata `jwdid` (measured 2026-07-26, `tests/test_etwfe_cs_stata_parity.py`):** the ATT(g,t) POINT estimates match `jwdid` exactly (~3e-8 on the `mpdta` panel, i.e. Stata's log-output rounding), but every `hc1` SE is SMALLER than `jwdid`'s by a factor that is **uniform across cells** (spread < 1e-6 within a fit) and shrinks as the cluster count grows: 1.0280 at G=20, 1.0132 at G=40, 1.0086 at G=60, 1.0046 at G=110, 1.00264 at G=191, 1.0010 at G=500. **Two ratios are PINNED by CI** (`tests/test_etwfe_cs_stata_parity.py`): G=500 on the full `mpdta` panel, and G=191 on the all-eventually-treated arm — each arm measures its own, because the constant does not transfer between cluster counts. The remaining smaller-G figures were measured ad hoc on subsampled panels during the #724 investigation and no committed artifact reproduces them, so treat them as the shape of the trend rather than as regression-gated constants. Committing that subsample ladder as a golden block is part of the derivation work tracked in `TODO.md` — the ladder is the instrument the derivation needs, not a separate chore. The library is therefore systematically **anti-conservative** relative to the reference - negligibly with many clusters (~0.1% at G=500) and materially in few-cluster designs (~2.8% at G=20). **The mechanism is NOT yet identified.** The gap tracks `sqrt(G/(G-1))` closely but lies consistently ABOVE it (by ~0.2% at G=20 down to ~0.001% at G=500), and `solve_ols` already applies the full CR1 `(G/(G-1)) * ((n-1)/(n-k))` - so a missing `G/(G-1)` is ruled out as the explanation. This is recorded as a measured deviation, not a diagnosed one; deriving the true factor from the within-transform `k` accounting against `hdfe`/`reghdfe`'s is tracked in `TODO.md`, and no correction should be applied until it reproduces `jwdid` exactly rather than approximately. This is distinct from the `lm + clubSandwich` deviation below, whose factor is `k`-based. `CallawaySantAnna` shows no such gap - its SEs match Stata `csdid` outright - which localizes this to the ETWFE path rather than a library-wide convention. **Deviation from R `lm + clubSandwich::vcovCR(type="CR1S")`:** the full-dummy `lm` SE differs by a factor of `sqrt((n - k_within) / (n - k_total))` because clubSandwich's `(n-1)/(n-p)` finite-sample correction counts ALL columns (intercept + treatment + unit dummies + time dummies = `k_total`) while WooldridgeDiD's `solve_ols` on the within-transformed design counts only the treatment-cell columns (`k_within`). On the 240-obs / 51-column R-parity fixture this is ~11%; on typical larger panels (n >> k_total) the gap shrinks to <2%. No public WooldridgeDiD code path exposes the `lm + CR1S` (CR1 cluster-robust on the full-dummy design) finite-sample correction — `vcov_type="hc2_bm"` routes to the CR2 Bell-McCaffrey sandwich on the full-dummy design (different variance estimator entirely), not CR1S. Users who need exact `lm + clubSandwich::vcovCR(type="CR1S")` parity must call `solve_ols` directly on a full-dummy design or fit via R. Same deviation pattern as SunAbraham PR #472 (`fixest::sunab` vs `lm + clubSandwich`). +- `hc1` (default) — CR1 Liang-Zeger cluster-robust on the within-transformed design. Bit-equal to prior behavior (FWL preserves the score). The natural R anchor is `fixest::feols(y ~ | unit + time, cluster=~unit)` or Stata `jwdid` (both within-transform). **Deviation from Stata `jwdid` (measured 2026-07-26, `tests/test_etwfe_cs_stata_parity.py`):** the ATT(g,t) POINT estimates match `jwdid` exactly (~3e-8 on the `mpdta` panel, i.e. Stata's log-output rounding), but every `hc1` SE is SMALLER than `jwdid`'s by a factor that is **uniform across cells** (spread < 1e-6 within a fit) and shrinks as the cluster count grows: 1.0280 at G=20, 1.0132 at G=40, 1.0086 at G=60, 1.0046 at G=110, 1.00264 at G=191, 1.0010 at G=500. **Two ratios are PINNED by CI** (`tests/test_etwfe_cs_stata_parity.py`): G=500 on the full `mpdta` panel, and G=191 on the all-eventually-treated arm — each arm measures its own, because the constant does not transfer between cluster counts. The remaining smaller-G figures were measured ad hoc on subsampled panels during the #724 investigation and no committed artifact reproduces them, so treat them as the shape of the trend rather than as regression-gated constants. Committing that subsample ladder as a golden block is part of the derivation work tracked in `TODO.md` — the ladder is the instrument the derivation needs, not a separate chore. The library is therefore systematically **anti-conservative** relative to the reference - negligibly with many clusters (~0.1% at G=500) and materially in few-cluster designs (~2.8% at G=20). **The mechanism IS now derived (2026-07, `docs/methodology/variance-conventions.md` defect D2):** the clustered CR1 factor's `k` counts only the visible treatment-cell columns, omitting the absorbed FE not nested in the cluster; the closed form `K_reference = explicit columns + (1 if no intercept column) + rank(non-nested FE | nested)` reproduces `jwdid`/reghdfe to ~1e-15 on all three committed arms and retrodicts the ad-hoc G=20/G=40 rungs (predicted 1.028016 vs recorded 1.0280; 1.013210 vs 1.0132). The correction itself is deliberately NOT yet applied — it lands with the rest of the clustered-CR1 convergence (PR B of the 3.9 variance-consolidation program), together with the committed subsample ladder, so every affected surface moves once under one ledger row rather than piecemeal. This is distinct from the `lm + clubSandwich` deviation below, whose factor is `k`-based. `CallawaySantAnna` shows no such gap - its SEs match Stata `csdid` outright - which localizes this to the ETWFE path rather than a library-wide convention. **Deviation from R `lm + clubSandwich::vcovCR(type="CR1S")`:** the full-dummy `lm` SE differs by a factor of `sqrt((n - k_within) / (n - k_total))` because clubSandwich's `(n-1)/(n-p)` finite-sample correction counts ALL columns (intercept + treatment + unit dummies + time dummies = `k_total`) while WooldridgeDiD's `solve_ols` on the within-transformed design counts only the treatment-cell columns (`k_within`). On the 240-obs / 51-column R-parity fixture this is ~11%; on typical larger panels (n >> k_total) the gap shrinks to <2%. No public WooldridgeDiD code path exposes the `lm + CR1S` (CR1 cluster-robust on the full-dummy design) finite-sample correction — `vcov_type="hc2_bm"` routes to the CR2 Bell-McCaffrey sandwich on the full-dummy design (different variance estimator entirely), not CR1S. Users who need exact `lm + clubSandwich::vcovCR(type="CR1S")` parity must call `solve_ols` directly on a full-dummy design or fit via R. Same deviation pattern as SunAbraham PR #472 (`fixest::sunab` vs `lm + clubSandwich`). - `hc2_bm` — CR2 Bell-McCaffrey via auto-route to full-dummy design (`[intercept, X_design, unit_dummies, time_dummies]`), then `solve_ols(..., vcov_type="hc2_bm")` through the clubSandwich port (PR #475). FWL does NOT preserve the hat matrix; HC2 leverage + BM DOF require the full-projection design. Per-coefficient SE matches `clubSandwich::vcovCR(lm(...), cluster=~unit, type="CR2")` at atol=1e-10. Per-cell `(g, t)` inference fields use `coef_test()$df_Satt` Bell-McCaffrey DOF (pinned at atol=1e-6 from CI half-width inversion). Aggregated inference (overall ATT + `.aggregate("group" | "calendar" | "event")`) uses contrast-specific BM DOFs from `_compute_cr2_bm_contrast_dof` (matches R `Wald_test(constraints=matrix(w, 1), vcov=vcov_CR2, test="HTZ")$df_denom`); the overall ATT contrast DOF is computed at fit time, the other three aggregations lazily on each `.aggregate(...)` call from BM artifacts (the REDUCED kept-column `X` / `cluster_ids` / bread matrix + the reduced-space coef-index map) stored on the Results object — using the reduced design after rank-deficient drops keeps the bread non-singular and matches the subspace `solve_ols` actually estimated in. Fail-closed across all surfaces: when BM DOF is unavailable (helper raises or returns non-finite), the affected inference fields are NaN — not normal-theory fallback (per `feedback_bm_contrast_dof_fail_closed`). - `classical`, `hc2` — supported via auto-route to full-dummy AND auto-drop of the unit auto-cluster (one-way families don't compose with `cluster_ids` per the linalg validator). Set `self.cluster=None` (default) for these; explicit `cluster="state"` + one-way family raises at the linalg validator. SE matches `summary(lm(...))$coefficients` (classical) and `sandwich::vcovHC(type="HC2")` respectively. Per-cell + aggregate p-values/CIs use the residual DOF `n - rank(X)` (matches R `lm()` / `coef_test()` t-distribution under both classical OLS SE and `sandwich::vcovHC` defaults) — not normal-theory, so inference is correct under small samples. - `conley` (spatial-HAC, Conley 1999) — supported on the **OLS path** via the within-transform design (or the full-dummy design when `cohort_trends=True`, like the other full-dummy families — see the cohort-trends row below), threading the `conley_*` params through `solve_ols` / `conley.py` (`conley_lag_cutoff=0` = within-period spatial only; `>0` adds within-unit Bartlett serial — the panel-aware path, since `conley_time`/`conley_unit` are always supplied, not pooled cross-sectional). Reuses the already-`conleyreg`-validated machinery (no new variance code). The unit auto-cluster is dropped on the conley path (an explicit `cluster=` enables the spatial+cluster product kernel); `survey_design=` / `weights` / `n_bootstrap>0` are rejected, and `method ∈ {logit, poisson}` + conley remains rejected (the `method != "ols"` guard — a QMLE-on-pseudo-residuals Conley sandwich is a separate derivation). FWL-composability (the within-transform conley SE equals the full-dummy conley SE) is pinned in `tests/test_conley_vcov.py::TestConleyWooldridge::test_fwl_composability_vs_full_dummy`. diff --git a/docs/methodology/variance-conventions.md b/docs/methodology/variance-conventions.md new file mode 100644 index 000000000..06676b498 --- /dev/null +++ b/docs/methodology/variance-conventions.md @@ -0,0 +1,118 @@ +# Clustered-variance conventions: measured inventory + +> **Repo-internal methodology document** (excluded from the published docs +> build). This is the map behind the 3.9 variance-consolidation program. The +> table below is **generated** from the expected literals in +> tests/test_variance_conventions.py — run +> python -m tests.test_variance_conventions and paste the output here after +> any change that legitimately moves a cell. The parametrized tests assert +> those literals against live instrumentation (fast subset in the default +> suite, full sweep under -m slow), so a stale table fails CI on the +> expected-literal side, not silently. + +## The measured matrix + +Shared DGP: numpy.default_rng(7), 60 units x periods 1..6, +first_treat = [0, 3, 4, 5][unit % 4], treated = 1{ft > 0 and t >= ft}, +grp = unit % 2 (time-invariant), post = 1{t >= 4}, +y = N(0,1) + 0.3 unit + 0.2 t + 2 treated; n = 360, G = 60 under +cluster="unit". Each row's exact fit kwargs are in the test file — the +k values and SEs in any claim below come from the SAME fit (mixing panels +produced wrong figures three separate times while this inventory was drafted). + +| surface | CR1 `k` (multiset) | tail df (multiset) | status | reason | +|---|---|---|---|---| +| `did_absorb_hc1_cluster_unit` | 2 | 294 | **defect** | D2: CR1 k omits absorbed FE not nested in the cluster (time) | +| `did_fixed_effects_hc1_cluster_unit` | 66 | 294 | **defect** | D1: same model as did_absorb yet k=66 vs 2 -> SEs differ 10.35%; full-dummy k also counts the cluster-nested unit FE the references drop | +| `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` | 2 | 298 | **defect** | D2 (within-transform k_visible); tail df is residual n-K_full | +| `wooldridge_hc1_within` | 9 | None, None, None, None, None, None, None, None, None, None | **defect** | D2 (k_visible=cells only) + normal-theory tail df with no df_convention knob | +| `sun_abraham_hc1` | 15 | 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, 280, None, None, None, None, None, None, None, None | **defect** | D2 + D4: residual df per cohort-period cell but normal theory on aggregates | +| `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 | +| `imputation_default` | — | None | **legitimate** | L3: BJS imputation variance, not the shared CR1 sandwich | +| `imputation_pretrends_event_study` | unpinned | unpinned | **defect** | pretrends lead regression runs the shared clustered CR1 with k_visible | +| `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 | + +cr1_k is the sorted multiset of visible column counts reaching the shared +clustered CR1 denominator (linalg._compute_robust_vcov_numpy with +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 (last measured: +k=3, normal-theory tail df on all 8 event-study leads/lags). Captured under the canonical Python backend; Rust +and Python agree to <= 8e-15 on every surface because both implement the same +conventions today. + +## Defects (scheduled: 3.9 consolidation program) + +- **D1 — absorb= vs fixed_effects=: 10.35% SE split on the same model.** + Identical ATT, se 0.2414226781 (k=2) vs 0.2664071714 (k=66), ratio + 1.103489 = sqrt((360-2)/(360-66)) exactly. The absorb-side deviation from + fixest is documented (REGISTRY "Known limitation (deviation from fixest)", + DEFERRED.md); the *user-facing consequence* — two documented-equivalent kwargs + disagreeing by 10%, in opposite directions from the reference — was not, until + this inventory. Fix: PR B converges both on K_reference. +- **D2 — clustered CR1 never counts absorbed FE not nested in the cluster.** + k_visible in the denominator; _absorbed_fe_vcov_scale exists but is + gated on cluster_ids is None. Anti-conservative: SEs understated by 0.10% + (n=2500), 1.30% (n=200), 5.51% (n=60). The correction + K_reference = explicit cols + (1 if no intercept col) + rank(non-nested FE | nested) + reproduces Stata reghdfe 3.2.9 (via jwdid) to ~1e-15 on three arms and + R fixest 0.14.2 to ~1e-12 on two, and retrodicts the two measured-but-unpinned + subsample rungs (predicted 1.028016 vs recorded 1.0280 at G=20; 1.013210 vs + 1.0132 at G=40). Fix: PR B. +- **D3 — absorbed rank assumed independent, connected FE dimensions.** + sum(levels - 1) over-counted on disconnected panels (true rank + sum(levels) - C) and on hierarchical specs + (absorb=["state", "state_year"]: true 29, old count 34). **Fixed** by + diff_diff.utils.absorbed_fe_rank (this PR): two-way rank from the + bipartite level graph's connected components; N >= 3 keeps sum(levels-1) + with the limitation documented (over-counts for duplicated/nested triples — + tracked in TODO.md). **External anchor:** the exact rank matches + ``fixest::ssc(K.exact = TRUE)`` at machine precision on the committed + hierarchical golden (``benchmarks/data/fixest_kexact_golden.json``); fixest's + 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. + +## 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. +- **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. + **ImputationDiD is conditional**: its pretrends=True + + aggregate="event_study" lead regression DOES run the shared clustered CR1 + and inherits D2 there (its own matrix row). +- **L4 — hc2/hc2_bm** (leverage / Satterthwaite DOF — no CR1 factor), + **survey TSL** (n_PSU - n_strata over the full design), and **Wooldridge + cohort_trends full-dummy** (documented opt-in landing on the L1 + convention). **conley** is out of this matrix by decision: the spatial-HAC + 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) + +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. 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. diff --git a/docs/performance-plan.md b/docs/performance-plan.md index 98af43f57..86f0732ae 100644 --- a/docs/performance-plan.md +++ b/docs/performance-plan.md @@ -4,6 +4,23 @@ This document outlines the strategy for improving diff-diff's performance on lar --- +## Component-aware absorbed-FE rank on the absorb path (v3.9, 2026-07) + +`diff_diff.utils.absorbed_fe_rank` adds one O(n) pass per absorbed fit (two +`pd.factorize` calls + a sparse weak-connectivity `connected_components` on the +bipartite level graph; the directed weak-connectivity form skips materializing +`A + A.T`, ~2.4x over the naive build). Measured on the county-class shape +(3,100 units x 60 periods, 186k rows): **1.9 ms/call**, ~10 ns/row. Against the +Rust-served TWFE hc1 fit on that shape (24 ms) the share is **7.7%** — above the +2% materiality gate set in the consolidation plan, but the denominator is the +library's fastest fit configuration; on fits with covariates, bootstrap, or the +pure-Python backend the share falls well under 2%, and the absolute cost is flat +O(n) with no allocation cliffs. The remaining mitigation — reusing the factorized +codes `demean_by_groups`/`within_transform` already produce instead of +re-factorizing — requires threading demeaner internals through four call sites +and is deliberately deferred (tracked alongside the D3 N-way row in TODO.md); +correctness on disconnected/hierarchical panels ships first. + ## Opt-in solve_ols normal-equations Cholesky fast path (v3.7, 2026-07) `solve_ols` is the universal OLS entry point; both backends run an equilibrated diff --git a/tests/test_utils.py b/tests/test_utils.py index f2ac8e9f4..112f73a3b 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1519,11 +1519,33 @@ def test_len1_byte_identical_to_demean_by_group(self, weighted): np.testing.assert_array_equal(out_groups["y_dm"].values, out_single["y_dm"].values) np.testing.assert_array_equal(out_groups["x1_dm"].values, out_single["x1_dm"].values) - def test_n_effects_is_sum_nunique_minus_one(self): + def test_n_effects_is_the_raw_level_count_not_df(self): + """The count return is the documented RAW level count, sum(nunique-1). + + It is deliberately NOT the absorbed dummy-space rank: on disconnected + or nested dimensions the two differ, and no library caller consumes + this value as a df adjustment any more (that is absorbed_fe_rank's + job). Pin both the raw contract and the distinction.""" + from diff_diff.utils import absorbed_fe_rank + df = _unbalanced_2way_panel(seed=2) _, n_eff = demean_by_groups(df, ["y"], ["unit", "period"], suffix="_dm") expected = (df["unit"].nunique() - 1) + (df["period"].nunique() - 1) assert n_eff == expected + # Hierarchical dims: the raw count and the component-aware rank split. + h = pd.DataFrame( + [ + {"state": s, "state_year": s * 100 + y, "y": float(s + y)} + for s in range(6) + for y in range(5) + for _ in range(2) + ] + ) + _, raw = demean_by_groups(h, ["y"], ["state", "state_year"], suffix="_dm") + assert raw == (6 - 1) + (30 - 1) == 34 # raw count + assert ( + absorbed_fe_rank(h, ["state", "state_year"], has_intercept_col=True) == 29 + ) # the df-valid rank @pytest.mark.parametrize("weighted", [False, True]) def test_unbalanced_2way_matches_full_dummy_ols(self, weighted): diff --git a/tests/test_variance_conventions.py b/tests/test_variance_conventions.py new file mode 100644 index 000000000..849a746b6 --- /dev/null +++ b/tests/test_variance_conventions.py @@ -0,0 +1,812 @@ +"""Audit matrix: which clustered-variance convention each estimator surface gets. + +This is the CI-enforced inventory behind ``docs/methodology/variance-conventions.md``. +Each row pins, for one (estimator, fit configuration) cell: + +- ``cr1_k`` — the sorted multiset of visible column counts ``k`` reaching the + shared clustered CR1 denominator (``linalg._compute_robust_vcov_numpy``), or + ``()`` when the surface's contract is that it makes NO shared-CR1 call. +- ``tail_df`` — the sorted multiset of ``df`` values passed to + ``safe_inference`` / ``safe_inference_batch`` (``None`` = normal theory). + +The point is visibility, not endorsement: several pinned values are DOCUMENTED +DEFECTS (anti-conservative k accounting, mixed tail-df conventions) scheduled to +change in the 3.9 consolidation program. Every row carries ``status`` and, for +legitimate differences, a ``reason``. When a later PR changes a convention, the +expected literal changes HERE, in one reviewable table. + +Instrumentation notes (each guards against a failure mode that produced wrong +inventory numbers during planning): + +- ``safe_inference`` is bound at import by ~20 modules AND imported dynamically + inside ``LinearRegression.get_inference``. Binder modules are discovered by + scanning ``diff_diff`` at runtime — the hand-maintained list was extended four + times and twice wrongly declared exhaustive. +- The Rust backend resolves ``DIFF_DIFF_BACKEND`` once at import, so the Python + lane is forced by nulling each module's own ``HAS_RUST_BACKEND`` binding — + also discovered programmatically (10 modules hold one). +- Wrappers bind via ``inspect.signature`` so positional and keyword call styles + both register; a naive wrapper silently captured zero calls. +- Rows whose contract is "no shared-CR1 call" assert EXACTLY zero CR1 captures + while still capturing tail df (the zero/non-zero rule is about CR1 only). +""" + +import importlib +import inspect +import pkgutil +import warnings + +import numpy as np +import pandas as pd +import pytest + +import diff_diff + +# --------------------------------------------------------------------------- +# Shared DGP (documented in docs/methodology/variance-conventions.md — keep in +# sync; the doc's table is generated from THIS fixture via +# `python -m tests.test_variance_conventions`). +# --------------------------------------------------------------------------- + + +def make_panel() -> pd.DataFrame: + rng = np.random.default_rng(7) + rows = [] + for unit in range(60): + first_treat = [0, 3, 4, 5][unit % 4] + for time in range(1, 7): + treated = int(first_treat > 0 and time >= first_treat) + grp = unit % 2 + post = int(time >= 4) + y = rng.normal() + 0.3 * unit + 0.2 * time + 2.0 * treated + rows.append( + dict( + unit=unit, + time=time, + first_treat=first_treat, + treated=treated, + grp=grp, + post=post, + y=y, + ) + ) + return pd.DataFrame(rows) + + +# --------------------------------------------------------------------------- +# Programmatic discovery of instrumentation targets. +# --------------------------------------------------------------------------- + + +def _diff_diff_modules(): + """Import and yield every diff_diff submodule (idempotent).""" + for info in pkgutil.iter_modules(diff_diff.__path__, prefix="diff_diff."): + try: + yield importlib.import_module(info.name) + except ImportError: + continue + + +def _modules_binding(attr: str): + """Modules holding their own module-level binding of ``attr``.""" + for mod in _diff_diff_modules(): + if attr in vars(mod): + yield mod + + +class Capture: + """Wrap the inference helpers + the shared CR1 kernel across every lane.""" + + def __init__(self, monkeypatch): + self.cr1_k: list = [] + self.tail_df: list = [] + # Clustered calls in a NON-hc1 family: (vcov_type, k). The CR1 factor + # this matrix audits is the clustered hc1 denominator, so a surface + # that silently switches clustered family must fail its row even when + # the design width is unchanged. + self.unexpected_clustered: list = [] + import diff_diff.linalg as lmod + import diff_diff.utils as umod + + # --- CR1 kernel (numpy lane; the Rust lane is disabled below) --- + orig_vcov = lmod._compute_robust_vcov_numpy + sig_vcov = inspect.signature(orig_vcov) + + def spy_vcov(*a, **k): + b = sig_vcov.bind(*a, **k) + b.apply_defaults() + if b.arguments.get("cluster_ids") is not None: + if b.arguments.get("vcov_type") == "hc1": + self.cr1_k.append(int(b.arguments["X"].shape[1])) + else: + self.unexpected_clustered.append( + (str(b.arguments.get("vcov_type")), int(b.arguments["X"].shape[1])) + ) + return orig_vcov(*a, **k) + + monkeypatch.setattr(lmod, "_compute_robust_vcov_numpy", spy_vcov) + + # --- inference helpers, canonical + every import-time binder --- + for name in ("safe_inference", "safe_inference_batch"): + orig = getattr(umod, name) + sig = inspect.signature(orig) + + def make_spy(orig=orig, sig=sig): + def spy(*a, **k): + b = sig.bind(*a, **k) + b.apply_defaults() + df = b.arguments.get("df") + self.tail_df.append(None if df is None else float(df)) + return orig(*a, **k) + + return spy + + spy = make_spy() + monkeypatch.setattr(umod, name, spy) # dynamic-import lane + for mod in _modules_binding(name): + if mod is not umod: + monkeypatch.setattr(mod, name, spy) + + # --- force the canonical numpy backend in EVERY module holding a flag --- + for mod in _modules_binding("HAS_RUST_BACKEND"): + monkeypatch.setattr(mod, "HAS_RUST_BACKEND", False) + for mod in _diff_diff_modules(): + for rust_name in [n for n in vars(mod) if n.startswith("_rust_")]: + monkeypatch.setattr(mod, rust_name, None) + + def snapshot(self): + return ( + tuple(sorted(self.cr1_k)), + tuple(sorted(self.tail_df, key=lambda v: (v is None, v))), + ) + + +# --------------------------------------------------------------------------- +# The matrix. One entry per (surface, configuration) cell. +# +# status: "defect" — scheduled to change in the 3.9 consolidation program +# "legitimate" — a declared exception with its reason +# cr1_k = () means the row's CONTRACT is "no shared clustered-CR1 call". +# --------------------------------------------------------------------------- + +ROWS = [ + dict( + key="did_absorb_hc1_cluster_unit", + fit=lambda df: diff_diff.DifferenceInDifferences(cluster="unit").fit( + df, outcome="y", treatment="grp", time="post", absorb=["unit", "time"] + ), + cr1_k=(2,), + tail_df=(294.0,), + status="defect", + reason="D2: CR1 k omits absorbed FE not nested in the cluster (time)", + ), + dict( + key="did_fixed_effects_hc1_cluster_unit", + fit=lambda df: diff_diff.DifferenceInDifferences(cluster="unit").fit( + df, outcome="y", treatment="grp", time="post", fixed_effects=["unit", "time"] + ), + cr1_k=(66,), + tail_df=(294.0,), + status="defect", + reason=( + "D1: same model as did_absorb yet k=66 vs 2 -> SEs differ 10.35%; " + "full-dummy k also counts the cluster-nested unit FE the references drop" + ), + ), + dict( + key="did_plain_hc1_cluster_unit", + fit=lambda df: diff_diff.DifferenceInDifferences(cluster="unit").fit( + df, outcome="y", treatment="grp", time="post" + ), + cr1_k=(4,), + tail_df=(356.0,), + status="legitimate", + reason="no absorbed FE: visible k is the whole design; nothing is omitted", + ), + dict( + key="twfe_hc1_cluster_unit_time_post", + fit=lambda df: diff_diff.TwoWayFixedEffects(vcov_type="hc1", cluster="unit").fit( + df, outcome="y", treatment="grp", time="post", unit="unit" + ), + cr1_k=(2,), + tail_df=(298.0,), + status="defect", + reason="D2 (within-transform k_visible); tail df is residual n-K_full", + ), + dict( + key="wooldridge_hc1_within", + fit=lambda df: diff_diff.WooldridgeDiD(method="ols").fit( + df, outcome="y", unit="unit", time="time", cohort="first_treat" + ), + cr1_k=(9,), + tail_df=(None,) * 10, + status="defect", + reason="D2 (k_visible=cells only) + normal-theory tail df with no df_convention knob", + ), + dict( + key="sun_abraham_hc1", + fit=lambda df: diff_diff.SunAbraham().fit( + df, outcome="y", unit="unit", time="time", first_treat="first_treat" + ), + cr1_k=(15,), + tail_df=(280.0,) * 15 + (None,) * 8, + status="defect", + reason="D2 + D4: residual df per cohort-period cell but normal theory on aggregates", + ), + dict( + key="stacked_did_hc1", + fit=lambda df: diff_diff.StackedDiD().fit( + df, outcome="y", unit="unit", time="time", first_treat="first_treat" + ), + cr1_k=(6,), + tail_df=(None,) * 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" + ), + ), + dict( + key="lpdid_pre2_post2", + fit=lambda df: diff_diff.LPDiD(pre_window=2, post_window=2).fit( + df, outcome="y", unit="unit", time="time", treatment="treated" + ), + 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", + ), + dict( + key="imputation_default", + fit=lambda df: diff_diff.ImputationDiD().fit( + df, outcome="y", unit="unit", time="time", first_treat="first_treat" + ), + cr1_k=(), + tail_df=(None,) * 1, + status="legitimate", + reason="L3: BJS imputation variance, not the shared CR1 sandwich", + ), + dict( + key="imputation_pretrends_event_study", + fit=lambda df: diff_diff.ImputationDiD(pretrends=True).fit( + df, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ), + cr1_k=None, # resolved at collection: must be NON-empty (D2 applies here) + tail_df=None, + status="defect", + reason="pretrends lead regression runs the shared clustered CR1 with k_visible", + ), + dict( + key="two_stage_default", + fit=lambda df: diff_diff.TwoStageDiD().fit( + df, outcome="y", unit="unit", time="time", first_treat="first_treat" + ), + cr1_k=(), + tail_df=(None,) * 1, + status="legitimate", + reason="L3: Gardner two-stage variance, not the shared CR1 sandwich", + ), + dict( + key="callaway_santanna_default", + fit=lambda df: diff_diff.CallawaySantAnna().fit( + df, outcome="y", unit="unit", time="time", first_treat="first_treat" + ), + cr1_k=(), + tail_df=(None,) * 16, + status="legitimate", + reason="L3: influence-function variance anchored to Stata csdid", + ), +] + +_FAST_KEYS = { + "did_absorb_hc1_cluster_unit", + "did_fixed_effects_hc1_cluster_unit", + "twfe_hc1_cluster_unit_time_post", + "wooldridge_hc1_within", + "lpdid_pre2_post2", + "callaway_santanna_default", +} + + +def _run_row(row, monkeypatch): + df = make_panel() + cap = Capture(monkeypatch) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + row["fit"](df) + return cap + + +def _assert_row(row, monkeypatch): + cap = _run_row(row, monkeypatch) + cr1_k, tail_df = cap.snapshot() + assert cap.unexpected_clustered == [], ( + f"{row['key']}: clustered vcov calls in a non-hc1 family " + f"{cap.unexpected_clustered} — the audit classifies only vcov_type='hc1' " + "as the shared clustered CR1; this surface switched clustered family" + ) + if row["cr1_k"] is None: + # Contract row: the shared CR1 must be REACHED (the exact k is + # configuration-detail); used where a literal would be brittle. + assert len(cr1_k) > 0, f"{row['key']}: expected shared-CR1 calls, saw none" + elif row["cr1_k"] == (): + assert cr1_k == (), ( + f"{row['key']}: contract is NO shared-CR1 call, captured k={cr1_k} — " + "a refactor routed this estimator through the shared sandwich" + ) + else: + assert cr1_k == row["cr1_k"], f"{row['key']}: cr1_k {cr1_k} != {row['cr1_k']}" + if row["tail_df"] is not None: + assert tail_df == row["tail_df"], f"{row['key']}: tail_df {tail_df} != {row['tail_df']}" + # Self-check: a row claiming tail-df expectations must actually capture some + # (guards the silent-zero-capture failure that produced wrong inventory + # numbers during planning). + if row["tail_df"] not in (None, ()): + assert len(tail_df) > 0, f"{row['key']}: instrumentation captured nothing" + + +@pytest.mark.parametrize("row", [r for r in ROWS if r["key"] in _FAST_KEYS], ids=lambda r: r["key"]) +def test_variance_convention_fast(row, monkeypatch): + """One config per major surface — runs in the default suite.""" + _assert_row(row, monkeypatch) + + +@pytest.mark.slow +@pytest.mark.parametrize( + "row", [r for r in ROWS if r["key"] not in _FAST_KEYS], ids=lambda r: r["key"] +) +def test_variance_convention_full(row, monkeypatch): + """The remaining rows — excluded from the default suite via addopts.""" + _assert_row(row, monkeypatch) + + +def test_every_row_declares_status_and_reason(): + for row in ROWS: + assert row["status"] in ("defect", "legitimate"), row["key"] + assert row["reason"], row["key"] + + +def test_capture_flags_non_hc1_clustered_family(monkeypatch): + """The CR1 classification is family-aware: a clustered fit in another + vcov family lands in unexpected_clustered, never in cr1_k (measured: a + clustered hc2_bm DiD reaches the same kernel with k=4 — the old + cluster_ids-only spy recorded it indistinguishably from the plain-hc1 + row's literal).""" + df = make_panel() + cap = Capture(monkeypatch) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + diff_diff.DifferenceInDifferences(vcov_type="hc2_bm", cluster="unit").fit( + df, outcome="y", treatment="grp", time="post" + ) + assert cap.cr1_k == [] + assert cap.unexpected_clustered, "non-hc1 clustered call was not flagged" + assert all(fam == "hc2_bm" for fam, _ in cap.unexpected_clustered) + + +def test_d1_divergence_is_pinned(): + """The absorb-vs-fixed_effects SE split: same model, same ATT, k=2 vs 66.""" + df = make_panel() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + a = diff_diff.DifferenceInDifferences(cluster="unit").fit( + df, outcome="y", treatment="grp", time="post", absorb=["unit", "time"] + ) + f = diff_diff.DifferenceInDifferences(cluster="unit").fit( + df, outcome="y", treatment="grp", time="post", fixed_effects=["unit", "time"] + ) + np.testing.assert_allclose(a.att, f.att, rtol=0, atol=1e-10) + ratio = f.se / a.se + predicted = np.sqrt((360 - 2) / (360 - 66)) + np.testing.assert_allclose(ratio, predicted, rtol=1e-9) + + +# --------------------------------------------------------------------------- +# Generator: emit the markdown table for docs/methodology/variance-conventions.md +# (run: python -m tests.test_variance_conventions) +# --------------------------------------------------------------------------- + + +def emit_markdown_table() -> str: + """Render the committed inventory table from the rows' EXPECTED literals. + + Deliberately no fitting here: the parametrized tests above assert these + literals against live instrumentation (fast subset in the default suite, + full sweep under ``-m slow``), so the doc-sync gate stays O(formatting) + and the slow split keeps bounding the default suite as the matrix grows. + ``None`` literals (contract rows) render as ``unpinned``. + """ + lines = [ + "| surface | CR1 `k` (multiset) | tail df (multiset) | status | reason |", + "|---|---|---|---|---|", + ] + + def fmt(t): + if t is None: + return "unpinned" + if not t: + return "—" + return ", ".join("None" if v is None else f"{v:g}" for v in t) + + for row in ROWS: + lines.append( + f"| `{row['key']}` | {fmt(row['cr1_k'])} | {fmt(row['tail_df'])} " + f"| **{row['status']}** | {row['reason']} |" + ) + return "\n".join(lines) + + +if __name__ == "__main__": + print(emit_markdown_table()) + + +# --------------------------------------------------------------------------- +# absorbed_fe_rank (D3 fix): unit + end-to-end verification arms. +# --------------------------------------------------------------------------- + + +class TestAbsorbedFeRank: + """Component-aware absorbed-FE rank: bit-identity where C=1, correctness + where the old ``sum(levels - 1)`` count over-stated rank.""" + + @staticmethod + def _connected(): + return pd.DataFrame( + {"unit": np.repeat(np.arange(60), 6), "time": np.tile(np.arange(6), 60)} + ) + + @staticmethod + def _disconnected(seed=11, effect=0.0): + """Two period-disjoint unit blocks: {1,2} and {3,4,5,6} (C=2). + + The block boundary is deliberately NOT the ``post`` cut — a symmetric + 3/3 split makes ``post`` unit-constant and the treatment column gets + rank-dropped, leaving no finite SE to assert. + """ + rng = np.random.default_rng(seed) + rows = [] + for u in range(30): + for t in (1, 2): + rows.append( + dict( + unit=u, + time=t, + grp=int(u % 2), + post=int(t == 2), + y=rng.normal() + 0.1 * u + effect * (u % 2) * (t == 2), + ) + ) + for u in range(30, 60): + for t in (3, 4, 5, 6): + rows.append( + dict( + unit=u, + time=t, + grp=int(u % 2), + post=int(t >= 5), + y=rng.normal() + 0.1 * u + effect * (u % 2) * (t >= 5), + ) + ) + return pd.DataFrame(rows) + + # ---- unit cases ------------------------------------------------------- + + def test_connected_two_way_reproduces_historical_counts(self): + from diff_diff.utils import absorbed_fe_rank + + d = self._connected() + # intercept-bearing designs (TWFE, DiD/MPD absorb): n_units + n_times - 2 + assert absorbed_fe_rank(d, ["unit", "time"], has_intercept_col=True) == 64 + # no-intercept design (SunAbraham): n_units + n_times - 1 + assert absorbed_fe_rank(d, ["unit", "time"], has_intercept_col=False) == 65 + + def test_hierarchical_dims_counted_as_components(self): + """absorb=["state", "state_year"]: each state is its own component, so + the absorbed rank is that of the finer dimension alone (measured true + rank 29 vs the old count 34).""" + from diff_diff.utils import absorbed_fe_rank + + h = pd.DataFrame( + [ + {"state": s, "state_year": s * 100 + y} + for s in range(6) + for y in range(5) + for _ in range(4) + ] + ) + assert absorbed_fe_rank(h, ["state", "state_year"], has_intercept_col=True) == 29 + + def test_disconnected_two_way(self): + from diff_diff.utils import absorbed_fe_rank + + d = self._disconnected() + assert absorbed_fe_rank(d, ["unit", "time"], has_intercept_col=True) == 63 + + def test_single_and_three_way_unchanged(self): + from diff_diff.utils import absorbed_fe_rank + + d = self._connected() + assert absorbed_fe_rank(d, ["unit"], has_intercept_col=True) == 59 + rng = np.random.default_rng(0) + t3 = pd.DataFrame( + { + "a": rng.integers(0, 5, 400), + "b": rng.integers(0, 4, 400), + "c": rng.integers(0, 4, 400), + } + ) + assert absorbed_fe_rank(t3, ["a", "b", "c"], has_intercept_col=True) == 10 + + def test_three_way_duplicated_dim_is_a_known_over_count(self): + """Pinned LIMITATION, not correct behavior: for N>=3 the helper keeps + sum(levels-1), which over-counts when dimensions are duplicated/nested + (true rank beyond intercept here is 7). Tracked in TODO.md.""" + from diff_diff.utils import absorbed_fe_rank + + rng = np.random.default_rng(0) + a = rng.integers(0, 5, 400) + b = rng.integers(0, 4, 400) + t3 = pd.DataFrame({"a": a, "b": b, "c": b.copy()}) + assert absorbed_fe_rank(t3, ["a", "b", "c"], has_intercept_col=True) == 10 + + def test_zero_weight_rows_do_not_contribute_levels_or_edges(self): + """REGISTRY guarantee: zero-weight padding is inference-invariant on the + generic paths — a level carried only by zero-weight rows adds no df.""" + from diff_diff.utils import absorbed_fe_rank + + d = self._connected() + padded = pd.concat( + [d, pd.DataFrame({"unit": [999, 999], "time": [0, 1]})], + ignore_index=True, + ) + w = np.r_[np.ones(len(d)), np.zeros(2)] + assert absorbed_fe_rank(padded, ["unit", "time"], has_intercept_col=True, weights=w) == 64 + + # ---- end-to-end: the df change reaches user-visible inference --------- + + def test_disconnected_end_to_end_did_absorb(self): + """On the C=2 panel the absorbed df is 63 (old count: 64). The + non-clustered hc1 rescale uses it, so the SE discriminates between the + two counts; the clustered lane and the full inference tuple stay + finite.""" + df = self._disconnected(effect=0.5) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + r_un = diff_diff.DifferenceInDifferences(vcov_type="hc1").fit( + df, outcome="y", treatment="grp", time="post", absorb=["unit", "time"] + ) + r_cl = diff_diff.DifferenceInDifferences(cluster="unit").fit( + df, outcome="y", treatment="grp", time="post", absorb=["unit", "time"] + ) + # new count (adj=63): measured; old count (adj=64) would be ~0.44% larger + np.testing.assert_allclose(r_un.se, 0.2811268249, rtol=1e-8) + old_count_se = r_un.se * np.sqrt( + ((180 - 2) / (180 - 2 - 64)) / ((180 - 2) / (180 - 2 - 63)) + ) + assert abs(r_un.se - old_count_se) / r_un.se > 0.003 + for v in (r_cl.se, r_cl.t_stat, r_cl.p_value, *r_cl.conf_int): + assert np.isfinite(v) + + def test_disconnected_end_to_end_twfe_and_mpd(self): + df = self._disconnected(effect=0.5) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + tw = diff_diff.TwoWayFixedEffects(vcov_type="hc1", cluster="unit").fit( + df, outcome="y", treatment="grp", time="time", unit="unit" + ) + mpd = diff_diff.MultiPeriodDiD(cluster="unit").fit( + df, outcome="y", treatment="grp", time="time", absorb=["unit", "time"] + ) + assert np.isfinite(tw.se) + # Discriminating df pins: the C=2 rank lowers the absorbed adjustment + # by 1, so each caller's reported residual df moves off the legacy + # value - reverting either caller to sum(levels-1) reads 114.0 (TWFE) + # / 111.0 (MPD) on this fixture and fails here. + assert tw.inference_df == 115.0 + assert mpd.inference_df == 112.0 + assert mpd.event_study_df is not None + assert set(mpd.event_study_df.values()) == {112.0} + # at least one period effect identified and finite on each component + finite = [e for e in mpd.period_effects.values() if np.isfinite(e.se)] + assert finite, "no finite MPD period effect on the disconnected panel" + + def test_fail_closed_boundary_moves_both_directions(self): + """Lowering the absorbed df by C-1 moves the n - k - adj boundary. + + Newly-finite direction: 8 obs, 2 components of 2x2 -> old adj 6 gives + n-k-adj = 0 (NaN vcov); new adj 5 gives 1 (finite SE). Still-NaN + direction: drop one row -> n-k-adj = 0 under the NEW count too. + """ + rng = np.random.default_rng(5) + rows = [] + for us, ts in [((0, 1), (1, 2)), ((2, 3), (3, 4))]: + for u in us: + for t in ts: + rows.append( + dict( + unit=u, + time=t, + grp=int(u % 2), + post=int(t in (2, 4)), + y=rng.normal() + 0.3 * u, + ) + ) + tiny = pd.DataFrame(rows) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + r = diff_diff.DifferenceInDifferences(vcov_type="hc1").fit( + tiny, outcome="y", treatment="grp", time="post", absorb=["unit", "time"] + ) + assert np.isfinite(r.se), "newly-finite direction: SE must be finite now" + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + r2 = diff_diff.DifferenceInDifferences(vcov_type="hc1").fit( + tiny.iloc[:-1], + outcome="y", + treatment="grp", + time="post", + absorb=["unit", "time"], + ) + assert not np.isfinite(r2.se), "saturated design must stay fail-closed NaN" + assert not np.isfinite(r2.p_value) and not np.isfinite(r2.conf_int[0]) + + def test_disconnected_end_to_end_sun_abraham(self): + """SunAbraham has no public df surface, so capture LinearRegression.df_ + at the fit boundary: on the C=2 panel the no-intercept raw rank is + levels - C = 43, giving df_ = 100 - 3 - 43 = 54 (old count: 53).""" + import diff_diff.linalg as lmod + + rng = np.random.default_rng(13) + rows = [] + for u in range(20): + ft = 0 if u < 10 else 2 + for t in (1, 2): + rows.append( + dict( + unit=u, + time=t, + first_treat=ft, + y=rng.normal() + 0.2 * u + 1.0 * (ft > 0 and t >= ft), + ) + ) + for u in range(20, 40): + ft = 0 if u < 30 else 5 + for t in (4, 5, 6): + rows.append( + dict( + unit=u, + time=t, + first_treat=ft, + y=rng.normal() + 0.2 * u + 1.0 * (ft > 0 and t >= ft), + ) + ) + df = pd.DataFrame(rows) + captured = [] + orig = lmod.LinearRegression.fit + + def spy(self, X, y, **kw): + out = orig(self, X, y, **kw) + captured.append(self.df_) + return out + + lmod.LinearRegression.fit = spy + try: + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + r = diff_diff.SunAbraham().fit( + df, outcome="y", unit="unit", time="time", first_treat="first_treat" + ) + finally: + lmod.LinearRegression.fit = orig + assert ( + captured and captured[-1] == 54 + ), f"SA df_ {captured} != 54 (the old sum(levels-1) count gives 53)" + assert np.isfinite(r.overall_se) + + def test_nan_group_key_raises_actionable_error(self): + """A NaN absorb key must name the offending column, not surface as + scipy's 'negative axis 0 index: -1'.""" + from diff_diff.utils import absorbed_fe_rank + + d = self._connected().astype({"unit": float}) + d.loc[3, "unit"] = np.nan + with pytest.raises(ValueError, match="'unit' contains NaN group keys"): + absorbed_fe_rank(d, ["unit", "time"], has_intercept_col=True) + + +def test_committed_doc_table_matches_generator(): + """The inventory doc's table is GENERATED — assert the committed markdown + block is byte-equal to ``emit_markdown_table()`` so the doc cannot drift + from the fixture (regenerate: ``python -m tests.test_variance_conventions``).""" + import pathlib + + doc = pathlib.Path(__file__).parent.parent / "docs" / "methodology" / "variance-conventions.md" + if not doc.exists(): + # CI's wheel-install legs run the test suite from a temp copy outside + # the repo checkout, where docs/ does not exist (same convention as + # the benchmarks/data golden skip guards). The checkout-based legs + # (pure Python fallback, local runs) still enforce the byte-equality. + pytest.skip("variance-conventions.md not present (run outside the repo checkout)") + text = doc.read_text() + generated = emit_markdown_table() + header = generated.splitlines()[0] + start = text.find(header) + assert start >= 0, "generated table header not found in variance-conventions.md" + committed = text[start : start + len(generated)] + assert committed == generated, ( + "committed inventory table is stale — regenerate with " + "`python -m tests.test_variance_conventions` and paste into the doc" + ) + + +# --------------------------------------------------------------------------- +# R parity: the component-aware rank matches fixest ssc(K.exact = TRUE). +# --------------------------------------------------------------------------- + +_KEXACT_GOLDEN = ( + __import__("pathlib").Path(__file__).parent.parent + / "benchmarks" + / "data" + / "fixest_kexact_golden.json" +) + + +@pytest.mark.skipif( + not _KEXACT_GOLDEN.exists(), + reason=( + "fixest_kexact_golden.json not present; regenerate via " + "`Rscript benchmarks/R/generate_fixest_kexact_golden.R`." + ), +) +class TestFixestKExactParity: + """On a hierarchical two-way design (state_year nested in state, C=6) the + component-aware absorbed rank matches ``fixest::ssc(K.exact = TRUE)`` at + machine precision — a documented deviation from the R DEFAULT + (``K.exact = FALSE``), whose naive count reproduces the library's OLD + ``sum(levels - 1)`` behavior.""" + + @staticmethod + def _load(): + import json + + with open(_KEXACT_GOLDEN) as fh: + return json.load(fh) + + def test_classical_se_matches_k_exact_not_default(self): + from diff_diff.linalg import LinearRegression + from diff_diff.utils import absorbed_fe_rank, demean_by_groups + + g = self._load() + d = pd.DataFrame(g["data"]) + # +1: the demeaned design carries no intercept column, but fixest's K + # counts one; absorbed_fe_rank(has_intercept_col=True) returns the + # rank BEYOND the intercept. + adj = absorbed_fe_rank(d, ["state", "state_year"], has_intercept_col=True) + 1 + dm, _ = demean_by_groups(d.copy(), ["out", "x"], ["state", "state_year"]) + reg = LinearRegression(include_intercept=False, robust=False).fit( + dm[["x"]].values, dm["out"].values, df_adjustment=adj + ) + se = float(np.sqrt(reg.vcov_[0, 0])) + np.testing.assert_allclose(reg.coefficients_[0], g["coef"], rtol=0, atol=1e-12) + np.testing.assert_allclose(se, g["iid_k_exact"]["se"], rtol=0, atol=1e-12) + assert reg.df_ == g["n_obs"] - g["iid_k_exact"]["df_k"] + # Discriminating: the R DEFAULT (naive count) must NOT match. + assert abs(se - g["iid_default"]["se"]) / se > 0.02 + + def test_rank_matches_fixest_exact_k(self): + from diff_diff.utils import absorbed_fe_rank + + g = self._load() + d = pd.DataFrame(g["data"]) + rank_beyond_intercept = absorbed_fe_rank(d, ["state", "state_year"], has_intercept_col=True) + # fixest df.K = x (1) + intercept (1) + absorbed rank + assert rank_beyond_intercept + 2 == g["iid_k_exact"]["df_k"] + # and the naive count reproduces the R DEFAULT + naive = (d["state"].nunique() - 1) + (d["state_year"].nunique() - 1) + assert naive + 2 == g["iid_default"]["df_k"]