From b6a64492a638f9f0d734b4d17dbe7131e51c70e5 Mon Sep 17 00:00:00 2001 From: igerber Date: Tue, 7 Jul 2026 18:34:09 -0400 Subject: [PATCH 1/3] test(survey): true half-sample BRR replicate regressions per estimator family MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 (test_survey_phase6). Adds a Hadamard-balanced half-sample generator (paired 2-PSU pseudo-strata; selected PSU w*2, the other exactly 0 — the survey::brrweights full-BRR convention) and a per-family regression class (DiD, DiD-absorb, MultiPeriodDiD, TWFE, SunAbraham, StackedDiD, ImputationDiD, TwoStageDiD) asserting finite positive replicate SEs under genuine half-samples plus base-weights point-estimate invariance (replicate columns drive only the variance). Construction sanity asserted (half the paired PSUs zeroed per replicate; multipliers in {0,2}). Genuine half-samples CAN lose identification inside a replicate refit on some designs (staggered binary-interaction TWFE) — TWFE fails loudly there; noted in-test with the family's dedicated panel. Co-Authored-By: Claude Fable 5 --- CHANGELOG.md | 13 ++ TODO.md | 1 - tests/test_replicate_weight_expansion.py | 201 +++++++++++++++++++++++ 3 files changed, 214 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7a5a82b4..5ed6a002 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -119,6 +119,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 by seven new adapter tests (wrapper + sanitization, zero-extracted, aggregate truncation, byte-vs-char cap parity, trusted-extractor sentinel, Git-quoted filename discovery, end-to-end prompt assembly order). +- **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 (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; diff --git a/TODO.md b/TODO.md index c9fe7934..850e2268 100644 --- a/TODO.md +++ b/TODO.md @@ -57,7 +57,6 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| | 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 | -| 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 | --- diff --git a/tests/test_replicate_weight_expansion.py b/tests/test_replicate_weight_expansion.py index 1a4fe229..20af63a4 100644 --- a/tests/test_replicate_weight_expansion.py +++ b/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) where checked, 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,144 @@ 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), + ) + assert np.isfinite(res.att) + assert np.isfinite(res.se) and res.se > 0 + + 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), + ) + assert np.isfinite(res.avg_att) + assert np.isfinite(res.avg_se) and res.avg_se > 0 + + 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) + ) + assert np.isfinite(res.overall_att) + assert np.isfinite(res.overall_se) and res.overall_se > 0 + + 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) + ) + assert np.isfinite(res.overall_att) + assert np.isfinite(res.overall_se) and res.overall_se > 0 From c6811cf2ce7a5b0e9d78c1e5bf423a6487e881de Mon Sep 17 00:00:00 2001 From: igerber Date: Wed, 8 Jul 2026 07:15:09 -0400 Subject: [PATCH 2/3] test(survey): base-weights invariance asserted on every true-BRR family (review P3) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The reviewer flagged CHANGELOG/coverage drift: invariance was asserted for DiD/TWFE/SunAbraham/ImputationDiD only. Rather than qualifying the prose, the assertion is added to the remaining families (DiD-absorb, MultiPeriodDiD, StackedDiD, TwoStageDiD) — all pass, so the invariance contract genuinely holds family-wide; docstring hedge removed. Co-Authored-By: Claude Fable 5 --- CHANGELOG.md | 4 +-- tests/test_replicate_weight_expansion.py | 42 +++++++++++++++++++++--- 2 files changed, 39 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5ed6a002..d2608cc9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -126,8 +126,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 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 (replicate columns - drive only the variance). Construction sanity is itself asserted (exactly half the + 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 diff --git a/tests/test_replicate_weight_expansion.py b/tests/test_replicate_weight_expansion.py index 20af63a4..441f0e47 100644 --- a/tests/test_replicate_weight_expansion.py +++ b/tests/test_replicate_weight_expansion.py @@ -583,7 +583,7 @@ class TestTrueBRRHalfSample: Each test asserts (a) finite positive replicate SE under genuine half-samples — half the paired PSUs at weight 0 per replicate — and - (b) where checked, the point estimate is IDENTICAL to the same fit + (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). """ @@ -757,8 +757,16 @@ def test_did_absorb_true_brr(self): absorb=["group"], survey_design=self._sd(rep_cols), ) - assert np.isfinite(res.att) + 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() @@ -771,8 +779,16 @@ def test_multiperiod_true_brr(self): post_periods=[1], survey_design=self._sd(rep_cols), ) - assert np.isfinite(res.avg_att) + 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 @@ -819,8 +835,16 @@ def test_stacked_true_brr(self): res = StackedDiD().fit( data, "outcome", "unit", "time", "first_treat", survey_design=self._sd(rep_cols) ) - assert np.isfinite(res.overall_att) + 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() @@ -845,5 +869,13 @@ def test_two_stage_true_brr(self): res = TwoStageDiD(n_bootstrap=0).fit( data, "outcome", "unit", "time", "first_treat", survey_design=self._sd(rep_cols) ) - assert np.isfinite(res.overall_att) + 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) From 9131400726ba5db1fb1cc40c823c3e783a344d19 Mon Sep 17 00:00:00 2001 From: igerber Date: Thu, 9 Jul 2026 06:24:16 -0400 Subject: [PATCH 3/3] test: repair class structure mangled by the union merge (review P2) The append-append conflict resolution had spliced the true-BRR banner + helper INSIDE TestReplicateVcovTypeWarn (leaving it an empty class) and absorbed its six methods into TestTrueBRRHalfSample, making class-level targeting silently skip them and the BRR docstring inaccurate. Blocks reordered: each class owns its methods (+ its own _sd helper); both class selectors collect (6 + 9), file total 37 passes. Co-Authored-By: Claude Fable 5 --- tests/test_replicate_weight_expansion.py | 127 ++++++++++++----------- 1 file changed, 67 insertions(+), 60 deletions(-) diff --git a/tests/test_replicate_weight_expansion.py b/tests/test_replicate_weight_expansion.py index 441f0e47..ca4f9a4f 100644 --- a/tests/test_replicate_weight_expansion.py +++ b/tests/test_replicate_weight_expansion.py @@ -527,66 +527,6 @@ 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): @@ -706,6 +646,73 @@ 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) + + +# --------------------------------------------------------------------------- +# 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): + return SurveyDesign(weights="weight", replicate_weights=rep_cols, replicate_method="BRR") + def test_construction_is_true_half_sample(self): data = _make_simple_panel() rep_cols = _add_true_brr_replicates(data)