From 05aebbe03f314eeb12f36cbc292130c704ca7634 Mon Sep 17 00:00:00 2001 From: igerber Date: Fri, 26 Jun 2026 16:30:10 -0400 Subject: [PATCH] feat(diagnostics): PlaceboTests methodology validation + exact RI p-value PR-B of the PlaceboTests methodology promotion (Bertrand-Duflo-Mullainathan 2004). Source validation + methodology tests + R parity + tracker -> Complete. Source fixes (diff_diff/diagnostics.py): - permutation_test: exact randomization-inference p-value (1 + count)/(B + 1) (Phipson & Smyth 2010), replacing count/B floored at 1/(B+1); the +1 includes the observed statistic, making the docstring's "exact" claim faithful (now qualified to "exact under the sharp null", BDM fn 12). - placebo_group_test: optional `treatment` parameter that drops ever-treated units for an uncontaminated placebo; degenerate-design ValueError guards (replacing a cryptic LinAlgError) + misuse UserWarning; docstring corrected for both modes. The run_placebo_test dispatcher's fake_group path now filters ever-treated units by default. Validation: - tests/test_methodology_placebo.py (19 tests, 5 BDM-anchored classes). - R parity via base-R combn exact enumeration: benchmarks/R/generate_placebo_golden.R -> benchmarks/data/placebo_golden.json (+ placebo_test_panel.csv). - REGISTRY.md ## PlaceboTests full entry; METHODOLOGY_REVIEW.md row -> Complete. - references.rst: Phipson & Smyth (2010), Rosenbaum (1996); doc-deps mapping; CHANGELOG. Co-Authored-By: Claude Opus 4.8 (1M context) --- CHANGELOG.md | 26 ++ METHODOLOGY_REVIEW.md | 63 +-- benchmarks/R/generate_placebo_golden.R | 133 +++++++ benchmarks/data/placebo_golden.json | 35 ++ benchmarks/data/placebo_test_panel.csv | 17 + diff_diff/diagnostics.py | 107 +++++- docs/doc-deps.yaml | 5 +- docs/methodology/REGISTRY.md | 49 ++- docs/references.rst | 8 + tests/test_methodology_placebo.py | 507 +++++++++++++++++++++++++ 10 files changed, 905 insertions(+), 45 deletions(-) create mode 100644 benchmarks/R/generate_placebo_golden.R create mode 100644 benchmarks/data/placebo_golden.json create mode 100644 benchmarks/data/placebo_test_panel.csv create mode 100644 tests/test_methodology_placebo.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 61140c38..1c8aa403 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,7 +7,23 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- **`placebo_group_test` gained an optional `treatment` parameter.** When supplied, units + that are ever real-treated are dropped before the placebo so it runs on never-treated units + only (the uncontaminated design); without it, behavior is unchanged and the caller must pass + control-only data. Degenerate designs (all fake-treated units dropped, or no controls + remaining) now raise a clear `ValueError` instead of a cryptic `LinAlgError`, and a + fake-treated unit that is itself real-treated emits a `UserWarning`. +- **PlaceboTests methodology validation:** `tests/test_methodology_placebo.py` (paper-anchored + to Bertrand-Duflo-Mullainathan 2004) plus base-R exact-enumeration R parity + (`benchmarks/R/generate_placebo_golden.R` → `benchmarks/data/placebo_golden.json`). The + `PlaceboTests` methodology-review row is promoted to **Complete**. + ### Changed +- **`run_placebo_test`'s `fake_group` path now filters ever-treated units by default.** The + dispatcher threads its `treatment` column into `placebo_group_test`, so the fake-group + placebo runs on never-treated units only (a more-correct placebo). Calling + `placebo_group_test` directly without `treatment` retains the previous behavior. - **Bumped the Rust backend's `blas-src` crate `0.10` → `0.14`.** `blas-src` is a linker-only crate pulled in **only by the `accelerate` (macOS) feature**; the Linux `openblas` path links system OpenBLAS via `build.rs` and the default/Windows builds use the @@ -17,6 +33,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 `maturin develop --features accelerate` against the pinned `ndarray 0.17`, the Rust unit tests, and the full Python⇄Rust equivalence suite (`tests/test_rust_backend.py`). +### Fixed +- **`permutation_test` now reports the randomization-inference p-value + `(1 + count) / (B + 1)`** (Phipson & Smyth 2010), replacing `count / B` floored at + `1/(B+1)`. The `+1` includes the observed statistic in both numerator and denominator + (the floor is now intrinsic). Because assignments are sampled with replacement, this is a + valid but slightly conservative Monte-Carlo randomization-inference p-value (not an exact + finite-sample value); it converges to the exact full-enumeration value `count/total` as the + number of permutations grows. (Permutation p-values shift by a small amount, at most + `~1/(B+1)`.) + ### Security - **Bumped the Rust backend's `pyo3` and `numpy` crates 0.28 → 0.29.** Resolves two RustSec advisories in `pyo3 < 0.29` — RUSTSEC-2026-0176 (out-of-bounds read in `PyList`/`PyTuple` diff --git a/METHODOLOGY_REVIEW.md b/METHODOLOGY_REVIEW.md index 82bbb9e9..bd52583b 100644 --- a/METHODOLOGY_REVIEW.md +++ b/METHODOLOGY_REVIEW.md @@ -24,7 +24,7 @@ A **Complete** entry has a documented review pass against the primary academic s The catalog grew incrementally over several quarters, so formats vary across the existing Complete entries; the consistent invariant is that someone walked through the implementation against the academic source and captured the result here. New reviews going forward should aim for the fuller structure (Verified Components + Corrections Made + Deviations + dedicated methodology test file) used by the more recent entries. -**In Progress** entries have a REGISTRY.md section and unit-test coverage, but no formal walk-through has been captured here yet. The In Progress band is wide — some entries also have some combination of a paper review (primary or companion), a dedicated methodology test file, and R parity fixtures; others have only the REGISTRY entry and unit tests (e.g., PlaceboTests). The "Documentation in place" sub-section enumerates what each entry already has; the "Outstanding for promotion" sub-section enumerates what's still needed to flip it to Complete. +**In Progress** entries have a REGISTRY.md section and unit-test coverage, but no formal walk-through has been captured here yet. The In Progress band is wide — some entries also have some combination of a paper review (primary or companion), a dedicated methodology test file, and R parity fixtures; others have only the REGISTRY entry and unit tests. The "Documentation in place" sub-section enumerates what each entry already has; the "Outstanding for promotion" sub-section enumerates what's still needed to flip it to Complete. **Not Started** entries have neither a tracker walk-through nor an REGISTRY.md section. This tracker no longer carries any Not Started rows; new estimators are expected to enter as In Progress when their REGISTRY entry lands. @@ -82,7 +82,7 @@ The catalog grew incrementally over several quarters, so formats vary across the | HonestDiD | `honest_did.py` | `HonestDiD` package | **Complete** | 2026-04-01 | | PreTrendsPower | `pretrends.py` | `pretrends` package | **Complete** | 2026-05-19 | | PowerAnalysis | `power.py` | `pwr` / `DeclareDesign` | **Complete** | 2026-05-31 | -| PlaceboTests | `diagnostics.py` | Bertrand-Duflo-Mullainathan (2004) (placebo laws); no canonical R | **In Progress** | — | +| PlaceboTests | `diagnostics.py` | Bertrand-Duflo-Mullainathan (2004); base-R exact-enumeration parity | **Complete** | 2026-06-26 | ### Cross-Cutting Inference Features @@ -1314,20 +1314,43 @@ CI and extending covariate-adjusted R parity are tracked follow-ups in `TODO.md` | Field | Value | |-------|-------| | Module | `diagnostics.py` | -| Primary Reference | Bertrand, Duflo & Mullainathan (2004), QJE 119(1):249-275 (placebo laws / randomization inference). Paper review on file: `docs/methodology/papers/bertrand-duflo-mullainathan-2004-review.md`. | -| R Reference | None canonical (no R package ships a generic placebo battery) | -| Status | **In Progress** | -| Last Review | — | +| Primary Reference | Bertrand, Duflo & Mullainathan (2004), QJE 119(1):249-275 (placebo laws / randomization inference). Paper review: `docs/methodology/papers/bertrand-duflo-mullainathan-2004-review.md`. | +| R Reference | No canonical R package for the DiD placebo battery; parity via base-R `combn` exact enumeration (`benchmarks/R/generate_placebo_golden.R`), optional `ri2`/`coin` convention check. | +| Status | **Complete** | +| Last Review | 2026-06-26 | -**Documentation in place:** -- REGISTRY.md section: `## PlaceboTests` (NaN-inference edge cases for `permutation_test` and `leave_one_out_test`) -- Paper review: `docs/methodology/papers/bertrand-duflo-mullainathan-2004-review.md` (BDM 2004 placebo-law / serial-correlation grounding; proposes a `## PlaceboTests` REGISTRY entry, not yet integrated) -- Implementation: tests embedded in `tests/test_diagnostics.py` +**Verified Components:** +- [x] Sampled randomization-inference p-value `(1 + count)/(B + 1)` — valid but slightly conservative (with-replacement Monte Carlo; Phipson & Smyth 2010; BDM fn 11), converging to the exact `count/total` full enumeration — `tests/test_methodology_placebo.py::TestPlaceboRandomizationInference` +- [x] R parity: Python exhaustive enumeration matches base-R `combn` exact p-value + observed ATT at `atol=1e-12`; deterministic `leave_one_out_test` / `placebo_group_test` match R at `atol≈1e-10`; sampled `permutation_test` within Monte-Carlo tolerance — `TestPlaceboParityR` (skip-guarded; golden `benchmarks/data/placebo_golden.json`) +- [x] `placebo_timing_test` detects differential pre-trends (significant under violated trends, null under parallel) and restricts to pre-treatment data — `TestPlaceboFakeTiming` +- [x] `placebo_group_test` never-treated `treatment` filter (drops ever-treated), degenerate-design `ValueError`, misuse `UserWarning`, backward-compatible without `treatment` — `TestPlaceboFakeGroup` +- [x] Permutation NaN-decoupling contract (RI p-value finite when `se` degenerate) + fail-closed `RuntimeError` on all-fail — `TestPlaceboInferenceContracts` +- [x] Functional coverage (dispatch routing, zero-SE / `<2`-LOO NaN-inference) — `tests/test_diagnostics.py` -**Outstanding for promotion:** -- Standalone-vs-absorb decision: **resolved — standalone** (`diagnostics.py` is an exported public surface distinct from per-estimator placebo/LOO) -- Integrate the proposed `## PlaceboTests` entry into REGISTRY.md (cite BDM 2004 + scope) and flip this row to Complete -- Dedicated `tests/test_methodology_placebo.py` with BDM-anchored Verified Components (empirical permutation p-value per fn 12; p-value floor; LOO; fake-timing/fake-group) + Deviations block (permutation path's deliberate non-`safe_inference` + percentile CI; the NaN-inference convention). R parity is N/A (no canonical R placebo battery) → self-consistency / analytic anchors +**Test Coverage:** +- `tests/test_methodology_placebo.py` (24 tests across 5 paper-anchored classes; 2 `@pytest.mark.slow`) +- `tests/test_diagnostics.py` (32 functional / edge-case tests) + +**R Comparison Results** (base-R `combn` exact enumeration; golden committed at `benchmarks/data/placebo_golden.json`): + +| Quantity | Tolerance | Result | +|----------|-----------|--------| +| Observed DiD ATT | `atol=1e-12` | match | +| Exact RI p-value (`count/total`) | `atol=1e-12` | match | +| Leave-one-out mean / se / CI / per-drop ATTs | `atol=1e-10` / `1e-9` | match | +| Fake-group ATT (never-treated filtered) | `atol=1e-10` | match | +| Sampled permutation p-value | Monte-Carlo | converges to exact | + +**Corrections Made:** +- **Permutation p-value (this PR):** replaced `count/B` floored at `1/(B+1)` with the Phipson-Smyth (2010) randomization-inference value `(1 + count)/(B + 1)` — a valid but slightly conservative estimator for with-replacement Monte-Carlo draws (floor now intrinsic); "exact" is reserved for the full enumeration (`diff_diff/diagnostics.py`). +- **`placebo_group_test` (this PR):** added an optional `treatment` parameter that drops ever-treated units so the placebo runs on never-treated data only (uncontaminated); added degenerate-design `ValueError` guards (replacing a cryptic `LinAlgError`) and a misuse `UserWarning`; corrected the docstring to describe both modes. The `run_placebo_test` dispatcher's `fake_group` path now filters ever-treated units by default (a more-correct placebo) — documented in CHANGELOG + REGISTRY. +- **Docstring (this PR):** reworded `permutation_test`'s docstring — dropped the "exact ... valid with any sample size" overclaim; the sampled value is the valid/conservative with-replacement RI p-value, with "exact" reserved for full enumeration (BDM fn 12). + +**Deviations:** (documented in REGISTRY `## PlaceboTests`) +- Permutation inference deliberately bypasses `safe_inference`: RI p-value + null-distribution percentile interval (not an effect CI) + null-mean `placebo_effect`. +- `leave_one_out_test` reports the dispersion of per-drop ATTs (a sensitivity spread), not a design-based jackknife SE. +- Permutation NaN-decoupling: the count-based RI p-value stays valid when the permutation `se` is degenerate (intentional departure from the bootstrap-NaN contract). +- BDM's serial-correlation SE corrections (parametric AR, block bootstrap, cluster VCV, aggregation) are out of scope for this diagnostic surface. --- @@ -1459,17 +1482,11 @@ whereas R's `did::att_gt` would error. This is a defensive enhancement that prov more graceful handling of edge cases while still signaling invalid inference to users. ``` -### Priority Order (2026-05-26) - -Promotion priority for the **In Progress** entries, ordered by what's blocked on substantive review work (top of list = needs review next) vs. consolidation pass (bottom of list = mostly tracker walk-through): - -**Substantive-review-blocked (each still missing one or more of: a methodology test file, R parity, or a paper review):** - -1. **PlaceboTests** — standalone-vs-absorb decision resolved (standalone) and the BDM (2004) paper review is now on file (`docs/methodology/papers/bertrand-duflo-mullainathan-2004-review.md`). Remaining for promotion: the dedicated methodology test file + REGISTRY integration (R parity N/A). Methodologically lightweight. +### Priority Order (updated 2026-06-26) -**Consolidation-pass-blocked (already has paper review or methodology file or R parity; mostly Verified Components walk-through):** +Only one **In Progress** entry remains: **Survey Data Support**. (PlaceboTests was promoted to Complete on 2026-06-26 — see its detail section above.) -4. **Survey Data Support** — cross-cutting feature; promotion requires the per-estimator integration paths to be locked down first. +- **Survey Data Support** — cross-cutting feature; consolidation-pass-blocked. Promotion requires the per-estimator integration paths to be locked down first, then a dedicated `tests/test_methodology_survey.py` (Binder-equation-numbered Verified Components), an R-parity table vs `survey::svyglm`/`svycontrast` wired into this tracker, a deviations block, and a consolidated cross-estimator `NotImplementedError`-gaps enumeration. --- diff --git a/benchmarks/R/generate_placebo_golden.R b/benchmarks/R/generate_placebo_golden.R new file mode 100644 index 00000000..8349a5cb --- /dev/null +++ b/benchmarks/R/generate_placebo_golden.R @@ -0,0 +1,133 @@ +#!/usr/bin/env Rscript +# Golden generator: PlaceboTests (diff_diff/diagnostics.py) R parity. +# +# Bertrand, Duflo & Mullainathan (2004) placebo-law / randomization-inference +# diagnostics. Writes a fixed, fully deterministic 2-period panel and the R +# reference values so tests/test_methodology_placebo.py::TestPlaceboParityR can +# pin Python output against R without requiring R at test time. +# +# Outputs (checked into the repo): +# benchmarks/data/placebo_test_panel.csv (unit, t, y, treatment) +# benchmarks/data/placebo_golden.json +# +# Usage: +# Rscript benchmarks/R/generate_placebo_golden.R +# +# Notes: +# - The panel is HARDCODED (not RNG-generated) so R and Python consume bit- +# identical data; no cross-language RNG matching is needed. +# - Permutation p-value uses EXACT enumeration of all C(8, 3) = 56 treated-group +# assignments (the observed assignment is one of them): exact p = +# #{|ATT*| >= |ATT_obs|} / total (observed included; min 1/total). This is the +# ground truth the library's SAMPLED (1 + count)/(B + 1) value converges to. +# - n_treated = 3 != N/2 = 4, so no assignment's complement shares its |ATT| +# (avoids exact-tie pairing); the panel is chosen with a clear boundary gap +# so the 1e-12 exact-parity comparison is not tie-flip fragile. +# - leave-one-out se is the dispersion (sd, ddof=1) of the per-drop ATTs (NOT a +# design-based jackknife SE), with a t-distribution (df = n_valid - 1), exactly +# matching diff_diff.leave_one_out_test via safe_inference. +# - Optional ri2/coin convention cross-check is guarded by requireNamespace and +# is NOT a committed dependency (base-R combn enumeration is the anchor). + +suppressMessages(library(jsonlite)) + +# ---- fixed panel (8 units x 2 periods; real treated = units 0,1,2) ---- +panel <- data.frame( + unit = rep(0:7, each = 2), + t = rep(c(0, 1), times = 8), + y = c( + -1.639137, -0.623634, 0.051834, 1.622805, 0.261434, 0.82986, + 0.337559, 1.580412, -1.055892, -1.067745, 1.062855, 1.478681, + 0.139217, 0.8575, -1.253286, -0.560034 + ) +) +panel$treatment <- as.integer(panel$unit %in% c(0, 1, 2)) +real_treated <- c(0, 1, 2) +n_treated <- 3L +units <- 0:7 + +# 2x2 DiD ATT = double difference of group means (post = t == 1). +did_att <- function(df, treated) { + is_t <- df$unit %in% treated + post <- df$t == 1 + (mean(df$y[is_t & post]) - mean(df$y[is_t & !post])) - + (mean(df$y[!is_t & post]) - mean(df$y[!is_t & !post])) +} + +att_obs <- did_att(panel, real_treated) + +# ---- permutation: EXACT randomization-inference p-value ---- +combos <- combn(units, n_treated, simplify = FALSE) +atts <- vapply(combos, function(s) did_att(panel, s), numeric(1)) +count <- sum(abs(atts) >= abs(att_obs) - 1e-12) +total <- length(atts) +p_exact <- count / total +# boundary gap: nearest distinct |ATT*| to |ATT_obs| (excluding the observed) +gap <- sort(abs(abs(atts) - abs(att_obs)))[2] + +# ---- leave-one-out (deterministic jackknife over treated units) ---- +loo_units <- real_treated +loo_atts <- sapply(loo_units, function(u) { + remaining <- panel[panel$unit != u, ] + treated_rem <- setdiff(real_treated, u) + did_att(remaining, treated_rem) +}) +loo_mean <- mean(loo_atts) +loo_se <- sd(loo_atts) # ddof = 1, the dispersion of LOO ATTs (not an SE-of-mean) +loo_df <- length(loo_atts) - 1L +loo_t <- loo_mean / loo_se +loo_p <- 2 * pt(-abs(loo_t), df = loo_df) +loo_crit <- qt(0.975, df = loo_df) +loo_ci <- c(loo_mean - loo_crit * loo_se, loo_mean + loo_crit * loo_se) + +# ---- fake-group (deterministic; drop ever-treated, fake-treat controls 3,4) ---- +fg_fake_treated <- c(3, 4) +fg_panel <- panel[!(panel$unit %in% real_treated), ] # never-treated only +fg_att <- did_att(fg_panel, fg_fake_treated) + +# ---- optional convention cross-check (NOT a committed dependency) ---- +ri2_ok <- requireNamespace("ri2", quietly = TRUE) + +golden <- list( + description = "PlaceboTests R parity (BDM 2004): exact RI permutation p-value + deterministic LOO + fake-group, on a fixed 2-period panel.", + panel_csv = "benchmarks/data/placebo_test_panel.csv", + real_treated = real_treated, + n_treated = n_treated, + observed_att = att_obs, + permutation = list( + convention = "exact enumeration: p = #{|ATT*| >= |ATT_obs|} / total (observed included)", + count = count, + total = total, + p_exact = p_exact, + boundary_gap = gap + ), + leave_one_out = list( + dropped_units = loo_units, + per_drop_att = as.list(setNames(loo_atts, as.character(loo_units))), + mean = loo_mean, + se = loo_se, + df = loo_df, + t_stat = loo_t, + p_value = loo_p, + ci_lower = loo_ci[1], + ci_upper = loo_ci[2] + ), + fake_group = list( + fake_treated_units = fg_fake_treated, + note = "ever-treated units dropped (treatment filter); ATT is the double-difference", + att = fg_att + ), + ri2_convention_checked = ri2_ok +) + +write.csv(panel, "benchmarks/data/placebo_test_panel.csv", row.names = FALSE) +write_json(golden, "benchmarks/data/placebo_golden.json", + auto_unbox = TRUE, pretty = TRUE, digits = 12 +) + +cat(sprintf("observed ATT = %.12f\n", att_obs)) +cat(sprintf("exact RI: count=%d total=%d p_exact=%.12f gap=%.4f\n", count, total, p_exact, gap)) +cat(sprintf("LOO: mean=%.12f se=%.12f df=%d p=%.6f\n", loo_mean, loo_se, loo_df, loo_p)) +cat(sprintf("fake_group ATT = %.12f\n", fg_att)) +cat(sprintf("ri2 convention cross-check available: %s\n", ri2_ok)) +cat("Wrote benchmarks/data/placebo_test_panel.csv + placebo_golden.json\n") diff --git a/benchmarks/data/placebo_golden.json b/benchmarks/data/placebo_golden.json new file mode 100644 index 00000000..90658320 --- /dev/null +++ b/benchmarks/data/placebo_golden.json @@ -0,0 +1,35 @@ +{ + "description": "PlaceboTests R parity (BDM 2004): exact RI permutation p-value + deterministic LOO + fake-group, on a fixed 2-period panel.", + "panel_csv": "benchmarks/data/placebo_test_panel.csv", + "real_treated": [0, 1, 2], + "n_treated": 3, + "observed_att": 0.4399611333333, + "permutation": { + "convention": "exact enumeration: p = #{|ATT*| >= |ATT_obs|} / total (observed included)", + "count": 15, + "total": 56, + "p_exact": 0.2678571428571, + "boundary_gap": 0.03574946666667 + }, + "leave_one_out": { + "dropped_units": [0, 1, 2], + "per_drop_att": { + "0": 0.4580263, + "1": 0.1802923, + "2": 0.6815648 + }, + "mean": 0.4399611333333, + "se": 0.2511240579855, + "df": 2, + "t_stat": 1.751967282079, + "p_value": 0.2218771530393, + "ci_lower": -0.6405384802636, + "ci_upper": 1.52046074693 + }, + "fake_group": { + "fake_treated_units": [3, 4], + "note": "ever-treated units dropped (treatment filter); ATT is the double-difference", + "att": 0.006379666666667 + }, + "ri2_convention_checked": false +} diff --git a/benchmarks/data/placebo_test_panel.csv b/benchmarks/data/placebo_test_panel.csv new file mode 100644 index 00000000..f5955687 --- /dev/null +++ b/benchmarks/data/placebo_test_panel.csv @@ -0,0 +1,17 @@ +"unit","t","y","treatment" +0,0,-1.639137,1 +0,1,-0.623634,1 +1,0,0.051834,1 +1,1,1.622805,1 +2,0,0.261434,1 +2,1,0.82986,1 +3,0,0.337559,0 +3,1,1.580412,0 +4,0,-1.055892,0 +4,1,-1.067745,0 +5,0,1.062855,0 +5,1,1.478681,0 +6,0,0.139217,0 +6,1,0.8575,0 +7,0,-1.253286,0 +7,1,-0.560034,0 diff --git a/diff_diff/diagnostics.py b/diff_diff/diagnostics.py index 509634c4..97be9bb3 100644 --- a/diff_diff/diagnostics.py +++ b/diff_diff/diagnostics.py @@ -19,7 +19,7 @@ from diff_diff.estimators import DifferenceInDifferences from diff_diff.results import _get_significance_stars -from diff_diff.utils import safe_inference +from diff_diff.utils import safe_inference, validate_binary @dataclass @@ -228,7 +228,7 @@ def run_placebo_test( test_type : str, default="fake_timing" Type of placebo test: - "fake_timing": Assign treatment at a fake (earlier) time period - - "fake_group": Run DiD designating some control units as "fake treated" + - "fake_group": Designate control units as "fake treated" (real-treated units, per the ``treatment`` column, are dropped first) - "permutation": Randomly reassign treatment and compute distribution - "leave_one_out": Drop each treated unit and re-estimate fake_treatment_period : any, optional @@ -313,6 +313,7 @@ def run_placebo_test( fake_treated_units=fake_treatment_group, post_periods=post_periods, alpha=alpha, + treatment=treatment, **estimator_kwargs, ) @@ -445,14 +446,20 @@ def placebo_group_test( fake_treated_units: List[Any], post_periods: Optional[List[Any]] = None, alpha: float = 0.05, + treatment: Optional[str] = None, **estimator_kwargs, ) -> PlaceboTestResults: """ - Test for differential trends among never-treated units. + Test for differential trends by designating control units as "fake treated". + + Designates ``fake_treated_units`` as fake-treated and estimates a DiD on the + resulting panel. A significant effect suggests heterogeneous trends in the + control group (a parallel-trends red flag). - Assigns some never-treated units as "fake treated" and estimates a - DiD model using only never-treated data. A significant effect suggests - heterogeneous trends in the control group. + If ``treatment`` is provided, units that are *ever* really treated are dropped + first, so the placebo runs on never-treated units only (the recommended, + uncontaminated design). If ``treatment`` is ``None``, the test runs on whatever + data is supplied, so the caller must pass control-only data for a valid placebo. Parameters ---------- @@ -470,6 +477,11 @@ def placebo_group_test( List of post-treatment period values. alpha : float, default=0.05 Significance level. + treatment : str, optional + Real treatment-indicator column. When given, units that are ever + real-treated (``data.groupby(unit)[treatment].max() == 1``) are dropped + before the placebo, so it runs on never-treated units only. When ``None`` + (default), no filtering is done and the caller must pass control-only data. **estimator_kwargs Arguments passed to DifferenceInDifferences. @@ -481,7 +493,35 @@ def placebo_group_test( if fake_treated_units is None or len(fake_treated_units) == 0: raise ValueError("fake_treated_units must be a non-empty list") - all_periods = sorted(data[time].unique()) + fake_data = data.copy() + + # Optionally restrict to never-treated units so the placebo is not contaminated + # by the real treatment effect (the BDM 2004 placebo-law design on controls). + if treatment is not None: + # Fail closed: a missing column or non-0/1 values would otherwise silently + # skip the ever-treated filter (groupby().max() drops NaN), running the + # placebo on contaminated data. + if treatment not in fake_data.columns: + raise ValueError(f"treatment column '{treatment}' not found in data") + if fake_data[treatment].isna().any(): + raise ValueError(f"treatment column '{treatment}' contains missing values") + validate_binary(fake_data[treatment].to_numpy(), "treatment") + ever_treated = fake_data.groupby(unit)[treatment].max() + ever_treated_units = set(ever_treated[ever_treated == 1].index) + misused = [u for u in fake_treated_units if u in ever_treated_units] + if misused: + import warnings + + warnings.warn( + f"{len(misused)} of fake_treated_units are themselves ever real-treated " + f"and will be dropped with the other real-treated units: {misused}. " + f"Pass only never-treated units as fake_treated_units for a valid placebo.", + UserWarning, + stacklevel=2, + ) + fake_data = fake_data[~fake_data[unit].isin(ever_treated_units)].copy() + + all_periods = sorted(fake_data[time].unique()) # Infer post periods if not provided if post_periods is None: @@ -489,14 +529,31 @@ def placebo_group_test( post_periods = all_periods[mid:] # Create fake treatment indicator - fake_data = data.copy() fake_data["_fake_treated"] = fake_data[unit].isin(fake_treated_units).astype(int) fake_data["_post"] = fake_data[time].isin(post_periods).astype(int) + # Guard degenerate designs (e.g., all fake_treated_units were dropped as + # real-treated, or no controls remain) before they surface as a cryptic + # LinAlgError inside the estimator. + if fake_data["_fake_treated"].sum() == 0: + raise ValueError( + "No fake-treated observations remain (all fake_treated_units were " + "dropped as real-treated, or are absent from the data). Pass " + "never-treated units as fake_treated_units." + ) + if (fake_data["_fake_treated"] == 0).sum() == 0: + raise ValueError("No control (non-fake-treated) units remain for the placebo comparison.") + # Fit DiD did = DifferenceInDifferences(**estimator_kwargs) results = did.fit(fake_data, outcome=outcome, treatment="_fake_treated", time="_post") + # Record the fake-treated units actually used (after any never-treated + # filtering), not just the originally requested list, to avoid metadata drift. + # Preserve the caller's order (sorting could raise TypeError on mixed-type IDs). + retained = set(fake_data.loc[fake_data["_fake_treated"] == 1, unit].unique()) + used_fake_treated = [u for u in fake_treated_units if u in retained] + return PlaceboTestResults( test_type="fake_group", placebo_effect=results.att, @@ -507,7 +564,7 @@ def placebo_group_test( n_obs=results.n_obs, is_significant=bool(results.p_value < alpha), alpha=alpha, - fake_group=list(fake_treated_units), + fake_group=used_fake_treated, ) @@ -526,8 +583,12 @@ def permutation_test( Compute permutation-based p-value for DiD estimate. Randomly reassigns treatment status at the unit level and computes the - DiD estimate for each permutation. The p-value is the proportion of - permuted estimates at least as extreme as the original. + DiD estimate for each permutation. The p-value is the randomization-inference + value ``(1 + count) / (B + 1)`` (Phipson & Smyth 2010), where ``count`` is the + number of permuted estimates at least as extreme as the observed and ``B`` is + the number of valid permutations. With ``B`` sampled permutations this is a + Monte-Carlo approximation that converges to the exact full-enumeration value + ``count / total`` as ``B`` grows. Parameters ---------- @@ -557,8 +618,17 @@ def permutation_test( Notes ----- - The permutation test is exact and does not rely on asymptotic - approximations, making it valid with any sample size. + This is a randomization-inference (permutation) test of the sharp null of no + effect for any unit; it does not rely on asymptotic approximations. Treatment + assignments are drawn independently each iteration (Monte-Carlo sampling *with + replacement* from the assignment space), so the reported p-value + ``(1 + count) / (B + 1)`` (Phipson & Smyth 2010) is a **valid but slightly + conservative** estimator -- the ``+1`` adds the observed assignment and + prevents a zero p-value. Here ``count`` is the number of permutations at least + as extreme as the observed estimate and ``B`` is the number of valid + permutations. As ``B`` grows it converges to the *exact* p-value obtained by + full enumeration of all assignments (the R-parity reference). "Exact" is + reserved for that full enumeration; the sampled value approximates it. """ rng = np.random.default_rng(seed) @@ -620,11 +690,12 @@ def permutation_test( stacklevel=2, ) - # Compute p-value: proportion of |permuted| >= |original| - p_value = np.mean(np.abs(valid_effects) >= np.abs(original_att)) - - # Ensure p-value is at least 1/(n_permutations + 1) - p_value = max(p_value, 1 / (len(valid_effects) + 1)) + # Randomization-inference p-value (Phipson & Smyth 2010): include the observed + # statistic in both numerator and denominator. The 1/(B+1) floor is intrinsic + # (count == 0 -> 1/(B+1)), so no separate clamp is needed. With sampled + # permutations this converges to the exact full-enumeration value count/total. + count = int(np.sum(np.abs(valid_effects) >= np.abs(original_att))) + p_value = (1 + count) / (len(valid_effects) + 1) # Compute SE and CI from permutation distribution se = np.std(valid_effects, ddof=1) diff --git a/docs/doc-deps.yaml b/docs/doc-deps.yaml index 1ed50392..7fd25a50 100644 --- a/docs/doc-deps.yaml +++ b/docs/doc-deps.yaml @@ -732,8 +732,11 @@ sources: drift_risk: low docs: - path: docs/methodology/REGISTRY.md - section: "Diagnostics" + section: "PlaceboTests" type: methodology + - path: docs/methodology/papers/bertrand-duflo-mullainathan-2004-review.md + type: methodology + note: "BDM (2004) placebo-law / randomization-inference grounding for PlaceboTests" - path: docs/api/diagnostics.rst type: api_reference - path: docs/tutorials/04_parallel_trends.ipynb diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index ad48999d..1b243da1 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -3097,13 +3097,56 @@ Shipped in `diff_diff/had_pretests.py` as `stute_joint_pretest()` (residuals-in ## PlaceboTests +**Primary source:** [Bertrand, M., Duflo, E., & Mullainathan, S. (2004). How Much Should We Trust Differences-in-Differences Estimates? *The Quarterly Journal of Economics*, 119(1), 249-275.](https://doi.org/10.1162/003355304772839588). Paper review on file: `docs/methodology/papers/bertrand-duflo-mullainathan-2004-review.md`. + **Module:** `diff_diff/diagnostics.py` +**Scope:** A battery of placebo / randomization-inference diagnostics for the parallel-trends assumption, built on the base 2×2 `DifferenceInDifferences`. BDM (2004) introduce the **placebo-law experiment** — randomly assign a fake treatment date and/or fake treated group, estimate the DD, and check whether a (necessarily spurious) "effect" is significant more than ~5% of the time. These are **diagnostics**, not estimators, and are distinct from the per-estimator placebo/LOO surfaces documented elsewhere (SyntheticDiD donor leave-one-out per ADH 2015 §4, HAD pre-tests, DCDH placebo). Public API: `run_placebo_test` (dispatcher) → `placebo_timing_test` / `placebo_group_test` / `permutation_test` / `leave_one_out_test`, plus `run_all_placebo_tests` and `PlaceboTestResults`. + +**Key implementation requirements:** + +*The four diagnostics:* +- **`placebo_timing_test` (fake timing):** restrict to pre-treatment periods, assign a fake post indicator (`time >= fake_treatment_period`), and run the DD on the real treated vs control units. A significant effect flags differential pre-trends. Requires ≥2 pre-periods; a `fake_treatment_period` inside `post_periods` raises `ValueError`. +- **`placebo_group_test` (fake group):** designate control units as fake-treated. With the optional `treatment=` column, units that are *ever* real-treated are dropped first (placebo on never-treated units only, uncontaminated by the real effect); without it, the caller must pass control-only data. Degenerate designs (all fake-treated dropped, or no controls remain) raise `ValueError`; a fake-treated unit that is itself real-treated emits a `UserWarning`. Via the `run_placebo_test` dispatcher (which always has the `treatment` column) the `fake_group` path filters ever-treated units by default. +- **`permutation_test` (randomization inference):** randomly reassign the treated-group label across units (BDM placebo-law over a fixed outcome set; the randomization-inference link is BDM fn 11, citing Rosenbaum 1996) and form the null distribution of the DD estimate. +- **`leave_one_out_test`:** drop each treated unit, refit, and report the per-drop ATTs (single-unit sensitivity). + +*Randomization-inference p-value (Phipson & Smyth 2010):* + +``` +p = (1 + #{ |ATT*_b| >= |ATT_obs| }) / (B + 1) +``` + +where `ATT*_b` are the `B` valid permutation estimates and the `+1` includes the observed statistic in both numerator and denominator (intrinsic floor `1/(B+1)`). Assignments are drawn independently each iteration (Monte-Carlo sampling **with replacement** from the assignment space), so this is the Phipson & Smyth (2010) **valid but slightly conservative** RI p-value — *not* an exact finite-sample value. The term **exact** is reserved for the full enumeration of all `C(N, n_treated)` assignments (observed included), `#{|ATT*| >= |ATT_obs|} / total`, to which the sampled value converges as `B → ∞` (this enumeration is the R-parity reference). + +*Standard errors / inference (stated honestly):* +- `placebo_timing_test` / `placebo_group_test` surface the underlying `DifferenceInDifferences` estimator's own jointly-computed inference (HC1 default via `safe_inference`). +- `permutation_test` does **not** use `safe_inference`: the p-value is the count-based RI value above, the confidence interval is the **percentile interval of the null (permutation) distribution — not an effect CI**, and `t_stat = original_att / se` is computed individually (`se` = std of the null distribution). `placebo_effect` reports the **mean of the null distribution** (≈0), with `original_effect` holding the observed ATT. +- `leave_one_out_test` uses `safe_inference` (t-distribution, `df = n_valid − 1`), but its `se` is the **dispersion (sample std) of the leave-one-out ATTs — a sensitivity spread, not a design-based jackknife SE**; the per-unit `leave_one_out_effects` dict is the primary output and `t_stat`/`p_value`/CI are heuristic. + *Edge cases:* +- **Permutation NaN-decoupling (deliberate):** the RI p-value is count-based and stays finite even when the permutation `se` is degenerate (`se == 0` for identical permutations, or `se` NaN at `n_valid == 1`), in which case `t_stat` is NaN. This intentionally departs from the bootstrap-NaN contract (non-finite SE → full NaN tuple), because the RI p-value does not depend on `se` (BDM fn 12: the placebo reference distribution is not standard normal, so the count-based RI p-value — not an `se`-based statistic — is the valid one). **Note:** intentional contract; see `tests/test_methodology_placebo.py::TestPlaceboInferenceContracts`. - NaN inference for undefined statistics: - - `permutation_test`: t_stat is NaN when permutation SE is zero (all permutations produce identical estimates) - - `leave_one_out_test`: t_stat, p_value, CI are NaN when LOO SE is zero (all LOO effects identical) - - **Note**: Defensive enhancement matching CallawaySantAnna NaN convention + - `permutation_test`: t_stat is NaN when the permutation SE is zero/degenerate (the p-value remains valid — see decoupling above). + - `leave_one_out_test`: t_stat, p_value, CI are NaN when the LOO SE is zero (all LOO effects identical) or `< 2` valid effects. + - **Note**: Defensive enhancement matching CallawaySantAnna NaN convention. +- Fail-closed: `permutation_test` and `leave_one_out_test` raise `RuntimeError` if **all** refits fail; `permutation_test` emits a `UserWarning` when `> 10%` of permutations fail. + +**Reference implementation(s):** +- No single canonical R/Stata package implements the DiD placebo *battery* as one command (BDM's own code is custom). R parity is anchored by **exact enumeration** in base R (`combn`) of the full randomization distribution, with an optional `ri2`/`coin` convention cross-check (guarded by `requireNamespace`, not a committed dependency). Generator: `benchmarks/R/generate_placebo_golden.R`; golden: `benchmarks/data/placebo_golden.json` (+ `placebo_test_panel.csv`). The deterministic `leave_one_out_test` / `placebo_group_test` / observed ATT match R exactly (`atol≈1e-10`); the sampled `permutation_test` matches the exact value within Monte-Carlo tolerance. + +**Deviations:** +- **Note (permutation inference is not `safe_inference`):** the permutation path uses the RI p-value + a null-distribution percentile interval (not an effect CI) and a null-mean `placebo_effect`. This is the standard randomization-inference convention and deliberately differs from the project-wide `safe_inference` t-based inference (which assumes a sampling-distribution SE the permutation test does not have). +- **Note (leave-one-out spread):** `leave_one_out_test`'s reported `se`/`t_stat`/`p_value`/CI summarize the dispersion of the per-drop ATTs (a robustness signal), not a design-based jackknife standard error; the per-unit effects are the primary output. +- **Note (BDM scope):** BDM (2004) is primarily about serial-correlation-robust *standard errors* (parametric AR, block bootstrap, cluster/arbitrary VCV, time aggregation). Those inference corrections are out of scope for this diagnostic surface — diff-diff's cluster-robust SE and bootstrap paths cover them under the relevant estimators. This entry covers only the placebo-law / randomization-inference diagnostics that BDM motivate. + +**Requirements checklist:** +- [x] RI p-value `(1+count)/(B+1)` (sampled) converges to the exact enumeration `count/total`: `tests/test_methodology_placebo.py::TestPlaceboRandomizationInference` +- [x] R parity: exact enumeration + observed ATT at `atol=1e-12`; deterministic LOO / fake-group at `atol=1e-10`; sampled permutation within MC tolerance: `TestPlaceboParityR` (skip-guarded) +- [x] Fake-timing detects differential pre-trends; null under parallel trends; pre-data only: `TestPlaceboFakeTiming` +- [x] Fake-group never-treated `treatment` filter + degenerate `ValueError` + misuse `UserWarning`; backward-compatible without `treatment`: `TestPlaceboFakeGroup` +- [x] Permutation NaN-decoupling + fail-closed `RuntimeError`: `TestPlaceboInferenceContracts` +- [x] Functional coverage (dispatch routing, zero-SE / `<2`-LOO NaN-inference): `tests/test_diagnostics.py` --- diff --git a/docs/references.rst b/docs/references.rst index 847ee7d5..05357a76 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -99,6 +99,14 @@ Placebo Tests and DiD Diagnostics - **Bertrand, M., Duflo, E., & Mullainathan, S. (2004).** "How Much Should We Trust Differences-in-Differences Estimates?" *The Quarterly Journal of Economics*, 119(1), 249-275. https://doi.org/10.1162/003355304772839588 +- **Rosenbaum, P. R. (1996).** "Observational Studies and Nonrandomized Experiments." In *Handbook of Statistics*, Vol. 13, S. Ghosh & C. R. Rao, eds., 181-197. Elsevier. https://doi.org/10.1016/S0169-7161(96)13009-0 + + Randomization-inference foundation for the ``PlaceboTests`` permutation test; cited by Bertrand-Duflo-Mullainathan (2004) footnote 11 as the basis for forming a randomization-inference test from the placebo-law distribution. + +- **Phipson, B., & Smyth, G. K. (2010).** "Permutation P-values Should Never Be Zero: Calculating Exact P-values When Permutations Are Randomly Drawn." *Statistical Applications in Genetics and Molecular Biology*, 9(1), Article 39. https://doi.org/10.2202/1544-6115.1585 + + Valid (slightly conservative) sampled randomization-inference p-value convention ``(1 + count) / (B + 1)`` used by ``permutation_test`` for randomly-drawn (with-replacement) permutations (the ``+1`` includes the observed statistic, giving an intrinsic ``1/(B+1)`` floor); the exact value is the full enumeration ``count/total``. + Synthetic Control Method ------------------------ diff --git a/tests/test_methodology_placebo.py b/tests/test_methodology_placebo.py new file mode 100644 index 00000000..2250bb1a --- /dev/null +++ b/tests/test_methodology_placebo.py @@ -0,0 +1,507 @@ +"""Methodology verification for PlaceboTests (``diff_diff/diagnostics.py``). + +Anchored to Bertrand, Duflo & Mullainathan (2004), "How Much Should We Trust +Differences-in-Differences Estimates?", QJE 119(1):249-275 — the placebo-law / +randomization-inference framework. Paper review on file: +``docs/methodology/papers/bertrand-duflo-mullainathan-2004-review.md``. + +These tests COMPLEMENT ``tests/test_diagnostics.py`` (which covers functional +behavior, dispatch routing, and the zero-SE / <2-LOO NaN-inference edge cases). +Here we verify the methodology itself: + +- the randomization-inference p-value convention ``(1 + count)/(B + 1)`` + (Phipson & Smyth 2010; BDM fn 11), converging to the exact full enumeration, +- R parity via exact enumeration + deterministic leave-one-out / fake-group, +- the fake-timing pre-trends logic, +- the never-treated ``treatment`` filter on ``placebo_group_test``, +- the deliberate permutation NaN-decoupling contract (a valid RI p-value even + when the permutation SE is degenerate; BDM fn 12). +""" + +import itertools +import json +import warnings +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from diff_diff.diagnostics import ( + leave_one_out_test, + permutation_test, + placebo_group_test, + placebo_timing_test, + run_placebo_test, +) + +_DATA_DIR = Path(__file__).resolve().parent.parent / "benchmarks" / "data" +_GOLDEN_PATH = _DATA_DIR / "placebo_golden.json" +_PANEL_PATH = _DATA_DIR / "placebo_test_panel.csv" + + +# --------------------------------------------------------------------------- # +# Helpers +# --------------------------------------------------------------------------- # +def _fixed_panel() -> pd.DataFrame: + """The committed 8-unit x 2-period golden panel, hardcoded (no fixture dep).""" + y = [ + -1.639137, -0.623634, 0.051834, 1.622805, 0.261434, 0.82986, + 0.337559, 1.580412, -1.055892, -1.067745, 1.062855, 1.478681, + 0.139217, 0.8575, -1.253286, -0.560034, + ] # fmt: skip + return pd.DataFrame( + { + "unit": np.repeat(range(8), 2), + "t": list(range(2)) * 8, + "y": y, + "treatment": np.repeat([1, 1, 1, 0, 0, 0, 0, 0], 2), + } + ) + + +def _did_att(df, treated, outcome="y", unit="unit", time="t") -> float: + """Closed-form 2-period 2x2 DiD ATT (double difference of group means).""" + is_t = df[unit].isin(treated) + post = df[time] == 1 + return float( + (df.loc[is_t & post, outcome].mean() - df.loc[is_t & ~post, outcome].mean()) + - (df.loc[~is_t & post, outcome].mean() - df.loc[~is_t & ~post, outcome].mean()) + ) + + +def _exact_ri(df, treated, n_treated, unit="unit"): + """Exact RI p-value over ALL C(N, n_treated) assignments (observed included). + + This is the ground truth the sampled ``(1 + count)/(B + 1)`` converges to; + here the observed assignment is one of the enumerated assignments, so the + convention is ``count / total`` (min ``1/total``). + """ + att_obs = _did_att(df, treated) + units = sorted(df[unit].unique()) + atts = np.array([_did_att(df, list(c)) for c in itertools.combinations(units, n_treated)]) + count = int(np.sum(np.abs(atts) >= np.abs(att_obs) - 1e-12)) + return att_obs, count, len(atts), count / len(atts) + + +def _make_timing_panel(violated: bool, seed: int = 7, n_per_group: int = 15): + """4-period panel; ``violated`` adds a differential pre-trend to treated units.""" + rng = np.random.default_rng(seed) + rows = [] + for u in range(2 * n_per_group): + treated = u < n_per_group + a = rng.normal(0, 1.0) + for t in range(4): + y = a + 0.4 * t + rng.normal(0, 0.1) + if violated and treated: + y += 1.2 * t # treated trend up faster, even pre-treatment + rows.append({"unit": u, "t": t, "y": y, "treatment": int(treated)}) + return pd.DataFrame(rows) + + +# --------------------------------------------------------------------------- # +# 1. Randomization-inference p-value convention (no R required) +# --------------------------------------------------------------------------- # +class TestPlaceboRandomizationInference: + """The sampled p-value is the Phipson-Smyth (2010) RI value ``(1 + count)/(B + 1)``. + + BDM (2004) fn 11 frames placebo laws as randomization inference; the p-value + convention is Phipson & Smyth (2010). Assignments are sampled with replacement, + so ``(1 + count)/(B + 1)`` is a *valid but slightly conservative* Monte-Carlo + p-value (not an exact finite-sample value); it converges to the exact + full-enumeration value ``count/total`` (see ``TestPlaceboParityR``). + """ + + def test_sampled_pvalue_matches_phipson_smyth_formula(self): + panel = _fixed_panel() + res = permutation_test( + panel, + outcome="y", + treatment="treatment", + time="t", + unit="unit", + n_permutations=200, + seed=123, + ) + dist = np.asarray(res.permutation_distribution) + b = len(dist) + assert res.original_effect is not None + count = int(np.sum(np.abs(dist) >= np.abs(res.original_effect))) + expected = (1 + count) / (b + 1) + assert res.p_value == pytest.approx(expected, abs=1e-12) + # With-replacement Monte Carlo: (1+count)/(B+1) is conservative vs the + # naive count/B (the +1/+1 only ever raises the p-value), and never zero. + assert res.p_value >= count / b - 1e-12 + assert res.p_value > 0.0 + + def test_with_replacement_sampling_tolerates_duplicate_assignments(self): + """Assignments are drawn independently (with replacement), so duplicates + (incl. the observed assignment) can recur across the B draws; the p-value + stays valid. This documents the conservative with-replacement contract.""" + # Tiny universe: C(4, 2) = 6 assignments, so B=50 draws must repeat some. + rows = [ + { + "unit": u, + "t": t, + "y": float(u) + 0.5 * t + 0.3 * (u % 2) * t, + "treatment": int(u < 2), + } + for u in range(4) + for t in range(2) + ] + panel = pd.DataFrame(rows) + res = permutation_test( + panel, + outcome="y", + treatment="treatment", + time="t", + unit="unit", + n_permutations=50, + seed=7, + ) + dist = np.asarray(res.permutation_distribution) + # duplicates present (6 distinct assignments, 50 draws) -> with replacement + assert len(np.unique(np.round(dist, 9))) < len(dist) + assert 0.0 < res.p_value <= 1.0 + + def test_pvalue_bounded_and_floored(self): + panel = _fixed_panel() + res = permutation_test( + panel, + outcome="y", + treatment="treatment", + time="t", + unit="unit", + n_permutations=200, + seed=1, + ) + b = len(np.asarray(res.permutation_distribution)) + assert 0.0 < res.p_value <= 1.0 + assert res.p_value >= 1.0 / (b + 1) - 1e-12 + + def test_exact_enumeration_is_count_over_total(self): + """From-scratch exhaustive enumeration uses ``count/total`` (observed incl.).""" + panel = _fixed_panel() + _, count, total, p_exact = _exact_ri(panel, [0, 1, 2], 3) + assert total == 56 # C(8, 3) + assert p_exact == pytest.approx(count / total, abs=1e-15) + # the observed assignment is always at least as extreme as itself + assert count >= 1 + + @pytest.mark.slow + def test_sampled_converges_to_exact(self, ci_params): + panel = _fixed_panel() + _, _, _, p_exact = _exact_ri(panel, [0, 1, 2], 3) + b = ci_params.bootstrap(20000, min_n=2000) + res = permutation_test( + panel, + outcome="y", + treatment="treatment", + time="t", + unit="unit", + n_permutations=b, + seed=20240101, + ) + assert res.p_value == pytest.approx(p_exact, abs=0.05) + + +# --------------------------------------------------------------------------- # +# 2. R parity (skip if golden JSON not committed) +# --------------------------------------------------------------------------- # +def _load_golden() -> dict: + if not _GOLDEN_PATH.exists() or not _PANEL_PATH.exists(): + pytest.skip( + f"Placebo R-parity goldens missing at {_DATA_DIR}. To regenerate: " + "`Rscript benchmarks/R/generate_placebo_golden.R`. The goldens are " + "committed by default; this skip covers partial checkouts." + ) + return json.loads(_GOLDEN_PATH.read_text()) + + +class TestPlaceboParityR: + """Parity with the R reference (base-R ``combn`` exact enumeration).""" + + @pytest.fixture + def golden(self) -> dict: + return _load_golden() + + @pytest.fixture + def panel(self) -> pd.DataFrame: + _load_golden() # trigger skip if missing + return pd.read_csv(_PANEL_PATH) + + def test_exact_enumeration_matches_r(self, golden, panel): + att_obs, count, total, p_exact = _exact_ri( + panel, golden["real_treated"], golden["n_treated"] + ) + assert att_obs == pytest.approx(golden["observed_att"], abs=1e-12) + assert count == golden["permutation"]["count"] + assert total == golden["permutation"]["total"] + assert p_exact == pytest.approx(golden["permutation"]["p_exact"], abs=1e-12) + + def test_leave_one_out_matches_r(self, golden, panel): + res = leave_one_out_test(panel, outcome="y", treatment="treatment", time="t", unit="unit") + gl = golden["leave_one_out"] + assert res.placebo_effect == pytest.approx(gl["mean"], abs=1e-10) + assert res.se == pytest.approx(gl["se"], abs=1e-10) + assert res.conf_int[0] == pytest.approx(gl["ci_lower"], abs=1e-9) + assert res.conf_int[1] == pytest.approx(gl["ci_upper"], abs=1e-9) + assert res.leave_one_out_effects is not None + for u, att in res.leave_one_out_effects.items(): + assert att == pytest.approx(gl["per_drop_att"][str(int(u))], abs=1e-10) + + def test_fake_group_matches_r(self, golden, panel): + fg = golden["fake_group"] + res = placebo_group_test( + panel, + outcome="y", + time="t", + unit="unit", + fake_treated_units=fg["fake_treated_units"], + treatment="treatment", + ) + assert res.placebo_effect == pytest.approx(fg["att"], abs=1e-10) + + @pytest.mark.slow + def test_sampled_permutation_matches_r_exact(self, golden, panel, ci_params): + b = ci_params.bootstrap(20000, min_n=2000) + res = permutation_test( + panel, + outcome="y", + treatment="treatment", + time="t", + unit="unit", + n_permutations=b, + seed=99, + ) + assert res.p_value == pytest.approx(golden["permutation"]["p_exact"], abs=0.05) + + +# --------------------------------------------------------------------------- # +# 3. Fake-timing pre-trends logic +# --------------------------------------------------------------------------- # +class TestPlaceboFakeTiming: + """``placebo_timing_test`` detects differential pre-trends (BDM placebo-law).""" + + def test_null_under_parallel_trends(self): + panel = _make_timing_panel(violated=False) + res = placebo_timing_test( + panel, + outcome="y", + treatment="treatment", + time="t", + fake_treatment_period=1, + post_periods=[2, 3], + ) + assert not res.is_significant + assert abs(res.placebo_effect) < 0.25 + + def test_detects_violated_pre_trends(self): + panel = _make_timing_panel(violated=True) + res = placebo_timing_test( + panel, + outcome="y", + treatment="treatment", + time="t", + fake_treatment_period=1, + post_periods=[2, 3], + ) + assert res.is_significant + assert res.placebo_effect > 0.8 # treated rose faster pre-treatment + + def test_uses_only_pre_treatment_data(self): + panel = _make_timing_panel(violated=False) + res = placebo_timing_test( + panel, + outcome="y", + treatment="treatment", + time="t", + fake_treatment_period=1, + post_periods=[2, 3], + ) + # pre-periods are {0, 1}: 30 units x 2 periods = 60 rows + assert res.n_obs == 60 + + def test_post_period_fake_date_rejected(self): + panel = _make_timing_panel(violated=False) + with pytest.raises(ValueError, match="pre-treatment"): + placebo_timing_test( + panel, + outcome="y", + treatment="treatment", + time="t", + fake_treatment_period=3, + post_periods=[2, 3], + ) + + +# --------------------------------------------------------------------------- # +# 4. Fake-group test + never-treated filter +# --------------------------------------------------------------------------- # +class TestPlaceboFakeGroup: + """``placebo_group_test`` and its optional never-treated ``treatment`` filter.""" + + def _control_panel(self, seed=3, n=10): + rng = np.random.default_rng(seed) + rows = [] + for u in range(n): + a = rng.normal(0, 1.0) + for t in range(2): + rows.append({"unit": u, "t": t, "y": a + 0.3 * t + rng.normal(0, 0.1)}) + return pd.DataFrame(rows) + + def test_null_on_valid_control_only_design(self): + panel = self._control_panel() + res = placebo_group_test( + panel, outcome="y", time="t", unit="unit", fake_treated_units=[0, 1, 2, 3, 4] + ) + assert not res.is_significant + + def test_treatment_filter_drops_ever_treated(self): + panel = _fixed_panel() # units 0,1,2 are real-treated + res = placebo_group_test( + panel, + outcome="y", + time="t", + unit="unit", + fake_treated_units=[3, 4], + treatment="treatment", + ) + # 5 never-treated units x 2 periods = 10 rows (3 real-treated dropped) + assert res.n_obs == 10 + # fake_group records the units actually used, not just the requested list + assert res.fake_group == [3, 4] + + def test_dispatcher_threads_treatment_into_fake_group(self): + """run_placebo_test(test_type='fake_group') filters ever-treated by default.""" + panel = _fixed_panel() # units 0,1,2 are real-treated + res = run_placebo_test( + panel, + outcome="y", + treatment="treatment", + time="t", + unit="unit", + test_type="fake_group", + fake_treatment_group=[3, 4], + ) + # dispatcher passes treatment -> ever-treated dropped -> 5 controls x 2 = 10 rows + assert res.n_obs == 10 + assert res.fake_group == [3, 4] + + def test_backward_compatible_without_treatment(self): + panel = _fixed_panel() + res = placebo_group_test( + panel, outcome="y", time="t", unit="unit", fake_treated_units=[3, 4] + ) + assert res.n_obs == 16 # all units retained (old behavior) + + def test_degenerate_filter_raises_valueerror(self): + panel = _fixed_panel() + # all requested fake-treated units are themselves real-treated -> dropped + # (the misuse UserWarning is intrinsic here and asserted separately below) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + with pytest.raises(ValueError, match="No fake-treated observations remain"): + placebo_group_test( + panel, + outcome="y", + time="t", + unit="unit", + fake_treated_units=[0, 1], + treatment="treatment", + ) + + def test_treatment_filter_rejects_non_binary(self): + """Fail closed: a non-0/1 treatment column must not silently skip filtering.""" + panel = _fixed_panel() + panel["treatment"] = panel["treatment"] * 2 # {0, 2} -> non-binary + with pytest.raises(ValueError, match="binary"): + placebo_group_test( + panel, + outcome="y", + time="t", + unit="unit", + fake_treated_units=[3, 4], + treatment="treatment", + ) + + def test_treatment_filter_rejects_missing_column(self): + panel = _fixed_panel() + with pytest.raises(ValueError, match="not found"): + placebo_group_test( + panel, + outcome="y", + time="t", + unit="unit", + fake_treated_units=[3, 4], + treatment="nonexistent", + ) + + def test_treatment_filter_rejects_missing_values(self): + """Fail closed: NaN treatment values (which groupby().max() silently skips).""" + panel = _fixed_panel() + panel.loc[0, "treatment"] = np.nan + with pytest.raises(ValueError, match="contains missing values"): + placebo_group_test( + panel, + outcome="y", + time="t", + unit="unit", + fake_treated_units=[3, 4], + treatment="treatment", + ) + + def test_misuse_warns_when_fake_unit_is_real_treated(self): + panel = _fixed_panel() + with pytest.warns(UserWarning, match="real-treated"): + placebo_group_test( + panel, + outcome="y", + time="t", + unit="unit", + fake_treated_units=[2, 3, 4], + treatment="treatment", + ) + + +# --------------------------------------------------------------------------- # +# 5. Inference-honesty contracts +# --------------------------------------------------------------------------- # +class TestPlaceboInferenceContracts: + """The deliberate permutation NaN-decoupling + fail-closed RuntimeErrors.""" + + def test_permutation_nan_decoupling(self): + """n_valid == 1 -> se is NaN and t_stat NaN, but the RI p-value stays finite. + + BDM fn 12: the RI p-value is count-based and valid even when the + permutation SE is degenerate, intentionally departing from the + bootstrap-NaN contract (non-finite SE -> full NaN tuple). With a single + permutation, ``se = std([x], ddof=1)`` is NaN, so ``t_stat`` is NaN, yet + the count-based p-value remains finite. + """ + panel = _fixed_panel() + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) # ddof=1 on a 1-element array + res = permutation_test( + panel, + outcome="y", + treatment="treatment", + time="t", + unit="unit", + n_permutations=1, + seed=5, + ) + assert not np.isfinite(res.se) # NaN at n_valid == 1 + assert np.isnan(res.t_stat) + assert np.isfinite(res.p_value) # the contract: p-value survives + + def test_leave_one_out_all_fail_raises(self): + """A single treated unit -> no LOO estimate possible -> fail-closed.""" + rows = [ + {"unit": u, "t": t, "y": float(u) + 0.5 * t + 0.1 * u * t, "treatment": int(u == 0)} + for u in range(5) + for t in range(2) + ] + panel = pd.DataFrame(rows) + with pytest.raises(RuntimeError, match="leave-one-out"): + leave_one_out_test(panel, outcome="y", treatment="treatment", time="t", unit="unit")