test(survey): true half-sample BRR replicate regressions per estimator family

CHANGELOG.md

@@ -80,6 +80,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   per-horizon event-study ATT and SE. Observed agreement on the reference platform:
   SE ~2e-10 (the covariate-augmented untreated `v_it` projection + clustering machinery),
   ATT ~2e-7; asserted at abs=1e-6/1e-7 for cross-platform robustness.
+- **True half-sample BRR replicate regressions per estimator family.** The
+  replicate-weight expansion tests used Fay-like 0.5/1.5 perturbations, under which every
+  unit keeps positive weight; true BRR was covered only at the vcov-helper level. A new
+  Hadamard-balanced half-sample generator (paired 2-PSU pseudo-strata; selected PSU w*2,
+  the other exactly 0 — the `survey::brrweights` full-BRR convention) now backs a
+  per-family regression class (DiD, DiD-absorb, MultiPeriodDiD, TWFE, SunAbraham,
+  StackedDiD, ImputationDiD, TwoStageDiD): finite positive replicate SEs under genuine
+  half-samples plus a base-weights point-estimate invariance check on EVERY family
+  (replicate columns drive only the variance). Construction sanity is itself asserted (exactly half the
+  paired PSUs zeroed per replicate, all multipliers in {0, 2}). Notably, genuine
+  half-samples CAN lose identification inside a replicate refit on some designs (the
+  staggered binary-interaction TWFE parameterization) — TWFE fails loudly there; the
+  family test uses its dedicated 2-period panel, with the behavior noted in-test.
 - **`CallawaySantAnna` ipw R-parity yardsticks folded into the golden fixture + no-covariate
   ipw structural-parity decision recorded.** `csdid_golden_values.json` regenerated (R 4.5.2,
   did 2.5.1, DRDID 1.3.0): all pre-existing data and result blocks reproduced byte-identically;

TODO.md

@@ -60,7 +60,6 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m
 |-------|----------|--------|--------|----------|
 | SE-audit CI-lock — remaining "fiddly bits" after the second coverage batch landed (that batch pinned C2 `dof_hc2_bm`/`dof_per_coef` via CI-inversion, C3 LOO `df`, C4 estimatr HC1/CR1 intercept SE, C5 Yatchew `p`/`sigma2_lin`/`sigma2_diff`, and the G2 fixest cluster-SE band). Still deferred, each needing a golden regeneration, new computation, or a documented-deviation call — **(a) G2 hetero/cluster is DONE (2026-07-07)**: the unbalanced/heteroskedastic-DGP regen landed, the DiD-path hetero AND cluster CR1 SEs are machine-precision-locked (the plain-OLS CR1 matches fixest exactly — the DOF-convention deviation is absorbed-FE-only), and the TWFE cluster band-pin is retained for the documented non-nested-FE ssc deviation (tracked under "Needs external reference"); TWFE has no public unclustered-hetero surface (auto-cluster convention). Remaining: **(b) PlaceboTests `boundary_gap`** — a permutation randomization-inference margin NOT computed anywhere in code (a new feature + result field, not a coverage lock); **(c) StackedDiD intercept SEs** (`se_cr1/cr2_intercept`, C1) — MEASURED to diverge ~0.3% from R: a nuisance-parameter reference-cell/parameterization gap, NOT machine-precision lockable (the event-study interaction SEs already match ~2e-13; surfacing it would add an unasserted, R-divergent public field); **(d) estimatr `classical` intercept SE** — same documented `O(1/n)` projection/DOF deviation as the slope (reference-only, excluded from parity). C6 CLOSED 2026-07-08 (dr remainder measured ~6e-11 ATT / ~2e-11 rel SE on the no-cov fixture and tightened to 1e-8 with the previously-missing SE assertion added; reg/ipw ATT bands 0.02/0.05 also tightened to 1e-8; the covariate-DR 2e-3 band stays as documented small-sample numerics). C7/C8 definitions were lost with the audit session scratchpad — re-derive from a fresh loose-tolerance inventory before acting. | `benchmarks/R/generate_fixest_did_twfe_golden.R`, `tests/test_fixest_did_twfe_parity.py`, `tests/test_methodology_stacked_did.py`, `tests/test_methodology_placebo.py` | SE-audit | Mid | Low |
 | Render `docs/methodology/REPORTING.md` and `REGISTRY.md` as in-site Sphinx pages so cross-refs can use `:doc:` instead of off-site `blob/main` URLs (stable-docs readers can otherwise land on a different revision than their package version). Two paths: (a) add `myst-parser` to `conf.py` + docs extras and link with `:doc:`, or (b) convert both to `.rst`. **Note:** REGISTRY.md is ~4.5k lines of LaTeX-heavy markdown — high risk under the `-W` (warnings-as-errors) Sphinx build; budget multiple rounds. | `docs/conf.py`, `docs/api/business_report.rst`, `docs/api/diagnostic_report.rst`, tutorials 18 & 19 | follow-up | Mid | Low |
-| Add true half-sample BRR replicate-weight regressions per estimator family (current tests use Fay-like 0.5/1.5 perturbations; `test_survey_phase6.py` covers true BRR at the helper level). | `tests/test_replicate_weight_expansion.py` | #253 | Mid | Low |
 | Port the CI `<notebook-prose>` extraction into the reviewer-eval harness so `docs/tutorials/*.ipynb` cases (currently guarded out of `verify-corpus`/`run`) can be reviewed with CI-equivalent context. | `tools/reviewer-eval/adapters/ci_prompt.py` | local-review | Mid | Low |
 
 ---

tests/test_replicate_weight_expansion.py

@@ -527,6 +527,66 @@ class TestReplicateVcovTypeWarn:
     point estimates only). Previously DiD silently ignored the kwarg and
     TWFE raised NotImplementedError — the warn-and-remap contract unifies
     the twins."""
+# ---------------------------------------------------------------------------
+# TRUE half-sample BRR (Hadamard-balanced) — per-estimator-family regressions
+# ---------------------------------------------------------------------------
+
+
+def _add_true_brr_replicates(data, unit_col="unit"):
+    """Add TRUE half-sample BRR replicate columns (Hadamard-balanced).
+
+    Pairs consecutive units into 2-PSU pseudo-strata and assigns
+    half-samples from a Sylvester-Hadamard matrix: within each pair the
+    selected PSU gets ``w*2`` and the other ``0`` (Wolter 2007 ch. 3; the
+    R ``survey::brrweights`` full-BRR convention with
+    ``combined_weights=True``). Unlike :func:`_add_brr_replicates`'
+    Fay-like 0.5/1.5 perturbation — under which every unit keeps positive
+    weight — every replicate here is a genuine half-sample: half the
+    paired PSUs carry exactly zero weight, exercising the zero-weight-PSU
+    code paths (FE identification drops, zero-mass cells) through each
+    estimator's replicate refit. An odd trailing unit (if any) is kept at
+    its base weight in every replicate (certainty PSU).
+    """
+    from scipy.linalg import hadamard
+
+    units = sorted(data[unit_col].unique())
+    n_paired = len(units) - (len(units) % 2)
+    pairs = [(units[i], units[i + 1]) for i in range(0, n_paired, 2)]
+    n_strata = len(pairs)
+    # Sylvester-Hadamard column 0 is all +1, which would leave the first
+    # pseudo-stratum permanently unbalanced (the same PSU selected in every
+    # replicate); survey::brrweights skips the constant column, so strata
+    # map to columns 1..n_strata — hence R >= n_strata + 1.
+    n_rep = 4
+    while n_rep < n_strata + 1:
+        n_rep *= 2
+    H = hadamard(n_rep)
+    base = data["weight"].to_numpy(dtype=float)
+    unit_vals = data[unit_col].to_numpy()
+    rep_cols = []
+    for r in range(n_rep):
+        w_r = base.copy()
+        for h, (u1, u2) in enumerate(pairs):
+            selected, dropped = (u1, u2) if H[r, h + 1] == 1 else (u2, u1)
+            w_r[unit_vals == selected] *= 2.0
+            w_r[unit_vals == dropped] = 0.0
+        col = f"tbrr_{r}"
+        data[col] = w_r
+        rep_cols.append(col)
+    return rep_cols
+
+
+class TestTrueBRRHalfSample:
+    """TRUE half-sample BRR per estimator family (TODO row: the smoke tests
+    above use Fay-like 0.5/1.5 perturbations, which never zero a unit;
+    ``test_survey_phase6.py`` covers true BRR only at the vcov-helper level).
+
+    Each test asserts (a) finite positive replicate SE under genuine
+    half-samples — half the paired PSUs at weight 0 per replicate — and
+    (b) the point estimate is IDENTICAL to the same fit
+    without replicate columns: replicate weights drive only the variance,
+    never the point estimate (base-weights invariance contract).
+    """
 
     @staticmethod
     def _sd(rep_cols):
@@ -646,3 +706,176 @@ def test_wild_bootstrap_rejection_fires_before_vcov_warning(self):
                     data, "outcome", "treated", "post", survey_design=self._sd(rep_cols)
                 )
         assert not any("has no effect with replicate-weight" in str(x.message) for x in w)
+    def test_construction_is_true_half_sample(self):
+        data = _make_simple_panel()
+        rep_cols = _add_true_brr_replicates(data)
+        units = sorted(data["unit"].unique())
+        n_paired = len(units) - (len(units) % 2)
+        base_per_unit = data.groupby("unit")["weight"].first()
+        for col in rep_cols:
+            per_unit = data.groupby("unit")[col].first()
+            zeroed = int((per_unit.loc[units[:n_paired]] == 0.0).sum())
+            assert zeroed == n_paired // 2, f"{col}: {zeroed} != {n_paired // 2}"
+            mult = per_unit.loc[units[:n_paired]] / base_per_unit.loc[units[:n_paired]]
+            assert set(np.round(mult.unique(), 12)) == {0.0, 2.0}
+        # Hadamard BALANCE: every paired PSU is selected (w*2) in exactly
+        # half the replicates — the property the constant Hadamard column
+        # would break for the first pseudo-stratum (each column 1..H of a
+        # Sylvester matrix has zero sum, so mean multiplier is exactly 1).
+        rep_mat = np.column_stack(
+            [data.groupby("unit")[c].first().loc[units[:n_paired]] for c in rep_cols]
+        )
+        base_vec = base_per_unit.loc[units[:n_paired]].to_numpy()
+        mean_mult = rep_mat.mean(axis=1) / base_vec
+        np.testing.assert_allclose(mean_mult, 1.0, rtol=0, atol=1e-12)
+
+    @staticmethod
+    def _assert_point_invariance(att_rep, att_base):
+        np.testing.assert_allclose(att_rep, att_base, rtol=0, atol=1e-12)
+
+    def test_did_true_brr(self):
+        data = _make_simple_panel()
+        rep_cols = _add_true_brr_replicates(data)
+        res = DifferenceInDifferences().fit(
+            data, "outcome", "treated", "post", survey_design=self._sd(rep_cols)
+        )
+        base = DifferenceInDifferences().fit(
+            data, "outcome", "treated", "post", survey_design=SurveyDesign(weights="weight")
+        )
+        assert np.isfinite(res.se) and res.se > 0
+        self._assert_point_invariance(res.att, base.att)
+
+    def test_did_absorb_true_brr(self):
+        data = _make_simple_panel()
+        data["group"] = (data["unit"] % 4).astype(str)
+        rep_cols = _add_true_brr_replicates(data)
+        res = DifferenceInDifferences().fit(
+            data,
+            "outcome",
+            "treated",
+            "post",
+            absorb=["group"],
+            survey_design=self._sd(rep_cols),
+        )
+        base = DifferenceInDifferences().fit(
+            data,
+            "outcome",
+            "treated",
+            "post",
+            absorb=["group"],
+            survey_design=SurveyDesign(weights="weight"),
+        )
+        assert np.isfinite(res.se) and res.se > 0
+        self._assert_point_invariance(res.att, base.att)
+
+    def test_multiperiod_true_brr(self):
+        data = _make_simple_panel()
+        rep_cols = _add_true_brr_replicates(data)
+        res = MultiPeriodDiD().fit(
+            data,
+            "outcome",
+            "treated",
+            "time",
+            post_periods=[1],
+            survey_design=self._sd(rep_cols),
+        )
+        base = MultiPeriodDiD().fit(
+            data,
+            "outcome",
+            "treated",
+            "time",
+            post_periods=[1],
+            survey_design=SurveyDesign(weights="weight"),
+        )
+        assert np.isfinite(res.avg_se) and res.avg_se > 0
+        self._assert_point_invariance(res.avg_att, base.avg_att)
+
+    def test_twfe_true_brr(self):
+        # The family's dedicated 2-period panel (the staggered fixture's
+        # binary treated x post parameterization can lose identification
+        # inside a genuine half-sample replicate refit — a real behavior
+        # of true BRR that the Fay-like perturbation never exercised;
+        # TWFE fails loudly there rather than silently degrading).
+        data = TestTWFEReplicate._make_twfe_panel()
+        rep_cols = _add_true_brr_replicates(data)
+        res = TwoWayFixedEffects().fit(
+            data, "outcome", "treated", "post", "unit", survey_design=self._sd(rep_cols)
+        )
+        base = TwoWayFixedEffects().fit(
+            data,
+            "outcome",
+            "treated",
+            "post",
+            "unit",
+            survey_design=SurveyDesign(weights="weight"),
+        )
+        assert np.isfinite(res.se) and res.se > 0
+        self._assert_point_invariance(res.att, base.att)
+
+    def test_sun_abraham_true_brr(self):
+        data = _make_staggered_panel()
+        rep_cols = _add_true_brr_replicates(data)
+        res = SunAbraham(n_bootstrap=0).fit(
+            data, "outcome", "unit", "time", "first_treat", survey_design=self._sd(rep_cols)
+        )
+        base = SunAbraham(n_bootstrap=0).fit(
+            data,
+            "outcome",
+            "unit",
+            "time",
+            "first_treat",
+            survey_design=SurveyDesign(weights="weight"),
+        )
+        assert np.isfinite(res.overall_se) and res.overall_se > 0
+        self._assert_point_invariance(res.overall_att, base.overall_att)
+
+    def test_stacked_true_brr(self):
+        data = _make_staggered_panel()
+        rep_cols = _add_true_brr_replicates(data)
+        res = StackedDiD().fit(
+            data, "outcome", "unit", "time", "first_treat", survey_design=self._sd(rep_cols)
+        )
+        base = StackedDiD().fit(
+            data,
+            "outcome",
+            "unit",
+            "time",
+            "first_treat",
+            survey_design=SurveyDesign(weights="weight"),
+        )
+        assert np.isfinite(res.overall_se) and res.overall_se > 0
+        self._assert_point_invariance(res.overall_att, base.overall_att)
+
+    def test_imputation_true_brr(self):
+        data = _make_staggered_panel()
+        rep_cols = _add_true_brr_replicates(data)
+        res = ImputationDiD(n_bootstrap=0).fit(
+            data, "outcome", "unit", "time", "first_treat", survey_design=self._sd(rep_cols)
+        )
+        base = ImputationDiD(n_bootstrap=0).fit(
+            data,
+            "outcome",
+            "unit",
+            "time",
+            "first_treat",
+            survey_design=SurveyDesign(weights="weight"),
+        )
+        assert np.isfinite(res.overall_se) and res.overall_se > 0
+        self._assert_point_invariance(res.overall_att, base.overall_att)
+
+    def test_two_stage_true_brr(self):
+        data = _make_staggered_panel()
+        rep_cols = _add_true_brr_replicates(data)
+        res = TwoStageDiD(n_bootstrap=0).fit(
+            data, "outcome", "unit", "time", "first_treat", survey_design=self._sd(rep_cols)
+        )
+        base = TwoStageDiD(n_bootstrap=0).fit(
+            data,
+            "outcome",
+            "unit",
+            "time",
+            "first_treat",
+            survey_design=SurveyDesign(weights="weight"),
+        )
+        assert np.isfinite(res.overall_se) and res.overall_se > 0
+        self._assert_point_invariance(res.overall_att, base.overall_att)