iRBCrd: IMPROVED ROBUST BIAS CORRECTION FOR RD DESIGNS Version 2.1 -- October 29, 2026 Website: https://irbc-packages.github.io/ AUTHORS AND CONTACTS -------------------- Giuseppe Cavaliere University of Bologna and Exeter Business School giuseppe.cavaliere@unibo.it Silvia Goncalves McGill University silvia.goncalves@mcgill.ca Morten Orregaard Nielsen Aarhus University mon@econ.au.dk Edoardo Zanelli Aarhus University ez@econ.au.dk 1. OVERVIEW ----------- iRBCrd provides bias-corrected point estimates, analytic confidence intervals, and bootstrap confidence intervals for: - Sharp regression discontinuity (RD). - Fuzzy RD. - Covariate-adjusted sharp RD. - Covariate-adjusted fuzzy RD, adjusting both outcome and treatment equations. - Sharp and fuzzy kink RD, with or without covariate adjustment. Use iRBCrd() for every design. Supply fuzzy for a treatment variable and covs for covariate adjustment. Add deriv = 1 for a kink design. The function selects the appropriate implementation from these inputs; deriv = NULL or 0 selects level RD. This is a distribution of R functions loaded with source(), not an installable R package or a CRAN release. No license file is included. Bootstrap inference uses a fixed-regressor wild bootstrap. 2. FILES AND LOADING -------------------- iRBC_functions.R Kernels, local-polynomial fits, variance calculations, and bandwidth regularization helpers. iRBCrd.R Unified entry point, level RD, and kink RD implementations. README.txt This manual. examples.R Reproducible level and kink examples. CHANGELOG.txt Release overview. smoke_test.R Standalone interface and bootstrap checks. kink_test.R Independent kink formulas, bootstrap refits, and option checks. output_test.R Returned-field and bandwidth-reporting checks for all designs. bootstrap_weights_test.R Independent refits for bootstrap weighting choices. VERIFICATION.txt Summary of the development regression checks. CONTACT_SOURCES.txt Sources used to verify the institutional contacts. Extract the ZIP into a folder. With R's working directory set to that folder: source("iRBC_functions.R") source("iRBCrd.R") Both files are needed, in that order. No additional R packages are required for estimation. rdrobust is optional for bandwidth selection or comparisons. Loading these files does not attach or install any additional packages. To run the examples or standalone checks after setting the working directory: source("examples.R") source("smoke_test.R") source("kink_test.R") source("output_test.R") source("bootstrap_weights_test.R") The functions can be used in R, RStudio, or via Rscript. They do not require RStudio or a compiler. This release was tested in R 4.5.2 on macOS. 3. QUICK START -------------- set.seed(17041998) n <- 400 x <- runif(n, -1, 1) z <- 0.3 * x + rnorm(n) d <- rbinom(n, 1, 0.2 + 0.6 * (x >= 0)) y <- 1 + x + x^2 + 0.8 * d + 0.5 * z + rnorm(n) Z <- cbind(z = z) sharp <- iRBCrd(y, x, h = 0.6) fuzzy_fit <- iRBCrd(y, x, h = 0.6, fuzzy = d) adjusted_sharp <- iRBCrd(y, x, h = 0.6, covs = Z) adjusted_fuzzy <- iRBCrd(y, x, h = 0.6, fuzzy = d, covs = Z) The bandwidth above is illustrative. It is not an estimated optimal bandwidth. The bootstrap and analytic reportable estimates are identical. Their intervals can differ, and the bootstrap interval need not be symmetric around that bias-corrected estimate. 4. FUNCTION INTERFACE --------------------- iRBCrd(y, x, h, c = NULL, p = NULL, kernel = NULL, fuzzy = NULL, deriv = NULL, covs = NULL, res = NULL, level = 95, bwcheck = .bwcheck_default, fast = NULL, boot.ci = NULL, boot.B = NULL, boot.seed = NULL, boot.wild = NULL, boot.keep = NULL, q.method = NULL, fuzzy.boot.weights = NULL, gamma.boot.weights = NULL) Design and estimation: y, x Numeric vectors in the same observation order. Prepare finite, complete-case data before calling the function. h Required positive bandwidth, scalar or a pair ordered left, right. c Cutoff; default 0. The right side includes x = c. p Local-polynomial order; default 1 for levels and 2 for kinks. Use a nonnegative integer for levels, and p >= 1 for kinks. kernel "tri" (default), "epa", or "uni". Long aliases "triangular", "epanechnikov", and "uniform" are accepted. A pair is allowed, ordered left, right, or named left/right or minus/plus. fuzzy Treatment vector. NULL selects sharp RD. This is not a Boolean. Binary, discrete, and continuous treatments are accepted. deriv Derivative order: 0 for discontinuities in levels, 1 for slope changes. NULL defaults to 0. Only numeric 0 and 1 are accepted. covs Numeric covariate vector, matrix, or data frame; default NULL. One row per observation. Both equations are adjusted in fuzzy RD. res "irbc-hc0" through "irbc-hc3", or "cct-hc0" through "cct-hc3". The default is pointwise iRBC-HC3 ("irbc-hc3"). level Confidence level in percent; default 95 gives a 95% interval. Supply a finite number strictly between 0 and 100. q.method "sample" or "asymptotic". Default "sample" when fast = FALSE; "asymptotic" when fast = TRUE. fast Standard-error shortcut; default FALSE. See Section 8. bwcheck Minimum-neighbor regularization rule; default 21. NULL disables it. Bootstrap: boot.ci TRUE to add a bootstrap interval; default FALSE. boot.B Number of draws; default 999 when boot.ci = TRUE. A supplied value must be a positive integer. At least two finite draws are needed for an interval; use substantially more for inference. boot.seed Optional numeric seed. Default NULL uses the current R RNG state. A supplied seed is passed to set.seed() after integer conversion. The prior random-number state is not restored. boot.wild "rademacher" (default), "mammen", or "normal". boot.keep Retain bootstrap statistic draws and weight diagnostics; default FALSE. fuzzy.boot.weights "fixed" (default) or "bootstrap" in every fuzzy design, including kink and covariate-adjusted designs. gamma.boot.weights "fixed" (default) or "bootstrap" in every covariate-adjusted design, including sharp, fuzzy, and kink designs. 5. BANDWIDTHS AND DATA ---------------------- The function does not select an optimal bandwidth. Supply h, for example from rdrobust::rdbwselect(), with the appropriate fuzzy and covariate inputs in the selector. For a kink, also set deriv = 1 and the same p in the selector. An unnamed bandwidth pair is left, right. Named pairs accept left/right or minus/plus. A scalar applies to both sides. Supply bandwidths only through h. Kernels follow the same convention through kernel. fit <- iRBCrd(y, x, h = c(left = 0.5, right = 0.7), kernel = c(left = "tri", right = "epa")) When using rdrobust for a comparison with rho = 1, explicitly use the same selected h and set b = h there. For an exact bandwidth comparison, consider bwcheck = NULL here, provided all needed local fits remain nonsingular. With bwcheck = k, a local bandwidth is enlarged if necessary to reach the kth closest observation, or the farthest observation if fewer than k are present. This counts observations, not distinct running-variable values. The initial cutoff regularization uses pooled observations for scalar h, and side-specific observations for paired bandwidth inputs. Auxiliary fits can enlarge their own bandwidths as well. For kink RD, cutoff regularization is side-specific even for scalar h. The reported h.minus and h.plus are the actual slope-fitting bandwidths. These effective bandwidths also start the side's auxiliary fits. pilot.h records any further bandwidth enlargement in the level pilots used for the bias. Kink bias, variance, and bootstrap calculations use the same coefficient maps, including when a pilot neighborhood is enlarged. Returned h is a named list of requested minus/plus bandwidths. h.minus and h.plus are output fields, not inputs; they report the effective cutoff bandwidths, not the bandwidth at every auxiliary fitting point. This rule does not guarantee full rank with ties, collinearity, or zero endpoint kernel weights. Do not trim the data to the cutoff h window before calling iRBCrd(): auxiliary fits need observations outside that window. Bias and variance calculations can involve observations up to 2h from the cutoff when regularization is disabled; pointwise residual fits use their own neighbors. 6. EXTENSIONS ------------- All extensions use iRBCrd() and support analytic and bootstrap inference. The design inputs can be combined; bandwidth and kernel conventions remain as in Section 5. Kink RD: Set deriv = 1 to estimate the right slope minus the left slope at the cutoff. The default is local quadratic estimation (p = 2); p = 1 or a higher order can also be supplied. Sharp kink reports the slope difference in units of y/x, without dividing by a known policy-schedule slope change. For that target, rescale y by the known change before estimation. Bias correction uses slope-specific Q adjustments; residuals and bootstrap outcomes remain in levels. kink_fit <- iRBCrd(y, x, h, deriv = 1) Fuzzy RD: Supply the treatment vector through fuzzy. The target is the outcome discontinuity divided by the treatment discontinuity, or the corresponding ratio of slope changes when deriv = 1. iRBC corrects the linearized ratio and retains the covariance between the outcome and treatment components. A near-zero first stage prevents ordinary ratio inference; a nonzero estimate alone does not establish a strong first stage. Bootstrap weighting choices are described in Section 7. fuzzy_fit <- iRBCrd(y, x, h, fuzzy = d) Covariate adjustment: Supply covs to add common covariate slopes to the side-specific order-p weighted regressions. The outcome discontinuity is adjusted by the estimated covariate coefficients times the covariate discontinuities. Fuzzy designs adjust both outcome and treatment equations, using separate coefficient vectors and retaining their covariance. In kink designs, the adjustment uses covariate slope changes. Preserving the unadjusted target requires covariate balance at the cutoff: continuity for level RD and no slope change for kink RD, alongside the design's identification conditions. adjusted_fit <- iRBCrd(y, x, h, covs = Z) 7. BOOTSTRAP ------------- The running variable, bandwidths, and Q adjustments are held fixed. For local observations, bootstrap outcomes equal pointwise local-polynomial fitted values plus HC-adjusted residuals multiplied by wild weights. The same draw is used for all components of an observation in joint designs, preserving their dependence. Observations outside the cutoff working window are left unchanged; the bandwidth is not re-estimated in each bootstrap draw. For each sharp component j, the bootstrap statistic is T_j* = Q_plus [ghat_j,plus*(c) - ghat_j,plus(c)] - Q_minus [ghat_j,minus*(c) - ghat_j,minus(c)]. For kink designs, replace the cutoff level fits in this expression with cutoff slope fits. Bootstrap outcomes and residuals are still constructed in levels, and Q uses the slope-specific adjustment. The sharp, fuzzy, and covariate-adjusted statistics use the relevant component linear combination. By default, fuzzy and covariate weights remain fixed. The two choices can be changed independently in level and kink designs: fuzzy.boot.weights = "bootstrap" uses the bootstrap first-stage discontinuity (or slope change) in the denominator of the fuzzy weights. The original outcome-to-treatment ratio remains fixed, as in the linearized statistic. gamma.boot.weights = "bootstrap" refits the common covariate coefficients using bootstrap outcomes and covariates. Fuzzy designs also refit the treatment equation. The bootstrap first stage is adjusted using whichever covariate coefficients are selected by gamma.boot.weights. Thus, covariate-adjusted fuzzy designs allow all four fixed/bootstrap combinations. These options change the bootstrap statistic's combination weights, not the original analytic estimate, variance, or Q adjustments. Draws with singular covariate refits or undefined fuzzy weights are omitted from the interval calculation. Let se_us denote the standard error of the uncorrected Q-weighted component combination, and mhat = se_irbc/se_us. The interval uses empirical quantiles of the bootstrap statistic at pnorm(mhat * qnorm(c((100-level)/200, (100+level)/200))). Quantiles use R's type = 7 interpolation and only finite draws. They are inverted around the corresponding uncorrected discontinuity or ratio estimate. The REPORTABLE bootstrap point estimate is nevertheless the same bias-corrected estimate as for analytic iRBC. This distinction is intentional. Rademacher weights take -1 or 1 with equal probability. Mammen weights use the usual mean-zero, variance-one asymmetric two-point distribution. Normal weights are also available. Bootstrap quantiles have simulation error, especially at high confidence levels. B = 999 is a default, not a guarantee of accurate tail quantiles. Check stability across seeds or larger B when differences between intervals are small. 8. OTHER DESIGN CHOICES ----------------------- Residuals: iRBC-HCk fits a p-th order local polynomial separately on each cutoff side, centered at each X_i. The raw residual is Y_i - ghat(X_i). There is no additional Q-weighted bias subtraction in this fitted value. The HC adjustment uses the pointwise fit's leverage and local sample size. CCT-HCk fits a (p+1)-th order local polynomial at the cutoff on each side. Evaluate it at X_i, subtract it from Y_i, and apply the HC adjustment using that cutoff fit's leverage and local sample size. Both residual choices use fits in levels, including for kink designs. For raw residual e_i and leverage ell_i: HC0: e_i HC1: e_i / sqrt((n_local - df) / n_local) HC2: e_i / sqrt(1 - ell_i) HC3: e_i / (1 - ell_i) Here df = p+1 for iRBC residuals and df = p+2 for CCT residuals; n_local counts nonzero-kernel observations. The same residuals are used in the variance estimator and bootstrap data-generating process. In joint designs, residuals are evaluated componentwise and their cross-products enter the covariance. "hc0" through "hc3" select iRBC-HCk; "rbc-hc0" through "rbc-hc3" select CCT-HCk. The returned res field uses "irbc-hck" or "cct-hck", respectively, with k replaced by 0, 1, 2, or 3. res = "loo" uses unscaled residuals from the ordinary pointwise auxiliary fits, without an HC adjustment. Q and analytic inference: q.method = "sample" uses finite-sample Q_n values (the default). q.method = "asymptotic" uses kernel-specific population values; keep fast = FALSE to retain the full variance calculation. Asymptotic Q requires p = 1 for level RD and also supports higher orders for kink RD. fast = TRUE uses asymptotic Q and a standard-error shortcut for sharp designs, with slope-specific constants for kink RD; joint designs retain full covariance matrices. It cannot use sample Q. The analytic interval is estimate_irbc +/- qnorm((100+level)/200) * se_irbc. 9. RESULTS ----------- The primary iRBCrd() result fields are: estimate_irbc Bias-corrected point estimate. se_irbc Analytic standard error. ci_irbc Analytic interval, lower then upper. level Confidence level in percent. deriv 0 for level RD; 1 for kink RD. estimate_irbc_boot Same point estimate, when boot.ci = TRUE. ci_irbc_boot Bootstrap interval, lower then upper. h Requested bandwidths, as a named minus/plus list. h.minus, h.plus Effective cutoff bandwidths. kernel, kernel.minus, kernel.plus Kernel specifications used on the two sides. res Residual specification. boot.B, boot.wild, boot.t Bootstrap details; boot.t is omitted by default. Set boot.keep = TRUE to retain bootstrap statistics and, where applicable, component statistics, covariate coefficients, and combination weights for individual draws. This does not change point estimates or confidence intervals. Design-specific details include uncorrected discontinuities, bias components, first-stage estimates, covariate slopes, covariance matrices, and sharp.* component objects. Residual arrays are computational inputs over their local support, not a full-sample diagnostic for fitted values. Kink results also report deriv = 1, p, design, linear.weights, kink.weights, and pilot.h. Component estimates are slope changes. Vectors in kink.weights and the residuals follow the original observation order within each side. Sharp kink returns tau.hat; fuzzy kink returns zeta.hat. 10. REFERENCES --------------- Cavaliere, G., Goncalves, S., Nielsen, M. O., and Zanelli, E. (2026): Improved Inference for Nonparametric Regression and Regression-Discontinuity Designs. arXiv:2512.00566. https://arxiv.org/abs/2512.00566. Calonico, S., Cattaneo, M. D., Farrell, M. H., and Titiunik, R. (2019): Regression Discontinuity Designs Using Covariates, Review of Economics and Statistics, 101(3), 442-451. The code is supplied as R functions. No install.packages("iRBCrd") or library(iRBCrd) step applies to this distribution. See CHANGELOG.txt for the release overview and VERIFICATION.txt for tests performed before preparing the ZIP.