feat(inference): opt-in df_convention="cluster" knob for clustered t/p/CI

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

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_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_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
 )
 ```