Regression Discontinuity#

Regression discontinuity estimation - sharp and fuzzy, with optional covariate adjustment - and robust bias-corrected inference, parity-targeting R rdrobust 4.0.0.

Sharp (default): treatment is assigned by a known threshold of an observed running variable (running >= cutoff; units exactly at the cutoff are treated, matching rdrobust). The effect is the jump in the conditional expectation of the outcome at the cutoff. Fuzzy (pass the observed take-up column via fit(..., treatment_col=...)): crossing the cutoff shifts take-up instead of determining it, and the estimand is the local Wald ratio - for binary take-up under monotonicity, the LATE for compliers at the cutoff; for non-binary take-up, the ratio of jumps (the estimand field says which) - with the first stage exposed as a full first_stage* block. Both designs use kernel-weighted local polynomials on each side with data-driven MSE/CER-optimal bandwidths (all 10 rdrobust selectors) and robust bias-corrected inference per Calonico, Cattaneo & Titiunik (2014). Covariate adjustment (fit(..., covariates=[...]); R’s covs=, per Calonico, Cattaneo, Farrell & Titiunik 2019): additive common-coefficient adjustment that leaves the estimand UNCHANGED - unlike the DiD estimators’ conditional-parallel-trends role, RD covariates buy precision only, and require covariate balance at the cutoff (testable by fitting each covariate as the outcome). Bandwidths are covariate-aware; collinear covariates are dropped with a warning naming them (covs_drop=True, R’s default).

Note

Canonical inference binding. The result’s att, se, t_stat, p_value, and conf_int are one internally coherent row - the ROBUST bias-corrected row (att is the bias-corrected estimate - the linearized bias-corrected ratio on fuzzy fits; conf_int is centered on it; t_stat == att/se). The estimand field names what att measures for the fit at hand. R’s rdrobust prints the conventional estimate as its headline coefficient while taking inference from the robust row; that estimate is available as att_conventional with its own full inference row, and summary() prints the familiar three-row rdrobust table (plus a first-stage block on fuzzy fits, as R does).

Note

Scope of this release. Sharp, fuzzy, and covariate-adjusted designs with the nearest-neighbor variance estimator (rdrobust’s default), plus the RDPlot diagnostic below. Cluster-robust variance, weights, kink estimands, weak-IV-robust fuzzy inference, a packaged covariate-balance helper, and the density-test diagnostic are documented follow-ups - see the methodology registry for the full deviations and seams list.

RD plots (RDPlot; Calonico, Cattaneo & Titiunik 2015, parity with R’s rdplot()): the standard exploratory companion - per-side global polynomial fits plus binned local means with a data-driven number of bins (all 8 binselect selectors: evenly/quantile spaced, IMSE-optimal or mimicking-variance, spacings or polynomial-regression variance estimators), implied-scale/WIMSE-weight reporting, optional covariate adjustment, and an optional matplotlib rendering via RDPlotResult.plot() (matplotlib is not a dependency; the numbers are the parity surface).

RegressionDiscontinuity#

class diff_diff.RegressionDiscontinuity[source]

Bases: object

Regression discontinuity estimator - sharp and fuzzy, with optional covariate adjustment (rdrobust 4.0.0 parity).

SHARP (default): treatment is defined by the running variable crossing a known cutoff (running >= cutoff treated, matching rdrobust: units exactly at the cutoff are treated). FUZZY: pass the observed take-up column via fit(..., treatment_col=...) - the estimand becomes the local Wald ratio (complier LATE at the cutoff for binary take-up under monotonicity; the estimand results field says which reading applies) and the results gain a first-stage block. COVARIATE ADJUSTMENT: pass fit(..., covariates=[...]) (R’s covs=) for the CCFT 2019 additive common-coefficient adjustment - the estimand is UNCHANGED (precision only; requires covariate balance at the cutoff, see the module docstring), bandwidths become covariate-aware, and collinear columns are dropped with a warning under covs_drop=True. Point estimation uses kernel-weighted local polynomials of order p on each side; inference is robust bias-corrected per Calonico, Cattaneo & Titiunik (2014). Defaults reproduce rdrobust(y, x) / rdrobust(y, x, fuzzy=t) / rdrobust(y, x, covs=Z): p=1, q=2, triangular kernel, bwselect="mserd", nearest-neighbor variance with 3 matches, masspoints="adjust", covs_drop=True.

Parameters:
  • cutoff (float, default 0.0) – The known threshold c of the running variable.

  • p (int, default 1) – Local-polynomial order for point estimation; integer in 0..20 (mirroring rdrobust’s accepted surface; p=0 is the local-constant fit).

  • q (int or None, default None) – Order for the bias regression; an explicit q must satisfy p < q <= 20. None resolves to p + 1 WITHOUT re-validation, exactly as R does (rdrobust.R:53-57) - so p=20 with the default q yields q=21 in both implementations while an explicit q=21 is rejected.

  • kernel (str, default "triangular") – “triangular”, “epanechnikov”, or “uniform” (R spellings “tri”/”epa”/”uni” accepted).

  • bwselect (str, default "mserd") – Data-driven bandwidth selector; one of the 10 rdrobust options (mserd, msetwo, msesum, msecomb1, msecomb2, cerrd, certwo, cersum, cercomb1, cercomb2). Ignored when h is supplied.

  • h (float or None) – Manual main / bias bandwidths (both sides). h alone implies b = h; h with rho implies b = h/rho (overriding a supplied b, as in R); b without h is ignored with a warning (R ignores it silently - documented deviation).

  • b (float or None) – Manual main / bias bandwidths (both sides). h alone implies b = h; h with rho implies b = h/rho (overriding a supplied b, as in R); b without h is ignored with a warning (R ignores it silently - documented deviation).

  • rho (float or None) – Bandwidth ratio h/b. Without h, applies to the SELECTED bandwidths (b = h_selected/rho), mirroring rdrobust.

  • vcov_type (str, default "nn") – Variance estimator. Only “nn” (same-side nearest-neighbor, rdrobust’s default) is implemented in this release; “hc0”-“hc3” and cluster modes raise NotImplementedError.

  • nnmatch (int, default 3) – Minimum number of nearest neighbors for the NN variance.

  • masspoints (str, default "adjust") – Mass-point handling: “adjust” (rdrobust default), “check”, “off”.

  • bwcheck (int or None, default None) – Minimum unique support points forced inside the bandwidth window.

  • bwrestrict (bool, default True) – Clamp bandwidths to the running variable’s observed range.

  • scaleregul (float, default 1.0) – Scale of the IK-style regularization in bandwidth selection (0 removes it).

  • sharpbw (bool, default False) – Fuzzy fits only (fit(..., treatment_col=...)): when True, bandwidths are selected for the SHARP reduced-form estimator on the outcome (rdrobust’s “approach 1”) instead of the default fuzzy-ratio objective. Automatically in effect - regardless of this flag - under one-sided perfect compliance (zero take-up variance on either side), exactly as in R. On sharp fits the flag has no effect and a warning is emitted (R ignores it silently - documented deviation). Never drops covariates from selection - with covariates it selects on the covariate-adjusted sharp objective, as in R.

  • covs_drop (bool, default True) – Covariate-adjusted fits only (fit(..., covariates=[...])): when True (R’s default), redundant (collinear) covariate columns are dropped with a warning naming them before fitting, and the covariate projection uses a pseudo-inverse; when False the solve is strict and collinear covariates raise a clear error. Without covariates the flag has no effect and setting it to False emits a warning (same pattern as sharpbw on sharp fits).

  • alpha (float, default 0.05) – Significance level (rdrobust level = 100*(1-alpha)).

Examples

>>> rd = RegressionDiscontinuity(cutoff=0.0)
>>> results = rd.fit(df, outcome_col="y", running_col="x")
>>> results.att, results.conf_int  # robust bias-corrected inference
>>> fuzzy = rd.fit(df, "y", "x", treatment_col="takeup")  # fuzzy RD
>>> fuzzy.att, fuzzy.first_stage  # local Wald ratio + take-up jump
__init__(cutoff=0.0, p=1, q=None, kernel='triangular', bwselect='mserd', h=None, b=None, rho=None, vcov_type='nn', nnmatch=3, masspoints='adjust', bwcheck=None, bwrestrict=True, scaleregul=1.0, sharpbw=False, covs_drop=True, alpha=0.05)[source]
Parameters:
get_params(deep=True)[source]

Return raw constructor parameters (sklearn-compatible).

Parameters:

deep (bool)

Return type:

Dict[str, Any]

set_params(**params)[source]

Transactionally update parameters (validate before mutating).

Parameters:

params (Any)

Return type:

RegressionDiscontinuity

fit(data, outcome_col, running_col, treatment_col=None, covariates=None)[source]

Estimate the RD effect at the cutoff (sharp or fuzzy, optionally covariate-adjusted).

Parameters:
  • data (pd.DataFrame) – Cross-sectional data.

  • outcome_col (str) – Column names of the outcome and the running variable.

  • running_col (str) – Column names of the outcome and the running variable.

  • treatment_col (str or None, default None) – None (sharp design): treatment is derived as running >= cutoff; no treatment column is needed. A column name activates the FUZZY design: the column holds the OBSERVED treatment take-up (typically binary, any numeric accepted, matching R’s fuzzy=), the estimand becomes the local Wald ratio, and the results gain the first_stage* block. The estimand label is data-dependent: for BINARY take-up it reads “fuzzy (LATE for compliers at the cutoff)” (the monotonicity-based complier reading); for non-binary (dose) take-up it reads “fuzzy (local Wald ratio at the cutoff; non-binary take-up)” - the complier-LATE interpretation does not apply there. A take-up column that is deterministic in the running variable reproduces the sharp fit exactly (first stage == 1).

  • covariates (list of str or None, default None) – Column names of pre-determined covariates for the additive common-coefficient adjustment of CCFT (2019) (R’s covs=). The estimand is UNCHANGED - unlike the DiD estimators’ conditional-parallel-trends role, RD covariates buy precision only, and require covariate BALANCE at the cutoff (zero RD effect on each covariate; testable by fitting each covariate as the outcome - imbalanced covariates make the adjusted estimator inconsistent). Continuous, discrete, or mixed columns are accepted; covariates propagate into bandwidth selection (covariate-aware, as in R). Collinear columns are dropped with a warning under covs_drop=True; see the covariates* results fields for the echo and the fitted projection coefficients.

Return type:

RegressionDiscontinuityResults

RegressionDiscontinuityResults#

class diff_diff.RegressionDiscontinuityResults[source]

Bases: BaseResults

Results of a regression discontinuity fit (sharp or fuzzy; the estimand field names which one, and first_stage* fields are populated on fuzzy fits only).

Canonical inference fields (att, se, t_stat, p_value, conf_int) all describe the ROBUST bias-corrected row: att is the bias-corrected point estimate and conf_int is centered on it (see module docstring for the binding rationale and the deviation from rdrobust’s printed headline). The conventional row is exposed as explicit *_conventional fields. The (rarely used) middle “Bias-Corrected” row shares its coefficient with att (both are tau_bc) and its standard error with se_conventional - only its inference triple carries the *_bias_corrected suffix (t_stat_bias_corrected, p_value_bias_corrected, conf_int_bias_corrected); there are deliberately no redundant att_bias_corrected / se_bias_corrected fields. Together the three rows mirror rdrobust’s output exactly.

att: float
se: float
t_stat: float
p_value: float
conf_int: Tuple[float, float]
alpha: float
att_conventional: float
se_conventional: float
t_stat_conventional: float
p_value_conventional: float
conf_int_conventional: Tuple[float, float]
t_stat_bias_corrected: float
p_value_bias_corrected: float
conf_int_bias_corrected: Tuple[float, float]
se_robust: float
h_left: float
h_right: float
b_left: float
b_right: float
n_obs: int
n_left: int
n_right: int
n_h_left: int
n_h_right: int
n_b_left: int
n_b_right: int
n_unique_left: int
n_unique_right: int
n_dropped: int
cutoff: float
p: int
q: int
kernel: str
bwselect: str
vcov_type: str
nnmatch: int
masspoints: str
bwcheck: int | None
bwrestrict: bool
scaleregul: float
h_input: float | None
b_input: float | None
rho_input: float | None
estimand: str
sharpbw: bool
treatment_col: str | None
covs_drop: bool
first_stage: float | None = None
first_stage_se: float | None = None
first_stage_t_stat: float | None = None
first_stage_p_value: float | None = None
first_stage_conf_int: Tuple[float, float] | None = None
first_stage_conventional: float | None = None
first_stage_se_conventional: float | None = None
first_stage_t_stat_conventional: float | None = None
first_stage_p_value_conventional: float | None = None
first_stage_conf_int_conventional: Tuple[float, float] | None = None
first_stage_t_stat_bias_corrected: float | None = None
first_stage_p_value_bias_corrected: float | None = None
first_stage_conf_int_bias_corrected: Tuple[float, float] | None = None
covariates: List[str] | None = None
covariates_dropped: List[str] | None = None
covariate_coefficients: Dict[str, float] | None = None
first_stage_covariate_coefficients: Dict[str, float] | None = None
beta_p_left: ndarray = None
beta_p_right: ndarray = None
beta_t_p_left: ndarray | None = None
beta_t_p_right: ndarray | None = None
summary()[source]

Human-readable summary with the three-row rdrobust table.

Return type:

str

print_summary()[source]
Return type:

None

to_dict()[source]

Flat scalar dict; confidence intervals split into lower/upper.

Return type:

Dict[str, Any]

to_dataframe()[source]
Return type:

DataFrame

__init__(att, se, t_stat, p_value, conf_int, alpha, att_conventional, se_conventional, t_stat_conventional, p_value_conventional, conf_int_conventional, t_stat_bias_corrected, p_value_bias_corrected, conf_int_bias_corrected, se_robust, h_left, h_right, b_left, b_right, n_obs, n_left, n_right, n_h_left, n_h_right, n_b_left, n_b_right, n_unique_left, n_unique_right, n_dropped, cutoff, p, q, kernel, bwselect, vcov_type, nnmatch, masspoints, bwcheck, bwrestrict, scaleregul, h_input, b_input, rho_input, estimand, sharpbw, treatment_col, covs_drop, first_stage=None, first_stage_se=None, first_stage_t_stat=None, first_stage_p_value=None, first_stage_conf_int=None, first_stage_conventional=None, first_stage_se_conventional=None, first_stage_t_stat_conventional=None, first_stage_p_value_conventional=None, first_stage_conf_int_conventional=None, first_stage_t_stat_bias_corrected=None, first_stage_p_value_bias_corrected=None, first_stage_conf_int_bias_corrected=None, covariates=None, covariates_dropped=None, covariate_coefficients=None, first_stage_covariate_coefficients=None, beta_p_left=None, beta_p_right=None, beta_t_p_left=None, beta_t_p_right=None)
Parameters:
Return type:

None

RDPlot#

class diff_diff.RDPlot[source]

Bases: object

Data-driven RD plot builder (CCT 2015; rdrobust 4.0.0 rdplot() parity).

Parameters mirror R’s rdplot() and keep its defaults - including kernel="uniform" and p=4, which deliberately differ from RegressionDiscontinuity’s estimation defaults (triangular, p=1): the global fit is a descriptive overlay, not the local RD estimator.

Parameters:
  • cutoff (float, default 0.0) – Threshold in the running variable.

  • p (int, default 4) – Order of the global polynomial fits (R accepts any integer >= 0).

  • nbins (int or (int, int), optional) – Manual number of bins per side; overrides the data-driven selector. Any positive count is accepted (R permissiveness); note that very large bin counts allocate proportional jump grids and, with covariates, a dense per-side bin-dummy matrix (documented performance seam).

  • binselect (str, default "esmv") – One of es, espr, esmv, esmvpr, qs, qspr, qsmv, qsmvpr - partition scheme (evenly/quantile spaced) x target (IMSE-optimal / mimicking variance) x variance estimator (spacings / polynomial regression). The spacings variants assume a continuously distributed OUTCOME (CCT 2015 Theorems 3-4); heavy outcome ties trigger a warning recommending the *pr sibling. QS cutpoints follow R: quantile type-7 interpolation, not the paper’s inf-form empirical inverse (documented deviation-from-paper inherited from rdrobust).

  • scale (float or (float, float), optional) – Multiplies the selected number of bins (undersmoothing knob). A fractional product is resolved with the ceiling of CCT 2015 Equation 2. Deviation from R: rdrobust 4.0.0 crashes on fractional products (vector-indexing accident) and rejects a length-2 scale outright on R >= 4.2 (vectorized-if error); both surfaces work here.

  • kernel (str, default "uniform") – Weighting for the global polynomial fit within h.

  • h (float or (float, float), optional) – Window for the global fit; default is each side’s full range.

  • support ((float, float), optional) – Extended support for bin construction; only WIDENS the observed range (R semantics).

  • masspoints (str, default "adjust") – "check" warns when either side has >= 20% duplicate running values; "adjust" additionally switches a spacings binselect to its polynomial-regression sibling; "off" disables detection.

  • ci (float, optional) – Confidence level (percent) for the per-bin intervals drawn by RDPlotResult.plot(). CI columns are ALWAYS computed (at 95 when ci is None, matching R); ci only sets the level and turns the error bars on.

  • covs_drop (bool, default True) – Drop redundant covariates (R’s pipeline: stable name-length sort + LINPACK-style QR check) with a warning naming the dropped columns; False raises a deterministic error on collinear covariates.

__init__(cutoff=0.0, p=4, nbins=None, binselect='esmv', scale=None, kernel='uniform', h=None, support=None, masspoints='adjust', ci=None, covs_drop=True)[source]
Parameters:
get_params()[source]

Get parameters of this plot builder.

Return type:

Dict[str, Any]

set_params(**params)[source]

Set parameters (transactional: validates the full candidate configuration before mutating this instance).

Parameters:

params (Any)

Return type:

RDPlot

fit(data, outcome_col, running_col, covariates=None)[source]

Build the RD plot quantities.

Parameters:
  • data (pd.DataFrame) – Cross-sectional data.

  • outcome_col (str) – Column names of the outcome and the running variable.

  • running_col (str) – Column names of the outcome and the running variable.

  • covariates (list of str, optional) – Covariate-adjusted plot per rdplot’s covs= (global fit partials the covariates out with a pooled-across-sides gamma; bin means use the fitted values of per-side regressions of the outcome on covariates and bin indicators, R’s lm(y ~ z + factor(bin)) with covs_eval="mean"). Note R’s caveat: covariate adjustment is meant for plotting alongside covariate-adjusted rdrobust ESTIMATES; the adjusted global fit may not be visually compatible with unadjusted binned means. covariate_coefficients is keyed in R’s processing order: name-length-sorted under covs_drop=True (the sort decides which of a collinear set survives), the user-given order under covs_drop=False (R sorts only inside its drop pipeline; same contract as RegressionDiscontinuity). Per-name values are identical either way on full-rank covariates.

Return type:

RDPlotResult

RDPlotResult#

class diff_diff.RDPlotResult[source]

Bases: Diagnostic

Quantities of a data-driven RD plot (field names follow R’s rdplot return list; vars_bins columns carry R’s rdplot_* names).

J is always integral: a fractional scale * J product takes the ceiling per CCT 2015 Equation 2 (see the Deviation note in RDPlot - R 4.0.0 crashes on that surface instead). vars_bins contains only NON-EMPTY bins (R’s tapply semantics) while bin_avg/bin_med summarize ALL bin lengths of the partition. CI columns are always computed (at ci_level); ci_requested records whether ci= was passed and is the default for drawing error bars in plot().

coef: DataFrame
vars_bins: DataFrame
vars_poly: DataFrame
J: Tuple[float, float]
J_IMSE: Tuple[float, float]
J_MV: Tuple[float, float]
scale: Tuple[float, float]
rscale: Tuple[float, float]
bin_avg: Tuple[float, float]
bin_med: Tuple[float, float]
p: int
cutoff: float
h: Tuple[float, float]
N: Tuple[int, int]
N_h: Tuple[int, int]
binselect: str
kernel_type: str
ci_level: float
ci_requested: bool
covariate_coefficients: Dict[str, float] | None = None
covariates_dropped: List[str] | None = None
property wimse_variance_weight: Tuple[float, float]

Implied WIMSE variance weights 1/(1 + rscale^3) per side (CCT 2015 Supplement S.1; printed by R’s summary.rdplot).

property wimse_bias_weight: Tuple[float, float]

Implied WIMSE bias weights rscale^3/(1 + rscale^3) per side.

to_dataframe()[source]

The non-empty-bin table (vars_bins), copied.

Return type:

DataFrame

summary()[source]

Formatted summary mirroring R’s summary.rdplot layout (rdplot.R:610-654), including the implied scale and WIMSE weights.

Return type:

str

plot(ax=None, title=None, xlabel=None, ylabel=None, dots_color=None, line_color=None, show_ci=None)[source]

Render the RD plot (binned means, global polynomial curves, cutoff line; error bars when CIs were requested at fit time or forced via show_ci=True). Requires matplotlib (an optional dependency).

Returns the matplotlib axes. Rendering is cosmetic - the parity surface of this port is the NUMBERS on the result object.

Parameters:
  • ax (Any | None)

  • title (str | None)

  • xlabel (str | None)

  • ylabel (str | None)

  • dots_color (str | None)

  • line_color (str | None)

  • show_ci (bool | None)

Return type:

Any

__init__(coef, vars_bins, vars_poly, J, J_IMSE, J_MV, scale, rscale, bin_avg, bin_med, p, cutoff, h, N, N_h, binselect, kernel_type, ci_level, ci_requested, covariate_coefficients=None, covariates_dropped=None)
Parameters:
Return type:

None