Add DiD OVB sensitivity analysis and reproduction benchmarks

README.md

@@ -122,6 +122,7 @@ Full guide: `diff_diff.get_llm_guide("practitioner")`.
 - [ChangesInChanges](https://diff-diff.readthedocs.io/en/stable/api/changes_in_changes.html) - Athey & Imbens (2006) nonlinear/distributional DiD for the 2x2 design: full counterfactual distribution and quantile treatment effects via CDF transformation, plus the QDiD comparison estimator via `method="qdid"`; bootstrap inference; R qte parity. Alias `CiC`
 - [LWDiD](https://diff-diff.readthedocs.io/en/stable/api/lwdid.html) - Lee & Wooldridge (2025, 2026) rolling-transformation DiD: unit-specific demean/detrend converts panel to cross-section, staggered adoption, `estimation_method` in `reg`/`ipw`/`dr`/`psm` (the papers' RA/IPW/IPWRA plus propensity-score matching), exact small-N inference on the classical collapsed regression
 - [DMLDiD](https://diff-diff.readthedocs.io/en/stable/api/dml_did.html) - Chang (2020) double/debiased machine learning DiD: staggered ATT(g,t) with cross-fitted ML nuisance learners (DML2) and Neyman-orthogonal scores, for flexible/high-dimensional covariate adjustment under conditional parallel trends; panel or declared repeated cross sections (`panel=False`); survey/cluster support on both lanes (bad-control lane: panel only, `cluster=` only); Caetano, Callaway, Payne & Sant'Anna (2026) bad-control score via `fit(bad_control=, bad_control_covariates=)`
+- [DIDOVBSensitivity](https://diff-diff.readthedocs.io/en/stable/api/did_ovb.html) - Wang, Sant'Anna, Chernozhukov & Cinelli (2026) canonical two-period DiD omitted-variable-bias sensitivity analysis with cross-fitted nuisance learners, restricted bias bounds, robustness values, and extreme robustness values
 - [DurationDiD](https://diff-diff.readthedocs.io/en/stable/api/duration_did.html) - Deaner & Ku (2026) causal duration DiD for a binary absorbing outcome (spell ended) in a two-group common-timing design: restricts the groups' untreated hazards (`method="cd"` additive gap or `method="ph"` ratio) instead of outcome levels, imputes the treated counterfactual survival, reports the per-date absorption ATT with whole-individual bootstrap pointwise and simultaneous bands plus a fixed-anchor pre-treatment specification test
 - [BaconDecomposition](https://diff-diff.readthedocs.io/en/stable/api/bacon.html) - Goodman-Bacon (2021) decomposition for diagnosing TWFE bias in staggered settings
 

benchmarks/R/benchmark_did_ovb_minwage.R

@@ -0,0 +1,134 @@
+#!/usr/bin/env Rscript
+
+# R application oracle for Wang et al.'s minimum-wage example.
+# This file uses ranger only as an executed reference implementation; no
+# ranger source is copied into diff-diff.
+
+suppressPackageStartupMessages({
+  library(jsonlite)
+  library(ranger)
+})
+
+args <- commandArgs(trailingOnly = TRUE)
+input <- if (length(args) >= 1) args[[1]] else "/private/tmp/CS_RR/data/min_wage_CS.rds"
+output <- if (length(args) >= 2) args[[2]] else "/private/tmp/did_ovb_minwage_r.json"
+fold_output <- if (length(args) >= 3) args[[3]] else "/private/tmp/did_ovb_minwage_folds.csv"
+seed <- if (length(args) >= 4) as.integer(args[[4]]) else 42L
+n_folds <- if (length(args) >= 5) as.integer(args[[5]]) else 5L
+num_trees <- if (length(args) >= 6) as.integer(args[[6]]) else 1000L
+mtry <- if (length(args) >= 7) as.integer(args[[7]]) else 2L
+min_node_size <- if (length(args) >= 8) as.integer(args[[8]]) else 10L
+splitrule <- if (length(args) >= 9) args[[9]] else "variance"
+
+if (grepl("\\.rds$", tolower(input))) {
+  raw <- readRDS(input)
+} else {
+  raw <- read.csv(input, stringsAsFactors = TRUE, check.names = FALSE)
+}
+
+dat <- raw[raw$year %in% c(2006, 2007) & raw$first.treat %in% c(0, 2007), , drop = FALSE]
+dat$treated <- as.integer(dat$first.treat == 2007)
+dat <- dat[order(dat$countyreal, dat$year), , drop = FALSE]
+
+units <- unique(dat$countyreal)
+if (length(units) != 1961L || sum(dat$treated[dat$year == 2006]) != 584L) {
+  stop("unexpected minimum-wage application sample")
+}
+
+wide <- dat[dat$year == 2006, c("countyreal", "treated", "region", "white", "hs", "pov", "lpop", "lmedinc"), drop = FALSE]
+post <- dat[dat$year == 2007, c("countyreal", "lemp"), drop = FALSE]
+pre <- dat[dat$year == 2006, c("countyreal", "lemp"), drop = FALSE]
+wide <- merge(wide, pre, by = "countyreal", suffixes = c("", "_pre"), sort = FALSE)
+wide <- merge(wide, post, by = "countyreal", suffixes = c("", "_post"), sort = FALSE)
+wide$delta_y <- wide$lemp_post - wide$lemp
+
+set.seed(seed)
+folds <- integer(nrow(wide))
+for (d in c(0L, 1L)) {
+  idx <- which(wide$treated == d)
+  folds[idx] <- sample(rep(0:(n_folds - 1L), length.out = length(idx)))
+}
+write.csv(data.frame(countyreal = wide$countyreal, fold_id = folds), fold_output, row.names = FALSE)
+
+x_names <- c("region", "white", "hs", "pov", "lpop", "lmedinc")
+x <- wide[, x_names, drop = FALSE]
+# The Python public API currently accepts numeric covariates. Encode region by
+# its public integer levels so both application runners consume identical X.
+x$region <- as.numeric(x$region)
+x <- as.matrix(x)
+d <- wide$treated
+y <- wide$delta_y
+ps <- numeric(nrow(wide))
+m <- numeric(nrow(wide))
+
+for (fold in 0:(n_folds - 1L)) {
+  test <- folds == fold
+  train <- !test
+  ps_fit <- ranger(
+    x = as.data.frame(x[train, , drop = FALSE]),
+    y = d[train],
+    num.trees = num_trees,
+    mtry = mtry,
+    min.node.size = min_node_size,
+    splitrule = splitrule,
+    seed = seed + fold
+  )
+  out_fit <- ranger(
+    x = as.data.frame(x[train & d == 0, , drop = FALSE]),
+    y = y[train & d == 0],
+    num.trees = num_trees,
+    mtry = mtry,
+    min.node.size = min_node_size,
+    splitrule = splitrule,
+    seed = seed + fold
+  )
+  ps[test] <- predict(ps_fit, data = as.data.frame(x[test, , drop = FALSE]))$predictions
+  m[test] <- predict(out_fit, data = as.data.frame(x[test, , drop = FALSE]))$predictions
+}
+
+ps <- pmin(pmax(ps, 0.01), 0.99)
+p <- mean(d)
+omega <- (ps / (1 - ps)) / (p / (1 - p))
+residual <- y - m
+score <- (d / p - (1 - d) / (1 - p) * omega) * residual
+short_att <- mean(score)
+sigma2 <- mean(residual[d == 0]^2)
+nu2 <- mean(omega[d == 0]^2)
+scale <- sqrt(sigma2 * nu2)
+theta_if <- score - short_att
+short_se <- sqrt(mean(theta_if^2) / length(y))
+
+# This is the same plug-in/IF root search used by the Python result object.
+sigma_if <- (1 - d) / (1 - p) * (residual^2 - sigma2)
+nu_if <- (1 - d) / (1 - p) * (omega^2 - nu2)
+scale_if <- (nu2 * sigma_if + sigma2 * nu_if) / (2 * scale)
+contains <- function(strength, kind, alpha = 0.05) {
+  multiplier <- if (kind == "rv") strength / sqrt(1 - strength) else sqrt(strength / (1 - strength))
+  radius <- multiplier * scale
+  lo_if <- theta_if - multiplier * scale_if
+  hi_if <- theta_if + multiplier * scale_if
+  z <- qnorm(1 - alpha / 2)
+  lo_se <- sqrt(mean(lo_if^2) / length(y))
+  hi_se <- sqrt(mean(hi_if^2) / length(y))
+  short_att - radius - z * lo_se <= 0 && 0 <= short_att + radius + z * hi_se
+}
+find_rv <- function(kind) {
+  if (!contains(1 - 1e-12, kind)) return(NaN)
+  lo <- 0
+  hi <- 1 - 1e-12
+  for (i in seq_len(70)) {
+    mid <- (lo + hi) / 2
+    if (contains(mid, kind)) hi <- mid else lo <- mid
+  }
+  hi
+}
+
+result <- list(
+  n = nrow(wide), treated = sum(d), control = sum(d == 0), n_folds = n_folds,
+  seed = seed, num_trees = num_trees, mtry = mtry,
+  min_node_size = min_node_size, splitrule = splitrule,
+  short_att = short_att, short_se = short_se, sigma2_control = sigma2,
+  nu2_selection = nu2, scale = scale, rv = find_rv("rv"), xrv = find_rv("xrv")
+)
+write_json(result, output, auto_unbox = TRUE, digits = 17, pretty = TRUE)
+cat(toJSON(result, auto_unbox = TRUE, digits = 8, pretty = TRUE), "\n")

benchmarks/R/generate_did_ovb_parity.R

@@ -0,0 +1,129 @@
+#!/usr/bin/env Rscript
+
+# Independent parity oracle for the canonical DiD OVB implementation.
+# The calculations below are written from the paper's displayed formulas;
+# this script does not import or call dml.sensemakr.
+
+suppressPackageStartupMessages(library(jsonlite))
+
+args <- commandArgs(trailingOnly = TRUE)
+output_path <- if (length(args) >= 1) args[[1]] else "benchmarks/data/did_ovb_r_results.json"
+panel_path <- if (length(args) >= 2) args[[2]] else "benchmarks/data/real/mpdta.csv"
+
+zcrit <- qnorm(0.975)
+data <- read.csv(panel_path, check.names = FALSE)
+data <- data[data$year %in% c(2006, 2007) & data$`first.treat` %in% c(0, 2007), ]
+data <- data[order(data$countyreal, data$year), ]
+
+pre <- data[data$year == 2006, ]
+post <- data[data$year == 2007, ]
+stopifnot(nrow(pre) == nrow(post), all(pre$countyreal == post$countyreal))
+
+d <- as.numeric(pre$`first.treat` == 2007)
+dy <- post$lemp - pre$lemp
+x <- pre$lpop
+n <- length(d)
+n_folds <- 2L
+fold_ids <- (seq_len(n) - 1L) %% n_folds
+p <- mean(d)
+
+ps <- numeric(n)
+m <- numeric(n)
+for (fold in 0:(n_folds - 1L)) {
+  test <- fold_ids == fold
+  train <- !test
+  ps_fit <- glm(d ~ x, family = binomial(), subset = train)
+  m_fit <- lm(dy ~ x, subset = train & d == 0)
+  ps[test] <- predict(ps_fit, newdata = data.frame(x = x[test]), type = "response")
+  m[test] <- predict(m_fit, newdata = data.frame(x = x[test]))
+}
+
+ps <- pmin(pmax(ps, 0.01), 0.99)
+omega <- (ps / (1 - ps)) / (p / (1 - p))
+residual <- dy - m
+score <- (d / p - (1 - d) / (1 - p) * omega) * residual
+short_att <- mean(score)
+sigma2 <- mean(residual[d == 0]^2)
+nu2 <- mean(omega[d == 0]^2)
+scale <- sqrt(sigma2 * nu2)
+theta_if <- score - short_att
+sigma_if <- (1 - d) / (1 - p) * (residual^2 - sigma2)
+nu_if <- (1 - d) / (1 - p) * (omega^2 - nu2)
+scale_if <- (nu2 * sigma_if + sigma2 * nu_if) / (2 * scale)
+short_se <- sqrt(mean(theta_if^2) / n)
+
+bounds <- function(trend_r2, selection_r2, rho_max = 1, alpha = 0.05) {
+  multiplier <- rho_max * sqrt(trend_r2) * sqrt(selection_r2 / (1 - selection_r2))
+  radius <- multiplier * scale
+  lower_if <- theta_if - multiplier * scale_if
+  upper_if <- theta_if + multiplier * scale_if
+  lower_se <- sqrt(mean(lower_if^2) / n)
+  upper_se <- sqrt(mean(upper_if^2) / n)
+  z <- qnorm(1 - alpha / 2)
+  list(
+    lower = short_att - radius,
+    upper = short_att + radius,
+    radius = radius,
+    lower_se = lower_se,
+    upper_se = upper_se,
+    lower_ci = short_att - radius - z * lower_se,
+    upper_ci = short_att + radius + z * upper_se,
+    trend_r2 = trend_r2,
+    selection_r2 = selection_r2,
+    rho_max = rho_max,
+    alpha = alpha
+  )
+}
+
+robustness <- function(null_value = 0, alpha = 0.05) {
+  base <- bounds(0, 0, alpha = alpha)
+  if (base$lower_ci <= null_value && null_value <= base$upper_ci) {
+    return(list(rv = 0, xrv = 0, rv_bounds = base, xrv_bounds = base))
+  }
+  search <- function(kind) {
+    contains <- function(s) {
+      b <- if (kind == "rv") bounds(s, s, alpha = alpha) else bounds(1, s, alpha = alpha)
+      b$lower_ci <= null_value && null_value <= b$upper_ci
+    }
+    lo <- 0
+    hi <- 1 - 1e-12
+    if (!contains(hi)) {
+      b <- if (kind == "rv") bounds(hi, hi, alpha = alpha) else bounds(1, hi, alpha = alpha)
+      return(list(value = NaN, bounds = b))
+    }
+    for (i in seq_len(80)) {
+      mid <- (lo + hi) / 2
+      if (contains(mid)) hi <- mid else lo <- mid
+    }
+    b <- if (kind == "rv") bounds(hi, hi, alpha = alpha) else bounds(1, hi, alpha = alpha)
+    list(value = hi, bounds = b)
+  }
+  rv <- search("rv")
+  xrv <- search("xrv")
+  list(rv = rv$value, xrv = xrv$value, rv_bounds = rv$bounds, xrv_bounds = xrv$bounds)
+}
+
+result <- list(
+  settings = list(
+    panel = panel_path,
+    periods = c(2006, 2007),
+    treatment_cohort = 2007,
+    n_folds = n_folds,
+    fold_ids = fold_ids,
+    pscore_trim = 0.01,
+    alpha = 0.05,
+    null_value = -0.1
+  ),
+  n_obs = n,
+  n_treated = sum(d == 1),
+  n_control = sum(d == 0),
+  short_att = short_att,
+  short_se = short_se,
+  sigma2_control = sigma2,
+  nu2_selection = nu2,
+  scale = scale,
+  bounds = bounds(1, 0.5),
+  robustness = robustness(-0.1, 0.05)
+)
+
+write_json(result, output_path, auto_unbox = TRUE, digits = 17, pretty = TRUE)

benchmarks/R/generate_did_ovb_simulation.R

@@ -0,0 +1,58 @@
+#!/usr/bin/env Rscript
+
+# Generate one draw from Appendix E.1 of Wang et al. and estimate the
+# short DiD OVB components using the correctly specified parametric nuisances.
+# The CSV is intentionally shared with the Python runner for cross-language
+# comparison; it contains the latent U only because this is a simulation.
+
+suppressPackageStartupMessages(library(jsonlite))
+
+args <- commandArgs(trailingOnly = TRUE)
+csv_path <- if (length(args) >= 1) args[[1]] else "benchmarks/data/did_ovb_simulation.csv"
+json_path <- if (length(args) >= 2) args[[2]] else "benchmarks/data/did_ovb_simulation_r.json"
+n <- if (length(args) >= 3) as.integer(args[[3]]) else 500L
+p <- if (length(args) >= 4) as.numeric(args[[4]]) else 0.5
+seed <- if (length(args) >= 5) as.integer(args[[5]]) else 20260928L
+
+set.seed(seed)
+d <- rbinom(n, 1L, p)
+x <- rnorm(n, mean = ifelse(d == 0, 0.3, 0), sd = ifelse(d == 0, sqrt(6), sqrt(3)))
+u <- rnorm(n, mean = ifelse(d == 0, 0.3, 0), sd = ifelse(d == 0, sqrt(6), sqrt(3)))
+delta_y <- 1 + x + u + 2 * d + rnorm(n, sd = sqrt(2))
+
+fold_ids <- (seq_len(n) - 1L) %% 10L
+ps <- numeric(n)
+m <- numeric(n)
+for (fold in 0:9) {
+  test <- fold_ids == fold
+  train <- !test
+  ps_fit <- glm(d ~ x + I(x^2), family = binomial(), subset = train)
+  m_fit <- lm(delta_y ~ x, subset = train & d == 0)
+  ps[test] <- predict(ps_fit, newdata = data.frame(x = x[test]), type = "response")
+  m[test] <- predict(m_fit, newdata = data.frame(x = x[test]))
+}
+
+ps <- pmin(pmax(ps, 0.01), 0.99)
+p_hat <- mean(d)
+omega <- (ps / (1 - ps)) / (p_hat / (1 - p_hat))
+residual <- delta_y - m
+score <- (d / p_hat - (1 - d) / (1 - p_hat) * omega) * residual
+short_att <- mean(score)
+sigma2 <- mean(residual[d == 0]^2)
+nu2 <- mean(omega[d == 0]^2)
+scale <- sqrt(sigma2 * nu2)
+
+write.csv(
+  data.frame(unit = seq_len(n), pre = 0, post = delta_y, treated = d, x = x, u = u),
+  csv_path,
+  row.names = FALSE
+)
+write_json(
+  list(
+    n = n, p = p, seed = seed, n_folds = 10L,
+    short_att = short_att, sigma2_control = sigma2,
+    nu2_selection = nu2, scale = scale,
+    treatment_effect = 2, omitted_bias = -mean(u[d == 1]) + mean(u[d == 0])
+  ),
+  json_path, auto_unbox = TRUE, digits = 17, pretty = TRUE
+)