fix(survey): unify explicit-vcov_type contract under replicate designs (warn + hc1 remap)

CHANGELOG.md

@@ -309,6 +309,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   supersedes it.
 
 ### Changed
+- **Explicit `vcov_type` with replicate-weight survey designs now warns and proceeds
+  (was: inconsistent raise/silent-ignore).** With replicate variance the analytical
+  sandwich is replaced wholesale — per-replicate refits return point estimates only,
+  identical across vcov families (FWL), so the requested vcov family cannot influence any
+  reported number. `TwoWayFixedEffects(vcov_type="hc2"/"hc2_bm")` previously raised
+  `NotImplementedError` while `DifferenceInDifferences` silently ignored the kwarg; both
+  (and `MultiPeriodDiD`) now share one contract: explicit non-`hc1` analytical
+  `vcov_type` (`hc2`/`hc2_bm`/`classical`) emits a `UserWarning` and the discarded base fit remaps to `hc1` (avoids wasted CR2-BM work,
+  one-way validator rejections, and the TWFE full-dummy route, which does not compose
+  with per-replicate re-demeaning). Explicit `hc1` — the old error message's own
+  workaround guidance — stays silent; `conley` keeps its own survey-design support contract and validators
+  (excluded from the remap). A per-replicate full-dummy HC2 implementation
+  (the TODO row) was investigated and rejected as a costly no-op: it cannot change the
+  replicate variance. Tests lock warn+bit-identity-to-hc1 on all three estimators.
 - **Per-cell solver fast paths for covariate fits (CallawaySantAnna and every
   estimator routing through the shared solvers).** Two pure-Python changes in
   `diff_diff/linalg.py`: (1) `solve_logit`'s IRLS inner step — previously a full

TODO.md

@@ -30,7 +30,6 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo
 |-------|----------|--------|--------|----------|
 | `SyntheticControl` conformal (CWZ 2021) extensions: (a) one-sided / signed-`t` variants (§7); (b) covariates in the conformal proxy (`X_jt`, eqs 4/6 — current proxy is outcomes-only); (c) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies. The joint test, pointwise CIs, and average-effect CI have landed. | `conformal.py`, `synthetic_control_results.py` | CWZ-2021 | Heavy | Low |
 | `ContinuousDiD` CGBS-2024 extensions. (a) `covariates=` kwarg — **DONE (reg/dr)**; (b) discrete-treatment saturated regression (`treatment_type="discrete"`) — **DONE**; (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0` (`control_group="lowest_dose"`) — **DONE** (discrete + continuous mass-point, single-cohort; estimand `ATT(d)−ATT(d_L)`; see REGISTRY Note #7). Remaining (all deferred `NotImplementedError`, documented): `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate); `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF); multi-cohort **heterogeneous-support** discrete aggregation (support-aware: average each dose only over the cohorts that observe it); **multi-cohort `lowest_dose`** (within-cohort `d_L` reference + support-aware cross-cohort aggregation); and **`covariates=` × `lowest_dose`** (conditional-PT-relative-to-`d_L` estimand). Single-cohort / 2-period / shared-support multi-cohort are supported. | `continuous_did.py` | CGBS-2024 | Heavy | Low |
-| `TwoWayFixedEffects(vcov_type in {hc2, hc2_bm})` with replicate-weight designs raises `NotImplementedError` (`twfe.py:~233`). The replicate path re-demeans per replicate, which doesn't compose with the full-dummy HC2/HC2-BM build — a correct impl needs per-replicate full-dummy refit. Workaround: `hc1` for replicate-weight CR1. | `twfe.py::fit` | follow-up | Heavy | Low |
 | TWFE's HC2/HC2-BM inline full-dummy build (`twfe.py:280-315`) duplicates the dummy-construction logic in `DifferenceInDifferences(fixed_effects=...)` (`estimators.py:478-486`). Extract a shared helper, or delegate TWFE's HC2/HC2-BM path to DiD's `fixed_effects=` branch (with TWFE-specific cluster-default threading), to reduce drift risk on FE naming / survey behavior / result-surface conventions. Substantive refactor — touches both estimators. | `twfe.py::fit`, `estimators.py::DifferenceInDifferences.fit` | follow-up | Heavy | Low |
 | `HonestDiD` Δ^SD optimal-FLCI center parity — **LANDED (SE-audit B2b)**: replaced the flat Nelder-Mead affine-estimator optimizer with a faithful port of R `HonestDiD::findOptimalFLCI`'s nested convex program (inner min-worst-case-bias at fixed estimator SD `h` via SLSQP QCQP; outer grid-zoom over `h`). Matches R's center + half-length + optimalVec to ~1e-3 (median ~1e-5) across a stress grid; the prior ~9% intermediate-M center drift is removed (widths always matched, coverage unaffected). Analytical folded-normal cv is more accurate than R's MC `.qfoldednormal`. Golden `honest_flci_golden.json` + `TestHonestFLCIParityR`. | `honest_did.py::_flci_solve` | SE-audit | Done | Low |
 

diff_diff/estimators.py

@@ -360,6 +360,7 @@ def fit(
                 "survey designs. Replicate weights provide their own variance "
                 "estimation."
             )
+        _replicate_vcov_remap = _uses_replicate and self._warn_replicate_vcov_ignored()
 
         # Handle absorbed fixed effects (within-transformation)
         working_data = data.copy()
@@ -411,7 +412,12 @@ def fit(
         # and solves WLS directly, with no within-transform step. R2 review
         # surfaced the scope mismatch (REGISTRY/CHANGELOG said "SUPPORTED" but
         # the survey guard fired first on weighted multi-absorb fits).
-        if absorb and self.vcov_type in ("hc2", "hc2_bm"):
+        # Route on the EFFECTIVE vcov family: under a replicate design the
+        # remap to hc1 must also disable this full-dummy swap, or an
+        # explicit hc2 request would still change the result surface
+        # (full-dummy coefficients vs absorbed reduced fit) despite the
+        # "has no effect" warning.
+        if absorb and not _replicate_vcov_remap and self.vcov_type in ("hc2", "hc2_bm"):
             fixed_effects = list(fixed_effects or []) + list(absorb)
             absorb = None
             absorbed_vars = []
@@ -576,7 +582,11 @@ def fit(
 
         # Remap implicit "classical" + cluster to CR1 for legacy-alias
         # backward compatibility (see `_resolve_effective_vcov_type`).
-        _fit_vcov_type = self._resolve_effective_vcov_type(effective_cluster_ids)
+        _fit_vcov_type = (
+            "hc1"
+            if _replicate_vcov_remap
+            else self._resolve_effective_vcov_type(effective_cluster_ids)
+        )
 
         # Build Conley coord/time/unit arrays when applicable. CRITICAL:
         # read from the ORIGINAL `data` frame, NOT `working_data` — `absorb`
@@ -1134,6 +1144,45 @@ def set_params(self, **params) -> "DifferenceInDifferences":
         self._vcov_type_explicit = self._vcov_type_arg is not None
         return self
 
+    def _warn_replicate_vcov_ignored(self) -> bool:
+        """Warn that an explicit ``vcov_type`` has no effect under a
+        replicate-weight survey design, and tell the caller to remap the
+        fit-time vcov to ``"hc1"``.
+
+        With ``uses_replicate_variance`` the analytical sandwich is replaced
+        wholesale by the replicate-refit variance (the per-replicate refits
+        return point estimates only), so the requested vcov family cannot
+        influence any reported number. Silently honoring the kwarg would
+        report hc1-identical output under an ``hc2``/``hc2_bm``/``classical``
+        label; remapping the (discarded) base-fit vcov to ``"hc1"`` also
+        avoids wasted CR2-BM work and one-way-only validator rejections.
+        ``conley`` is excluded — it carries its own survey-design support
+        contract and validators, which must keep firing unchanged. Returns
+        True when a remap should be applied (explicit non-hc1, non-conley
+        vcov_type).
+        """
+        if not self._vcov_type_explicit or self.vcov_type == "hc1":
+            # Explicit hc1 is exactly the remap target (and the value the
+            # old TwoWayFixedEffects NotImplementedError guidance told
+            # users to pass) — nothing is being overridden, stay quiet.
+            return False
+        if self.vcov_type == "conley":
+            # Conley keeps its own survey-design support contract (TSL via
+            # the stratified-Conley sandwich; dedicated per-design
+            # validators in diff_diff.conley) — do not warn-and-remap past
+            # those gates; let the conley validation speak for itself.
+            return False
+        warnings.warn(
+            f"vcov_type={self.vcov_type!r} has no effect with replicate-weight "
+            "survey designs: the replicate-refit variance replaces the "
+            "analytical vcov entirely (per-replicate refits return point "
+            "estimates only, identical across vcov families). Proceeding "
+            "with replicate variance; the base fit uses 'hc1'.",
+            UserWarning,
+            stacklevel=3,
+        )
+        return True
+
     def _resolve_effective_vcov_type(self, effective_cluster_ids) -> str:
         """Pick the ``vcov_type`` to use for a given fit given cluster context.
 
@@ -1555,6 +1604,7 @@ def fit(  # type: ignore[override]
                 "survey designs. Replicate weights provide their own variance "
                 "estimation."
             )
+        _replicate_vcov_remap_mp = _uses_replicate_mp and self._warn_replicate_vcov_ignored()
 
         # Handle absorbed fixed effects (within-transformation)
         working_data = data.copy()
@@ -1605,7 +1655,8 @@ def fit(  # type: ignore[override]
         # dimensions") doesn't apply when we're about to swap absorb for
         # fixed_effects: the fixed_effects= path builds the full-dummy
         # design and solves WLS directly, with no within-transform step.
-        if absorb and self.vcov_type in ("hc2", "hc2_bm"):
+        # Route on the EFFECTIVE vcov family (see DifferenceInDifferences).
+        if absorb and not _replicate_vcov_remap_mp and self.vcov_type in ("hc2", "hc2_bm"):
             fixed_effects = list(fixed_effects or []) + list(absorb)
             absorb = None
             n_absorbed_effects = 0
@@ -1803,7 +1854,11 @@ def fit(  # type: ignore[override]
         _use_survey_vcov = resolved_survey is not None and resolved_survey.needs_survey_vcov
 
         # Remap implicit "classical" + cluster to CR1 (legacy backward compat).
-        _fit_vcov_type = self._resolve_effective_vcov_type(effective_cluster_ids)
+        _fit_vcov_type = (
+            "hc1"
+            if _replicate_vcov_remap_mp
+            else self._resolve_effective_vcov_type(effective_cluster_ids)
+        )
 
         # Cluster + CR2 Bell-McCaffrey (non-survey, unweighted) shares the SAME
         # expensive CR2 precomputes (per-cluster A_g eigendecompositions, S_W,

diff_diff/twfe.py

@@ -242,23 +242,20 @@ def fit(  # type: ignore[override]
                 "survey designs. Replicate weights provide their own variance "
                 "estimation."
             )
-        # Replicate weights + HC2 / HC2-BM is incompatible with the
-        # full-dummy auto-route: the replicate path re-demeans per
-        # replicate (re-demeaning depends on the per-replicate weight
-        # vector), which doesn't compose with the full-dummy design
-        # build. A correct implementation would need to re-build the
-        # full-dummy X per replicate and recompute the HC2 leverage,
-        # which is deferred. Mirrors the
-        # ``linalg.py::_validate_vcov_args`` ``hc2_bm + weights`` gate.
-        if _uses_replicate_twfe and self.vcov_type in ("hc2", "hc2_bm"):
-            raise NotImplementedError(
-                f"TwoWayFixedEffects(vcov_type={self.vcov_type!r}) with "
-                "replicate-weight survey designs is not yet supported: the "
-                "replicate path re-demeans per replicate, which does not "
-                "compose with the full-dummy HC2/HC2-BM build (would need "
-                "per-replicate full-dummy refit). Use vcov_type='hc1' for "
-                "replicate-weight CR1, or drop to analytical inference."
-            )
+        # Replicate designs replace the analytical vcov wholesale: the
+        # per-replicate refits return point estimates only, so vcov_type
+        # cannot influence any reported number (FWL: the full-dummy and
+        # within-transformed fits give identical coefficients, hence an
+        # identical replicate variance). An explicit vcov_type therefore
+        # warns and the (discarded) base fit remaps to hc1 on the
+        # within-transformed design — this also disables the full-dummy
+        # HC2/HC2-BM auto-route, which does not compose with the
+        # per-replicate re-demeaning. Previously hc2/hc2_bm raised
+        # NotImplementedError here; the warn-and-proceed contract matches
+        # DifferenceInDifferences/MultiPeriodDiD (shared helper).
+        _replicate_vcov_remap_twfe = _uses_replicate_twfe and self._warn_replicate_vcov_ignored()
+        if _replicate_vcov_remap_twfe:
+            use_full_dummy = False
 
         # Unit-level clustering is the TWFE default when `cluster` is not
         # explicitly provided. But the one-way ``classical`` and ``hc2``
@@ -472,7 +469,11 @@ def fit(  # type: ignore[override]
         # Don't forward `robust=self.robust` to LinearRegression when the
         # remapped vcov_type disagrees; the remapped `vcov_type` is the
         # single source of truth.
-        _fit_vcov_type = self._resolve_effective_vcov_type(survey_cluster_ids)
+        _fit_vcov_type = (
+            "hc1"
+            if _replicate_vcov_remap_twfe
+            else self._resolve_effective_vcov_type(survey_cluster_ids)
+        )
 
         # Phase 2 panel-Conley: build coord array + row-aligned time/unit
         # vectors from the original (un-demeaned) data. FWL composability:

docs/methodology/REGISTRY.md

@@ -4790,6 +4790,7 @@ variance from the distribution of replicate estimates.
   - JK1: `V = (R-1)/R * sum_r (theta_r - theta)^2`
   - JKn: `V = sum_h ((n_h-1)/n_h) * sum_{r in h} (theta_r - theta)^2`
 - **Note:** SDR (Successive Difference Replication) uses variance factor 4/R, following Fay & Train (1995). Used by ACS PUMS (80 replicate columns). Treated identically to BRR for scaling purposes — no fay_rho, no replicate_strata, custom scale/rscales ignored.
+- **Note (vcov_type has no effect under replicate variance, 2026-07):** with `uses_replicate_variance`, the analytical sandwich is replaced wholesale by the replicate-refit variance — the per-replicate refits return point estimates only, which are identical across vcov families (FWL: full-dummy and within-transformed fits give the same coefficients). An explicit non-`hc1` analytical `vcov_type` (`hc2`, `hc2_bm`, `classical`) on `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` therefore emits a `UserWarning` and the (discarded) base fit remaps to `hc1` — avoiding wasted CR2-BM work, one-way-only validator rejections, and the TWFE full-dummy auto-route (which does not compose with per-replicate re-demeaning). Explicit `hc1` stays silent (it is the remap target). `conley` is excluded from the remap: it carries its own survey-design support contract (TSL stratified-Conley sandwich; dedicated per-design validators), which keeps firing unchanged. This replaces two prior inconsistent behaviors: `TwoWayFixedEffects(vcov_type="hc2"/"hc2_bm")` raised `NotImplementedError`, while `DifferenceInDifferences` silently ignored the kwarg — a per-replicate full-dummy HC2 implementation was considered and rejected as a costly no-op (it cannot change the replicate variance).
 - **IF-based replicate variance**: For influence-function estimators (CS
   aggregation, ContinuousDiD, EfficientDiD, TripleDifference), replicate
   contrasts are formed via weight-ratio rescaling:

tests/test_estimators_vcov_type.py

@@ -971,13 +971,17 @@ def test_twfe_hc2_wild_bootstrap_keeps_auto_cluster(self):
         assert res.n_clusters == 20
 
     @pytest.mark.parametrize("vcov", ["hc2", "hc2_bm"])
-    def test_twfe_rejects_replicate_weights_under_hc2(self, vcov):
-        """TWFE + hc2/hc2_bm + replicate-weight survey design raises
-        NotImplementedError.
-
-        The replicate path re-demeans per replicate (re-demeaning depends
-        on the per-replicate weight vector), which doesn't compose with
-        the full-dummy build. Documented scope limit; tracked in TODO.md.
+    def test_twfe_replicate_weights_under_hc2_warns_and_matches_hc1(self, vcov):
+        """TWFE + hc2/hc2_bm + replicate-weight survey design warns and
+        proceeds with replicate variance, bit-identical to the hc1 request.
+
+        Replicate designs replace the analytical vcov wholesale (the
+        per-replicate refits return point estimates only, identical across
+        vcov families by FWL), so vcov_type cannot influence any reported
+        number. The previous NotImplementedError is replaced by the
+        warn-and-remap contract shared with DifferenceInDifferences /
+        MultiPeriodDiD; explicit hc1 (the old error's own guidance) stays
+        silent.
         """
         data = _make_did_panel(n_units=20).copy()
         # Attach full-sample weight + 4 BRR replicate-weight columns.
@@ -992,11 +996,9 @@ def test_twfe_rejects_replicate_weights_under_hc2(self, vcov):
             replicate_method="BRR",
             weight_type="pweight",
         )
-        with pytest.raises(
-            NotImplementedError,
-            match=r"replicate-weight.*not yet supported",
-        ):
-            TwoWayFixedEffects(vcov_type=vcov).fit(
+
+        def _fit(vc):
+            return TwoWayFixedEffects(vcov_type=vc).fit(
                 data,
                 outcome="y",
                 treatment="treated",
@@ -1005,6 +1007,15 @@ def test_twfe_rejects_replicate_weights_under_hc2(self, vcov):
                 survey_design=sd,
             )
 
+        with pytest.warns(UserWarning, match="has no effect with replicate-weight"):
+            res = _fit(vcov)
+        with warnings.catch_warnings():
+            warnings.simplefilter("error", UserWarning)
+            base = _fit("hc1")  # explicit hc1 must NOT warn
+        assert res.att == base.att
+        assert res.se == base.se
+        assert np.isfinite(res.se) and res.se > 0
+
     def test_twfe_hc2_always_treated_unit_finite_att(self):
         """Always-treated unit (D=1 in all periods) doesn't poison the ATT
         on the full-dummy HC2 path.