diff --git a/CHANGELOG.md b/CHANGELOG.md index ed8ab360..f86d5967 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,6 +40,41 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 `survey_design=SurveyDesign(psu=)`. No behavior change for unclustered fits. ### Fixed +- **`CallawaySantAnna` base period is now selected positionally (sorted-index), matching R + `did::att_gt` on gapped panels.** Previously the base period used literal calendar arithmetic + (`t-1` pre-treatment, `g-1-anticipation` post/universal), so on non-consecutive period grids + (e.g. biennial surveys, skipped years) cells whose calendar base was unobserved were NaN'd as + `missing_period` even though R estimates them from the nearest observed period — which also + corrupted the aggregated event-study / group SEs. The base is now the nearest observed period + (largest observed `p < t` pre-treatment; largest observed `p` with `p + anticipation < g` + post/universal), and the pre/post split is on `t < g` (independent of anticipation, matching a + deparse of `did` 2.5.1 `compute.att_gt`) — resolving the internal inconsistency with the + library's own dCDH estimator, which already used positional neighbors. On consecutive grids + this is byte-identical to the old rule (all existing goldens unchanged); on gapped panels we now + reproduce every R cell (e.g. `fewer_periods` {1,3,4,6} 7/15 previously-NaN cells, `reg` to + ~1e-11 / `dr` to ~1e-4). A single shared `_select_base_period` helper is used across all + estimation paths (panel fast / vectorized / covariate-reg, and repeated cross-sections). For + `base_period="universal"`, each cohort's positional base is now materialized as a zero reference + cell (`att=0`, `se=NaN`) in `group_time_effects` / `to_dataframe("group_time")` at its positional + base event time (`e = base - g`, which can be `-2`, `-3`, … on gapped grids), matching R's + `att_gt` table and `aggte(type="dynamic")` — including the overlapping-reference case where the + zero base dilutes another cohort's estimated pre-trend at the same event time (analytical AND + multiplier-bootstrap paths verified vs `did` 2.5.1 to ~1e-5). The `group` / `simple` aggregations + use post-treatment cells only and exclude the reference. +- **Absorbed-FE (`absorb=`) non-clustered classical/hetero standard errors now use the full-K + finite-sample scale, matching `fixest`.** The within-transform variance scale (`sse/(n-k)` for + `classical`, `n/(n-k)` for `hc1`) counted only the *visible* regressors `k`, excluding the + absorbed fixed effects, so `TwoWayFixedEffects(vcov_type="classical")`, + `DifferenceInDifferences(absorb=..., vcov_type in {classical,hc1})`, and + `MultiPeriodDiD(absorb=..., vcov_type in {classical,hc1})` non-clustered SEs sat ~6.5% below + `fixest::feols(vcov="iid"/"hetero")` — even though the reported t-`df` already counted the + absorbed FE (an internal inconsistency). A single scalar rescale + (`(n-k)/(n-k-df_adjustment)`, fail-closed when the full-K residual dof is non-positive) now + aligns the SE's `k` with `K_full`, so the absorb path equals the explicit full-dummy + (`fixed_effects=`) path and fixest. **Clustered** SEs are unchanged (fixest's nested-FE `ssc` + convention already matches with `k_visible`); `hc2`/`hc2_bm` (leverage / Satterthwaite DOF) and + survey vcov are unaffected; `SunAbraham` `hc1` auto-clusters at unit and so keeps its documented + deviation. - **Non-finite degrees of freedom now fail closed to all-NaN inference.** `safe_inference()` previously rejected `df <= 0` but let a non-finite `df` (NaN) through, producing an inconsistent tuple (finite t-stat, NaN p-value/CI). It now returns all-NaN for diff --git a/TODO.md b/TODO.md index 09ef8439..8cd58779 100644 --- a/TODO.md +++ b/TODO.md @@ -111,6 +111,7 @@ 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 | | `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 | | 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 | @@ -128,6 +129,8 @@ Doable in principle, but no current caller and/or explicitly out of paper scope. | Issue | Location | PR | Priority | |-------|----------|----|----------| +| CallawaySantAnna **unbalanced-panel event-study weighting**: a dynamic horizon pooling >1 (g,t) cell weights each cell by its per-cell aggregation weight (valid `n_treated` / `agg_weight`), but R `did::aggte` weights by the fixed cohort probability `pg = n_g/N` (or survey mass) from the group column. Balanced panels coincide exactly (verified ~1e-5, incl. universal zero-reference dilution); unbalanced panels can differ at multi-cell horizons (e.g. a two-real-cell horizon `-0.065` vs R `-0.136`). Pre-existing + general (independent of universal reference cells; documented Deviation from R in REGISTRY). Fix = weight every event-study cell by fixed cohort mass (touches `_aggregate_event_study` + the bootstrap bucket + WIF `pg` consistency); gate on balanced byte-identity + fresh unbalanced R parity. | `staggered_aggregation.py`, `staggered_bootstrap.py` | SE-audit D3 | Low | +| CallawaySantAnna event-study bucket/weight construction is duplicated between the analytical aggregator (`staggered_aggregation.py::_aggregate_event_study`) and the multiplier bootstrap (`staggered_bootstrap.py`): both group (g,t) by `e = t - g`, apply the finite/NaN/reference masks, and read cohort weights. Both already consume the same source-materialized universal reference cells (so they agree), but the bucket logic is copy-pasted. Extract one shared helper returning per-event-time buckets (finite cells, NaN cells, reference flags, cohort weights, combined-IF inputs) used by both. Pure refactor; gate on byte-identical analytical + bootstrap output. | `staggered_aggregation.py`, `staggered_bootstrap.py` | SE-audit D3 | Low | | `StackedDiD` survey re-resolution intra-file dedup (raw-weight extraction ×3, compose-normalize ×3, resolve-on-stacked ×2). The cross-estimator ContinuousDiD/EfficientDiD panel-to-unit collapse consolidation LANDED (#226 shared helpers `ResolvedSurveyDesign.subset_to_units_by_row_idx` / `build_unit_first_row_index`); StackedDiD is deliberately NOT on that path (control units are duplicated across sub-experiments, so it re-resolves at stacked granularity rather than collapsing to one row per unit). The residual is stacked-specific, low value, and touches the numerically-sensitive composed-weight renormalization. Post-filter re-resolution / metadata-recompute unification across the three estimators was assessed and is not warranted — they use genuinely different mechanisms and already delegate to shared `_resolve_survey_for_fit` / `compute_survey_metadata`. | `stacked_did.py` | #226 | Low | | `solve_ols` residual memory floor is inherent to SVD-based lstsq after the PR-E marshalling slim (rust-side allocator high-water 15.32 -> 7.81 GB on the 2.4M x 130 clustered solve; remaining = one fused equilibrated input copy + thin-SVD U transient + vcov scores block + LAPACK gelsd internals on the python path). Further reduction requires the tall-skinny-QR algorithm change - same trade-off as the parked QR-reuse row below (library-wide last-digit perturbation + loses the second independent collinearity detector); the two land together if ever reopened. Diagnosis + measurements in docs/performance-plan.md. | `rust/src/linalg.rs::solve_ols`, `diff_diff/linalg.py` | PR-E | Low | | SunAbraham-family fits are SOLVER-bound after the v3.6.x demeaning speedups: on the county-class shape (3.1k units x 60 months, ~100 interaction columns) `solve_ols` = ~60% of fit time (faer SVD 1.24s + pivoted-QR rank detection 0.70s of 3.27s), and the share grows further under the Rust demean kernel; the pivoted QR alone is ~20% of the whole fit (measured attribution in docs/performance-plan.md). Candidate fix is QR-reuse (one factorization for rank detection + solve). **Parked per the 2026-07 correctness-first decision**: the change perturbs every coefficient in the library at the last digits (golden/parity/bit-identity recapture) and removes the SVD's independent singular-value truncation - a second, different near-collinearity detector that is deliberate defense-in-depth for the no-silent-failures contract (the FE-span guard episode showed borderline cases are real). Re-open only on practitioner demand for wide event studies; any implementation gates on fixest/R end-to-end parity, not just internal consistency. | `linalg.py::solve_ols`, `linalg.py::_detect_rank_deficiency` | PR-C attribution | Medium | diff --git a/diff_diff/estimators.py b/diff_diff/estimators.py index 1aba3041..ad20a483 100644 --- a/diff_diff/estimators.py +++ b/diff_diff/estimators.py @@ -21,6 +21,7 @@ from diff_diff.linalg import ( LinearRegression, + _absorbed_fe_vcov_scale, _expand_vcov_with_nan, compute_r_squared, solve_ols, @@ -1991,6 +1992,29 @@ def _refit_mp_absorb(w_r): if survey_weights is not None and survey_weight_type == "fweight": n_eff_df = int(round(np.sum(survey_weights))) df = n_eff_df - k_effective - n_absorbed_effects + + # Absorbed-FE variance scale (fixest full-K convention): the within- + # transform solve_ols above scales the non-clustered classical/hc1 vcov + # by k_visible, but the correct finite-sample count is + # K_full = k_effective + n_absorbed_effects (matching `df` just above and + # fixest feols(vcov="iid"/"hetero")). Rescale so the SE's k agrees with + # the t-df's. Gated exactly as LinearRegression.fit: clustered SEs keep + # k_visible (fixest ssc nested-FE convention), hc2/hc2_bm use + # leverage/Satterthwaite DOF, survey has its own df. When the full-K + # residual dof is non-positive the helper returns NaN and we void the + # vcov -> NaN inference (fail-closed, per the non-finite-df contract). + if ( + n_absorbed_effects > 0 + and effective_cluster_ids is None + and not _use_survey_vcov + and _fit_vcov_type in ("classical", "hc1") + ): + _fe_scale_mp = _absorbed_fe_vcov_scale(n_eff_df, k_effective, n_absorbed_effects) + if np.isnan(_fe_scale_mp): + vcov = np.full_like(vcov, np.nan) + elif _fe_scale_mp != 1.0: + vcov = vcov * _fe_scale_mp + if resolved_survey is not None and resolved_survey.df_survey is not None: df = resolved_survey.df_survey # Replicate df: rank-deficient → NaN inference; dropped replicates → n_valid-1 diff --git a/diff_diff/linalg.py b/diff_diff/linalg.py index 33183803..6647d370 100644 --- a/diff_diff/linalg.py +++ b/diff_diff/linalg.py @@ -267,6 +267,47 @@ def _expand_coefficients_with_nan( return coef_full +def _absorbed_fe_vcov_scale(n_eff: float, k_eff: int, df_adjustment: int) -> float: + """Finite-sample vcov rescale for absorbed fixed effects (fixest full-K). + + Within-transform (``absorb=``) fits residualize the FE out of the design, so + the classical ``sse/(n-k)`` and HC1 ``n/(n-k)`` variance factors use only the + *visible* regressor count ``k_eff = k_visible``. fixest (and the reported + t-``df``) count the absorbed FE too, i.e. ``K_full = k_visible + + df_adjustment``. This returns the scalar that maps the ``k_visible`` vcov to + the ``K_full`` one:: + + vcov_full = vcov_visible * (n_eff - k_eff) / (n_eff - k_eff - df_adjustment) + + which is algebraically exact for both the classical and HC1 factors (both are + pure scalars in ``1/(n-k)``; the single-coefficient robust variance is + FWL-invariant between the demeaned and full-dummy designs, so only the + ``1/(n-k)`` scalar differs). + + Returns: + - ``1.0`` when ``df_adjustment <= 0`` (no absorbed FE -- a no-op). + - the finite scale when the full-K residual dof + ``(n_eff - k_eff - df_adjustment)`` is positive. + - ``nan`` (fail-closed) when that full-K residual dof (or the visible dof) is + non-positive: the full-K variance is undefined for such a saturated + within-transform design, so callers void the vcov to NaN -> NaN inference + (per the non-finite-df fail-closed contract) rather than leaving a + misleading ``k_visible`` SE in place. + + Callers must gate on non-clustered ``classical``/``hc1``: clustered SEs + follow fixest's ``ssc`` nested-FE convention (FE nested in the cluster are + not counted, so ``k_visible`` already matches for the nested case) and + ``hc2``/``hc2_bm`` use leverage / Satterthwaite DOF -- none must be rescaled. + """ + denom_visible = n_eff - k_eff + denom_full = n_eff - k_eff - df_adjustment + if df_adjustment <= 0: + return 1.0 + if denom_full <= 0 or denom_visible <= 0: + return float("nan") + return denom_visible / denom_full + + def _expand_vcov_with_nan( vcov_reduced: np.ndarray, k_full: int, @@ -3897,6 +3938,29 @@ def fit( n_eff_df = int(round(np.sum(_fit_weights))) elif np.any(_fit_weights == 0): n_eff_df = int(np.count_nonzero(_fit_weights > 0)) + # Absorbed-FE variance scale (fixest full-K convention): the classical + # sse/(n-k) and HC1 n/(n-k) factors computed above use k_visible, but + # with absorbed FE the correct finite-sample count is + # K_full = n_params_effective_ + df_adjustment (the t-df below already + # uses it). Rescale the NON-CLUSTERED iid/hetero vcov so the SE's k + # agrees with the t-df's and with fixest feols(vcov="iid"/"hetero"). + # Clustered SEs keep k_visible (fixest ssc nested-FE convention already + # matches); hc2/hc2_bm use leverage/Satterthwaite DOF; survey has its + # own df; full-dummy fits carry df_adjustment == 0. When the full-K + # residual dof is non-positive the helper returns NaN and we void the + # vcov -> NaN inference (fail-closed, per the non-finite-df contract). + if ( + df_adjustment > 0 + and effective_cluster_ids is None + and not _use_survey_vcov + and _fit_vcov_type in ("classical", "hc1") + ): + _fe_scale = _absorbed_fe_vcov_scale(n_eff_df, self.n_params_effective_, df_adjustment) + if np.isnan(_fe_scale): + self.vcov_ = np.full_like(self.vcov_, np.nan) + elif _fe_scale != 1.0: + self.vcov_ = self.vcov_ * _fe_scale + self.df_ = n_eff_df - self.n_params_effective_ - df_adjustment # Survey degrees of freedom: n_PSU - n_strata (overrides standard df) diff --git a/diff_diff/staggered.py b/diff_diff/staggered.py index 1d4ef247..225476e7 100644 --- a/diff_diff/staggered.py +++ b/diff_diff/staggered.py @@ -5,6 +5,7 @@ including the Callaway-Sant'Anna (2021) estimator. """ +import bisect import warnings from typing import Any, Dict, List, Optional, Tuple @@ -346,15 +347,25 @@ class CallawaySantAnna( - "silent": Drop columns silently without warning base_period : str, default="varying" Method for selecting the base (reference) period for computing - ATT(g,t). Options: - - - "varying": For pre-treatment periods (t < g - anticipation), use - t-1 as base (consecutive comparisons). For post-treatment, use - g-1-anticipation. Requires t-1 to exist in data. - - "universal": Always use g-1-anticipation as base period. - + ATT(g,t). Base periods are selected *positionally* (by the nearest + observed period in the sorted panel), matching R ``did::att_gt`` -- so + on gapped (non-consecutive) grids the base is the nearest observed + period, not literal ``t-1`` / ``g-1``. The pre/post split is on the + current period vs the cohort (``t < g`` -> pre), independent of + anticipation; anticipation only shifts the post/universal base. Options: + + - "varying": pre-treatment (``t < g``) uses the immediately-preceding + observed period as base; post-treatment uses the last observed + pre-treatment period (largest observed ``p`` with + ``p + anticipation < g``). + - "universal": always uses that last observed pre-treatment period as + base. + + On consecutive grids these reduce to ``t-1`` / ``g-1-anticipation``. Both produce identical post-treatment effects. Matches R's - did::att_gt() base_period parameter. + ``did::att_gt()`` on gapped panels (base selection, estimable ATT/SE + cells, the ``"universal"`` zero reference cells, and all aggregations). + See :func:`_select_base_period`. cband : bool, default=True Whether to compute simultaneous confidence bands (sup-t) for event study aggregation. Requires ``n_bootstrap > 0``. @@ -460,8 +471,11 @@ class CallawaySantAnna( where G=g indicates treatment cohort g and C=1 indicates control units. This uses g-1 as the base period, which applies to post-treatment (t >= g). - With base_period="varying" (default), pre-treatment uses t-1 as base for - consecutive comparisons useful in parallel trends diagnostics. + With base_period="varying" (default), pre-treatment uses the immediately- + preceding observed period as base for the consecutive comparisons useful in + parallel trends diagnostics. Base periods are selected positionally (nearest + observed period), matching R did::att_gt on gapped grids (see + ``_select_base_period``). References ---------- @@ -801,6 +815,7 @@ def _precompute_structures( "unit_cohorts": unit_cohorts, "outcome_matrix": outcome_matrix, "period_to_col": period_to_col, + "observed_sorted": sorted(period_to_col), "cohort_masks": cohort_masks, "never_treated_mask": never_treated_mask, "covariate_by_period": covariate_by_period, @@ -816,6 +831,71 @@ def _precompute_structures( ), } + def _select_base_period(self, g: Any, t: Any, observed_sorted: List) -> Optional[Any]: + """Select the base period for cell ``(g, t)``, matching R ``did::att_gt``. + + R selects the base period by *position* in the sorted list of observed + periods, not by literal calendar arithmetic: on gapped (non-consecutive) + period grids the base is the nearest observed period, not ``t-1`` / + ``g-1``. This reproduces ``did`` 2.5.1 ``compute.att_gt`` exactly and is + byte-identical to the old ``t-1`` / ``g-1-anticipation`` rule on + consecutive grids (where positional == calendar). + + Parameters + ---------- + g, t : cohort and evaluation period (calendar values). + observed_sorted : ascending list of the unique observed periods. + + Returns + ------- + The base period value, or ``None`` when no valid earlier observed + period exists (a non-estimable cell, materialized by callers as a + ``missing_period`` NaN; R cannot estimate it either). + + Notes + ----- + Following R ``compute.att_gt``, the pre/post split is on the current + period vs the cohort (``t < g`` -> pre), **independent of + anticipation**; anticipation only enters the post/universal base. + - ``universal``, or post-treatment (``t >= g``): base is the last + pre-treatment observed period, i.e. the largest observed ``p`` with + ``p + anticipation < g``. + - ``varying`` pre-treatment (``t < g``): base is the immediately + preceding observed period, i.e. the largest observed ``p < t``. + """ + if self.base_period == "universal" or t >= g: + threshold = g - self.anticipation + else: # varying pre-treatment + threshold = t + # Largest observed period strictly below `threshold` (positional + # neighbor). bisect_left gives the insertion index of `threshold`; the + # element just before it is the largest observed value < threshold. + idx = bisect.bisect_left(observed_sorted, threshold) + return observed_sorted[idx - 1] if idx > 0 else None + + def _valid_periods_for_group(self, g: Any, time_periods: List, observed_sorted: List) -> List: + """Evaluation periods ``t`` to attempt as ``ATT(g, t)`` for cohort ``g``. + + Centralizes the per-group period filter used by every estimation path + (single source of truth, so the positional-base contract cannot drift): + + - ``universal``: all observed periods except the (positional) base + reference period (which is the trivial ``ATT = 0`` cell); the base is + the last observed pre-treatment period, matching R -- NOT literal + ``g-1-anticipation`` (which on gapped grids would leave the real base + in and materialize a fake zero cell). + - ``varying``: all periods except the earliest observed one (which can + never be a current period -- it has no earlier base). + + Cells that remain non-estimable (``_select_base_period`` returns + ``None``) are pruned downstream as ``missing_period`` NaNs. + """ + if self.base_period == "universal": + universal_base = self._select_base_period(g, g, observed_sorted) + return [t for t in time_periods if t != universal_base] + min_period = observed_sorted[0] + return [t for t in time_periods if t >= g - self.anticipation or t > min_period] + def _compute_att_gt_fast( self, precomputed: PrecomputedData, @@ -855,24 +935,12 @@ def _compute_att_gt_fast( unit_cohorts = precomputed["unit_cohorts"] covariate_by_period = precomputed["covariate_by_period"] - # Base period selection based on mode - if self.base_period == "universal": - # Universal: always use g - 1 - anticipation - base_period_val = g - 1 - self.anticipation - else: # varying - if t < g - self.anticipation: - # Pre-treatment: use t - 1 (consecutive comparison) - base_period_val = t - 1 - else: - # Post-treatment: use g - 1 - anticipation - base_period_val = g - 1 - self.anticipation - - if base_period_val not in period_to_col: - # Base period must exist; no fallback to maintain methodological consistency - return None, 0.0, 0, 0, None, None, "missing_period" + # Base period selection: positional (sorted-index) neighbor, matching + # R did::att_gt (see _select_base_period). Returns None for a + # non-estimable cell (no earlier observed period). + base_period_val = self._select_base_period(g, t, precomputed["observed_sorted"]) - # Check if periods exist in the data - if base_period_val not in period_to_col or t not in period_to_col: + if base_period_val is None or t not in period_to_col: return None, 0.0, 0, 0, None, None, "missing_period" base_col = period_to_col[base_period_val] @@ -1063,25 +1131,16 @@ def _compute_all_att_gt_vectorized( # Collect all valid (g, t, base_col, post_col) tuples tasks = [] for g in treatment_groups: - if self.base_period == "universal": - universal_base = g - 1 - self.anticipation - valid_periods = [t for t in time_periods if t != universal_base] - else: - valid_periods = [ - t for t in time_periods if t >= g - self.anticipation or t > min_period - ] + valid_periods = self._valid_periods_for_group( + g, time_periods, precomputed["observed_sorted"] + ) for t in valid_periods: - # Base period selection - if self.base_period == "universal": - base_period_val = g - 1 - self.anticipation - else: - if t < g - self.anticipation: - base_period_val = t - 1 - else: - base_period_val = g - 1 - self.anticipation + # Base period selection: positional (sorted-index) neighbor, + # matching R did::att_gt (see _select_base_period). + base_period_val = self._select_base_period(g, t, precomputed["observed_sorted"]) - if base_period_val not in period_to_col or t not in period_to_col: + if base_period_val is None or t not in period_to_col: skipped_missing_period.append((g, t)) group_time_effects[(g, t)] = _nan_gt_entry(skip_reason="missing_period") continue @@ -1312,24 +1371,16 @@ def _compute_all_att_gt_covariate_reg( # Collect all valid (g, t) tasks with their base periods tasks_by_group = {} # control_key -> list of (g, t, base_period_val, base_col, post_col) for g in treatment_groups: - if self.base_period == "universal": - universal_base = g - 1 - self.anticipation - valid_periods = [t for t in time_periods if t != universal_base] - else: - valid_periods = [ - t for t in time_periods if t >= g - self.anticipation or t > min_period - ] + valid_periods = self._valid_periods_for_group( + g, time_periods, precomputed["observed_sorted"] + ) for t in valid_periods: - if self.base_period == "universal": - base_period_val = g - 1 - self.anticipation - else: - if t < g - self.anticipation: - base_period_val = t - 1 - else: - base_period_val = g - 1 - self.anticipation + # Base period selection: positional (sorted-index) neighbor, + # matching R did::att_gt (see _select_base_period). + base_period_val = self._select_base_period(g, t, precomputed["observed_sorted"]) - if base_period_val not in period_to_col or t not in period_to_col: + if base_period_val is None or t not in period_to_col: skipped_missing_period.append((g, t)) group_time_effects[(g, t)] = _nan_gt_entry(skip_reason="missing_period") continue @@ -2085,13 +2136,9 @@ def fit( ) for g in treatment_groups: - if self.base_period == "universal": - universal_base = g - 1 - self.anticipation - valid_periods = [t for t in time_periods if t != universal_base] - else: - valid_periods = [ - t for t in time_periods if t >= g - self.anticipation or t > min_period - ] + valid_periods = self._valid_periods_for_group( + g, time_periods, precomputed["observed_sorted"] + ) for t in valid_periods: rc_result = self._compute_att_gt_rc( @@ -2207,13 +2254,9 @@ def fit( ) for g in treatment_groups: - if self.base_period == "universal": - universal_base = g - 1 - self.anticipation - valid_periods = [t for t in time_periods if t != universal_base] - else: - valid_periods = [ - t for t in time_periods if t >= g - self.anticipation or t > min_period - ] + valid_periods = self._valid_periods_for_group( + g, time_periods, precomputed["observed_sorted"] + ) for t in valid_periods: att_gt, se_gt, n_treat, n_ctrl, inf_info, sw_sum, skip_reason = ( @@ -2336,6 +2379,69 @@ def fit( stacklevel=2, ) + # Universal base period: materialize each cohort's positional base as a + # zero reference cell (att=0, se=NaN, zero influence function) so it + # appears in the group-time table AND is weighted into the dynamic + # event-study aggregation and the multiplier bootstrap exactly like R + # `did` (which lists the base as a zero att_gt cell and includes it in + # aggte(type="dynamic")). Injected once here — before every aggregation + # and the bootstrap — so all consumers see the same cells uniformly. The + # `group` / `simple` aggregations filter to post-treatment cells + # (t >= g - anticipation) and so exclude the reference (t = base is the + # last observed period with base + anticipation < g). + if self.base_period == "universal": + _obs_sorted = precomputed["observed_sorted"] + # Each cohort's reference weight is the FIXED cohort mass = R's aggte + # `pg` numerator: the survey-weighted sum (else the unit / observation + # count) over the WHOLE cohort, read straight from the cohort column. + # NOT a single cell's valid-treated count (which under-weights the + # reference on unbalanced panels). Because it needs no estimable cell, + # the reference is materialized for every cohort with a valid base -- + # matching R, which materializes the base zero cell before estimating + # the other cells, even for a cohort whose non-reference cells are all + # non-estimable. RCS cells carry `agg_weight` (fixed mass); panel cells + # read `n_treated`; survey aggregation overrides both with the cohort + # survey mass. + _unit_cohorts = precomputed.get("unit_cohorts") + _survey_w = precomputed.get("survey_weights") + _is_panel = precomputed.get("is_panel", True) + for _g in treatment_groups: + _base = self._select_base_period(_g, _g, _obs_sorted) + if _base is None or (_g, _base) in group_time_effects or _unit_cohorts is None: + continue + _cmask = np.asarray(_unit_cohorts) == _g + _mass = ( + float(np.sum(np.asarray(_survey_w)[_cmask])) + if _survey_w is not None + else float(np.count_nonzero(_cmask)) + ) + if _mass <= 0: + continue # cohort has no units -> nothing to weight + _ref: Dict[str, Any] = { + "effect": 0.0, + "se": np.nan, + "t_stat": np.nan, + "p_value": np.nan, + "conf_int": (np.nan, np.nan), + "n_treated": int(round(_mass)), + "n_control": 0, + "skip_reason": None, + "is_reference": True, + } + if not _is_panel: + _ref["agg_weight"] = _mass + if _survey_w is not None: + _ref["survey_weight_sum"] = _mass + group_time_effects[(_g, _base)] = _ref + influence_func_info[(_g, _base)] = { + "treated_idx": np.array([], dtype=int), + "control_idx": np.array([], dtype=int), + "treated_inf": np.array([]), + "control_inf": np.array([]), + "treated_units": np.array([]), + "control_units": np.array([]), + } + # Compute overall ATT (simple aggregation) overall_att, overall_se, overall_effective_df = self._aggregate_simple( group_time_effects, influence_func_info, df, unit, precomputed @@ -3447,6 +3553,7 @@ def _precompute_structures_rc( "never_treated_mask": never_treated_mask, "time_periods": time_periods, "period_to_col": period_to_col, + "observed_sorted": sorted(period_to_col), "is_balanced": False, "survey_weights": survey_weights_arr, "resolved_survey": resolved_survey, @@ -3506,16 +3613,11 @@ def _compute_att_gt_rc( obs_outcome = precomputed["obs_outcome"] period_to_col = precomputed["period_to_col"] - # Base period selection (same logic as panel) - if self.base_period == "universal": - base_period_val = g - 1 - self.anticipation - else: # varying - if t < g - self.anticipation: - base_period_val = t - 1 - else: - base_period_val = g - 1 - self.anticipation + # Base period selection: positional (sorted-index) neighbor, matching + # R did::att_gt (see _select_base_period). Same logic as the panel path. + base_period_val = self._select_base_period(g, t, precomputed["observed_sorted"]) - if base_period_val not in period_to_col or t not in period_to_col: + if base_period_val is None or t not in period_to_col: return None, 0.0, 0, 0, None, None, None, "missing_period" # Treated mask = cohort g diff --git a/diff_diff/staggered_aggregation.py b/diff_diff/staggered_aggregation.py index b52a7c98..4d3ebbb7 100644 --- a/diff_diff/staggered_aggregation.py +++ b/diff_diff/staggered_aggregation.py @@ -848,6 +848,18 @@ def _aggregate_event_study( ) effects_by_e = balanced_effects + # Universal base period: each cohort's positional base is materialized in + # `group_time_effects` / `influence_func_info` (with a zero effect and a + # zero influence function) by `fit()` before aggregation, so it is already + # grouped into `effects_by_e` above and weighted into the dynamic horizon + # exactly like R `did::aggte(type="dynamic")` (a reference cell dilutes the + # real cells at an overlapping negative horizon). We only flag which cells + # are references so a reference-only horizon reports NaN (not a spurious + # se=0) and does not count toward `n_groups`. + reference_cells: Set[Tuple[Any, Any]] = { + (g, t) for (g, t), data in group_time_effects.items() if data.get("is_reference") + } + # Compute aggregated effects and SEs for all relative periods sorted_periods = sorted(effects_by_e.items()) agg_effects_list = [] @@ -877,10 +889,27 @@ def _aggregate_event_study( # which already drops all-NaN groups. continue + # Reference-only horizon (universal base): every cell is a zero + # reference (att=0, no influence function), so there is no estimated + # effect. Report att=0, se=NaN — matching R `did` (base rows carry + # `se = NA`) — instead of a spurious se=0 from the all-zero IF. + if reference_cells and all(gt in reference_cells for gt in gt_pairs): + agg_effects_list.append(0.0) + agg_ses_list.append(np.nan) + agg_n_groups.append(0) + agg_effective_dfs.append(None) + agg_periods.append(e) + # No influence-function column for a reference-only horizon (it + # carries no estimated effect); leaving it out of _psi_event_times + # keeps the VCV index aligned with the VCV columns (valid_psi). + continue + weights = ns / np.sum(ns) agg_effect = np.sum(weights * effs) - # Compute SE with WIF adjustment (matching R's did::aggte) + # Compute SE with WIF adjustment (matching R's did::aggte). Zero-IF + # reference cells contribute nothing to the variance but their cohort + # weight dilutes the real cells, matching R's dynamic aggregation. groups_for_gt = np.array([g for (g, t) in gt_pairs]) agg_se, psi_e, eff_df = self._compute_aggregated_se_with_wif( gt_pairs, @@ -896,10 +925,10 @@ def _aggregate_event_study( agg_effects_list.append(agg_effect) agg_ses_list.append(agg_se) - # Count only finite-contributing cells (gt_pairs is finite-filtered - # above) so materialized NaN cells don't inflate n_groups — matches - # the all-NaN early-return which already reports 0. - agg_n_groups.append(len(gt_pairs)) + # Count only finite-contributing NON-reference cells so materialized + # NaN cells and zero references don't inflate n_groups — matches the + # all-NaN early-return which already reports 0. + agg_n_groups.append(sum(1 for gt in gt_pairs if gt not in reference_cells)) agg_effective_dfs.append(eff_df) agg_periods.append(e) _psi_vectors.append(psi_e) @@ -947,23 +976,19 @@ def _aggregate_event_study( "n_groups": agg_n_groups[idx], } - # Add reference period for universal base period mode (matches R did package) - if getattr(self, "base_period", "varying") == "universal": - ref_period = -1 - self.anticipation - if event_study_effects and ref_period not in event_study_effects: - event_study_effects[ref_period] = { - "effect": 0.0, - "se": np.nan, - "t_stat": np.nan, - "p_value": np.nan, - "conf_int": (np.nan, np.nan), - "n_groups": 0, - } + # (Universal-mode zero reference rows are now materialized per cohort at + # their positional base event time e = base - g during the aggregation + # above — matching R `did::aggte(type="dynamic")` — rather than as a + # single fixed e = -1-anticipation display row.) # Compute full event-study VCV from per-event-time IF vectors (Phase 7d) # This enables HonestDiD to use the full covariance structure event_study_vcov = None - valid_psi = [p for p in _psi_vectors if len(p) > 0] + # Pair event times with their IF vectors and keep only non-empty psi, so + # the stored VCV index (below) always aligns 1:1 with the VCV columns. + _valid_pairs = [(et, p) for et, p in zip(_psi_event_times, _psi_vectors) if len(p) > 0] + valid_psi = [p for _, p in _valid_pairs] + valid_event_times = [et for et, _ in _valid_pairs] if valid_psi: try: Psi = np.column_stack(valid_psi) # (n_units, n_event_times) @@ -1001,8 +1026,10 @@ def _aggregate_event_study( pass # Fall back to diagonal (None) # Store the event-time index that matches VCV columns (for subsetting - # in HonestDiD when some event times are filtered out) - self._event_study_vcov_index = _psi_event_times if event_study_vcov is not None else None + # in HonestDiD when some event times are filtered out). Uses the + # non-empty-psi event times so the index aligns 1:1 with the VCV columns + # (reference-only and empty-IF horizons never get a column). + self._event_study_vcov_index = valid_event_times if event_study_vcov is not None else None # Attach VCV to self for CallawaySantAnna to pick up self._event_study_vcov = event_study_vcov diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 39a80873..9b32fe98 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -339,6 +339,26 @@ 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` +- **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`). + 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})`, + and `MultiPeriodDiD(absorb=..., vcov_type in {classical,hc1})`. **Clustered** SEs are unchanged + (this fix is gated on `cluster_ids is None`): the clustered CR1 `k_visible` scale matches fixest + for absorbed FE **nested** in the cluster (e.g. unit FE with unit clustering, per fixest's `ssc` + nested-FE convention, which does not count nested FE). **Known limitation (deviation from fixest):** + when an absorbed FE is *not* nested in the cluster (e.g. `absorb=["unit","time"]` clustered by + `unit`, where the time FE are non-nested), fixest counts the non-nested FE in the CR1 + finite-sample denominator, but the current clustered path uses only `k_visible` -- a small, + pre-existing deviation left out of this D4 (non-clustered) scope and tracked in `TODO.md`. This + is distinct from the SunAbraham / Wooldridge `hc1` deviations below (whose event-study / + aggregation paths auto-cluster or use a different k-convention). `hc2`/`hc2_bm` use leverage / + Satterthwaite DOF and are unaffected. The full-dummy (`fixed_effects=`) idiom carries + `df_adjustment == 0` and is unchanged (it already matched fixest). *Edge cases:* - Singleton units/periods are automatically dropped @@ -485,7 +505,7 @@ ATT(g,t) = E[Y_t - Y_{g-1} | G_g=1] - E[Y_t - Y_{g-1} | C=1] ``` where G_g=1 indicates units first treated in period g, and C=1 indicates never-treated. -*Note:* This equation uses g-1 as the base period, which applies to post-treatment effects (t ≥ g) and `base_period="universal"`. With `base_period="varying"` (default), pre-treatment effects use t-1 as base for consecutive comparisons (see Base period selection in Edge cases). +*Note:* This equation uses g-1 as the base period, which applies to post-treatment effects (t ≥ g) and `base_period="universal"`. With `base_period="varying"` (default), pre-treatment effects use the immediately-preceding observed period as base. Base periods are selected *positionally* (nearest observed period), so on gapped grids the base is the nearest observed period rather than literal g-1 / t-1 (see Base period selection in Edge cases). With covariates (doubly robust): ``` @@ -537,6 +557,19 @@ Aggregations: adjustment, matching R's `did::aggte()`. The WIF accounts for uncertainty in estimating group-size aggregation weights. Group aggregation uses equal time weights (deterministic), so WIF is zero. + - **Deviation from R (unbalanced-panel event-study weighting):** within an event-study + (dynamic) horizon that pools MORE THAN ONE (g,t) cell, diff-diff weights each cell by its + per-cell aggregation weight (`agg_weight` / valid `n_treated`), whereas R `did::aggte` + weights by the fixed cohort probability `pg = n_g / N` (or survey mass) computed once from + the group column. On **balanced** panels these coincide exactly (a cell's valid count equals + the cohort mass), so event-study / group / simple aggregates match R to ~1e-5 -- verified, + including the universal zero-reference dilution case. On **unbalanced** panels a cell can + carry fewer valid units than the cohort mass, so a multi-cell horizon's relative weights (and + thus that horizon's ATT/SE) can differ from R -- e.g. on a dropped-unit panel a two-real-cell + horizon was diff-diff `-0.065` vs R `-0.136`. This is a **pre-existing, general** aggregation + convention (it predates and is independent of the universal reference cells -- it applies to + any multi-cell horizon, references or not) and is tracked in `TODO.md`. Single-cell horizons + (the majority) normalize to weight 1 and are unaffected. - Bootstrap: Multiplier bootstrap with Rademacher, Mammen, or Webb weights. Bootstrap perturbs the combined influence function (standard IF + WIF) directly, not just fixed-weight re-aggregation. This correctly propagates weight estimation uncertainty. @@ -618,20 +651,52 @@ The multiplier bootstrap uses random weights w_i with E[w]=0 and Var(w)=1: - Uses NaN when SE is non-finite or zero (matches per-effect and overall t_stat behavior) - Previous behavior (0.0 default) was inconsistent and misleading - Base period selection (`base_period` parameter): - - "varying" (default): Pre-treatment uses t-1 as base (consecutive comparisons). - Post-treatment uses long differences from g-1-anticipation. The parallel trends - assumption for treatment effects is about long-run trends, but pre-treatment tests - only check consecutive periods. - - "universal": ALL effects (pre and post) are long differences from g-1-anticipation. - The parallel trends assumption is about long-run trends. Pre-treatment coefficients - test cumulative divergence from the base period. + - **Positional (sorted-index) selection** (`_select_base_period`): base periods are + selected by *position* in the sorted list of observed periods, not by literal + calendar arithmetic. On consecutive grids this reduces to `t-1` / `g-1-anticipation`; + on **gapped (non-consecutive) grids** (e.g. biennial surveys, skipped years) the base + is the nearest observed period, so R-estimable cells that literal `t-1` / `g-1` would + NaN are estimated. The pre/post split is on the current period vs the cohort + (`t < g` -> pre), **independent of anticipation**; only the post/universal base uses + anticipation. This matches R `did::att_gt()` (verified against a deparse of `did` 2.5.1 + `compute.att_gt`) and resolves the prior internal inconsistency with the library's own + dCDH estimator, which already used positional neighbors. + - "varying" (default): pre-treatment (`t < g`) uses the immediately-preceding observed + period as base (consecutive comparisons); post-treatment uses the last observed + pre-treatment period (largest observed `p` with `p + anticipation < g`) as a long + difference. Pre-treatment tests check nearest-observed-period comparisons. + - "universal": ALL effects (pre and post) are long differences from the last observed + pre-treatment period. Pre-treatment coefficients test cumulative divergence from it. - Both produce identical post-treatment ATT(g,t); differ only pre-treatment - - `anticipation` shifts the base period to g-1-anticipation, which moves it further - from treatment and strengthens the parallel trends assumption. - - Matches R `did::att_gt()` base_period parameter - - **Event study output**: With "universal", includes reference period (e=-1-anticipation) - with effect=0, se=NaN, conf_int=(NaN, NaN). Inference fields are NaN since this is - a normalization constraint, not an estimated effect. Only added when real effects exist. + - `anticipation` shifts the post/universal base to the last observed `p` with + `p + anticipation < g`, moving it further from treatment and strengthening the + parallel trends assumption (it does not change the pre/post split). + - Matches R `did::att_gt()` base_period parameter, including on gapped panels (base + selection, estimable ATT/SE cells, zero reference cells, and all aggregations). + - **Event study output**: With "universal", each cohort's positional base contributes a + zero reference row at `e = base - g` (effect=0, se=NaN, conf_int=(NaN, NaN)); on a + consecutive grid these coincide at `e = -1-anticipation`, on a gapped grid they can fall + at `e = -2, -3, …`. Inference fields are NaN since this is a normalization constraint, not + an estimated effect (see the full-R-parity note below). + - **Universal-mode zero reference cells (full R parity):** with `base_period="universal"`, + R `did::att_gt()` materializes each cohort's base reference period as a zero cell + (`att = 0`, `se = NA`) in its `att_gt` table and includes it in `aggte(type="dynamic")`. + diff-diff now does the same: `fit()` materializes each cohort's positional base as a zero + reference cell (`att = 0`, `se = NaN`, zero influence function, flagged `is_reference`) in + `group_time_effects` / `to_dataframe("group_time")` **at its positional base event time** + (`e = base - g`, which on gapped grids can be `-2`, `-3`, … not just `-1`). Because the + reference is a real cell, every consumer weights it uniformly — the **event-study dynamic + aggregation**, the **multiplier bootstrap**, and `balance_e` — matching R exactly, including + the overlapping-reference case where a cohort's zero base shares an event time with another + cohort's estimated pre-trend cell (the reference correctly dilutes that horizon; verified vs + `did` 2.5.1 on gapped **balanced** universal panels to ~1e-5 for the analytical AND bootstrap + paths). The reference weight is the fixed cohort mass (R's `pg` numerator), so on unbalanced + panels the general per-cell multi-cell-horizon weighting deviation documented above applies to + the reference exactly as it does to any real cell (the reference is not special). + Reference-only horizons report `att = 0, se = NaN`. The zero reference carries no influence + function, so it adds nothing to any variance; the `group` / `simple` aggregations use + post-treatment cells only (`t >= g - anticipation`) and exclude it. Regression-guarded by + `tests/test_csdid_ported.py::TestCSDIDPositionalBasePeriod`. - Base period interaction with Sun-Abraham comparison: - CS with `base_period="varying"` produces different pre-treatment estimates than SA - This is expected: CS uses consecutive comparisons, SA uses fixed reference (e=-1-anticipation) diff --git a/tests/test_csdid_ported.py b/tests/test_csdid_ported.py index 129a00b4..4bbe9e34 100644 --- a/tests/test_csdid_ported.py +++ b/tests/test_csdid_ported.py @@ -1086,7 +1086,16 @@ def test_golden_dynamic_effects(self, golden_values): ), f"Dynamic e={e}: Py={py_att:.4f}, R={r_att:.4f}" def test_golden_non_consecutive_periods(self, golden_values): - """Non-consecutive periods event study matches R.""" + """Non-consecutive periods {1,2,5,7} match R at the per-cell level. + + R's ``did::att_gt`` selects the base period *positionally* (the nearest + observed period), so it estimates every (g,t) cell on a gapped grid. + Since the positional-base fix we reproduce all of them; the ``reg`` + method is linear, so per-cell att AND se match R to ~machine precision + (measured max |dATT|~3e-11, |dSE|~6e-12). Previously we NaN'd 4/9 cells + (calendar t-1 base not observed) and this test only checked the loose + event-study aggregate at abs<0.5. + """ if "non_consecutive_periods" not in golden_values: pytest.skip("Scenario not in golden values") scenario = golden_values["non_consecutive_periods"] @@ -1102,17 +1111,21 @@ def test_golden_non_consecutive_periods(self, golden_values): covariates=["X"], aggregate="event_study", ) - r_dyn = scenario["results"]["dynamic"] - if results.event_study_effects: - for i, e in enumerate(r_dyn["egt"]): - e = int(e) - if e in results.event_study_effects: - py_att = results.event_study_effects[e]["effect"] - r_att = r_dyn["att"][i] - # Non-consecutive periods with small n have higher variance - assert ( - abs(py_att - r_att) < 0.5 - ), f"Nonconsec e={e}: Py={py_att:.4f}, R={r_att:.4f}" + r_gt = scenario["results"]["group_time"] + for i, (g, t) in enumerate(zip(r_gt["group"], r_gt["time"])): + g, t = int(g), int(t) + assert (g, t) in results.group_time_effects + cell = results.group_time_effects[(g, t)] + py_att, py_se = cell["effect"], cell["se"] + assert np.isfinite(py_att) and np.isfinite( + py_se + ), f"Nonconsec (g={g},t={t}) should be estimable (positional base), got NaN" + assert ( + abs(py_att - r_gt["att"][i]) < 1e-7 + ), f"Nonconsec ATT(g={g},t={t}): Py={py_att:.8f}, R={r_gt['att'][i]:.8f}" + assert ( + abs(py_se - r_gt["se"][i]) < 1e-7 + ), f"Nonconsec SE(g={g},t={t}): Py={py_se:.8f}, R={r_gt['se'][i]:.8f}" def test_golden_anticipation(self, golden_values): """Anticipation=1 matches R.""" @@ -1211,7 +1224,18 @@ def test_golden_two_group(self, golden_values): ) def test_golden_fewer_periods(self, golden_values): - """Fewer-periods-than-groups matches R.""" + """Fewer-periods-than-groups (gapped panel {1,3,4,6}) matches R. + + R's ``did::att_gt`` selects the base period *positionally* (the nearest + observed period), so it estimates every (g,t) cell even when the + calendar t-1 / g-1 base is not observed. Since the positional-base fix + we reproduce all 15 cells (previously 7/15 were NaN'd because the + calendar base was absent, and this test skipped them). Assert att AND se + on every committed cell. The ``dr`` method's per-cell agreement with R + is nonlinear-propensity limited (~1e-4 on the base=1 comparisons; the + linear ``reg`` path matches to ~1e-11, pinned by + ``test_golden_non_consecutive_periods``). + """ if "fewer_periods" not in golden_values: pytest.skip("Scenario not in golden values") scenario = golden_values["fewer_periods"] @@ -1230,19 +1254,18 @@ def test_golden_fewer_periods(self, golden_values): r_gt = scenario["results"]["group_time"] for i, (g, t) in enumerate(zip(r_gt["group"], r_gt["time"])): g, t = int(g), int(t) - if (g, t) in results.group_time_effects: - py_att = results.group_time_effects[(g, t)]["effect"] - # Skip cells we materialize as non-estimable (e.g. a gapped panel - # where the base period g-1 is not observed -> missing_period). R - # falls back to an available base and reports a value where our - # impl does not; compare only cells both estimate (R-parity on the - # finite cells, which is what this golden test pins). - if not np.isfinite(py_att): - continue - r_att = r_gt["att"][i] - assert abs(py_att - r_att) < 0.05, ( - f"Fewer periods ATT(g={g},t={t}): " f"Py={py_att:.4f}, R={r_att:.4f}" - ) + assert (g, t) in results.group_time_effects + cell = results.group_time_effects[(g, t)] + py_att, py_se = cell["effect"], cell["se"] + assert np.isfinite(py_att) and np.isfinite( + py_se + ), f"Fewer periods (g={g},t={t}) should be estimable (positional base), got NaN" + assert ( + abs(py_att - r_gt["att"][i]) < 2e-3 + ), f"Fewer periods ATT(g={g},t={t}): Py={py_att:.6f}, R={r_gt['att'][i]:.6f}" + assert ( + abs(py_se - r_gt["se"][i]) < 5e-4 + ), f"Fewer periods SE(g={g},t={t}): Py={py_se:.6f}, R={r_gt['se'][i]:.6f}" def test_golden_zero_pretreatment(self, golden_values): """Zero pre-treatment outcomes match R. @@ -1510,3 +1533,385 @@ def test_golden_ipw_aggregation_se_vs_r_did_251(self, golden_values): grp = results.group_effects[g] assert abs(grp["effect"] - r_att) < 1e-6, f"ipw group ATT(g={g})" assert abs(grp["se"] - r_se) < 1e-6, f"ipw group SE(g={g})" + + +def _make_gapped_panel(seed: int = 0, periods=(1, 3, 4, 6)) -> pd.DataFrame: + """Deterministic gapped (non-consecutive) panel with a unit-level covariate. + + Cohorts: never-treated + treated-at-4 + treated-at-6. On {1, 3, 4, 6} the + calendar base ``g-1`` for cohort 6 (period 5) is unobserved, so the old + calendar rule NaN'd cohort 6's post cells; the positional rule uses the + nearest observed pre-period (4) and estimates them. + """ + rng = np.random.default_rng(seed) + periods = list(periods) + rows = [] + uid = 0 + for cohort, n in [(0, 60), (4, 40), (6, 40)]: + for _ in range(n): + x = float(rng.normal()) + ui = float(rng.normal()) + for p in periods: + treated_now = cohort != 0 and p >= cohort + y = ( + ui + + 0.3 * x + + 0.5 * p + + (1.5 if treated_now else 0.0) + + float(rng.normal(0, 0.5)) + ) + rows.append({"unit": uid, "period": p, "first_treat": cohort, "outcome": y, "X": x}) + uid += 1 + return pd.DataFrame(rows) + + +class TestCSDIDPositionalBasePeriod: + """Positional (sorted-index) base-period selection matching R ``did::att_gt``. + + R selects the base period by *position* in the sorted list of observed + periods, not by literal calendar arithmetic (``t-1`` / ``g-1``). On gapped + (non-consecutive) grids the two differ; on consecutive grids they coincide. + See ``CallawaySantAnna._select_base_period``. These cases were verified + against a deparse of ``did`` 2.5.1 ``compute.att_gt``. + """ + + @pytest.mark.parametrize( + "observed, g, t, base_period, anticipation, expected", + [ + # Varying, gapped {1,3,4,6}: pre = nearest observed < current; + # post = last observed before the cohort. + ([1, 3, 4, 6], 4, 3, "varying", 0, 1), # pre (3<4): nearest <3 -> 1 + ([1, 3, 4, 6], 4, 4, "varying", 0, 3), # post: last <4 -> 3 + ([1, 3, 4, 6], 4, 6, "varying", 0, 3), # post + ([1, 3, 4, 6], 6, 6, "varying", 0, 4), # post: last <6 -> 4 (calendar 5 absent) + # KEY anticipation>0 case: R splits pre/post on current 2 + ([2, 5, 6, 9], 7, 6, "varying", 2, 5), # pre (6<7): nearest <6 -> 5, NOT 2 + ([2, 5, 6, 9], 7, 9, "varying", 2, 2), # post: last < 7-2=5 -> 2 + # Universal: base is the last pre-treatment observed period, + # independent of the current period. + ([1, 3, 4, 6], 6, 1, "universal", 0, 4), + ([1, 3, 4, 6], 6, 6, "universal", 0, 4), + # Non-estimable: earliest observed period as current -> no base. + ([1, 3, 4, 6], 6, 1, "varying", 0, None), + # Cohort treated from the first observed period -> no pre-base. + ([2, 5, 6, 9], 2, 5, "varying", 0, None), + ], + ) + def test_select_base_period_matches_r( + self, observed, g, t, base_period, anticipation, expected + ): + est = CallawaySantAnna(base_period=base_period, anticipation=anticipation) + assert est._select_base_period(g, t, sorted(observed)) == expected + + @pytest.mark.parametrize("method", ["reg", "dr", "ipw"]) + def test_gapped_panel_all_cells_estimated_with_covariates(self, method): + """Every cell R estimates on a gapped grid is finite (no calendar NaN). + + ``reg`` + covariates routes through ``_compute_all_att_gt_covariate_reg``; + ``dr``/``ipw`` route through the fast/vectorized paths. All must use the + shared positional base. Cohort 6's post cell (t=6, calendar base 5 + absent) was ``missing_period`` NaN pre-fix. + """ + data = _make_gapped_panel() + cs = CallawaySantAnna(estimation_method=method) + res = cs.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X"], + ) + for g, t in [(6, 6), (6, 4), (6, 3), (4, 4), (4, 6), (4, 3)]: + cell = res.group_time_effects[(g, t)] + assert np.isfinite(cell["effect"]) and np.isfinite( + cell["se"] + ), f"{method}: (g={g},t={t}) should be estimable on the gapped grid" + assert ( + cell.get("skip_reason") is None + ), f"{method}: (g={g},t={t}) unexpectedly skipped: {cell.get('skip_reason')}" + + def test_gapped_panel_repeated_cross_section_path(self): + """panel=False routes through ``_compute_att_gt_rc`` (the RCS base site). + + Repeated cross-sections: each row is a distinct unit sampled fresh in + its period, so the positional base must be selected there too. + """ + rng = np.random.default_rng(1) + rows = [] + uid = 0 + for p in (1, 3, 4, 6): + for cohort, n in [(0, 50), (4, 40), (6, 40)]: + for _ in range(n): + x = float(rng.normal()) + treated_now = cohort != 0 and p >= cohort + y = ( + 0.3 * x + + 0.5 * p + + (1.5 if treated_now else 0.0) + + float(rng.normal(0, 0.5)) + ) + rows.append( + {"unit": uid, "period": p, "first_treat": cohort, "outcome": y, "X": x} + ) + uid += 1 + data = pd.DataFrame(rows) + cs = CallawaySantAnna(estimation_method="reg", panel=False) + res = cs.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X"], + ) + cell = res.group_time_effects[(6, 6)] + assert np.isfinite(cell["effect"]) and np.isfinite(cell["se"]) + assert cell.get("skip_reason") is None + + def test_universal_gapped_uses_positional_base(self): + """Universal base on {1,3,4,6} for cohort 6 is period 4 (last observed + < g), not calendar 5; cohort-6 cells are estimable.""" + data = _make_gapped_panel() + cs = CallawaySantAnna(estimation_method="dr", base_period="universal") + res = cs.fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X"], + ) + # Every non-base observed period for cohort 6 is estimable. + for t in [1, 3, 6]: + cell = res.group_time_effects[(6, t)] + assert np.isfinite(cell["effect"]), f"universal (g=6,t={t}) NaN" + # The positional base (period 4, not calendar 5) is materialized as a + # zero reference cell (att=0, se=NaN), matching R `did`'s att_gt table + # rather than the stale-calendar-base filter. + ref = res.group_time_effects[(6, 4)] + assert ref["effect"] == 0.0 and np.isnan(ref["se"]) and ref.get("is_reference") + + def test_universal_gapped_materializes_positional_base_cell(self): + """The universal *positional* base is materialized as a zero reference + cell (att=0, se=NaN) in the group-time table across every estimation path + -- the fast (dr/ipw), covariate-reg, and repeated-cross-section paths must + all place it at the positional base (period 4 for cohort 6 on {1,3,4,6}), + matching R `did::att_gt(base_period="universal")`, NOT the calendar base 5. + """ + data = _make_gapped_panel() # cohort 6 positional universal base = period 4 + for method in ["dr", "ipw", "reg"]: + res = CallawaySantAnna(estimation_method=method, base_period="universal").fit( + data, + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X"], + ) + ref = res.group_time_effects.get((6, 4)) + assert ( + ref is not None + and ref["effect"] == 0.0 + and np.isnan(ref["se"]) + and ref.get("is_reference") + ), f"{method}: (6,4) positional base must be a zero reference cell" + # The unobserved calendar base (period 5) is never materialized. + assert (6, 5) not in res.group_time_effects + + # Repeated cross-sections path (unique unit per row). + rng = np.random.default_rng(2) + rows = [] + uid = 0 + for p in (1, 3, 4, 6): + for cohort, n in [(0, 50), (4, 40), (6, 40)]: + for _ in range(n): + x = float(rng.normal()) + tn = cohort != 0 and p >= cohort + rows.append( + { + "unit": uid, + "period": p, + "first_treat": cohort, + "outcome": 0.3 * x + + 0.5 * p + + (1.5 if tn else 0.0) + + float(rng.normal(0, 0.5)), + "X": x, + } + ) + uid += 1 + res = CallawaySantAnna(estimation_method="reg", base_period="universal", panel=False).fit( + pd.DataFrame(rows), + outcome="outcome", + unit="unit", + time="period", + first_treat="first_treat", + covariates=["X"], + ) + ref = res.group_time_effects.get((6, 4)) + assert ( + ref is not None + and ref["effect"] == 0.0 + and np.isnan(ref["se"]) + and ref.get("is_reference") + ), "RCS: (6,4) positional base must be a zero reference cell" + + def test_universal_gapped_event_study_reference_dilution_matches_r(self): + """Universal event-study dynamic aggregation materializes each cohort's + positional zero reference at ``e = base - g`` and weights it into the + horizon, matching R ``did::aggte(type="dynamic")``. + + Periods {1,3,4,5,7}, equal cohorts 5 and 7: cohort 7's base is period 5 + (``e=-2``), which OVERLAPS cohort 5's estimated cell ``(5,3)`` at + ``e=-2``. R dilutes that horizon by averaging the real cell with the zero + reference; with equal cohort sizes the dynamic ATT is exactly half the + real cell. (Previously diff-diff added a single fixed ``e=-1`` display row + that never entered any horizon's weighting, so ``e=-2`` reported the + undiluted real cell.) Reference-only horizons stay ``att=0, se=NaN``, and + the reference cells are omitted from the group-time table. + """ + rng = np.random.default_rng(11) + rows = [] + uid = 0 + for cohort, n in [(0, 60), (5, 40), (7, 40)]: + for _ in range(n): + ui = float(rng.normal()) + for p in (1, 3, 4, 5, 7): + tr = cohort != 0 and p >= cohort + rows.append( + { + "unit": uid, + "period": p, + "first_treat": cohort, + "y": ui + 0.5 * p + (1.5 if tr else 0.0) + float(rng.normal(0, 0.5)), + } + ) + uid += 1 + res = CallawaySantAnna(estimation_method="reg", base_period="universal").fit( + pd.DataFrame(rows), + outcome="y", + unit="unit", + time="period", + first_treat="first_treat", + aggregate="event_study", + ) + es = res.event_study_effects + gt = res.group_time_effects + real = gt[(5, 3)]["effect"] # cohort 5 estimated cell at e=-2 + real_se = gt[(5, 3)]["se"] + # Equal cohorts -> the zero reference (7,5) halves the e=-2 horizon. + assert es[-2]["effect"] == pytest.approx(0.5 * real, rel=1e-6), ( + f"overlap horizon must be diluted by the zero reference: " + f"es[-2]={es[-2]['effect']:.6f} vs 0.5*real={0.5*real:.6f}" + ) + # SE is meaningfully reduced (would equal real_se if the reference were + # missing or placed at a non-overlapping fixed e=-1). + assert 0.4 * real_se < es[-2]["se"] < 0.7 * real_se + # Reference-only horizon (cohort 5's base at e=-1): att=0, se=NaN. + assert es[-1]["effect"] == 0.0 and np.isnan(es[-1]["se"]) + # Both cohorts' positional bases are materialized as zero reference cells + # in the group-time table (matching R's att_gt): (5,4) and (7,5). + for gt_key in [(5, 4), (7, 5)]: + r = gt[gt_key] + assert r["effect"] == 0.0 and np.isnan(r["se"]) and r.get("is_reference") + + def test_universal_reference_weight_ignores_skipped_cells(self): + """A cohort's zero reference is weighted by the FIXED cohort mass (R's + pg = n_g/N), read from the cohort column -- never from a skipped NaN cell + (which carries n_treated=0). Otherwise the reference would get a zero + weight and fail to dilute an overlapping horizon. + + Same overlap panel as the dilution test, but cohort 7's period-1 outcomes + are NaN'd so its first cell (7,1) is a `zero_treated_control` skip. The + (7,5) reference must still carry the cohort's positive fixed mass, so the + e=-2 horizon is still diluted to half the real cell -- for the analytical + AND bootstrap paths. + """ + rng = np.random.default_rng(11) + rows = [] + uid = 0 + for cohort, n in [(0, 60), (5, 40), (7, 40)]: + for _ in range(n): + ui = float(rng.normal()) + for p in (1, 3, 4, 5, 7): + tr = cohort != 0 and p >= cohort + y = ui + 0.5 * p + (1.5 if tr else 0.0) + float(rng.normal(0, 0.5)) + # NaN out cohort 7's earliest period -> (7,1) is a skipped cell. + if cohort == 7 and p == 1: + y = np.nan + rows.append({"unit": uid, "period": p, "first_treat": cohort, "y": y}) + uid += 1 + data = pd.DataFrame(rows) + for n_boot in (0, 199): + res = CallawaySantAnna( + estimation_method="reg", base_period="universal", n_bootstrap=n_boot, seed=5 + ).fit( + data, + outcome="y", + unit="unit", + time="period", + first_treat="first_treat", + aggregate="event_study", + ) + gt = res.group_time_effects + # (7,1) is skipped; (7,5) reference still carries a positive weight. + assert gt[(7, 1)]["skip_reason"] is not None + assert gt[(7, 5)].get("is_reference") and gt[(7, 5)]["n_treated"] > 0 + # e=-2 is still diluted to half the real cell (weight not zeroed by the + # skipped cell). Bootstrap SE differs from analytical but the POINT is + # the same diluted estimand. + assert res.event_study_effects[-2]["effect"] == pytest.approx( + 0.5 * gt[(5, 3)]["effect"], rel=1e-6 + ) + + def test_universal_reference_materialized_for_all_skipped_cohort(self): + """A cohort whose non-reference cells are ALL non-estimable still gets its + base zero reference cell (weighted by the fixed cohort mass), matching R + `did`, which materializes the base zero cell before estimating others. + + Checked end to end through BOTH the analytical and the multiplier-bootstrap + paths -- the latter is where the reference-weight regression originally + surfaced. + """ + rng = np.random.default_rng(11) + rows = [] + uid = 0 + for cohort, n in [(0, 60), (5, 40), (7, 40)]: + for _ in range(n): + ui = float(rng.normal()) + for p in (1, 3, 4, 5, 7): + tr = cohort != 0 and p >= cohort + y = ui + 0.5 * p + (1.5 if tr else 0.0) + float(rng.normal(0, 0.5)) + # NaN out ALL of cohort 7's outcomes -> every (7,t) cell skips. + if cohort == 7: + y = np.nan + rows.append({"unit": uid, "period": p, "first_treat": cohort, "y": y}) + uid += 1 + data = pd.DataFrame(rows) + for n_boot in (0, 199): + res = CallawaySantAnna( + estimation_method="reg", base_period="universal", n_bootstrap=n_boot, seed=5 + ).fit( + data, + outcome="y", + unit="unit", + time="period", + first_treat="first_treat", + aggregate="event_study", + ) + gt = res.group_time_effects + # Every non-reference cohort-7 cell is skipped, yet its base (7,5) + # reference is materialized with the fixed cohort mass (40), not + # dropped -- for the analytical AND bootstrap paths. + assert all(gt[(7, t)]["skip_reason"] is not None for t in (1, 3, 4, 7) if (7, t) in gt) + ref = gt.get((7, 5)) + assert ( + ref is not None + and ref.get("is_reference") + and ref["effect"] == 0.0 + and np.isnan(ref["se"]) + and ref["n_treated"] == 40 # cohort 7's fixed mass (40 units) + ) diff --git a/tests/test_estimators_vcov_type.py b/tests/test_estimators_vcov_type.py index 61053677..39c55473 100644 --- a/tests/test_estimators_vcov_type.py +++ b/tests/test_estimators_vcov_type.py @@ -20,8 +20,9 @@ import pandas as pd import pytest -from diff_diff import SurveyDesign +from diff_diff import SunAbraham, SurveyDesign from diff_diff.estimators import DifferenceInDifferences, MultiPeriodDiD +from diff_diff.linalg import _absorbed_fe_vcov_scale from diff_diff.twfe import TwoWayFixedEffects @@ -2490,3 +2491,225 @@ def test_fit_is_reproducible(self): for p in r1.period_effects: assert r.period_effects[p].effect == r1.period_effects[p].effect assert r.period_effects[p].se == r1.period_effects[p].se + + +def _make_absorb_panel(seed: int = 11, n_units: int = 50, n_periods: int = 6) -> pd.DataFrame: + """Heteroskedastic multi-period panel with a unit-level covariate. + + Heteroskedasticity makes hc1 differ from classical; absorbed unit + time FE + give a nonzero ``df_adjustment`` so the full-K vcov rescale is exercised. + """ + rng = np.random.default_rng(seed) + rows = [] + for u in range(n_units): + treated = int(u >= n_units // 2) + ui = float(rng.normal()) + x = float(rng.normal()) + for t in range(n_periods): + post = int(t >= n_periods // 2) + sd = 0.5 + 1.2 * treated + 0.25 * t + y = ( + ui + + 0.4 * x + + 0.5 * t + + (0.8 if treated and post else 0.0) + + float(rng.normal(0, sd)) + ) + rows.append({"unit": u, "time": t, "post": post, "treated": treated, "x": x, "y": y}) + return pd.DataFrame(rows) + + +class TestAbsorbedFEFullKParity: + """Absorbed-FE variance scale matches the full-dummy (fixest full-K) SE. + + Within-transform (``absorb=``) fits previously scaled the *non-clustered* + classical / hc1 vcov by ``k_visible`` (absorbed FE excluded), sitting ~6.5% + below ``fixest`` feols(vcov="iid"/"hetero") even though the reported t-df + already used ``K_full``. The absorbed-FE vcov rescale + (``linalg._absorbed_fe_vcov_scale``) aligns the SE's k with ``K_full`` so the + absorb path equals the explicit full-dummy path -- the in-repo oracle, + verified equal to ``fixest`` to 6 digits (iid 0.202158 / hetero 0.432290 on + the audit smoke panels). Clustered SEs (fixest nested-FE convention) and + hc2 / hc2_bm (leverage / Satterthwaite DOF) are intentionally unaffected. + """ + + @pytest.mark.parametrize("vcov", ["classical", "hc1"]) + def test_did_absorb_matches_full_dummy_oracle(self, vcov): + df = _make_absorb_panel() + ab = DifferenceInDifferences(vcov_type=vcov).fit( + df, + outcome="y", + treatment="treated", + time="post", + absorb=["unit", "time"], + covariates=["x"], + ) + fe = DifferenceInDifferences(vcov_type=vcov).fit( + df, + outcome="y", + treatment="treated", + time="post", + fixed_effects=["unit", "time"], + covariates=["x"], + ) + assert np.isfinite(ab.se) + np.testing.assert_allclose(ab.att, fe.att, rtol=1e-9) + np.testing.assert_allclose(ab.se, fe.se, rtol=1e-9) + + @pytest.mark.parametrize("vcov", ["classical", "hc1"]) + def test_mpd_absorb_matches_full_dummy_oracle(self, vcov): + df = _make_absorb_panel() + ab = MultiPeriodDiD(vcov_type=vcov).fit( + df, outcome="y", treatment="treated", time="time", absorb=["unit"], covariates=["x"] + ) + fe = MultiPeriodDiD(vcov_type=vcov).fit( + df, + outcome="y", + treatment="treated", + time="time", + fixed_effects=["unit"], + covariates=["x"], + ) + assert np.isfinite(ab.avg_se) + np.testing.assert_allclose(ab.avg_att, fe.avg_att, rtol=1e-9) + np.testing.assert_allclose(ab.avg_se, fe.avg_se, rtol=1e-9) + + def test_twfe_classical_matches_full_dummy_oracle(self): + """TWFE always within-transforms; the equivalent explicit full-dummy + DiD is the oracle. After the fix TWFE(classical) SE == that oracle.""" + df = _make_absorb_panel() + tw = TwoWayFixedEffects(vcov_type="classical").fit( + df, outcome="y", treatment="treated", time="post", unit="unit" + ) + fe = DifferenceInDifferences(vcov_type="classical").fit( + df, outcome="y", treatment="treated", time="post", fixed_effects=["unit", "post"] + ) + np.testing.assert_allclose(tw.att, fe.att, rtol=1e-9) + np.testing.assert_allclose(tw.se, fe.se, rtol=1e-9) + + def test_absorb_cluster_not_rescaled(self): + """The absorbed-FE full-K rescale must NOT touch clustered SEs. + + The rescale is gated on ``cluster_ids is None``, so the cluster-absorb + SE stays at ``k_visible`` and differs from the full-dummy path here. + (Full fixest cluster parity is a *separate*, out-of-scope matter: fixest + counts non-nested absorbed FE in the CR1 denominator, so for + ``absorb=["unit","time"], cluster="unit"`` the non-nested time FE would + need counting -- a documented pre-existing limitation, see REGISTRY; + this test only pins that D4 does not rescale the cluster path.) + """ + df = _make_absorb_panel() + ab = DifferenceInDifferences(vcov_type="hc1", cluster="unit").fit( + df, + outcome="y", + treatment="treated", + time="post", + absorb=["unit", "time"], + covariates=["x"], + ) + fe = DifferenceInDifferences(vcov_type="hc1", cluster="unit").fit( + df, + outcome="y", + treatment="treated", + time="post", + fixed_effects=["unit", "time"], + covariates=["x"], + ) + assert np.isfinite(ab.se) + assert not np.isclose(ab.se, fe.se, rtol=1e-6), ( + "cluster-absorb SE must keep k_visible (nested-FE), not be rescaled to " + "the full-dummy K_full value" + ) + + def test_absorb_hc2_bm_not_rescaled(self): + """hc2_bm auto-routes absorb -> full-dummy and uses Satterthwaite DOF; + the classical/hc1 rescale must not touch it (absorb == full-dummy).""" + df = _make_absorb_panel() + ab = DifferenceInDifferences(vcov_type="hc2_bm").fit( + df, + outcome="y", + treatment="treated", + time="post", + absorb=["unit", "time"], + covariates=["x"], + ) + fe = DifferenceInDifferences(vcov_type="hc2_bm").fit( + df, + outcome="y", + treatment="treated", + time="post", + fixed_effects=["unit", "time"], + covariates=["x"], + ) + np.testing.assert_allclose(ab.se, fe.se, rtol=1e-9) + + def test_sunab_hc1_autoclusters_so_gate_skips(self): + """SunAbraham hc1 auto-clusters at unit, so the D4 rescale (gated on + ``cluster_ids is None``) never fires -> its documented hc1 deviation is + preserved. Invariant: default hc1 == explicit unit cluster.""" + rng = np.random.default_rng(3) + rows = [] + for u in range(40): + cohort = 0 if u < 20 else (3 if u < 30 else 5) + ui = float(rng.normal()) + for t in range(6): + treated_now = cohort != 0 and t >= cohort + rows.append( + { + "unit": u, + "time": t, + "first_treat": cohort, + "y": ui + + 0.5 * t + + (1.2 if treated_now else 0.0) + + float(rng.normal(0, 0.7)), + } + ) + df = pd.DataFrame(rows) + r_def = SunAbraham(vcov_type="hc1").fit( + df, outcome="y", first_treat="first_treat", time="time", unit="unit" + ) + r_cl = SunAbraham(vcov_type="hc1", cluster="unit").fit( + df, outcome="y", first_treat="first_treat", time="time", unit="unit" + ) + np.testing.assert_allclose(r_def.overall_se, r_cl.overall_se, rtol=1e-12) + + def test_absorbed_fe_vcov_scale_fail_closed(self): + """The rescale helper is fail-closed on non-positive full-K residual dof.""" + # Normal: (n-k)/(n-k-adj) = 7/5. + assert _absorbed_fe_vcov_scale(10, 3, 2) == pytest.approx(1.4) + # No absorbed FE -> no-op scale. + assert _absorbed_fe_vcov_scale(10, 3, 0) == 1.0 + # Fail-closed: full-K residual dof <= 0 -> NaN, so the caller voids the + # vcov to NaN inference (rather than leaving a misleading k_visible SE). + assert np.isnan(_absorbed_fe_vcov_scale(10, 3, 8)) # full-K denom = -1 + assert np.isnan(_absorbed_fe_vcov_scale(10, 2, 8)) # full-K denom = 0 + assert np.isnan(_absorbed_fe_vcov_scale(10, 12, 2)) # visible denom < 0 + + def test_absorb_saturated_full_k_df_le_zero_nan_inference(self): + """A saturated within-transform design (full-K residual dof <= 0) yields + NaN SE/inference end-to-end, not a misleading finite k_visible SE. + + 2 units x 5 periods with ``absorb=["unit"]`` drives K_full past n, so the + classical full-K variance is undefined -> the rescale helper returns NaN + and the vcov is voided (fail-closed), even though the point estimate is + still computed. + """ + rng = np.random.default_rng(0) + rows = [] + for u in range(2): + tr = int(u >= 1) + for t in range(5): + rows.append( + { + "unit": u, + "time": t, + "treated": tr, + "y": float(rng.normal()) + 0.5 * tr * (t >= 1), + } + ) + df = pd.DataFrame(rows) + r = MultiPeriodDiD(vcov_type="classical").fit( + df, outcome="y", treatment="treated", time="time", absorb=["unit"] + ) + assert np.isnan(r.avg_se), "saturated full-K design must yield NaN SE (fail-closed)" diff --git a/tests/test_linalg.py b/tests/test_linalg.py index 14c75e18..c07c90d3 100644 --- a/tests/test_linalg.py +++ b/tests/test_linalg.py @@ -1095,7 +1095,15 @@ def test_cluster_ids_in_fit(self, clustered_data): np.testing.assert_allclose(reg1.get_se(1), reg2.get_se(1), rtol=1e-10) def test_df_adjustment(self, simple_data): - """Test degrees of freedom adjustment parameter.""" + """Test degrees of freedom adjustment parameter. + + ``df_adjustment`` accounts for parameters not present in ``X`` (e.g. + absorbed fixed effects). It reduces the reported ``df_`` AND, for the + non-clustered classical / hc1 variance families, rescales the finite- + sample variance to the full parameter count ``n - k - df_adjustment`` + (fixest full-K convention), so the SE's ``k`` agrees with the reported + t-``df`` instead of ignoring the absorbed parameters. + """ X, y, _ = simple_data reg = LinearRegression().fit(X, y) reg_adj = LinearRegression().fit(X, y, df_adjustment=10) @@ -1107,9 +1115,12 @@ def test_df_adjustment(self, simple_data): result = reg.get_inference(1) result_adj = reg_adj.get_inference(1) - # Same coefficient and SE + # Coefficient unchanged; SE is rescaled to the full parameter count, + # larger by sqrt((n-k) / (n-k-df_adjustment)) = sqrt(reg.df_ / reg_adj.df_). assert result.coefficient == result_adj.coefficient - assert result.se == result_adj.se + expected_se = result.se * np.sqrt(reg.df_ / reg_adj.df_) + np.testing.assert_allclose(result_adj.se, expected_se, rtol=1e-12) + assert result_adj.se > result.se # Different df affects p-value and CI (though often slightly) assert result.df != result_adj.df diff --git a/tests/test_methodology_callaway.py b/tests/test_methodology_callaway.py index 2d12dc68..64e5716d 100644 --- a/tests/test_methodology_callaway.py +++ b/tests/test_methodology_callaway.py @@ -493,13 +493,12 @@ def test_group_time_effects_match_r(self, require_r, benchmark_data): """Test individual ATT(g,t) values match R for post-treatment periods. Post-treatment effects (t >= g) should match closely since both - Python and R use g-1 as the base period for these. + Python and R use the last observed pre-treatment period as base. - Pre-treatment effects may differ due to base_period handling: - - Python varying: uses t-1 as base for pre-treatment - - R varying: may handle differently - - We focus on post-treatment where alignment is expected. + Pre-treatment effects also match: since the positional-base fix, Python + selects base periods positionally (nearest observed period) exactly like + R did::att_gt (varying uses the immediately-preceding observed period). + We focus on post-treatment where alignment has always been expected. """ data, csv_path = benchmark_data diff --git a/tests/test_staggered.py b/tests/test_staggered.py index fbe4cc69..b1cdb1b4 100644 --- a/tests/test_staggered.py +++ b/tests/test_staggered.py @@ -2035,9 +2035,9 @@ def test_underdetermined_control_cell_reg_no_crash(self): def test_universal_base_period_anticipation_reg_smoke(self): """reg+cov under base_period="universal" + anticipation=1: every - produced cell has finite inference and the reference period - t = g-1-anticipation is excluded from the grid (no degenerate - zero-IF cell exists to divide by).""" + estimated cell has finite inference, and each cohort's positional base + period is materialized as a zero reference cell (att=0, se=NaN, matching + R `did`'s att_gt table) rather than omitted.""" data = generate_staggered_data_with_covariates(seed=97) with warnings.catch_warnings(): warnings.simplefilter("ignore") @@ -2049,9 +2049,14 @@ def test_universal_base_period_anticipation_reg_smoke(self): ) assert len(res.group_time_effects) > 0 for (g, t), cell in res.group_time_effects.items(): - assert t != g - 2, "universal reference period must be excluded" + if cell.get("is_reference"): + # Zero reference cell: att=0, se=NaN by construction. + assert cell["effect"] == 0.0 and np.isnan(cell["se"]) + continue if cell["skip_reason"] is None: assert np.isfinite(cell["se"]), f"cell ({g},{t})" + # The reference period is now materialized as a zero cell. + assert any(c.get("is_reference") for c in res.group_time_effects.values()) assert np.isfinite(res.overall_se) def test_reg_constant_only_covariate_matches_no_covariate(self): @@ -4960,26 +4965,30 @@ def test_item8_no_warning_when_first_treat_zero(self): assert len(inf_warnings) == 0 def test_item4_consolidated_skip_warning(self): - """Item 4: Consolidated warning when (g,t) cells are skipped.""" + """Item 4: Consolidated warning when (g,t) cells are non-estimable. + + With positional base-period selection a cell is only non-estimable when + no earlier observed period exists. Cohort ``g=2`` treated at the earliest + observed period (periods ``{2,3,4,5}``) has no pre-treatment period, so + R cannot estimate it either -> ``missing_period`` skips + a consolidated + warning. Cohort ``g=4`` (base = observed period 3) is estimable. + """ import warnings - # Two cohorts: g=4 succeeds (base=3 exists), g=6 fails (base=5 - # exists but post periods 7 need base=5 which exists — so use - # g=8 whose base=7 exists). Actually simplest: periods [1,2,4,5] - # with cohort g=4 (base=3 missing → some skips) and g=2 (base=1 - # exists → succeeds). rng = np.random.default_rng(42) n_units = 40 rows = [] for u in range(n_units): - for t in [1, 2, 4, 5]: - # u < 10: never-treated; u < 25: cohort g=2; rest: cohort g=4 + for t in [2, 3, 4, 5]: + # u < 10: never-treated; u < 25: cohort g=2 (treated at the + # earliest observed period -> no pre-period -> skipped); + # rest: cohort g=4 (base = observed 3 -> succeeds) if u < 10: ft = 0 elif u < 25: - ft = 2 # base=1 exists → succeeds + ft = 2 # no earlier observed period -> non-estimable else: - ft = 4 # base=3 missing → skipped + ft = 4 # base=3 exists -> succeeds outcome = rng.standard_normal() + (2.0 if (ft > 0 and t >= ft) else 0.0) rows.append( {