feat: fuzzy regression discontinuity via treatment_col on RegressionDiscontinuity

CHANGELOG.md

@@ -8,31 +8,49 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
-- **`RegressionDiscontinuity` - sharp regression discontinuity estimation with robust
-  bias-corrected inference (alias `RDD`).** Local-polynomial sharp RD per Calonico,
-  Cattaneo & Titiunik (2014), parity-targeting R `rdrobust` 4.0.0 end-to-end: all 10
-  data-driven bandwidth selectors (`mserd` default, `msetwo`/`msesum`/comb and the
-  CER-optimal variants), triangular/epanechnikov/uniform kernels, `masspoints`
-  adjust/check/off, manual `h`/`b`/`rho` with R-exact resolution semantics (including
-  `rho`-without-`h` applying to selected bandwidths and the unconditional N<20
-  full-range fallback), and the three-row Conventional / Bias-Corrected / Robust
-  output. **Canonical binding:** `att`/`se`/`t_stat`/`p_value`/`conf_int` are ONE
+- **`RegressionDiscontinuity` - sharp AND fuzzy regression discontinuity estimation
+  with robust bias-corrected inference (alias `RDD`).** Local-polynomial RD per
+  Calonico, Cattaneo & Titiunik (2014), parity-targeting R `rdrobust` 4.0.0
+  end-to-end: all 10 data-driven bandwidth selectors (`mserd` default,
+  `msetwo`/`msesum`/comb and the CER-optimal variants),
+  triangular/epanechnikov/uniform kernels, `masspoints` adjust/check/off, manual
+  `h`/`b`/`rho` with R-exact resolution semantics (including `rho`-without-`h`
+  applying to selected bandwidths and the unconditional N<20 full-range fallback),
+  and the three-row Conventional / Bias-Corrected / Robust output. **Fuzzy designs**
+  via `fit(..., treatment_col=...)` (R's `fuzzy=`): the estimand is the local Wald
+  ratio (complier LATE at the cutoff for binary take-up under monotonicity;
+  ratio-of-jumps otherwise, with a data-dependent `estimand` label) with the LINEARIZED
+  bias correction and delta-method variance of rdrobust (T stacked as a second
+  response column; Y-T covariance via the `res @ s_Y` collapse); bandwidths select
+  on the fuzzy-ratio objective by default with R's exact `sharpbw`/one-sided
+  perfect-compliance auto-switch to the sharp reduced form; the first stage
+  (take-up jump) is a full three-row `first_stage*` mirror plus a first-stage
+  `summary()` block; a weak-first-stage `UserWarning` fires when the first-stage
+  robust CI contains zero (documented deviation - R is silent; CCT 2014 Theorem 3
+  "guard and warn", Feir-Lemieux-Marmer weak-IV inference a documented seam); R's
+  exact no-variation-no-jump identification error is raised on both entry points.
+  **Canonical binding:** `att`/`se`/`t_stat`/`p_value`/`conf_int` are ONE
   coherent row - the robust bias-corrected row (`att = tau_bc`, CI centered on it,
-  `t_stat == att/se`), preserving the library-wide field identities; rdrobust's
+  `t_stat == att/se`), preserving the library-wide field identities; the new
+  `estimand` results field names what `att` measures per fit; rdrobust's
   printed headline coefficient is exposed as `att_conventional` with a full inference
   row, and `summary()` prints the familiar three-row table. Estimation-path port
-  (`rdrobust_fit_sharp`: Q_q bias-correction score matrix, conventional/robust NN
-  sandwiches) validated against a new estimates golden
-  (`benchmarks/data/rdrobust_estimates_golden.json`, 16 configurations incl. the
-  Senate anchors) at rtol=1e-9 in `tests/test_rdd_parity.py`; R-free methodology
+  (`rdrobust_fit`: Q_q bias-correction score matrix, conventional/robust NN
+  sandwiches, fuzzy ratio/first-stage variances) validated against a new estimates
+  golden (`benchmarks/data/rdrobust_estimates_golden.json`, 23 configurations incl.
+  the Senate anchors and 7 fuzzy configs - default/sharpbw/manual-h/epa/msetwo/
+  one-sided-compliance/ties) at rtol=1e-9 in `tests/test_rdd_parity.py` (+ port-level
+  linearized-bias pins in `tests/test_rdrobust_port.py`); R-free methodology
   anchors (CCT 2014 Remark 7 bias-corrected == local-quadratic equivalence at rel
-  1e-10 across all kernels, invariances, joint-NaN degenerate contracts) in
+  1e-10 across all kernels, perfect-compliance == sharp reproduction, bandwidth
+  auto-switch locks, invariances, joint-NaN degenerate contracts) in
   `tests/test_rdd_methodology.py`; API/validation suite in `tests/test_rdd.py`.
   Deviations from R (each labeled in the REGISTRY section): warn-instead-of-silent
-  NaN drops, warn-and-ignore `b`-without-`h`, fail-closed targeted errors on
-  degenerate designs, and the canonical-binding note above. Sharp designs only:
-  fuzzy RD, covariates (CCFT 2019 - review on file), cluster-robust variance,
-  weights, kink estimands, and rdplot/density diagnostics are documented follow-ups.
+  NaN drops, warn-and-ignore `b`-without-`h`, the weak-first-stage warning,
+  warn-and-ignore `sharpbw` on sharp fits, fail-closed targeted errors on
+  degenerate designs, and the canonical-binding note above. Covariates (CCFT 2019 -
+  review on file), cluster-robust variance, weights, kink estimands, weak-IV-robust
+  fuzzy inference, and rdplot/density diagnostics are documented follow-ups.
 - **Internal: mypy enforced at zero errors.** Triaged the 184 pre-existing
   `mypy diff_diff` errors to an enforceable zero and added a blocking Mypy job to
   the Lint CI workflow (pinned `mypy==2.1.0` + pinned numpy/pandas/scipy for stub
@@ -91,8 +109,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 - `diff_diff/guides/llms-autonomous.txt` no longer lists regression discontinuity as
-  out of scope: sharp RD routes to `RegressionDiscontinuity`; only fuzzy RD is
-  referred to external tooling.
+  out of scope: sharp AND fuzzy RD route to `RegressionDiscontinuity`; only kink
+  designs and the covariate-adjusted / cluster-robust RD variants are referred to
+  external tooling.
 - **Internal: repo-wide lint normalization + pinned tooling.** black/ruff/mypy are now
   pinned exactly in the `dev` extra (`black==26.3.1`, `ruff==0.15.13`, `mypy==2.1.0`;
   the tools require Python >= 3.10 — the library floor stays 3.9); full `black` +

README.md

@@ -112,7 +112,7 @@ Full guide: `diff_diff.get_llm_guide("practitioner")`.
 - [TripleDifference](https://diff-diff.readthedocs.io/en/stable/api/triple_diff.html) - triple difference (DDD) estimator for designs requiring two criteria for treatment eligibility
 - [ContinuousDiD](https://diff-diff.readthedocs.io/en/stable/api/continuous_did.html) - Callaway, Goodman-Bacon & Sant'Anna (2024) continuous treatment DiD with dose-response curves
 - [HeterogeneousAdoptionDiD](https://diff-diff.readthedocs.io/en/stable/api/had.html) - de Chaisemartin, Ciccia, D'Haultfœuille & Knau (2026) for designs where **no unit remains untreated**; local-linear estimator at the dose support boundary returning Weighted Average Slope (WAS) on Design 1' (`d̲ = 0` / QUG) or `WAS_{d̲}` on Design 1 (`d̲ > 0`, continuous-near-d̲ or mass-point), with a multi-period event-study extension (last-treatment cohort, pointwise CIs). **Panel-only** in this release - repeated cross-sections rejected by the validator. Alias `HAD`.
-- [RegressionDiscontinuity](https://diff-diff.readthedocs.io/en/stable/api/regression_discontinuity.html) - Calonico, Cattaneo & Titiunik (2014) sharp regression discontinuity with robust bias-corrected inference and rdrobust-parity bandwidth selection (all 10 selectors, mass-point handling). **Sharp designs only** in this release; canonical `att` is the bias-corrected estimate with a coherent robust CI (rdrobust's printed headline is `att_conventional`). Alias `RDD`.
+- [RegressionDiscontinuity](https://diff-diff.readthedocs.io/en/stable/api/regression_discontinuity.html) - Calonico, Cattaneo & Titiunik (2014) sharp AND fuzzy regression discontinuity with robust bias-corrected inference and rdrobust-parity bandwidth selection (all 10 selectors, mass-point handling; fuzzy via `treatment_col=` with a first-stage block and weak-identification warning). Canonical `att` is the bias-corrected estimate with a coherent robust CI (rdrobust's printed headline is `att_conventional`). Alias `RDD`.
 - [StackedDiD](https://diff-diff.readthedocs.io/en/stable/api/stacked_did.html) - Wing, Freedman & Hollingsworth (2024) stacked DiD with Q-weights and sub-experiments; optional covariate balancing (Ustyuzhanin 2026)
 - [EfficientDiD](https://diff-diff.readthedocs.io/en/stable/api/efficient_did.html) - Chen, Sant'Anna & Xie (2025) efficient DiD with optimal weighting for tighter SEs
 - [TROP](https://diff-diff.readthedocs.io/en/stable/api/trop.html) - Triply Robust Panel estimator (Athey et al. 2025) with nuclear norm factor adjustment

benchmarks/R/README.md

@@ -137,9 +137,14 @@ flatten misaligns when reshaped row-major.
 `benchmarks/R/generate_rdrobust_golden.R` produces
 `benchmarks/data/rdrobust_golden.json`, consumed by
 `tests/test_rdrobust_port.py` to verify that
-`diff_diff._rdrobust_port.rdbwselect_sharp` matches R `rdrobust::rdbwselect`
-(Calonico, Cattaneo, Farrell & Titiunik) across all 10 bandwidth selectors at
-rtol ≤ 1e-9.
+`diff_diff._rdrobust_port.rdbwselect` matches R `rdrobust::rdbwselect`
+(Calonico, Cattaneo, Farrell & Titiunik) on SHARP bandwidth selection
+across all 10 selectors at rtol ≤ 1e-9 (17 configs; this fixture predates
+fuzzy support and is deliberately never regenerated). Fuzzy bandwidth and
+estimation parity lives in `benchmarks/data/rdrobust_estimates_golden.json`
+(generator `generate_rdrobust_estimates_golden.R`), pinned by
+`tests/test_rdrobust_port.py::TestFuzzyPortGoldenParity` and
+`tests/test_rdd_parity.py`.
 
 ## Version pin
 

benchmarks/R/generate_rdrobust_estimates_golden.R

@@ -1,6 +1,8 @@
-# Golden-value generator for the diff-diff sharp-RD ESTIMATION port
-# (diff_diff/_rdrobust_port.py::rdrobust_fit_sharp and the public
-# RegressionDiscontinuity estimator).
+# Golden-value generator for the diff-diff RD ESTIMATION port - sharp AND
+# fuzzy (diff_diff/_rdrobust_port.py::rdrobust_fit and the public
+# RegressionDiscontinuity estimator; 23 configs across four synthetic DGPs
+# + the Senate data, incl. 7 fuzzy configs with full first-stage
+# tau_T/se_T/z_T/pv_T/ci_T blocks and per-side take-up coefficients).
 #
 # Deliberately a SEPARATE file/JSON from generate_rdrobust_golden.R so the
 # bandwidth fixtures reviewed in the machinery PR are never regenerated.
@@ -24,16 +26,20 @@ TARBALL_SHA256 <- "78f0d6b4bdec4091cc8f42f6f1598704747f95926446d3aaee381ea1d613a
 
 run_estimate <- function(y, x, c = 0, masspoints = "adjust", kernel = "tri",
                          p = 1, q = 2, h = NULL, b = NULL, rho = NULL,
-                         level = 95, bwselect = "mserd") {
+                         level = 95, bwselect = "mserd",
+                         fuzzy = NULL, sharpbw = FALSE) {
   args <- list(y = y, x = x, c = c, masspoints = masspoints, kernel = kernel,
-               p = p, q = q, level = level, bwselect = bwselect)
+               p = p, q = q, level = level, bwselect = bwselect,
+               sharpbw = sharpbw)
   if (!is.null(h)) args$h <- h
   if (!is.null(b)) args$b <- b
   if (!is.null(rho)) args$rho <- rho
+  if (!is.null(fuzzy)) args$fuzzy <- fuzzy
   r <- suppressWarnings(do.call(rdrobust, args))
-  list(
+  out <- list(
     c = c, masspoints = masspoints, kernel = kernel, p = p, q = q,
     bwselect = bwselect,
+    fuzzy_in = !is.null(fuzzy), sharpbw = sharpbw,
     h_in = if (is.null(h)) NA else h,
     b_in = if (is.null(b)) NA else b,
     rho_in = if (is.null(rho)) NA else rho,
@@ -49,6 +55,18 @@ run_estimate <- function(y, x, c = 0, masspoints = "adjust", kernel = "tri",
     beta_p_l = unname(as.vector(r$beta_Y_p_l)),
     beta_p_r = unname(as.vector(r$beta_Y_p_r))
   )
+  if (!is.null(fuzzy)) {
+    # First-stage three-row block (Conventional / Bias-Corrected / Robust)
+    out$tau_T <- unname(r$tau_T)
+    out$se_T <- unname(r$se_T)
+    out$z_T <- unname(r$z_T)
+    out$pv_T <- unname(r$pv_T)
+    out$ci_T_lower <- unname(r$ci_T[, 1])
+    out$ci_T_upper <- unname(r$ci_T[, 2])
+    out$beta_t_p_l <- unname(as.vector(r$beta_T_p_l))
+    out$beta_t_p_r <- unname(as.vector(r$beta_T_p_r))
+  }
+  out
 }
 
 golden <- list()
@@ -57,12 +75,14 @@ golden$metadata <- list(
   rdrobust_version = as.character(packageVersion("rdrobust")),
   rdrobust_tarball_sha256 = TARBALL_SHA256,
   seeds = list(dgp_lee_smooth = 42L, dgp_ties_moderate = 123L,
-               dgp_asymmetric_scaled = 777L),
+               dgp_asymmetric_scaled = 777L, dgp_fuzzy = 314L),
   generator = "benchmarks/R/generate_rdrobust_estimates_golden.R",
   algorithm = paste(
-    "rdrobust() sharp-RD estimation blocks (three-row coef/se/z/pv/ci,",
-    "counts, per-side beta_p) for the vce='nn' no-covariate path,",
-    "complementing the bandwidth fixtures in rdrobust_golden.json."
+    "rdrobust() sharp AND fuzzy estimation blocks (three-row coef/se/z/pv/ci,",
+    "counts, per-side beta_p; fuzzy configs add the first-stage",
+    "tau_T/se_T/z_T/pv_T/ci_T rows and per-side beta_T_p) for the vce='nn'",
+    "no-covariate path, complementing the bandwidth fixtures in",
+    "rdrobust_golden.json."
   ),
   r_version = R.version.string
 )
@@ -116,6 +136,32 @@ golden$dgp_asymmetric_scaled <- list(
   )
 )
 
+# Fuzzy DGP: two-sided imperfect compliance (take-up jumps 0.15 -> 0.75).
+set.seed(314)
+n4 <- 1500
+x4 <- 2 * rbeta(n4, 2, 4) - 1
+t4 <- rbinom(n4, 1, ifelse(x4 >= 0, 0.75, 0.15))
+y4 <- 0.5 * x4 + 1.2 * t4 + rnorm(n4, sd = 0.3)
+# One-sided perfect compliance variant (T == 0 left of the cutoff):
+# exercises the perf_comp bandwidth auto-switch (rdbwselect.R:334-346).
+t4_one <- ifelse(x4 >= 0, t4, 0)
+# Tied running variable variant (2dp rounding; keeps all masspoints modes
+# runnable in R per the sharp-golden lesson).
+x4_ties <- round(x4, 2)
+
+golden$dgp_fuzzy <- list(
+  x = x4, y = y4, t = t4, t_one = t4_one, x_ties = x4_ties,
+  configs = list(
+    default      = run_estimate(y4, x4, fuzzy = t4),
+    sharpbw_true = run_estimate(y4, x4, fuzzy = t4, sharpbw = TRUE),
+    manual_h     = run_estimate(y4, x4, fuzzy = t4, h = 0.2),
+    epa          = run_estimate(y4, x4, fuzzy = t4, kernel = "epa"),
+    msetwo       = run_estimate(y4, x4, fuzzy = t4, bwselect = "msetwo"),
+    one_sided    = run_estimate(y4, x4, fuzzy = t4_one),
+    ties_adjust  = run_estimate(y4, x4_ties, fuzzy = t4)
+  )
+)
+
 senate_path <- "benchmarks/data/rdrobust_senate.csv"
 stopifnot(file.exists(senate_path))
 senate <- read.csv(senate_path)