|
| 1 | +# Golden-value generator for the diff-diff RD PLOT port |
| 2 | +# (diff_diff/rdplot.py::RDPlot). Parity target: CRAN rdrobust 4.0.0 |
| 3 | +# (tarball sha256 |
| 4 | +# 78f0d6b4bdec4091cc8f42f6f1598704747f95926446d3aaee381ea1d613a36f); |
| 5 | +# the GitHub 4.1.0-dev tree must NOT be installed when regenerating. |
| 6 | +# |
| 7 | +# Outputs benchmarks/data/rdplot_golden.json: per config the input |
| 8 | +# arguments, the selector outputs (J / J_IMSE / J_MV / rscale / bin_avg / |
| 9 | +# bin_med), counts, descriptive strings, per-side global-fit coefficients, |
| 10 | +# the FULL vars_bins table, and vars_poly at a fixed downsample of the |
| 11 | +# 1000 curve points. DGP inputs are embedded at 17 significant digits. |
| 12 | +# |
| 13 | +# Configs stay inside the surface where R 4.0.0's rdplot() actually runs: |
| 14 | +# fractional scale*J products and length-2 scale CRASH R (vector-indexing |
| 15 | +# accident / vectorized-if error on R >= 4.2) and are covered by |
| 16 | +# Deviation-locked Python unit tests instead of goldens. |
| 17 | +# |
| 18 | +# Run from the repo root: Rscript benchmarks/R/generate_rdplot_golden.R |
| 19 | + |
| 20 | +library(rdrobust) |
| 21 | +library(jsonlite) |
| 22 | + |
| 23 | +stopifnot(packageVersion("rdrobust") == "4.0.0") |
| 24 | + |
| 25 | +TARBALL_SHA256 <- "78f0d6b4bdec4091cc8f42f6f1598704747f95926446d3aaee381ea1d613a36f" |
| 26 | + |
| 27 | +POLY_IDX <- sort(unique(c(seq(1, 1000, by = 25), 500L, 501L, 1000L))) |
| 28 | + |
| 29 | +run_rdplot <- function(name, y, x, c = 0, p = 4, nbins = NULL, |
| 30 | + binselect = "esmv", scale = NULL, kernel = "uni", |
| 31 | + h = NULL, support = NULL, masspoints = "adjust", |
| 32 | + ci = NULL, covs = NULL, covs_drop = TRUE) { |
| 33 | + args <- list(y = y, x = x, c = c, p = p, binselect = binselect, |
| 34 | + kernel = kernel, masspoints = masspoints, hide = TRUE, |
| 35 | + covs_drop = covs_drop) |
| 36 | + if (!is.null(nbins)) args$nbins <- nbins |
| 37 | + if (!is.null(scale)) args$scale <- scale |
| 38 | + if (!is.null(h)) args$h <- h |
| 39 | + if (!is.null(support)) args$support <- support |
| 40 | + if (!is.null(ci)) args$ci <- ci |
| 41 | + if (!is.null(covs)) args$covs <- covs |
| 42 | + r <- tryCatch( |
| 43 | + suppressWarnings(suppressMessages(do.call(rdplot, args))), |
| 44 | + error = function(e) stop(sprintf("config '%s' failed in R: %s", |
| 45 | + name, conditionMessage(e))) |
| 46 | + ) |
| 47 | + vb <- r$vars_bins |
| 48 | + vp <- r$vars_poly |
| 49 | + out <- list( |
| 50 | + # input echoes (the Python test rebuilds RDPlot(...) from these) |
| 51 | + args = list( |
| 52 | + c = c, p = p, |
| 53 | + nbins = if (is.null(nbins)) NA else nbins, |
| 54 | + binselect = binselect, |
| 55 | + scale = if (is.null(scale)) NA else scale, |
| 56 | + kernel = kernel, |
| 57 | + h = if (is.null(h)) NA else h, |
| 58 | + support = if (is.null(support)) NA else support, |
| 59 | + masspoints = masspoints, |
| 60 | + ci = if (is.null(ci)) NA else ci, |
| 61 | + covs_in = !is.null(covs), |
| 62 | + covs_names = if (is.null(covs)) NA else colnames(covs), |
| 63 | + covs_drop = covs_drop |
| 64 | + ), |
| 65 | + J = unname(r$J), J_IMSE = unname(r$J_IMSE), J_MV = unname(r$J_MV), |
| 66 | + scale_out = unname(r$scale), rscale = unname(r$rscale), |
| 67 | + bin_avg = unname(r$bin_avg), bin_med = unname(r$bin_med), |
| 68 | + h_out = unname(r$h), N = unname(r$N), N_h = unname(r$N_h), |
| 69 | + binselect_type = r$binselect, kernel_type = r$kernel, |
| 70 | + coef_left = unname(r$coef[, 1]), coef_right = unname(r$coef[, 2]), |
| 71 | + vars_bins = list( |
| 72 | + mean_bin = vb$rdplot_mean_bin, mean_x = vb$rdplot_mean_x, |
| 73 | + mean_y = vb$rdplot_mean_y, min_bin = vb$rdplot_min_bin, |
| 74 | + max_bin = vb$rdplot_max_bin, se_y = vb$rdplot_se_y, |
| 75 | + N = vb$rdplot_N, ci_l = vb$rdplot_ci_l, ci_r = vb$rdplot_ci_r |
| 76 | + ), |
| 77 | + vars_poly_idx = POLY_IDX, |
| 78 | + vars_poly_x = vp$rdplot_x[POLY_IDX], |
| 79 | + vars_poly_y = vp$rdplot_y[POLY_IDX] |
| 80 | + ) |
| 81 | + if (!is.null(covs)) out$coef_covs <- unname(as.vector(r$coef_covs)) |
| 82 | + out |
| 83 | +} |
| 84 | + |
| 85 | +golden <- list() |
| 86 | + |
| 87 | +golden$metadata <- list( |
| 88 | + rdrobust_version = as.character(packageVersion("rdrobust")), |
| 89 | + rdrobust_tarball_sha256 = TARBALL_SHA256, |
| 90 | + seeds = list(dgp_ties = 123L, dgp_covs = 2718L, dgp_small = 555L, |
| 91 | + dgp_ladder = 999L), |
| 92 | + generator = "benchmarks/R/generate_rdplot_golden.R", |
| 93 | + algorithm = paste( |
| 94 | + "rdplot() bin selectors (all 8 binselect values), manual", |
| 95 | + "nbins/scale/h/support/ci knobs, masspoints adjust remap, covariate", |
| 96 | + "adjustment, and the k=4->3->2 selector-fit fallback ladder; per", |
| 97 | + "config the full vars_bins table and downsampled vars_poly curve." |
| 98 | + ), |
| 99 | + r_version = R.version.string |
| 100 | +) |
| 101 | + |
| 102 | +## --- Senate data (vendored CSV; same complete-case filter as rdplot) ---- |
| 103 | +senate_path <- "benchmarks/data/rdrobust_senate.csv" |
| 104 | +stopifnot(file.exists(senate_path)) |
| 105 | +senate <- read.csv(senate_path) |
| 106 | +ok <- complete.cases(senate$vote, senate$margin) |
| 107 | +sv <- senate$vote[ok] |
| 108 | +sm <- senate$margin[ok] |
| 109 | + |
| 110 | +golden$senate <- list( |
| 111 | + csv = "benchmarks/data/rdrobust_senate.csv", |
| 112 | + configs = list( |
| 113 | + default = run_rdplot("default", sv, sm), |
| 114 | + es = run_rdplot("es", sv, sm, binselect = "es"), |
| 115 | + espr = run_rdplot("espr", sv, sm, binselect = "espr"), |
| 116 | + esmvpr = run_rdplot("esmvpr", sv, sm, binselect = "esmvpr"), |
| 117 | + qs = run_rdplot("qs", sv, sm, binselect = "qs"), |
| 118 | + qspr = run_rdplot("qspr", sv, sm, binselect = "qspr"), |
| 119 | + qsmv = run_rdplot("qsmv", sv, sm, binselect = "qsmv"), |
| 120 | + qsmvpr = run_rdplot("qsmvpr", sv, sm, binselect = "qsmvpr"), |
| 121 | + p2 = run_rdplot("p2", sv, sm, p = 2), |
| 122 | + p0 = run_rdplot("p0", sv, sm, p = 0), |
| 123 | + p1_tri_h20 = run_rdplot("p1_tri_h20", sv, sm, p = 1, kernel = "tri", |
| 124 | + h = 20), |
| 125 | + epa_h_asym = run_rdplot("epa_h_asym", sv, sm, kernel = "epa", |
| 126 | + h = c(15, 25)), |
| 127 | + nbins_asym = run_rdplot("nbins_asym", sv, sm, nbins = c(10, 14)), |
| 128 | + scale2 = run_rdplot("scale2", sv, sm, scale = 2), |
| 129 | + support_w = run_rdplot("support_w", sv, sm, support = c(-110, 110)), |
| 130 | + ci90 = run_rdplot("ci90", sv, sm, ci = 90), |
| 131 | + qs_nbins = run_rdplot("qs_nbins", sv, sm, binselect = "qs", nbins = 8) |
| 132 | + ) |
| 133 | +) |
| 134 | + |
| 135 | +## --- Tied running variable (masspoints machinery) ----------------------- |
| 136 | +# Same seeded recipe as generate_rdrobust_estimates_golden.R dgp_ties. |
| 137 | +set.seed(123) |
| 138 | +n2 <- 800 |
| 139 | +x2 <- round(2 * rbeta(n2, 2, 4) - 1, 2) |
| 140 | +y2 <- 0.5 + 0.8 * x2 + (x2 >= 0) * 1.0 + rnorm(n2, sd = 0.3) |
| 141 | + |
| 142 | +golden$dgp_ties <- list( |
| 143 | + x = x2, y = y2, |
| 144 | + configs = list( |
| 145 | + # default esmv under adjust remaps to esmvpr (rdplot.R:130-135) |
| 146 | + adjust = run_rdplot("adjust", y2, x2, masspoints = "adjust"), |
| 147 | + off = run_rdplot("off", y2, x2, masspoints = "off"), |
| 148 | + # explicit esmvpr with detection off must EQUAL the adjust config |
| 149 | + # (the Python suite asserts the cross-config identity) |
| 150 | + esmvpr_off = run_rdplot("esmvpr_off", y2, x2, binselect = "esmvpr", |
| 151 | + masspoints = "off") |
| 152 | + ) |
| 153 | +) |
| 154 | + |
| 155 | +## --- Covariates (same seeded recipe as the estimates golden) ------------ |
| 156 | +set.seed(2718) |
| 157 | +n5 <- 1200 |
| 158 | +x5 <- 2 * rbeta(n5, 2, 4) - 1 |
| 159 | +zlong <- 0.5 * x5 + rnorm(n5, sd = 0.8) |
| 160 | +zb <- rbinom(n5, 1, 0.4) |
| 161 | +y5 <- 0.4 * x5 + 0.9 * (x5 >= 0) + 0.7 * zlong + 0.3 * zb + rnorm(n5, sd = 0.3) |
| 162 | +zdup <- 1.5 * zlong - 0.5 * zb |
| 163 | +covs2 <- cbind(zlong = zlong, zb = zb) # UNSORTED names: pins order(nchar) |
| 164 | +covs3 <- cbind(zlong = zlong, zb = zb, zdup = zdup) |
| 165 | + |
| 166 | +golden$dgp_covs <- list( |
| 167 | + x = x5, y = y5, zlong = zlong, zb = zb, zdup = zdup, |
| 168 | + configs = list( |
| 169 | + covs_default = run_rdplot("covs_default", y5, x5, covs = covs2), |
| 170 | + covs_collinear = run_rdplot("covs_collinear", y5, x5, covs = covs3) |
| 171 | + ) |
| 172 | +) |
| 173 | + |
| 174 | +## --- Small sample (n = 25; just above rdplot's n >= 20 hard floor) ------ |
| 175 | +set.seed(555) |
| 176 | +n6 <- 25 |
| 177 | +x6 <- runif(n6, -1, 1) |
| 178 | +y6 <- 0.3 + 0.5 * x6 + 0.8 * (x6 >= 0) + rnorm(n6, sd = 0.4) |
| 179 | + |
| 180 | +golden$dgp_small <- list( |
| 181 | + x = x6, y = y6, |
| 182 | + configs = list(default = run_rdplot("small_default", y6, x6)) |
| 183 | +) |
| 184 | + |
| 185 | +## --- Selector-fit fallback ladder (k = 4 -> 3) -------------------------- |
| 186 | +# |x| ~ 1e40 keeps every x^4 design entry finite but overflows the k=4 |
| 187 | +# Gram (entries ~x^8 -> Inf), so R's try(qrXXinv) errors and the ladder |
| 188 | +# drops to k=3 (rdplot.R:281-298). Extreme magnitudes make float surfaces |
| 189 | +# BLAS-sensitive, so the Python suite asserts only the integer/string |
| 190 | +# outputs for this DGP. |
| 191 | +set.seed(999) |
| 192 | +n7 <- 120 |
| 193 | +x7 <- c(-runif(50, 1e39, 2e39), runif(70, 1e39, 2e39)) |
| 194 | +y7 <- rnorm(n7) |
| 195 | + |
| 196 | +golden$dgp_ladder <- list( |
| 197 | + x = x7, y = y7, |
| 198 | + configs = list(default = run_rdplot("ladder_default", y7, x7)) |
| 199 | +) |
| 200 | + |
| 201 | +out_path <- "benchmarks/data/rdplot_golden.json" |
| 202 | +write_json(golden, out_path, auto_unbox = TRUE, pretty = TRUE, digits = I(17)) |
| 203 | +cat("Wrote", out_path, "\n") |
0 commit comments