diff --git a/CHANGELOG.md b/CHANGELOG.md index ae86a0c9..2c46655c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -223,6 +223,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 assert between the two (defense in depth). ### Added +- **Opt-in `df_convention="cluster"` inference-df knob (DiD / TWFE / MultiPeriodDiD + + `LinearRegression`).** Clustered analytical fits historically compute t-statistics, + p-values, and CIs at the fitted **residual df** (`n − K_full`), while `fixest`/Stata use + the **cluster df** `G − 1` (the documented clustered-CR1 inference-df deviation, REGISTRY + §TwoWayFixedEffects). `df_convention="cluster"` opts into the Stata/fixest convention: + only the reference t-distribution changes — point estimates, SEs, and t-statistics are + untouched; survey df and per-coefficient Bell-McCaffrey DOF (`hc2_bm`) keep precedence; + unclustered fits are unaffected. Default remains `"residual"`; **the default flips at + v4** (major-version change — it moves every clustered p-value/CI). sklearn-compatible: + `get_params`/`set_params` round-trip with transactional validation. Locked by + `TestDfConvention` (exact `t(G−1)` tail match, precedence ordering, no-op default). - **REGISTRY.md and REPORTING.md are now published on Read the Docs.** The two methodology markdown pages render as in-site Sphinx pages (MyST) under a new "Methodology" toctree section, so cross-references from the API docs use `:doc:` links diff --git a/TODO.md b/TODO.md index 460dc5e6..bda44ece 100644 --- a/TODO.md +++ b/TODO.md @@ -50,7 +50,6 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| -| Clustered-CR1 inference **df convention**: Python uses the residual df (`n − K_full`, e.g. t(148) on the TWFE live-R panel) for clustered t-stats/p-values/CIs while `fixest`/Stata use the cluster df `G − 1` (t(49) same panel) — the common applied recommendation for few-cluster inference. Both are t-distributions; at large \|t\| tails diverge by orders of magnitude (documented as a REGISTRY Note (deviation from R), 2026-07-09). Decide: adopt cluster-df library-wide (moves EVERY clustered p-value/CI — needs a full-suite impact sweep + R parity regen) or keep residual-df as the documented convention. | `diff_diff/linalg.py::LinearRegression.get_inference` + `safe_inference` df-passing callers | SE-audit C7/C8 | Heavy | Medium | --- @@ -113,6 +112,7 @@ Doable in principle, but no current caller and/or explicitly out of paper scope. | Issue | Location | PR | Priority | |-------|----------|----|----------| | Rust-backend CR2 Bell-McCaffrey port (`return_dof` in the Rust vcov dispatch + CR2 algebra) — **premise re-scoped 2026-07-09**: the scores-based DOF + low-rank factored `A_g` changes made the NumPy CR2-BM path BLAS-bound (`O(n_g k²)` per cluster; 4.1s→38ms at n=100k/k=40), so a Rust port buys ~nothing and adds a parity surface. Revisit only if profiling shows CR2-BM hot again. | `rust/src/linalg.rs` | — | Low | +| Clustered-CR1 inference df **default flip to `"cluster"` (G−1) at v4** — the opt-in `df_convention=` knob landed 2026-07 (DiD/TWFE/MPD + LinearRegression; REGISTRY §TwoWayFixedEffects deviation note); the remaining work is the major-version default change (moves every clustered p-value/CI) + migration note + flipping `TestDfConvention`/`test_moderate_t_pins_residual_df_convention` expectations. Also evaluate extending the knob to standalone estimators with CR1-t inference at that time. | `diff_diff/linalg.py::LinearRegression`, `diff_diff/estimators.py`, `diff_diff/twfe.py` | — | Medium | | CallawaySantAnna **unbalanced-panel R parity — LANDED** via `allow_unbalanced_panel=True` (matches R `did::att_gt(allow_unbalanced_panel=TRUE)` / `DRDID::reg_did_rc`: ATT bit-exact on cells AND dynamic aggregation via fixed unit-cohort-mass `pg` + a per-unit WIF; SE up to the documented CR1 `sqrt(G/(G-1))` factor). The earlier "weighting" framing was a mis-diagnosis — on unbalanced panels the dominant divergence from R is the *estimator* (within-cell differencing vs RC-on-pooled-obs), not only the weighting; both are resolved by the flag. The DEFAULT path keeps within-cell differencing as a documented design choice and now emits a `UserWarning` on unbalanced input (no-silent-failures). **Remaining deferred:** `survey_design=` × `allow_unbalanced_panel=` (per-obs vs per-unit weight resolution — currently fail-closed `NotImplementedError`); and covariate / ipw / dr × the flag R-parity verification (the RC path supports them; the committed golden covers `reg` no-cov). | `staggered.py`, `staggered_aggregation.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 | diff --git a/diff_diff/estimators.py b/diff_diff/estimators.py index 9b976ea2..2fb01ba6 100644 --- a/diff_diff/estimators.py +++ b/diff_diff/estimators.py @@ -133,6 +133,19 @@ class DifferenceInDifferences: path — absent ``cluster=``, pure Conley spatial HAC applies. ``survey_design=`` + Conley and ``inference='wild_bootstrap'`` + Conley both raise ``NotImplementedError``. + df_convention : str, default "residual" + Degrees-of-freedom convention for t-statistics, p-values, and CIs on + clustered analytical fits. ``"residual"`` (default) uses the fitted + residual df (``n − K_full``); ``"cluster"`` uses the Stata/fixest + cluster df ``G − 1``. Applies only at the fallback level of the df + resolution: survey df and per-coefficient Bell-McCaffrey DOF + (``vcov_type="hc2_bm"``) are more refined small-sample corrections + and always take precedence. Point estimates, SEs, and t-statistics + are unaffected — only the reference t-distribution changes. Has no + effect on unclustered fits or on ``vcov_type="conley"`` (the combined + Conley+cluster product kernel has no documented ``G − 1`` df + reference and keeps the residual df). The default flips to ``"cluster"`` at + v4 (see the REGISTRY clustered-CR1 inference-df deviation note). Attributes ---------- @@ -198,11 +211,17 @@ def __init__( conley_metric: str = "haversine", conley_kernel: str = "bartlett", conley_lag_cutoff: Optional[int] = None, + df_convention: str = "residual", ): # Resolve vcov_type from the legacy `robust` alias via the shared # helper so __init__ and set_params use identical validation logic. from diff_diff.linalg import resolve_vcov_type + if df_convention not in ("residual", "cluster"): + raise ValueError( + f"df_convention must be 'residual' or 'cluster', got {df_convention!r}" + ) + self.robust = robust self.cluster = cluster self.vcov_type = resolve_vcov_type(robust, vcov_type) @@ -233,6 +252,13 @@ def __init__( # arrays are auto-derived from data[time].values + data[unit].values # at fit-time (panel estimators already take time/unit as column names). self.conley_lag_cutoff = conley_lag_cutoff + # Inference df convention for clustered analytical fits: "residual" + # (default; t/p/CI at the fitted residual df) or "cluster" (the + # Stata/fixest G-1 convention). Survey df and per-coefficient + # Bell-McCaffrey DOF always take precedence over either. The default + # flips to "cluster" at v4 (REGISTRY clustered-CR1 inference-df + # deviation note). + self.df_convention = df_convention self.is_fitted_ = False self.results_ = None @@ -637,6 +663,7 @@ def fit( conley_time=_conley_time_arr, conley_unit=_conley_unit_arr, conley_lag_cutoff=self.conley_lag_cutoff, + df_convention=self.df_convention, ).fit(X, y, df_adjustment=n_absorbed_effects) coefficients = reg.coefficients_ @@ -714,8 +741,12 @@ def _refit_did_absorb(w_r): if survey_metadata is not None: survey_metadata.df_survey = _df_rep if _df_rep > 0 else None t_stat, p_value, conf_int = safe_inference(att, se, alpha=self.alpha, df=_df_rep) + _inference_df_used = float(_df_rep) if _df_rep is not None and _df_rep > 0 else None elif self.inference == "wild_bootstrap" and self.cluster is not None: - # Override with wild cluster bootstrap inference + # Override with wild cluster bootstrap inference (bootstrap + # test-inversion based; no reference t-distribution, so no + # effective inference df). + _inference_df_used = None se, p_value, conf_int, t_stat, vcov, _ = self._run_wild_bootstrap_inference( X, y, residuals, cluster_ids, att_idx ) @@ -728,6 +759,9 @@ def _refit_did_absorb(w_r): t_stat = inference.t_stat p_value = inference.p_value conf_int = inference.conf_int + _inference_df_used = ( + float(inference.df) if inference.df is not None and inference.df > 0 else None + ) r_squared = compute_r_squared(y, residuals) @@ -776,6 +810,8 @@ def _refit_did_absorb(w_r): vcov_type=_fit_vcov_type, cluster_name=self.cluster, conley_lag_cutoff=(self.conley_lag_cutoff if _fit_vcov_type == "conley" else None), + df_convention=self.df_convention, + inference_df=_inference_df_used, ) self._coefficients = coefficients @@ -1080,6 +1116,7 @@ def get_params(self) -> Dict[str, Any]: "conley_metric": self.conley_metric, "conley_kernel": self.conley_kernel, "conley_lag_cutoff": self.conley_lag_cutoff, + "df_convention": self.df_convention, } def set_params(self, **params) -> "DifferenceInDifferences": @@ -1110,6 +1147,11 @@ def set_params(self, **params) -> "DifferenceInDifferences": # `vcov_type` on local variables, then apply all mutations atomically. pending_robust = params.get("robust", self.robust) pending_vcov_type = params.get("vcov_type", self.vcov_type) + pending_df_convention = params.get("df_convention", self.df_convention) + if pending_df_convention not in ("residual", "cluster"): + raise ValueError( + "df_convention must be 'residual' or 'cluster', " f"got {pending_df_convention!r}" + ) # First pass: validate that every incoming key is a known attribute # so we don't partially apply a batch that ends in "Unknown parameter". @@ -2045,6 +2087,41 @@ 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 + _df_cluster_knob_invalid = False + # Opt-in Stata/fixest cluster-df convention (df_convention="cluster"): + # the shared analytical df becomes G - 1 on a clustered fit. Placed + # BEFORE the survey/replicate overrides below (which overwrite df, so + # survey df keeps precedence) and upstream of the per-period BM-DOF + # branch (which wins per coefficient on the hc2_bm path). Mirrors + # LinearRegression's resolution: only positive-weight clusters count + # on a weighted fit. + if ( + self.df_convention == "cluster" + and effective_cluster_ids is not None + and _fit_vcov_type != "conley" + ): + # conley is excluded: the combined Conley+cluster product kernel is + # a diff-diff convention with no documented G-1 df reference (see + # the REGISTRY Conley section); its inference keeps the residual df. + from diff_diff.linalg import effective_cluster_count + + _g_eff_mp = effective_cluster_count(effective_cluster_ids, survey_weights) + if _g_eff_mp <= 1: + # Cluster df G - 1 undefined: fail closed with NaN inference + # (df=0 forces NaN through safe_inference), mirroring + # LinearRegression.get_inference's guard. Unreachable via the + # CR1 vcov path (its validator now counts positive-weight + # clusters for all weight types and raises) — defense-in-depth. + warnings.warn( + "df_convention='cluster' requires at least 2 effective " + f"clusters; got {_g_eff_mp}. Inference fields will be NaN.", + UserWarning, + stacklevel=2, + ) + _df_cluster_knob_invalid = True + df = 0 + else: + df = _g_eff_mp - 1 # Absorbed-FE variance scale (fixest full-K convention): the within- # transform solve_ols above scales the non-clustered classical/hc1 vcov @@ -2081,7 +2158,7 @@ def _refit_mp_absorb(w_r): # Guard: fall back to normal distribution if df is non-positive # Skip for replicate designs — df=0 is intentional for NaN inference - if df is not None and df <= 0 and not _uses_replicate_mp: + if df is not None and df <= 0 and not _uses_replicate_mp and not _df_cluster_knob_invalid: warnings.warn( f"Degrees of freedom is non-positive (df={df}). " "Using normal distribution instead of t-distribution for inference.", @@ -2238,6 +2315,7 @@ def _refit_mp_absorb(w_r): # R-style NA propagation: if ANY post-period effect is NaN, average is undefined effect_arr = np.array(post_effect_values) + _avg_df = None if np.any(np.isnan(effect_arr)): # Some period effects are NaN (unidentified) - cannot compute valid average # This follows R's default behavior where mean(c(1, 2, NA)) returns NA @@ -2308,6 +2386,12 @@ def _refit_mp_absorb(w_r): len(np.unique(effective_cluster_ids)) if effective_cluster_ids is not None else None ), conley_lag_cutoff=(self.conley_lag_cutoff if _fit_vcov_type == "conley" else None), + df_convention=self.df_convention, + inference_df=( + float(_avg_df) + if _avg_df is not None and np.isfinite(_avg_df) and _avg_df > 0 + else None + ), ) self._coefficients = coefficients diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 05d1c7fe..25e037a5 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -60,6 +60,7 @@ DifferenceInDifferences( bootstrap_weights: str = "rademacher", # "rademacher", "webb", or "mammen" seed: int | None = None, # Random seed rank_deficient_action: str = "warn", # "warn", "error", or "silent" + df_convention: str = "residual", # Clustered t/p/CI df: "residual" (n-K, default) or "cluster" (Stata/fixest G-1); survey df + hc2_bm BM-DOF keep precedence; default flips at v4 ) ``` @@ -106,6 +107,7 @@ TwoWayFixedEffects( robust: bool = True, cluster: str | None = None, # Auto-clusters at unit level if None alpha: float = 0.05, + df_convention: str = "residual", # Clustered t/p/CI df: "residual" (default) or "cluster" (G-1); flips at v4 ) ``` @@ -145,6 +147,7 @@ MultiPeriodDiD( robust: bool = True, cluster: str | None = None, alpha: float = 0.05, + df_convention: str = "residual", # Clustered t/p/CI df: "residual" (default) or "cluster" (G-1); flips at v4 ) ``` diff --git a/diff_diff/linalg.py b/diff_diff/linalg.py index 84c4a873..9af39a1a 100644 --- a/diff_diff/linalg.py +++ b/diff_diff/linalg.py @@ -626,6 +626,7 @@ def _solve_ols_rust( """ # Convert cluster_ids to int64 for Rust (handles string/categorical IDs) if cluster_ids is not None: + _validate_cluster_ids(cluster_ids) cluster_ids = _factorize_cluster_ids(cluster_ids) # Call Rust backend with fallback on numerical instability @@ -1510,6 +1511,50 @@ def _validate_vcov_args( ) +def _validate_cluster_ids(cluster_ids: np.ndarray) -> None: + """Front-door check shared by every clustered vcov backend: missing + (NaN/None) cluster labels are rejected outright. The pandas groupby + aggregating cluster scores drops NaN-labelled rows while np.unique / + factorize-based counts keep (or sentinel) them, so no count can agree + with the meat's partition — silently wrong CR1 SEs. Matches R fixest / + Stata, which error on missing cluster values. Called by + ``compute_robust_vcov``, ``_solve_ols_rust``, and the CR2-BM shared + core so NumPy, Rust, and CR2 routes all fail closed identically.""" + if pd.isna(np.asarray(cluster_ids)).any(): + raise ValueError( + "cluster_ids contain missing values (NaN/None). Drop or " + "impute those rows before requesting cluster-robust SEs." + ) + + +def effective_cluster_count(cluster_ids: np.ndarray, weights: Optional[np.ndarray] = None) -> int: + """Effective cluster count for clustered inference metadata. + + Unweighted: the number of unique cluster labels. Weighted: only clusters + with positive total weight count (zero-weight rows are inert per the + linalg contract). The positive-total-weight definition applies to ALL + weight types (pweight/aweight/fweight alike — an all-zero-weight cluster + contributes nothing to any weighted sandwich), which can be STRICTER + than the raw-unique count some vcov validators use; consumers that need + a defined cluster df must fail closed when this count is < 2 (see + ``get_inference``'s df_convention="cluster" guard). Grouped reduction + via factorize + bincount, O(n). + + Missing (NaN/None) cluster labels raise ``ValueError`` — the CR1 meat + aggregation (a pandas groupby) drops NaN-labelled rows, so no count + can agree with the sandwich's partition; the vcov validation rejects + them for the same reason (matching R fixest / Stata, which error on + missing cluster values). + """ + arr = np.asarray(cluster_ids) + _validate_cluster_ids(arr) + codes, uniques = pd.factorize(arr, sort=False) + if weights is None: + return int(len(uniques)) + sums = np.bincount(codes, weights=np.asarray(weights, dtype=np.float64)) + return int(np.sum(sums > 0)) + + def resolve_vcov_type( robust: bool = True, vcov_type: Optional[str] = None, @@ -1784,6 +1829,7 @@ def compute_robust_vcov( cluster_ids_int = None if cluster_ids is not None: + _validate_cluster_ids(cluster_ids) cluster_ids_int = pd.factorize(cluster_ids)[0].astype(np.int64) try: @@ -1950,6 +1996,7 @@ def _compute_cr2_bm_vcov_and_dof( """ n, k = X.shape cluster_ids_arr = np.asarray(cluster_ids) + _validate_cluster_ids(cluster_ids_arr) unique_clusters = np.unique(cluster_ids_arr) # When weights are provided, enforce subpopulation invariance: zero-weight # rows must contribute nothing to the sandwich. The earlier "drop zero- @@ -3049,8 +3096,14 @@ def _compute_robust_vcov_numpy( # weighted groupby below cannot index-align a pandas Series grouper # against the freshly-created Series(weights) and miscount clusters. cluster_ids_arr = np.asarray(cluster_ids) + _validate_cluster_ids(cluster_ids_arr) n_clusters_check = len(np.unique(cluster_ids_arr)) - if weights is not None and weight_type != "fweight" and np.any(weights == 0): + # Zero-total-weight clusters are inert for ALL weight types (a + # zero-frequency fweight cluster contributes nothing to the sandwich, + # exactly like a subpopulation-zeroed pweight cluster) — the prior + # fweight carve-out let a one-effective-cluster fweight fit through + # to a degenerate CR1 SE. + if weights is not None and np.any(weights == 0): cluster_weight_sums = pd.Series(weights).groupby(cluster_ids_arr).sum() n_clusters_check = int((cluster_weight_sums > 0).sum()) if n_clusters_check < 2: @@ -3097,8 +3150,9 @@ def _compute_robust_vcov_numpy( unique_clusters = np.unique(cluster_ids) n_clusters = len(unique_clusters) - # Exclude clusters with zero total weight (subpopulation-zeroed) - if weights is not None and weight_type != "fweight" and np.any(weights == 0): + # Exclude clusters with zero total weight (subpopulation-zeroed + # pweight/aweight AND zero-frequency fweight alike — inert either way) + if weights is not None and np.any(weights == 0): cluster_weights = pd.Series(weights).groupby(cluster_ids).sum() n_clusters = int((cluster_weights > 0).sum()) @@ -3754,6 +3808,15 @@ class LinearRegression: neither is formally PSD-guaranteed in the radial pairwise form (Conley 1999's explicit PSD Bartlett formula is the 2-D separable product window, Eq 3.14, not the 1-D radial pairwise form). + df_convention : {"residual", "cluster"}, default "residual" + Degrees-of-freedom convention for ``get_inference`` t/p/CI on + clustered fits. ``"residual"`` uses the fitted residual df; + ``"cluster"`` uses the Stata/fixest cluster df ``G − 1`` (from + ``n_clusters_``). Fallback-level only: survey df and per-coefficient + Bell-McCaffrey DOF always take precedence. No effect on + coefficients, SEs, unclustered fits, or ``vcov_type="conley"`` (no + documented ``G − 1`` reference for the Conley+cluster product + kernel). Default flips at v4. Attributes ---------- @@ -3772,6 +3835,10 @@ class LinearRegression: n_params_effective_ : int Effective number of parameters after dropping linearly dependent columns. Equals n_params_ for full-rank matrices (available after fit). + n_clusters_ : int or None + Effective cluster count on a clustered fit (positive-weight clusters + only when weighted); None on unclustered fits. Feeds the + ``df_convention="cluster"`` inference df (available after fit). df_ : int Degrees of freedom (n - n_params_effective) (available after fit). @@ -3820,7 +3887,12 @@ def __init__( conley_time: Optional[np.ndarray] = None, conley_unit: Optional[np.ndarray] = None, conley_lag_cutoff: Optional[int] = None, + df_convention: str = "residual", ): + if df_convention not in ("residual", "cluster"): + raise ValueError( + f"df_convention must be 'residual' or 'cluster', got {df_convention!r}" + ) self.include_intercept = include_intercept self.robust = robust self.cluster_ids = cluster_ids @@ -3839,6 +3911,15 @@ def __init__( self.conley_time = conley_time self.conley_unit = conley_unit self.conley_lag_cutoff = conley_lag_cutoff + # Inference df convention for clustered analytical fits: + # "residual" (default) keeps the fitted residual df `n - K_full` for + # t/p/CI; "cluster" uses the Stata/fixest cluster df `G - 1` instead. + # Applies ONLY at the fallback level of the `get_inference` df + # resolution: survey df and per-coefficient Bell-McCaffrey DOF (which + # are more refined small-sample corrections) always take precedence. + # The default flips to "cluster" at v4 (see REGISTRY clustered-CR1 + # inference-df deviation note). + self.df_convention = df_convention # Resolve vcov_type from the legacy `robust` alias via the shared helper. self.vcov_type = resolve_vcov_type(robust, vcov_type) # Preserve the raw constructor arg (possibly None) so `fit()` can @@ -3864,6 +3945,9 @@ def __init__( self.n_params_effective_: Optional[int] = None self.df_: Optional[int] = None self.survey_df_: Optional[int] = None + # Effective cluster count on a clustered fit (None otherwise); + # feeds the df_convention="cluster" G-1 inference df. + self.n_clusters_: Optional[int] = None # Per-coefficient Bell-McCaffrey DOF vector when vcov_type="hc2_bm". # None for all other vcov_types; preserves df_ as the fallback. self._bm_dof: Optional[np.ndarray] = None @@ -4271,6 +4355,14 @@ def fit( self.df_ = n_eff_df - self.n_params_effective_ - df_adjustment + # Effective cluster count for the df_convention="cluster" inference + # df (G - 1). Mirrors the vcov path's convention: on a weighted fit + # only clusters with positive total weight count (zero-weight rows + # are inert per the linalg contract). + self.n_clusters_ = None + if effective_cluster_ids is not None: + self.n_clusters_ = effective_cluster_count(effective_cluster_ids, _fit_weights) + # Survey degrees of freedom: n_PSU - n_strata (overrides standard df) self.survey_df_ = None if _effective_survey_design is not None: @@ -4496,6 +4588,43 @@ def get_inference( stacklevel=2, ) effective_df = 0 # Forces NaN from t-distribution + elif ( + self.df_convention == "cluster" + and self.n_clusters_ is not None + and self.vcov_type != "conley" + ): + # Opt-in Stata/fixest cluster-df convention for clustered + # analytical fits: t/p/CI at df = G - 1 instead of the residual + # df. conley is excluded: the combined Conley+cluster product + # kernel is a diff-diff convention with no documented G-1 df + # reference (REGISTRY Conley section); it keeps the residual df. + # Deliberately the LAST branch before the residual fallback: + # survey df and per-coefficient Bell-McCaffrey DOF are more + # refined small-sample corrections and always win. Default flips + # at v4 (REGISTRY clustered-CR1 inference-df deviation note). + if self.n_clusters_ <= 1: + # Cluster df G - 1 is undefined for an effectively + # one-cluster fit (e.g. weighted fit where only one cluster + # carries positive weight). Fail closed with NaN inference + # (df=0 forces NaN through safe_inference) instead of + # silently degrading to normal-theory inference. + warnings.warn( + "df_convention='cluster' requires at least 2 effective " + f"clusters; got {self.n_clusters_}. Inference fields " + "will be NaN.", + UserWarning, + stacklevel=2, + ) + return InferenceResult( + coefficient=coef, + se=se, + t_stat=float("nan"), + p_value=float("nan"), + conf_int=(float("nan"), float("nan")), + df=None, + alpha=effective_alpha, + ) + effective_df = self.n_clusters_ - 1 else: effective_df = self.df_ diff --git a/diff_diff/results.py b/diff_diff/results.py index 2a611876..5be4fc37 100644 --- a/diff_diff/results.py +++ b/diff_diff/results.py @@ -150,6 +150,15 @@ class DiDResults: vcov_type: Optional[str] = field(default=None) cluster_name: Optional[str] = field(default=None) conley_lag_cutoff: Optional[int] = field(default=None) + # Inference df metadata: the df convention the estimator was configured + # with ("residual" | "cluster") and the effective df the reported t/p/CI + # actually used. Set only when finite and strictly positive — None under + # wild bootstrap (test-inversion), normal-distribution fallback, and NaN + # inference (incl. the df=0 replicate / invalid-cluster-df sentinels). + # Lets consumers audit which reference distribution produced the + # p-value/CI without refitting. + df_convention: Optional[str] = field(default=None) + inference_df: Optional[float] = field(default=None) def __repr__(self) -> str: """Concise string representation.""" @@ -308,6 +317,10 @@ def to_dict(self) -> Dict[str, Any]: result["n_clusters"] = self.n_clusters if self.p_val_type is not None: result["p_val_type"] = self.p_val_type + if self.df_convention is not None: + result["df_convention"] = self.df_convention + if self.inference_df is not None: + result["inference_df"] = self.inference_df if self.survey_metadata is not None: sm = self.survey_metadata result["weight_type"] = sm.weight_type @@ -689,6 +702,15 @@ class MultiPeriodDiDResults: vcov_type: Optional[str] = field(default=None) cluster_name: Optional[str] = field(default=None) conley_lag_cutoff: Optional[int] = field(default=None) + # Inference df metadata: the df convention the estimator was configured + # with ("residual" | "cluster") and the effective df the reported t/p/CI + # actually used. Set only when finite and strictly positive — None under + # wild bootstrap (test-inversion), normal-distribution fallback, and NaN + # inference (incl. the df=0 replicate / invalid-cluster-df sentinels). + # Lets consumers audit which reference distribution produced the + # p-value/CI without refitting. + df_convention: Optional[str] = field(default=None) + inference_df: Optional[float] = field(default=None) # --- Inference-field aliases (balance/external-adapter compatibility) --- @property @@ -949,6 +971,10 @@ def to_dict(self) -> Dict[str, Any]: result["cluster_name"] = self.cluster_name if self.conley_lag_cutoff is not None: result["conley_lag_cutoff"] = self.conley_lag_cutoff + if self.df_convention is not None: + result["df_convention"] = self.df_convention + if self.inference_df is not None: + result["inference_df"] = self.inference_df # Add period-specific effects for period, pe in self.period_effects.items(): diff --git a/diff_diff/twfe.py b/diff_diff/twfe.py index b98747a2..33c628f2 100644 --- a/diff_diff/twfe.py +++ b/diff_diff/twfe.py @@ -56,6 +56,12 @@ class TwoWayFixedEffects(DifferenceInDifferences): auto-cluster is also preserved (routes to CR2-BM at unit). alpha : float, default=0.05 Significance level for confidence intervals. + df_convention : str, default "residual" + Inherited from :class:`DifferenceInDifferences`: df convention for + t/p/CI on clustered analytical fits — ``"residual"`` (fitted + residual df, default) or ``"cluster"`` (Stata/fixest ``G − 1``). + Survey df and per-coefficient Bell-McCaffrey DOF keep precedence; + inert on unclustered and Conley fits. Default flips at v4. Notes ----- @@ -525,6 +531,7 @@ def fit( # type: ignore[override] conley_time=_conley_time_arr, conley_unit=_conley_unit_arr, conley_lag_cutoff=self.conley_lag_cutoff, + df_convention=self.df_convention, ).fit(X, y, df_adjustment=df_adjustment) else: # Suppress generic warning, TWFE provides context-specific messages below @@ -546,6 +553,7 @@ def fit( # type: ignore[override] conley_time=_conley_time_arr, conley_unit=_conley_unit_arr, conley_lag_cutoff=self.conley_lag_cutoff, + df_convention=self.df_convention, ).fit(X, y, df_adjustment=df_adjustment) coefficients = reg.coefficients_ @@ -662,8 +670,12 @@ def _refit_twfe(w_r): if survey_metadata is not None: survey_metadata.df_survey = _df_rep if _df_rep > 0 else None t_stat, p_value, conf_int = _safe_inf(att, se, alpha=self.alpha, df=_df_rep) + _inference_df_used = float(_df_rep) if _df_rep is not None and _df_rep > 0 else None elif self.inference == "wild_bootstrap": - # Override with wild cluster bootstrap inference + # Override with wild cluster bootstrap inference (bootstrap + # test-inversion based; no reference t-distribution, so no + # effective inference df). + _inference_df_used = None se, p_value, conf_int, t_stat, vcov, _ = self._run_wild_bootstrap_inference( X, y, residuals, cluster_ids, att_idx ) @@ -675,6 +687,9 @@ def _refit_twfe(w_r): t_stat = inference.t_stat p_value = inference.p_value conf_int = inference.conf_int + _inference_df_used = ( + float(inference.df) if inference.df is not None and inference.df > 0 else None + ) # Count observations treated_units = data[data[treatment] == 1][unit].unique() @@ -754,6 +769,8 @@ def _refit_twfe(w_r): vcov_type=_fit_vcov_type, cluster_name=_twfe_cluster_label, conley_lag_cutoff=(self.conley_lag_cutoff if _fit_vcov_type == "conley" else None), + df_convention=self.df_convention, + inference_df=_inference_df_used, ) self.is_fitted_ = True diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 743315a1..b9cbf7e1 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -369,9 +369,16 @@ This matches the behavior of R's `fixest::feols()` with absorbed FE. numerically zero, which is what `tests/test_methodology_twfe.py::test_pvalue_matches_r_twfe` pins; the convention itself is locked at a distinguishable moderate |t|~1.8 by `tests/test_methodology_twfe.py::test_moderate_t_pins_residual_df_convention`, where the - t(148) and t(49) tails differ by ~5%). The SE itself carries the separate ~0.25% non-nested-FE ssc band above. Whether to - adopt the cluster-df convention library-wide is tracked as an Actionable TODO row (it is the - common applied recommendation, but changing it moves every clustered p-value/CI). The full-dummy (`fixed_effects=`) idiom carries + t(148) and t(49) tails differ by ~5%). The SE itself carries the separate ~0.25% non-nested-FE ssc band above. **Opt-in knob (2026-07):** + `df_convention="cluster"` on `DifferenceInDifferences` / `TwoWayFixedEffects` / `MultiPeriodDiD` + (and the `LinearRegression` linalg surface) switches clustered analytical t/p/CI to the + Stata/fixest `G − 1` convention — fallback-level only (survey df and per-coefficient + Bell-McCaffrey DOF keep precedence), point estimates/SEs/t-statistics unchanged, inert on + unclustered fits. The default remains `"residual"` and flips at v4 (tracked in TODO Deferred → + Parked; flipping moves every clustered p-value/CI, so it is a major-version change). Standalone + estimators (CS, SA, imputation-family, etc.) carry their own inference stacks and are out of the + knob's scope. Locked by + `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). *Edge cases:* diff --git a/tests/test_estimators_vcov_type.py b/tests/test_estimators_vcov_type.py index 23be4857..d869cce1 100644 --- a/tests/test_estimators_vcov_type.py +++ b/tests/test_estimators_vcov_type.py @@ -2755,3 +2755,344 @@ def test_absorb_saturated_full_k_df_le_zero_nan_inference(self): 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)" + + +class TestDfConvention: + """`df_convention="cluster"` (opt-in Stata/fixest G-1 inference df on + clustered analytical fits; REGISTRY clustered-CR1 inference-df deviation + note). The knob changes ONLY the reference t-distribution for t/p/CI — + never the point estimate, SE, or t-statistic — and sits at the fallback + level of the df resolution (survey df and per-coefficient Bell-McCaffrey + DOF always win). Default stays "residual" until the v4 flip.""" + + @staticmethod + def _clustered_panel(n_units=50, seed=42): + rng = np.random.default_rng(seed) + rows = [] + for i in range(n_units): + group = int(i < n_units // 2) + for t in range(4): + post = int(t >= 2) + rows.append( + { + "unit": i, + "post": post, + "group": group, + "time": t, + "y": 1 + + 0.3 * i / n_units + + 0.2 * t + + 0.25 * (group * post) + + rng.normal(0, 1.0), + } + ) + return pd.DataFrame(rows) + + def test_default_residual_is_noop(self): + """Explicit df_convention="residual" is byte-identical to the default.""" + from scipy import stats + + data = self._clustered_panel() + r0 = DifferenceInDifferences(cluster="unit").fit( + data, outcome="y", treatment="group", time="post", unit="unit" + ) + r1 = DifferenceInDifferences(cluster="unit", df_convention="residual").fit( + data, outcome="y", treatment="group", time="post", unit="unit" + ) + assert (r0.att, r0.se, r0.t_stat, r0.p_value, r0.conf_int) == ( + r1.att, + r1.se, + r1.t_stat, + r1.p_value, + r1.conf_int, + ) + # and the default really is the residual-df tail, not G-1 + assert abs(r0.p_value - 2 * stats.t.sf(abs(r0.t_stat), 49)) > 1e-4 + + @pytest.mark.parametrize("estimator_cls", [DifferenceInDifferences]) + def test_cluster_convention_matches_g_minus_1(self, estimator_cls): + """Knob on: p/CI follow t(G-1) exactly; att/SE/t unchanged.""" + from scipy import stats + + data = self._clustered_panel() + kw = dict(outcome="y", treatment="group", time="post", unit="unit") + r0 = estimator_cls(cluster="unit").fit(data, **kw) + r1 = estimator_cls(cluster="unit", df_convention="cluster").fit(data, **kw) + assert r0.att == r1.att and r0.se == r1.se and r0.t_stat == r1.t_stat + G = data["unit"].nunique() + np.testing.assert_allclose(r1.p_value, 2 * stats.t.sf(abs(r1.t_stat), G - 1), rtol=1e-12) + crit = stats.t.ppf(0.975, G - 1) + np.testing.assert_allclose( + r1.conf_int, (r1.att - crit * r1.se, r1.att + crit * r1.se), rtol=1e-12 + ) + assert r0.p_value != r1.p_value + + def test_twfe_inherits_knob(self): + from scipy import stats + + from diff_diff import TwoWayFixedEffects + + data = self._clustered_panel() + data["treated"] = data["group"] * data["post"] + kw = dict(outcome="y", treatment="treated", time="post", unit="unit") + r0 = TwoWayFixedEffects(robust=True).fit(data, **kw) + r1 = TwoWayFixedEffects(robust=True, df_convention="cluster").fit(data, **kw) + assert r0.se == r1.se and r0.t_stat == r1.t_stat + G = data["unit"].nunique() + np.testing.assert_allclose(r1.p_value, 2 * stats.t.sf(abs(r1.t_stat), G - 1), rtol=1e-12) + + def test_mpd_period_effects_follow_convention(self): + from scipy import stats + + data = self._clustered_panel() + data["treated"] = data["group"] + kw = dict(outcome="y", treatment="treated", time="time", unit="unit") + m0 = MultiPeriodDiD(cluster="unit").fit(data, **kw) + m1 = MultiPeriodDiD(cluster="unit", df_convention="cluster").fit(data, **kw) + G = data["unit"].nunique() + for period, pe1 in m1.period_effects.items(): + pe0 = m0.period_effects[period] + assert pe0.se == pe1.se + if np.isfinite(pe1.p_value): + np.testing.assert_allclose( + pe1.p_value, + 2 * stats.t.sf(abs(pe1.effect / pe1.se), G - 1), + rtol=1e-12, + ) + + def test_mpd_avg_att_follows_convention(self): + """MPD's headline avg_att p-value/CI route through the shared df, + so the knob moves them to t(G-1); avg_att/avg_se/avg_t_stat are + unchanged. Also locks the new inference_df/df_convention metadata.""" + from scipy import stats + + data = self._clustered_panel() + data["treated"] = data["group"] + kw = dict(outcome="y", treatment="treated", time="time", unit="unit") + m0 = MultiPeriodDiD(cluster="unit").fit(data, **kw) + m1 = MultiPeriodDiD(cluster="unit", df_convention="cluster").fit(data, **kw) + G = data["unit"].nunique() + assert m0.avg_att == m1.avg_att + assert m0.avg_se == m1.avg_se + assert m0.avg_t_stat == m1.avg_t_stat + np.testing.assert_allclose( + m1.avg_p_value, 2 * stats.t.sf(abs(m1.avg_t_stat), G - 1), rtol=1e-12 + ) + crit = stats.t.ppf(0.975, G - 1) + np.testing.assert_allclose( + m1.avg_conf_int, + (m1.avg_att - crit * m1.avg_se, m1.avg_att + crit * m1.avg_se), + rtol=1e-12, + ) + assert m0.avg_p_value != m1.avg_p_value + # result metadata + assert m1.df_convention == "cluster" and m1.inference_df == G - 1 + assert m0.df_convention == "residual" and m0.inference_df > G - 1 + assert m1.to_dict()["inference_df"] == G - 1 + + def test_results_metadata_did_twfe(self): + """DiDResults carries df_convention + the effective inference_df, + included in to_dict() only when set.""" + data = self._clustered_panel() + r1 = DifferenceInDifferences(cluster="unit", df_convention="cluster").fit( + data, outcome="y", treatment="group", time="post", unit="unit" + ) + G = data["unit"].nunique() + assert r1.df_convention == "cluster" and r1.inference_df == G - 1 + d = r1.to_dict() + assert d["df_convention"] == "cluster" and d["inference_df"] == G - 1 + + def test_conley_fits_excluded_from_knob(self): + """vcov_type="conley" + explicit cluster keeps the residual df even + with the knob on (no documented G-1 reference for the product + kernel) — inference identical with the knob on and off.""" + rng = np.random.default_rng(9) + data = self._clustered_panel() + units = data["unit"].unique() + lat = {u: 40 + rng.uniform(-2, 2) for u in units} + lon = {u: -100 + rng.uniform(-2, 2) for u in units} + data["lat"] = data["unit"].map(lat) + data["lon"] = data["unit"].map(lon) + kw = dict(outcome="y", treatment="group", time="post", unit="unit") + common = dict( + vcov_type="conley", + cluster="unit", + conley_coords=("lat", "lon"), + conley_cutoff_km=100.0, + conley_lag_cutoff=0, + ) + r0 = DifferenceInDifferences(**common).fit(data, **kw) + r1 = DifferenceInDifferences(**common, df_convention="cluster").fit(data, **kw) + assert (r0.p_value, r0.conf_int) == (r1.p_value, r1.conf_int) + + def test_bm_dof_precedence_over_knob(self): + """hc2_bm's per-coefficient Satterthwaite DOF beats the knob: the + knob is fallback-level only, so hc2_bm inference is IDENTICAL with + the knob on and off.""" + data = self._clustered_panel(n_units=20) + kw = dict(outcome="y", treatment="group", time="post", unit="unit") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + r0 = DifferenceInDifferences(cluster="unit", vcov_type="hc2_bm").fit(data, **kw) + r1 = DifferenceInDifferences( + cluster="unit", vcov_type="hc2_bm", df_convention="cluster" + ).fit(data, **kw) + assert (r0.se, r0.p_value, r0.conf_int) == (r1.se, r1.p_value, r1.conf_int) + + def test_linalg_survey_df_precedence(self): + """LinearRegression level: a survey df on the fit wins over the knob.""" + from diff_diff.linalg import LinearRegression + + rng = np.random.default_rng(3) + X = np.column_stack([np.ones(120), rng.normal(size=(120, 2))]) + y = X @ np.array([1.0, 0.4, 2.0]) + rng.normal(size=120) + cl = np.repeat(np.arange(12), 10) + reg = LinearRegression( + cluster_ids=cl, include_intercept=False, df_convention="cluster" + ).fit(X, y) + assert reg.n_clusters_ == 12 + # simulate a survey fit: survey_df_ set -> must win over G-1 + reg.survey_df_ = 7 + inf = reg.get_inference(1) + assert inf.df == 7 + + def test_one_effective_cluster_fweight_fit_raises(self): + """Weighted CR1 fit where only ONE cluster carries positive weight + must fail at the vcov layer for ALL weight types: a zero-frequency + fweight cluster is as inert as a subpopulation-zeroed pweight + cluster, so the prior fweight carve-out let a degenerate + one-effective-cluster CR1 SE through (local-review round-3 P1).""" + from diff_diff.linalg import LinearRegression + + rng = np.random.default_rng(5) + n = 60 + X = np.column_stack([np.ones(n), rng.normal(size=(n, 2))]) + y = X @ np.array([1.0, 0.5, 2.0]) + rng.normal(size=n) + cl = np.repeat(np.arange(3), 20) + w = np.where(cl == 0, 2.0, 0.0) # only cluster 0 has positive weight + for wt in ("fweight", "pweight", "aweight"): + with pytest.raises(ValueError, match="at least 2 clusters"): + LinearRegression( + cluster_ids=cl, + include_intercept=False, + weights=w, + weight_type=wt, + df_convention="cluster", + ).fit(X, y) + + def test_get_inference_guard_fails_closed_at_one_cluster(self): + """Defense-in-depth: if a fit ever surfaces n_clusters_ <= 1 (all + supported routes now raise at the vcov layer), the knob's + get_inference branch returns all-NaN t/p/CI with a warning instead + of degrading to normal-theory inference via the non-positive-df + fallback (local-review round-2 P1).""" + from diff_diff.linalg import LinearRegression + + rng = np.random.default_rng(5) + n = 60 + X = np.column_stack([np.ones(n), rng.normal(size=(n, 2))]) + y = X @ np.array([1.0, 0.5, 2.0]) + rng.normal(size=n) + cl = np.repeat(np.arange(3), 20) + reg = LinearRegression( + cluster_ids=cl, include_intercept=False, df_convention="cluster" + ).fit(X, y) + reg.n_clusters_ = 1 # simulate the degenerate state directly + with pytest.warns(UserWarning, match="at least 2 effective"): + inf = reg.get_inference(1) + assert np.isnan(inf.t_stat) and np.isnan(inf.p_value) + assert np.isnan(inf.conf_int[0]) and np.isnan(inf.conf_int[1]) + assert np.isfinite(inf.se) # SE reported; only the reference- + # distribution-dependent fields are NaN + + def test_mpd_one_effective_cluster_fails_closed(self): + """MPD twin of the guard: a one-effective-cluster fit fails at the + CR1 vcov layer (the estimator surface has no raw weights= kwarg, so + the weighted one-positive-cluster route is exercised at the linalg + level above; here a single cluster label locks the MPD route).""" + data = self._clustered_panel(n_units=4) + data["treated"] = data["group"] + data["one_cluster"] = 0 + with pytest.raises(ValueError, match="at least 2 clusters"): + MultiPeriodDiD(cluster="one_cluster", df_convention="cluster").fit( + data, outcome="y", treatment="treated", time="time", unit="unit" + ) + + def test_nan_cluster_labels_rejected(self): + """Missing cluster labels are REJECTED before any clustered CR1 + computation (round-5 P0): the meat's pandas groupby drops + NaN-labelled rows, so any NaN-inclusive count (np.unique keeps a + NaN group) would disagree with the sandwich's partition — silently + wrong SEs. R fixest / Stata error on missing cluster values too. + Both the vcov layer and effective_cluster_count fail loud.""" + from diff_diff.linalg import LinearRegression, effective_cluster_count + + rng = np.random.default_rng(7) + n = 80 + cl = np.repeat(np.arange(4).astype(float), 20) + cl[cl == 3] = np.nan + with pytest.raises(ValueError, match="missing values"): + effective_cluster_count(cl) + with pytest.raises(ValueError, match="missing values"): + effective_cluster_count(cl, np.ones(n)) + X = np.column_stack([np.ones(n), rng.normal(size=(n, 2))]) + y = X @ np.array([1.0, 0.5, 2.0]) + rng.normal(size=n) + with pytest.raises(ValueError, match="missing values"): + LinearRegression(cluster_ids=cl, include_intercept=False, df_convention="cluster").fit( + X, y + ) + # default convention rejects too — the contract is the vcov layer's + with pytest.raises(ValueError, match="missing values"): + LinearRegression(cluster_ids=cl, include_intercept=False).fit(X, y) + + def test_nan_cluster_labels_rejected_all_backends(self): + """Round-6 P0: every clustered backend route shares the front-door + rejection — solve_ols(return_vcov=True) (the Rust HC1/CR1 dispatch + when available), CR2-BM (hc2_bm + cluster), and the MPD estimator + route (which calls solve_ols directly).""" + from diff_diff.linalg import compute_robust_vcov, solve_ols + + rng = np.random.default_rng(7) + n = 80 + cl = np.repeat(np.arange(4).astype(float), 20) + cl[cl == 3] = np.nan + X = np.column_stack([np.ones(n), rng.normal(size=(n, 2))]) + y = X @ np.array([1.0, 0.5, 2.0]) + rng.normal(size=n) + with pytest.raises(ValueError, match="missing values"): + solve_ols(X, y, cluster_ids=cl, return_vcov=True) + resid = y - X @ np.linalg.lstsq(X, y, rcond=None)[0] + with pytest.raises(ValueError, match="missing values"): + compute_robust_vcov(X, resid, cluster_ids=cl) # rust or numpy HC1 + with pytest.raises(ValueError, match="missing values"): + compute_robust_vcov(X, resid, cluster_ids=cl, vcov_type="hc2_bm") + data = self._clustered_panel(n_units=4) + data["treated"] = data["group"] + data.loc[data["unit"] == 0, "unit"] = np.nan + with pytest.raises(ValueError, match="missing values"): + MultiPeriodDiD(cluster="unit").fit( + data, outcome="y", treatment="treated", time="time", unit="unit" + ) + + def test_validation_and_transactional_set_params(self): + with pytest.raises(ValueError, match="df_convention"): + DifferenceInDifferences(df_convention="bogus") + est = DifferenceInDifferences() + with pytest.raises(ValueError, match="df_convention"): + est.set_params(df_convention="bogus") + assert est.df_convention == "residual" # unchanged after failed set + est.set_params(df_convention="cluster") + assert est.df_convention == "cluster" + + def test_get_params_roundtrip(self): + est = DifferenceInDifferences(cluster="unit", df_convention="cluster") + params = est.get_params() + assert params["df_convention"] == "cluster" + clone = DifferenceInDifferences(**params) + assert clone.df_convention == "cluster" + + def test_unclustered_fit_knob_is_inert(self): + """No cluster -> n_clusters_ is None -> knob has zero effect.""" + data = self._clustered_panel() + kw = dict(outcome="y", treatment="group", time="post") + r0 = DifferenceInDifferences().fit(data, **kw) + r1 = DifferenceInDifferences(df_convention="cluster").fit(data, **kw) + assert (r0.p_value, r0.conf_int) == (r1.p_value, r1.conf_int)