statspai.rd

rd

Regression Discontinuity (RD) module for StatsPAI.

Regression-discontinuity tools exposed through the StatsPAI API:

Core estimation: - Sharp, Fuzzy, and Kink RD estimation with robust bias-corrected inference (CCT 2014) - Covariate-adjusted local polynomial estimation (Calonico et al. 2019) - Donut-hole RD for manipulation near cutoff - Regression Discontinuity in Time (Hausman-Rapson 2018) - Multi-cutoff and multi-score RD (Cattaneo et al. 2024) - Boundary discontinuity / 2D RD designs (Cattaneo-Titiunik-Yu 2025)

Bandwidth selection: - MSE-optimal: mserd, msetwo, msecomb1, msecomb2 - CER-optimal: cerrd, certwo, cercomb1, cercomb2 (Calonico-Cattaneo-Farrell 2020) - Fuzzy-specific and covariate-adjusted bandwidth selection

Inference: - Honest confidence intervals (Armstrong-Kolesar 2018, 2020) - Local randomization inference with Fisher exact tests (Cattaneo-Titiunik-VB 2016) - Window selection and sensitivity analysis - Rosenbaum sensitivity bounds

Treatment effect heterogeneity: - CATE estimation via fully interacted local linear (Calonico et al. 2025) - ML + RD: Causal Forest, Gradient Boosting, LASSO-assisted RD

External validity: - Extrapolation away from cutoff (Angrist-Rokkanen 2015) - Multi-cutoff extrapolation (Cattaneo et al. 2024) - External validity diagnostics

Diagnostics & visualization: - Diagnostic dashboard (rdsummary) - IMSE-optimal binned scatter with pointwise CI bands - Density manipulation testing (CJM 2020) - Bandwidth sensitivity, covariate balance, placebo cutoff tests - Power analysis and sample size calculations

Examples:

sp.rd is a callable dispatcher: the default path is the CCT robust bias-corrected estimator (method="rdrobust"), and every other RD variant is one method= away.

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 1000
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.8 * (x >= 0) + 0.5 * x + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({"x": x, "y": y})
>>> res = sp.rd(df, y="y", x="x", c=0)          # CCT robust bias-corrected
>>> type(res).__name__
'CausalResult'
>>> round(float(res.estimate), 3)
0.738
>>> honest = sp.rd(df, y="y", x="x", c=0, method="honest")
>>> round(float(honest.estimate), 3)            # Armstrong-Kolesar honest CI
0.745

BayesRDHTEResult dataclass

Bases: ResultProtocolMixin

Bayesian RDD with HTE posterior summaries.

Returned by :func:rd_bayes_hte. posterior_mean / posterior_sd / posterior_ci summarise the average treatment effect at the cutoff, while cate / cate_sd give the per-unit covariate-conditional effect posteriors.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> run = rng.uniform(-1, 1, size=n)
>>> x = rng.normal(size=n)
>>> treat = (run >= 0).astype(float)
>>> y = (1.0 + 0.5 * run + (0.4 + 0.3 * x) * treat
...      + rng.normal(scale=0.2, size=n))
>>> df = pd.DataFrame({'y': y, 'run': run, 'x': x})
>>> res = sp.rd_bayes_hte(df, y='y', running='run', covariates=['x'],
...                       cutoff=0.0, n_draws=500, seed=0)
>>> res.cate.shape == (res.n_obs,)
True
>>> bool(res.posterior_sd >= 0)
True

DistRDResult dataclass

Bases: ResultProtocolMixin

RDD effect on the distribution of the outcome at the cutoff.

Attributes:

Name	Type	Description
quantiles	ndarray	Quantile levels at which the effects are reported.
qte	ndarray	Quantile treatment effect at the cutoff: the q-quantile of the outcome just right of the cutoff minus the one just left of it, in units of the outcome.
se	ndarray	Bootstrap standard error of each qte.
bandwidth	float	Bandwidth used in the local-linear fits.
n_obs	int	Number of complete observations.
ci_lower, ci_upper	ndarray	Percentile-bootstrap interval for each qte.
cdf_effect	ndarray	Jump at the cutoff in P(Y <= y_q), where y_q is the pooled q-quantile of the outcome (y_at_quantile). A probability, and negative where the treatment raises the outcome.
cdf_se	ndarray	Standard error of each cdf_effect.
y_at_quantile	ndarray	The thresholds y_q behind cdf_effect.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> x = rng.uniform(-1, 1, n)
>>> y = 1.0 + 0.5 * x + (x >= 0) * 0.8 + rng.normal(0, 1, n)
>>> data = pd.DataFrame({"y": y, "x": x})
>>> res = sp.rd_distribution(data, y="y", running="x")
>>> bool(res.qte.shape == res.quantiles.shape)
True

DDDResult dataclass

Bases: ResultProtocolMixin

Distributional discontinuity design output.

Attributes:

Name	Type	Description
quantiles	ndarray	Quantiles of y at which the effects are evaluated.
rdd_effect	ndarray	Level jump in the conditional CDF at the cutoff (RDD).
rkd_effect	ndarray	Slope jump in the conditional CDF at the cutoff (RKD).
bandwidth	float	Bandwidth used in the local-linear fit.
n_obs	int	Number of complete observations.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> x = rng.uniform(-1, 1, n)
>>> y = 1.0 + 0.5 * x + (x >= 0) * 0.8 + rng.normal(0, 1, n)
>>> data = pd.DataFrame({"y": y, "x": x})
>>> res = sp.rd_distributional_design(data, y="y", running="x")
>>> bool(res.rdd_effect.shape == res.quantiles.shape)
True

RDInterferenceResult dataclass

Bases: ResultProtocolMixin

Direct + spillover RDD effects under network interference.

Returned by :func:rd_interference. Carries the boundary direct effect and the neighbour-boundary spillover effect, each with a local-linear standard error, plus the sample size and bandwidth.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng(1)
>>> n = 600
>>> r = rng.uniform(-1, 1, size=n)
>>> rn = rng.uniform(-1, 1, size=n)
>>> y = (0.5 * (r >= 0) + 0.3 * (rn >= 0)
...      + 0.2 * r + 0.1 * rn + rng.normal(scale=0.3, size=n))
>>> df = pd.DataFrame({"y": y, "run": r, "nbr_run": rn})
>>> res = sp.rd_interference(df, y="y", running="run",
...                          neighbour_running="nbr_run")
>>> isinstance(res, sp.RDInterferenceResult)
True
>>> isinstance(res.summary(), str)
True

MultiScoreRDResult dataclass

Bases: ResultProtocolMixin

Multi-score RDD effect on the eligibility boundary.

Attributes:

Name	Type	Description
boundary_effect	float	Local average treatment effect on the boundary manifold.
se	float	Standard error of boundary_effect.
n_obs	int	Number of observations used (after dropping missing rows).
boundary_share	float	Fraction of observations inside the kernel window.
bandwidth	float	Bandwidth applied to the (max-) distance to the boundary.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 2000
>>> m = rng.uniform(400, 800, n)
>>> v = rng.uniform(400, 800, n)
>>> treat = ((m >= 600) & (v >= 600)).astype(int)
>>> y = 10.0 + 0.01 * m + 0.01 * v + 5.0 * treat + rng.normal(0, 2, n)
>>> df = pd.DataFrame({"y": y, "math": m, "verbal": v})
>>> res = sp.rd_multi_score(
...     df, y="y", running_vars=["math", "verbal"], cutoffs=[600.0, 600.0])
>>> isinstance(res, sp.MultiScoreRDResult)
True
>>> res.n_obs
2000
>>> bool(res.se >= 0.0)
True

RDMultiResult

Bases: ResultProtocolMixin

Results from multi-cutoff/multi-score RD.

Returned by :func:rdmc. Holds per-cutoff estimates in cutoff_results plus the pooled estimate / SE / CI, and exposes :meth:summary and a forest :meth:plot.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 900
>>> running_var = rng.uniform(0, 100, n)
>>> cutoffs = [30, 60]
>>> nearest = np.array([min(cutoffs, key=lambda c: abs(rv - c))
...                     for rv in running_var])
>>> treated = (running_var >= nearest).astype(float)
>>> score = 1.0 + 0.02 * running_var + 2.0 * treated + rng.normal(0, 1, n)
>>> df = pd.DataFrame({"score": score, "running_var": running_var})
>>> result = sp.rdmc(df, y="score", x="running_var", cutoffs=[30, 60])
>>> isinstance(result, sp.RDMultiResult)
True
>>> result.n_cutoffs
2

plot

plot(ax: Optional[Any] = None, **kwargs: Any) -> Any

Forest plot of cutoff-specific and pooled estimates.

RDPowerResult dataclass

Bases: ResultProtocolMixin

Power of an RD design against one alternative.

summary

summary() -> str

Plain-text table of the calculation.

RDSampSiResult dataclass

Bases: ResultProtocolMixin

Sample size needed to reach a target power.

summary

summary() -> str

One-line statement of the required sample.

boundary_rd

boundary_rd(*args: Any, **kwargs: Any) -> CausalResult

User-friendly alias for :func:sp.rd2d (boundary discontinuity design).

Cattaneo, Titiunik & Yu (2025) boundary discontinuity design for 2D running variables (lat/long, for example).

See :func:statspai.rd.rd2d.rd2d for the full signature.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 1500
>>> x1 = rng.uniform(-1, 1, n)
>>> x2 = rng.uniform(-1, 1, n)
>>> treat = (x1 >= 0).astype(int)
>>> y = 0.7 * treat + 0.4 * x1 + 0.2 * x2 + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({"x1": x1, "x2": x2, "treat": treat, "y": y})
>>> res = sp.boundary_rd(df, y="y", x1="x1", x2="x2", treatment="treat")
>>> type(res).__name__
'CausalResult'
>>> bool(hasattr(res, "estimate"))
True

References

cattaneo2025boundary

geographic_rd

geographic_rd(*args: Any, **kwargs: Any) -> CausalResult

User-friendly alias for :func:sp.rdms (multi-score RD).

Geographic RD is the most common multi-score special case: the running score is a signed distance to a political boundary. Cattaneo et al. (2024) dispatch it via rdms; this alias makes the intent explicit.

See :func:statspai.rd.rdmulti.rdms for the full signature.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 2000
>>> x1 = rng.uniform(-1, 1, n)
>>> x2 = rng.uniform(-1, 1, n)
>>> treat = ((x1 >= 0) & (x2 >= 0)).astype(float)
>>> y = 0.8 * treat + 0.5 * x1 + 0.3 * x2 + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({"x1": x1, "x2": x2, "y": y})
>>> res = sp.geographic_rd(df, y="y", x1="x1", x2="x2", bandwidth=0.5)
>>> type(res).__name__
'CausalResult'
>>> bool(hasattr(res, "estimate"))
True

References

keele2015geographic

multi_cutoff_rd

multi_cutoff_rd(*args: Any, **kwargs: Any) -> RDMultiResult

User-friendly alias for :func:sp.rdmc (multi-cutoff RD).

See :func:statspai.rd.rdmulti.rdmc for full documentation.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 600
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.8 * (x >= 0) + 0.5 * x + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({"x": x, "y": y})
>>> res = sp.multi_cutoff_rd(
...     df, y="y", x="x", cutoffs=[-0.4, 0.0, 0.4]
... )
>>> res.n_cutoffs
3
>>> len(res.cutoff_results)
3
>>> round(float(res.pooled_estimate), 3)
0.258

multi_score_rd

multi_score_rd(*args: Any, **kwargs: Any) -> MultiScoreRDResult

User-friendly alias for :func:sp.rd_multi_score.

Multi-score RD when eligibility depends on more than one discontinuous rule (e.g. income AND age thresholds). Separate from boundary RDD in that the rules are axis-aligned.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 800
>>> x1 = rng.uniform(-1, 1, n)
>>> x2 = rng.uniform(-1, 1, n)
>>> treat = (x1 >= 0) & (x2 >= 0)
>>> y = (0.8 * treat + 0.5 * x1 + 0.3 * x2
...      + rng.normal(0, 0.3, n))
>>> df = pd.DataFrame({"x1": x1, "x2": x2, "y": y})
>>> res = sp.multi_score_rd(
...     df, y="y", running_vars=["x1", "x2"], cutoffs=[0.0, 0.0]
... )
>>> res.n_obs
800
>>> round(float(res.boundary_effect), 3)  # true jump is 0.8
0.76

rdbwselect

rdbwselect(data: DataFrame, y: str, x: str, c: float = 0, fuzzy: Optional[str] = None, deriv: int = 0, p: int = 1, q: Optional[int] = None, covs: Optional[List[str]] = None, kernel: str = 'triangular', bwselect: str = 'mserd', cluster: Optional[str] = None, vce: str = 'nn', all: bool = False, masspoints: str = 'adjust') -> DataFrame

Bandwidth selection for local polynomial RD estimation.

Runs the Calonico-Cattaneo-Titiunik three-stage cascade and returns all ten MSE-optimal and CER-optimal selectors rdrobust offers. MSE-optimal bandwidths minimize the mean squared error of the RD point estimator, while CER-optimal bandwidths minimize the coverage error rate of robust bias-corrected confidence intervals.

This is the same code path sp.rdrobust uses to pick its own default bandwidth, so a value read here and passed back as h= reproduces sp.rdrobust's default fit exactly. Pinned against rdrobust::rdbwselect (R) and the rdbwselect ado (Stata) by Track A module 88_rdbwselect to 1.8e-12 and 3.7e-9 respectively, across all ten selectors, polynomial orders 1-3, three kernels, covariate adjustment, clustering and the regression-kink derivative.

.. versionchanged:: 1.27.0 Both properties above are new. This function previously ran a separate single-step rule of thumb and returned bandwidths 2.8x-4.8x too narrow, and the four comb selectors resolved to the plain rd cascade. See MIGRATION.md — results should be recomputed.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable column name.	required
x	str	Running variable column name.	required
c	float	RD cutoff value.	0
fuzzy	str	Treatment variable name for fuzzy RD. The MSE being minimised is then the one for the Wald ratio rather than for the reduced form, so every stage of the cascade changes. .. versionchanged:: 1.27.0 This argument was parsed and then discarded: it was never passed to the bandwidth cascade, so fuzzy= returned the sharp bandwidth while the docstring claimed otherwise. On a two-sided noncompliance replica of the Lee 2008 senate data that is a 9% to 16% error in h. Designs with one-sided noncompliance were unaffected, because rdbwselect itself falls back to the sharp bandwidth there (R's perf_comp) -- which is also why the defect survived: the fixture in the repository was one-sided.	None
deriv	int	Derivative order. 0 = standard RD (jump in level), 1 = regression kink design (change in slope).	0
p	int	Polynomial order for point estimation (1 = local linear).	1
q	int	Polynomial order for bias correction. Default is p + 1.	None
covs	list of str	Covariate column names. When provided, the variance estimates used in bandwidth selection account for covariate adjustment, typically yielding narrower bandwidths.	None
kernel	str	Kernel function: 'triangular', 'uniform', or 'epanechnikov'.	'triangular'
bwselect	str	Bandwidth selection method. One of: 'mserd' : MSE-optimal common bandwidth (default) 'msetwo' : MSE-optimal separate left/right bandwidths 'msesum' : MSE-optimal for the sum of the two intercepts 'msecomb1' : min(mserd, msesum), per side 'msecomb2' : median(msetwo, mserd, msesum), per side 'cerrd', 'certwo', 'cersum', 'cercomb1', 'cercomb2' : the CER-optimal counterparts of the above The combination rules are applied to the finished h and b of each cascade, element-wise per side -- not stage by stage.	'mserd'
cluster	str	Cluster variable name for cluster-robust variance estimation.	None
vce	(nn, hc0, hc1, hc2, hc3, cr1, cr2, cr3)	Variance estimator inside the bandwidth cascade, R rdbwselect(vce=): 'nn' uses nearest-neighbour residuals, the hc* family regression residuals. With cluster=, 'nn' / 'hc0' / 'hc1' become R's 'cr1' and 'hc2' / 'hc3' R 4.0.0's 'cr2' / 'cr3'; 'cr*' without cluster= fall back to 'hc*' (with a warning). The bandwidths equal sp.rdrobust's with the same vce. Stata's spellings 'cluster v' and 'nncluster v' (also inside 'vce(...)') are read as cluster='v' with that variance.	'nn'
all	bool	If True, compute and return all ten bandwidth types.	False
masspoints	(adjust, check, off)	R rdbwselect(masspoints=). When 20% or more of either side's observations are tied, 'adjust' uses the count of unique running-variable values in the reference bandwidth and floors the pilot bandwidth at the 10th unique value from the cutoff; 'check' only warns; 'off' treats the data as continuous.	'adjust'

Returns:

Type	Description
DataFrame	DataFrame with columns [method, h_left, h_right, b_left, b_right, n_left, n_right]. When all=False, contains a single row for the selected method. When all=True, contains ten rows, one per method.

Notes

The CER-optimal bandwidth is related to the MSE-optimal bandwidth by:

h_CER = h_MSE * n^{-1/((2p+3)(2p+5))}

For p=1 (local linear), this gives h_CER ~ h_MSE * n^{-1/35}, which is strictly narrower than h_MSE. The narrower bandwidth yields confidence intervals with better coverage properties at the cost of slightly wider intervals.

The MSE-optimal bandwidth has rate n^{-1/(2p+3)} while the CER rate is n^{-1/(2p+3) - 1/((2p+3)(2p+5))}. For large samples the difference is meaningful: CER bandwidths produce robust CIs that achieve their nominal coverage rate, whereas MSE bandwidths can exhibit substantial coverage distortion.

References

Calonico, S., Cattaneo, M.D. and Farrell, M.H. (2020). "Optimal Bandwidth Choice for Robust Bias-Corrected Inference in Regression Discontinuity Designs." Econometrics Journal, 23(2), 192-210. [@calonico2020optimal]

Calonico, S., Cattaneo, M.D. and Titiunik, R. (2014). "Robust Nonparametric Confidence Intervals for Regression-Discontinuity Designs." Econometrica, 82(6), 2295-2326. [@calonico2014robust]

Examples:

Basic MSE-optimal bandwidth:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> n = 2000
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.5 * x + 3.0 * (x >= 0) + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({'y': y, 'x': x})
>>> bw = sp.rdbwselect(df, y='y', x='x', c=0)
>>> bool('h_left' in bw.columns)  # DataFrame with bandwidth info
True

Compare all eight bandwidth methods:

>>> bw_all = sp.rdbwselect(df, y='y', x='x', c=0, all=True)
>>> bool(len(bw_all) >= 1)
True

CER-optimal bandwidth for better coverage:

>>> bw_cer = sp.rdbwselect(df, y='y', x='x', c=0, bwselect='cerrd')
>>> bool(bw_cer is not None)
True

Fuzzy RD with covariates:

>>> treated = ((x >= 0) & (rng.uniform(0, 1, n) < 0.85)).astype(float)
>>> df['treated'] = treated
>>> df['age'] = rng.normal(40, 10, n)
>>> df['gender'] = rng.integers(0, 2, n).astype(float)
>>> bw_fuzzy = sp.rdbwselect(df, y='y', x='x', c=0,
...                          fuzzy='treated', covs=['age', 'gender'])
>>> bool(bw_fuzzy is not None)
True

rd_bayes_hte

rd_bayes_hte(data: DataFrame, y: str, running: str, covariates: List[str], cutoff: float = 0.0, bandwidth: Optional[float] = None, kernel: str = 'triangular', alpha: float = 0.05, n_draws: int = 2000, seed: int = 0) -> BayesRDHTEResult

Bayesian RDD allowing CATE to depend on covariates.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
running	str		required
covariates	list of str		required
cutoff	float		0.0
bandwidth	float		None
kernel	str		'triangular'
alpha	float		0.05
n_draws	int		2000
seed	int		0

Returns:

Type	Description
BayesRDHTEResult

Examples:

Sharp RD with a covariate-modulated jump at the cutoff. Without PyMC installed this falls back to the conjugate closed-form posterior:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> run = rng.uniform(-1, 1, size=n)
>>> x = rng.normal(size=n)
>>> treat = (run >= 0).astype(float)
>>> y = (1.0 + 0.5 * run + (0.4 + 0.3 * x) * treat
...      + rng.normal(scale=0.2, size=n))
>>> df = pd.DataFrame({'y': y, 'run': run, 'x': x})
>>> res = sp.rd_bayes_hte(df, y='y', running='run', covariates=['x'],
...                       cutoff=0.0, n_draws=500, seed=0)
>>> res.n_obs
400
>>> len(res.cate)
400
>>> bool(np.isfinite(res.posterior_mean))
True
>>> lo, hi = res.posterior_ci
>>> bool(lo <= hi)
True
>>> print(res.summary())

rd_bias_aware_fuzzy

rd_bias_aware_fuzzy(data: DataFrame, y: str, x: str, fuzzy: str, c: float = 0.0, M_y: Optional[float] = None, M_d: Optional[float] = None, h: Optional[float] = None, kernel: str = 'triangular', alpha: float = 0.05, cluster: Optional[str] = None, n_grid: int = 401) -> CausalResult

Bias-aware confidence interval for fuzzy RD (Noack & Rothe 2024).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome.	required
x	str	Running variable; units with x >= c are assigned to treatment.	required
fuzzy	str	Treatment received.	required
c	float	RD cutoff.	0
M_y	float	Bounds on the second derivatives of the outcome and treatment means on either side. When omitted, the rule of thumb of Armstrong and Kolesar (global quartic on each side) is applied to each.	None
M_d	float	Bounds on the second derivatives of the outcome and treatment means on either side. When omitted, the rule of thumb of Armstrong and Kolesar (global quartic on each side) is applied to each.	None
h	float	Bandwidth. When omitted, the bandwidth minimising the worst-case MSE of Delta_Y_hat - tau_h Delta_D_hat (tau_h the Wald ratio at that bandwidth), with worst-case bias M_y B_1(h) (an initial effect of zero) and a preliminary covariance of (Y, D) that is homoskedastic on each side.	None
kernel	('triangular', 'uniform', 'epanechnikov')	Kernel of the local linear fits.	'triangular'
alpha	float	One minus the confidence level of the bias-aware set.	0.05
cluster	str	Cluster identifier. Standard errors become cluster-robust sums of the weighted local linear residuals, and the bandwidth search adds a common within-cluster covariance.	None
n_grid	int	Points of the diagnostic rejection grid returned in model_info['bias_aware']['rejection_grid']. The endpoints of the confidence set are solved exactly, not read off the grid.	401

Returns:

Type	Description
CausalResult	estimate is the Wald ratio, se its delta-method standard error, ci the Anderson-Rubin bias-aware set (with infinite ends when it is unbounded). model_info['bias_aware'] carries the jumps, the 2x2 variance of (Delta_Y_hat, Delta_D_hat), the bias per unit curvature, the first-stage F statistic, whether the set is a single interval, and the linearised interval under 'rdhonest'.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(-1, 1, 1500)
>>> d = (rng.uniform(size=1500) < 0.2 + 0.5 * (x >= 0)).astype(float)
>>> yy = 1 + 2 * d + 0.5 * x + rng.normal(0, 0.7, 1500)
>>> df = pd.DataFrame({"y": yy, "d": d, "x": x})
>>> r = sp.rd_bias_aware_fuzzy(df, y="y", x="x", fuzzy="d", M_y=2, M_d=0.5)

References

noack2024biasaware, armstrong2020simple

rd_compare

rd_compare(data: DataFrame, y: str, x: str, c: float = 0.0, methods: Sequence[str] = _DEFAULT_COMPARE_METHODS, fuzzy: Optional[str] = None, method_kwargs: Optional[Dict[str, Dict[str, Any]]] = None, alpha: float = 0.05) -> DataFrame

Compare multiple RD estimators on the same data.

Returns a tidy DataFrame: one row per method with point estimate, SE, p-value and CI. Useful for robustness tables across estimation families (local-polynomial vs honest vs local-randomisation vs flexible-covariate-adjusted).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
x	str		required
c	float		0.0
methods	sequence of str	Method aliases recognised by the :data:sp.rd._RD_METHOD_ALIASES dispatcher. Default: ('rdrobust', 'honest'). Add 'randinf' together with its window, e.g. method_kwargs={'randinf': {'wl': -0.1, 'wr': 0.1}}.	_DEFAULT_COMPARE_METHODS
fuzzy	str	Fuzzy treatment column passed through to all methods that accept it.	None
method_kwargs	dict of dict	Per-method extra kwargs, e.g. {'rdrobust': {'kernel': 'uniform'}}.	None
alpha	float	Confidence level forwarded to each estimator that accepts alpha.	0.05

Returns:

Type	Description
DataFrame	Columns: method, estimate, se, pvalue, ci_lower, ci_upper, n_obs, status.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 500
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.8 * (x >= 0) + 0.5 * x + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({"x": x, "y": y})
>>> tab = sp.rd_compare(df, y="y", x="x", c=0.0,
...                     methods=["rdrobust"])
>>> tab.shape
(1, 8)
>>> tab.loc[0, "method"], tab.loc[0, "status"]
('rdrobust', 'ok')
>>> round(float(tab.loc[0, "estimate"]), 3)
0.818

rd_dashboard

rd_dashboard(data: DataFrame, y: str, x: str, c: float = 0.0, covs: Optional[List[str]] = None, fuzzy: Optional[str] = None, bw_grid: Optional[Sequence[float]] = None, h: Optional[float] = None, figsize: Tuple[float, float] = (12, 9), title: Optional[str] = None, save: Optional[str] = None) -> Tuple[Any, Any]

Plot a four-panel RD diagnostic dashboard.

Combines the four checks every RD analysis should report:

1. RD plot (top-left): IMSE-binned scatter with polynomial fit.
2. Density discontinuity (top-right): :func:rdplotdensity output for a McCrary-style manipulation test.
3. Covariate balance (bottom-left): mean of each cov on each side of the cutoff with point ranges.
4. Bandwidth sensitivity (bottom-right): point estimate and CI across a grid of bandwidths.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
x	str		required
c	float	Cutoff.	0.0
covs	list of str	Covariates for the balance panel. If None, the panel shows the density's binomial near-cutoff test instead.	None
fuzzy	str	Fuzzy treatment column (passed through to RD plot/sensitivity).	None
bw_grid	sequence of float	Bandwidths to evaluate in the sensitivity panel. If None, uses [0.5, 0.75, 1.0, 1.25, 1.5, 2.0] × h_mse.	None
h	float	Reference bandwidth used for plotting and as the basis for bw_grid. If None, MSE-optimal bandwidth from rdrobust.	None
figsize	tuple		(12, 9)
title	str	Suptitle.	None
save	str	If a path is given, also save the figure (extension determines format).	None

Returns:

Type	Description
(fig, axes) where axes is a 2x2 numpy array.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 500
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.8 * (x >= 0) + 0.5 * x + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({"x": x, "y": y})
>>> fig, axes = sp.rd_dashboard(df, y="y", x="x", c=0.0)
>>> axes.shape
(2, 2)
>>> fig.savefig("rd_dashboard.png")

rd_robustness_table

rd_robustness_table(data: DataFrame, y: str, x: str, c: float = 0.0, fuzzy: Optional[str] = None, kernels: Sequence[str] = ('triangular', 'epanechnikov', 'uniform'), bwselects: Sequence[str] = ('mserd', 'cerrd', 'msetwo'), polynomials: Sequence[int] = (1, 2), donuts: Sequence[float] = (0.0,), covs: Optional[List[str]] = None, cluster: Optional[str] = None, vce: str = 'nn', alpha: float = 0.05) -> DataFrame

Sweep over (kernel, bwselect, polynomial, donut) and return a robustness table for reporting.

Each row is one specification. The DataFrame contains both "Conventional" and "Robust" point estimates and CIs from :func:rdrobust, plus the bandwidth used. Suitable for direct .to_latex(...) / .to_excel(...) export, or feed into :func:statspai.outreg2 for native multi-column tables.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
x	str		required
c	float		0.0
fuzzy	optional, forwarded to rdrobust		None
covs	optional, forwarded to rdrobust		None
cluster	optional, forwarded to rdrobust		None
vce	('nn', 'hc0', 'hc1', 'hc2', 'hc3')	Variance estimator of every rdrobust fit, as in sp.rdrobust (it also enters the bandwidth selection). With cluster=, 'nn' / 'hc0' / 'hc1' become R's 'cr1'.	'nn'
kernels	sequences	Specification grid.	('triangular', 'epanechnikov', 'uniform')
bwselects	sequences	Specification grid.	('triangular', 'epanechnikov', 'uniform')
polynomials	sequences	Specification grid.	('triangular', 'epanechnikov', 'uniform')
donuts	sequences	Specification grid.	('triangular', 'epanechnikov', 'uniform')
alpha	float		0.05

Returns:

Type	Description
pd.DataFrame with one row per specification and columns:	kernel, bwselect, p, donut, h, b, estimate_conv, se_conv, ci_conv_lo, ci_conv_hi, pvalue_conv, estimate_rbc, se_rbc, ci_rbc_lo, ci_rbc_hi, pvalue_rbc, n_left, n_right, status

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 500
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.8 * (x >= 0) + 0.5 * x + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({"x": x, "y": y})
>>> tab = sp.rd_robustness_table(
...     df, y="y", x="x", c=0.0,
...     kernels=["triangular", "uniform"],
...     bwselects=["mserd"], polynomials=[1],
... )
>>> tab.shape
(2, 19)
>>> tab["kernel"].tolist()
['triangular', 'uniform']
>>> round(float(tab.loc[0, "estimate_conv"]), 3)
0.806

rdbalance

rdbalance(data: DataFrame, x: str, c: float = 0, covs: Optional[List[str]] = None, p: int = 1, kernel: str = 'triangular', h: Optional[float] = None, alpha: float = 0.05, cluster: Optional[str] = None, vce: str = 'nn') -> DataFrame

Covariate balance test at the RD cutoff.

For each covariate, estimates the discontinuity at the cutoff using local polynomial regression. Pre-treatment covariates should show no significant jump at the cutoff if the RD design is valid.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
x	str	Running variable name.	required
c	float	Cutoff.	0
covs	list of str	Covariate names to test. If None, tests all numeric columns except x.	None
p	int		1
kernel	str		'triangular'
h	float	Manual bandwidth. If None, uses MSE-optimal per covariate.	None
alpha	float		0.05
cluster	str	Cluster column for the standard errors of every rdrobust fit (vce(cluster); also the bandwidth selection, as rdrobust does).	None
vce	(nn, hc0, hc1, hc2, hc3)	Variance estimator of every rdrobust fit, as in sp.rdrobust (it also enters the bandwidth selection). With cluster=, 'nn' / 'hc0' / 'hc1' become R's 'cr1'.	'nn'

Returns:

Type	Description
DataFrame	Columns: covariate, estimate, se, z, pvalue, significant.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> x = rng.uniform(-1, 1, 500)
>>> y = 0.5 * (x >= 0) + 0.8 * x + rng.normal(0, 0.3, 500)
>>> w1 = 0.3 * x + rng.normal(0, 0.5, 500)
>>> w2 = rng.normal(0, 1, 500)
>>> df = pd.DataFrame({"x": x, "y": y, "w1": w1, "w2": w2})
>>> bal = sp.rdbalance(df, x="x", covs=["w1", "w2"])
>>> bal["covariate"].tolist()
['w1', 'w2']
>>> list(bal.columns)
['covariate', 'estimate', 'se', 'z', 'pvalue', 'significant']
>>> bool((bal["pvalue"] > 0.01).all())    # no covariate jump at the 1% level
True

rdbwsensitivity

rdbwsensitivity(data: DataFrame, y: str, x: str, c: float = 0, fuzzy: Optional[str] = None, p: int = 1, kernel: str = 'triangular', bw_grid: Optional[List[float]] = None, n_grid: int = 15, bw_range: Tuple[float, float] = (0.5, 2.0), alpha: float = 0.05, cluster: Optional[str] = None, vce: str = 'nn', ax: Optional[Any] = None, figsize: Tuple[float, float] = (10, 6)) -> DataFrame

Bandwidth sensitivity analysis for RD estimates.

Re-estimates the RD effect across a grid of bandwidths to assess robustness. Plots point estimates with confidence intervals.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome and running variable names.	required
x	str	Outcome and running variable names.	required
c	float	Cutoff.	0
fuzzy	str	Treatment variable for fuzzy RD.	None
p	int	Polynomial order.	1
kernel	str		'triangular'
bw_grid	list of float	Explicit bandwidth values to evaluate. If None, auto-generates a grid as multiples of the MSE-optimal bandwidth.	None
n_grid	int	Number of grid points if bw_grid is None.	15
bw_range	tuple	Range of multipliers for the optimal bandwidth.	(0.5, 2.0)
alpha	float		0.05
ax	matplotlib Axes		None
figsize	tuple		(10, 6)
cluster	str	Cluster column for the standard errors of every rdrobust fit (vce(cluster); also the bandwidth selection, as rdrobust does).	None
vce	(nn, hc0, hc1, hc2, hc3)	Variance estimator of every rdrobust fit, as in sp.rdrobust (it also enters the bandwidth selection). With cluster=, 'nn' / 'hc0' / 'hc1' become R's 'cr1'.	'nn'

Returns:

Type	Description
DataFrame	Columns: bandwidth, estimate, se, ci_lower, ci_upper, pvalue.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 500
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.8 * (x >= 0) + 0.5 * x + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({"x": x, "y": y})
>>> tab = sp.rdbwsensitivity(df, y="y", x="x", c=0.0, n_grid=5)
>>> tab.shape
(5, 6)
>>> list(tab.columns)
['bandwidth', 'estimate', 'se', 'ci_lower', 'ci_upper', 'pvalue']
>>> [round(v, 3) for v in tab["estimate"].tolist()]
[0.805, 0.863, 0.825, 0.825, 0.819]

rdplacebo

rdplacebo(data: DataFrame, y: str, x: str, c: float = 0, placebo_cutoffs: Optional[List[float]] = None, n_placebo: int = 10, side: str = 'both', fuzzy: Optional[str] = None, p: int = 1, kernel: str = 'triangular', alpha: float = 0.05, cluster: Optional[str] = None, vce: str = 'nn', ax: Optional[Any] = None, figsize: Tuple[float, float] = (10, 6)) -> DataFrame

Placebo cutoff test for RD validity.

Estimates the RD effect at fake cutoff points where no treatment effect should exist. Significant effects at placebo cutoffs suggest the RD design may be invalid.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome and running variable names.	required
x	str	Outcome and running variable names.	required
c	float	True cutoff.	0
placebo_cutoffs	list of float	Explicit placebo cutoff values. If None, auto-generates from the data distribution.	None
n_placebo	int	Number of placebo cutoffs if auto-generating.	10
side	str	Which side of the true cutoff to place placebos: 'left', 'right', or 'both'.	'both'
fuzzy	str	Treatment variable for fuzzy RD.	None
p	int		1
kernel	str		'triangular'
alpha	float		0.05
ax	matplotlib Axes		None
figsize	tuple		(10, 6)
cluster	str	Cluster column for the standard errors of every rdrobust fit (vce(cluster); also the bandwidth selection, as rdrobust does).	None
vce	(nn, hc0, hc1, hc2, hc3)	Variance estimator of every rdrobust fit, as in sp.rdrobust (it also enters the bandwidth selection). With cluster=, 'nn' / 'hc0' / 'hc1' become R's 'cr1'.	'nn'

Returns:

Type	Description
DataFrame	Columns: cutoff, estimate, se, ci_lower, ci_upper, pvalue, is_true_cutoff.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> x = rng.uniform(-1, 1, 500)
>>> y = 0.5 * (x >= 0) + 0.8 * x + rng.normal(0, 0.3, 500)
>>> df = pd.DataFrame({"x": x, "y": y})
>>> pl = sp.rdplacebo(df, y="y", x="x", n_placebo=4)
>>> true_row = pl[pl["is_true_cutoff"]]
>>> round(float(true_row["estimate"].iloc[0]), 2)
0.52
>>> int((pl[~pl["is_true_cutoff"]]["pvalue"] < 0.05).sum())
0

rdsummary

rdsummary(data: DataFrame, y: str, x: str, c: float = 0, fuzzy: Optional[str] = None, covs: Optional[List[str]] = None, p: int = 1, kernel: str = 'triangular', alpha: float = 0.05, cluster: Optional[str] = None, vce: str = 'nn', verbose: bool = True, plot: bool = False, full: bool = False) -> Dict[str, Any]

One-stop RD diagnostic battery.

Runs the main RD estimate plus all standard validation checks in a single call, returning a structured summary.

Diagnostics run (standard): 1. Main RD estimate (conventional + robust) 2. Density manipulation test (CJM 2020) 3. Covariate balance at cutoff (if covs provided) 4. Bandwidth sensitivity

Extended diagnostics (full=True): 5. Honest CI (Armstrong-Kolesár 2020) 6. Power analysis (MDE at 80%) 7. Placebo cutoff tests (5 placebos per side) 8. Multiple bandwidth selectors comparison

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome and running variable.	required
x	str	Outcome and running variable.	required
c	float	Cutoff.	0
fuzzy	str	Treatment variable for fuzzy RD.	None
covs	list of str	Pre-treatment covariates for balance test.	None
p	int		1
kernel	str		'triangular'
alpha	float		0.05
verbose	bool	Print formatted summary to console.	True
plot	bool	Generate a multi-panel diagnostic plot.	False
full	bool	Run extended diagnostics (honest CI, power, placebos).	False
cluster	str	Cluster column for the standard errors of every rdrobust fit (vce(cluster); also the bandwidth selection, as rdrobust does).	None
vce	(nn, hc0, hc1, hc2, hc3)	Variance estimator of every rdrobust fit, as in sp.rdrobust (it also enters the bandwidth selection). With cluster=, 'nn' / 'hc0' / 'hc1' become R's 'cr1'.	'nn'

Returns:

Type	Description
dict with keys:	'estimate': CausalResult from rdrobust 'density_test': CausalResult from rddensity 'balance': pd.DataFrame (if covs given) 'bw_sensitivity': pd.DataFrame 'honest_ci': CausalResult (if full=True) 'power': RDPowerResult (if full=True) 'placebos': pd.DataFrame (if full=True) 'bandwidth_comparison': pd.DataFrame (if full=True) 'figure': matplotlib Figure (if plot=True)

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> x = rng.uniform(-1, 1, 500)
>>> y = 0.5 * (x >= 0) + 0.8 * x + rng.normal(0, 0.3, 500)
>>> w1 = 0.3 * x + rng.normal(0, 0.5, 500)
>>> df = pd.DataFrame({"x": x, "y": y, "w1": w1})
>>> res = sp.rdsummary(df, y="y", x="x", covs=["w1"],
...                    verbose=False)
>>> sorted(res.keys())
['balance', 'bw_sensitivity', 'density_test', 'estimate']
>>> round(float(res["estimate"].estimate), 2)
0.52

rd_diff_in_disc

rd_diff_in_disc(data: DataFrame, y: str, x: str, post: str, c: float = 0.0, h: Optional[float] = None, p: int = 1, kernel: str = 'triangular', cluster: Optional[str] = None, covariates: Optional[List[str]] = None, alpha: float = 0.05) -> CausalResult

Difference-in-discontinuities estimate with a pooled local polynomial.

Parameters:

Name	Type	Description	Default
data	DataFrame	Pooled pre- and post-period observations.	required
y	str	Outcome and running variable.	required
x	str	Outcome and running variable.	required
post	str	0/1 period indicator (1 = the period in which the cutoff policy is active).	required
c	float	Cutoff; observations with x >= c are above it.	0
h	float	Bandwidth, common to both sides and periods. Default: the smaller of the two periods' MSE-optimal rdbwselect bandwidths (mserd) -- a conservative choice; pass h to set it.	None
p	int	Polynomial order on each side and period.	1
kernel	('triangular', 'uniform', 'epanechnikov')		'triangular'
cluster	str	Cluster the standard error (e.g. by unit or site, which lets the two periods' errors correlate): CR1 as Stata regress [aw=], vce(cluster c). Otherwise HC1 (vce(robust)).	None
covariates	list of str	Additive covariates in the pooled regression (covs= is accepted, as in rdrobust).	None
alpha	float		0.05

Returns:

Type	Description
CausalResult	estimate is the change in the discontinuity; detail has the pre- and post-period discontinuities; model_info the bandwidth and sample sizes.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 2000
>>> d = pd.DataFrame({"x": rng.uniform(-1, 1, n), "post": rng.integers(0, 2, n)})
>>> d["y"] = d.x + 0.3 * (d.x >= 0) + 0.5 * (d.x >= 0) * d.post + rng.normal(0, 0.3, n)
>>> r = sp.rd_diff_in_disc(d, y="y", x="x", post="post", h=0.5)
>>> round(r.estimate, 1)
0.5

References

Grembi, V., Nannicini, T. and Troiano, U. (2016). "Do Fiscal Rules Matter?" American Economic Journal: Applied Economics, 8(3), 1-30.

rd_distribution

rd_distribution(data: DataFrame, y: str, running: str, cutoff: float = 0.0, quantiles: Optional[ndarray] = None, bandwidth: Optional[float] = None, kernel: str = 'triangular', alpha: float = 0.05, n_boot: int = 200, seed: int = 0) -> DistRDResult

Distribution-valued sharp RDD.

Estimates the distribution of the outcome just left and just right of the cutoff by local-linear regression of 1{Y <= y} on the running variable, and reports the discontinuity two ways: as a quantile treatment effect (qte, the horizontal distance between the two distributions, in units of the outcome) and as the jump in the CDF at the pooled quantiles (cdf_effect, the vertical distance, a probability).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
running	str		required
cutoff	float		0.0
quantiles	array - like	Defaults to (0.1, 0.25, 0.5, 0.75, 0.9).	None
bandwidth	float	Defaults to the interquartile range of the running variable.	None
kernel	str		'triangular'
alpha	float	Level of the percentile-bootstrap interval for qte.	0.05
n_boot	int	Bootstrap replications for the qte standard errors.	200
seed	int		0

Returns:

Type	Description
DistRDResult	qte with bootstrap se / ci_lower / ci_upper, and cdf_effect with cdf_se, at each requested quantile.

Notes

Until the 2026-10 fix qte held what is now cdf_effect: a change in probability, negative where the treatment raises the outcome. On a design that shifts the outcome up by 2.0 it reported -0.48 at the median.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> x = rng.uniform(-1, 1, n)
>>> y = 1.0 + 0.5 * x + (x >= 0) * 0.8 + rng.normal(0, 1, n)
>>> data = pd.DataFrame({"y": y, "x": x})
>>> res = sp.rd_distribution(data, y="y", running="x", cutoff=0.0)
>>> res.n_obs
400
>>> bool(res.qte.shape == res.se.shape == res.quantiles.shape)
True
>>> bool(res.ci_lower[2] < 0.8 < res.ci_upper[2])  # a 0.8 shift
True
>>> bool((res.cdf_effect < 0).all())  # mass moves to higher outcomes
True

rd_distributional_design

rd_distributional_design(data: DataFrame, y: str, running: str, cutoff: float = 0.0, quantiles: Optional[ndarray] = None, bandwidth: Optional[float] = None, kernel: str = 'triangular') -> DDDResult

Joint RDD + RKD on the conditional distribution of Y.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
running	str		required
cutoff	float		0.0
quantiles	array - like		None
bandwidth	float		None
kernel	str		'triangular'

Returns:

Type	Description
DDDResult	Quantile-by-quantile level jump (rdd_effect) and slope jump (rkd_effect) in the conditional CDF of y at the cutoff.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> x = rng.uniform(-1, 1, n)
>>> y = 1.0 + 0.5 * x + (x >= 0) * 0.8 + rng.normal(0, 1, n)
>>> data = pd.DataFrame({"y": y, "x": x})
>>> res = sp.rd_distributional_design(data, y="y", running="x", cutoff=0.0)
>>> res.n_obs
400
>>> len(res.rdd_effect)
5

rd_external_validity

rd_external_validity(data: DataFrame, y: str, x: str, c: float = 0, covs: Optional[List[str]] = None, target_x_range: Optional[tuple] = None, alpha: float = 0.05) -> Dict[str, Any]

Diagnostic assessment of RD external validity.

Compares covariate distributions between the cutoff neighborhood and a target population, runs the conditional independence test, and provides a recommendation on whether extrapolation is credible.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable name.	required
x	str	Running variable name.	required
c	float	RD cutoff value.	0
covs	list of str	Covariate names for overlap and CI testing.	None
target_x_range	tuple of (float, float)	Running variable range (x_low, x_high) defining the target population. Defaults to the full data range.	None
alpha	float	Significance level.	0.05

Returns:

Type	Description
dict	Keys: 'local_estimate': local RD effect at the cutoff (CausalResult). 'ci_test': conditional independence test results (dict or None). 'overlap': covariate overlap statistics (dict or None). 'extrapolated_estimate': extrapolated ATE for the target population if CI test passes (CausalResult or None). 'recommendation': human-readable guidance string.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> n = 2000
>>> Z = rng.normal(0, 1, n)
>>> X = Z + rng.normal(0, 0.5, n)
>>> D = (X >= 0).astype(int)
>>> Y = 1.0 + 2.0 * Z + 3.0 * D + rng.normal(0, 0.5, n)
>>> df = pd.DataFrame({'y': Y, 'x': X, 'z': Z})
>>> diag = sp.rd_external_validity(df, y='y', x='x', c=0, covs=['z'])
>>> bool('recommendation' in diag)
True

rd_extrapolate

rd_extrapolate(data: DataFrame, y: str, x: str, c: float = 0, covs: Optional[List[str]] = None, treatment: Optional[str] = None, eval_points: Optional[ndarray] = None, n_eval: int = 20, method: str = 'ols', h_local: Optional[float] = None, alpha: float = 0.05) -> CausalResult

Angrist-Rokkanen (2015) extrapolation of RD effects away from the cutoff.

If Y(0), Y(1) are independent of X conditional on covariates Z, then treatment effects can be estimated at any value of the running variable, not just at the cutoff.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable name.	required
x	str	Running variable name.	required
c	float	RD cutoff value.	0
covs	list of str	Covariate names for conditional independence. Required.	None
treatment	str	Treatment variable for fuzzy RD. If None, sharp design assumed (D = 1{X >= c}).	None
eval_points	ndarray	Running variable values at which to extrapolate the CATE. If None, n_eval equally spaced points spanning the data range.	None
n_eval	int	Number of evaluation points when eval_points is None.	20
method	str	Estimation method: 'ols', 'ipw', or 'doubly_robust'.	'ols'
h_local	float	Bandwidth for local RD estimate at the cutoff (for comparison). If None, the MSE-optimal bandwidth from rdrobust is used.	None
alpha	float	Significance level.	0.05

Returns:

Type	Description
CausalResult	estimate is the average treatment effect across evaluation points. detail is a DataFrame with columns [x_value, cate, se, ci_lower, ci_upper]. model_info contains conditional independence test results and the local RD estimate for comparison.

Notes

The key identifying assumption is conditional independence: Y(0), Y(1) ⊥ X | Z. A partial F-test is run on each side of the cutoff to check whether X retains predictive power for Y after conditioning on Z. A rejection suggests the assumption may fail.

References

Angrist, J.D. and Rokkanen, M. (2015). "Wanna Get Away? Regression Discontinuity Estimation of Exam School Effects Away from the Cutoff." JASA, 110(512), 1331-1344. [@angrist2015wanna]

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> n = 2000
>>> Z = rng.normal(0, 1, n)
>>> X = Z + rng.normal(0, 0.5, n)
>>> D = (X >= 0).astype(int)
>>> Y = 1.0 + 2.0 * Z + 3.0 * D + rng.normal(0, 0.5, n)
>>> df = pd.DataFrame({'y': Y, 'x': X, 'z': Z})
>>> result = sp.rd_extrapolate(df, y='y', x='x', c=0, covs=['z'])
>>> bool(abs(result.estimate - 3.0) < 1.0)
True

rd_multi_extrapolate

rd_multi_extrapolate(data: DataFrame, y: str, x: str, cutoffs: List[float], eval_points: Optional[ndarray] = None, method: str = 'linear', alpha: float = 0.05) -> CausalResult

Multi-cutoff RD extrapolation (Cattaneo, Keele, Titiunik, Vazquez-Bare 2021).

Estimates local RD effects at each cutoff, then interpolates/extrapolates a treatment effect function tau(x) through those point estimates.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable name.	required
x	str	Running variable name.	required
cutoffs	list of float	Cutoff values. Must contain at least 2 cutoffs.	required
eval_points	ndarray	Running variable values at which to predict tau(x). Defaults to 30 equally spaced points spanning the data range.	None
method	str	Interpolation method: 'linear': tau(x) = a + b*x 'polynomial': polynomial of degree min(len(cutoffs)-1, 3) 'weighted': inverse-variance weighted local linear	'linear'
alpha	float	Significance level.	0.05

Returns:

Type	Description
CausalResult	estimate is the average extrapolated effect. detail is a DataFrame with columns [x_value, cate_extrapolated, se, ci_lower, ci_upper]. model_info contains cutoff-specific estimates and heterogeneity test results.

References

Cattaneo, M.D., Keele, L., Titiunik, R. and Vazquez-Bare, G. (2021). "Extrapolating Treatment Effects in Multi-Cutoff Regression Discontinuity Designs." Journal of the American Statistical Association, 116(536), 1941-1952. doi:10.1080/01621459.2020.1751646 [@cattaneo2021extrapolating]

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> n = 3000
>>> X = rng.uniform(-2, 4, n)
>>> tau_true = 2.0 + 0.5 * X
>>> D = ((X >= 1) | (X >= 3)).astype(int)
>>> Y = 0.5 * X + tau_true * D + rng.normal(0, 0.5, n)
>>> df = pd.DataFrame({'y': Y, 'x': X})
>>> result = sp.rd_multi_extrapolate(df, y='y', x='x', cutoffs=[1.0, 3.0])
>>> bool(result.estimate > 0)
True

rd_honest

rd_honest(data: DataFrame, y: str, x: str, c: float = 0.0, M: Optional[float] = None, kernel: str = 'triangular', h: Optional[float] = None, alpha: float = 0.05, opt_criterion: str = 'mse', sclass: str = 'H', cluster: Optional[str] = None) -> CausalResult

Honest confidence intervals for regression discontinuity designs.

Implements Armstrong & Kolesár (2018, 2020): CIs that are valid uniformly over the class of regression functions whose second derivative is bounded by M. These "honest" CIs account for the smoothing bias that standard local-polynomial CIs ignore, yielding correct coverage even in finite samples.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data.	required
y	str	Outcome variable name.	required
x	str	Running variable name.	required
c	float	RD cutoff.	0
M	float	Upper bound on |f''(c)|. If None, estimated from a local quadratic fit on each side of the cutoff.	None
kernel	str	Kernel function: "triangular", "epanechnikov", or "uniform".	"triangular"
h	float	Bandwidth. If None, chosen by the criterion in opt_criterion.	None
alpha	float	Significance level for the confidence interval.	0.05
opt_criterion	str	Bandwidth selection criterion when h is None: "mse" (MSE-optimal, Imbens-Kalyanaraman style), "flci" (minimises honest CI length), or "oci" (one-sided CI length).	"mse"
sclass	str	Smoothness class the bound M constrains, matching RDHonest. "H" (Holder — f' is M-Lipschitz) is the default and lets the two sides' curvature contributions cancel. "T" (Taylor) only requires |f(x) - f(0) - f'(0)x| <= M x^2 / 2 on each side, which permits no cancellation, so its worst-case bias — and hence the interval — is always at least as wide for the same M.	"H"
cluster	str	Cluster column, as RDHonest(clusterid=): the standard error sums w_i e_i within clusters, and the bandwidth search uses the preliminary variance plus a Moulton within-cluster correlation.	None

Returns:

Type	Description
CausalResult	Result object with model_info containing: - honest_ci : tuple – honest confidence interval - naive_ci : tuple – standard CI for comparison - M : float – smoothness bound used - bias_bound: float – estimated worst-case bias - bandwidth : float – bandwidth used

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 1000
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.5 * x + 2.0 * (x >= 0) + rng.normal(0, 0.4, n)
>>> df = pd.DataFrame({'y': y, 'x': x})
>>> result = sp.rd_honest(df, y='y', x='x', c=0)
>>> bool(abs(result.estimate - 2.0) < 0.3)  # true jump is 2.0
True
>>> lo, hi = result.model_info['honest_ci']
>>> nlo, nhi = result.model_info['naive_ci']
>>> bool((hi - lo) > (nhi - nlo))  # honest CI widens for worst-case bias
True

References

Armstrong, T. B. and Kolesár, M. (2018). Optimal Inference in a Class of Regression Models. Econometrica. [@armstrong2018optimal]

Armstrong, T. B. and Kolesár, M. (2020). Simple and honest confidence intervals in nonparametric regression. Quantitative Economics. [@armstrong2020simple]

rdbwhte

rdbwhte(data: DataFrame, y: str, x: str, z: Union[str, List[str]], c: float = 0, p: int = 1, kernel: str = 'triangular', q: Optional[int] = None, bwselect: str = 'mserd', vce: Optional[str] = None, cluster: Optional[str] = None) -> Any

Bandwidth selection for :func:rdhte.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome and running variable.	required
x	str	Outcome and running variable.	required
z	str or list of str	Moderator(s), as in :func:rdhte.	required
c	float		0
p	int		1
kernel	(triangular, uniform, epanechnikov)		'triangular'
q	int	Bias-correction order (default p + 1).	None
bwselect	str		'mserd'
vce	str	Variance used inside the selector (default 'hc3'; 'cr1' with cluster).	None
cluster	str		None

Returns:

Type	Description
float or DataFrame	One common bandwidth as a float; otherwise a frame with columns level, h_left, h_right (one row per subgroup, or one row with level=None for side-specific bandwidths).

Notes

Continuous moderators: the CCT selector on (y, x). Subgroups: the selector on each level's observations (one bandwidth pair per level).

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(-1, 1, 2000); z = rng.normal(size=2000)
>>> y = x + (2 + z) * (x >= 0) + rng.normal(0, 0.5, 2000)
>>> h = sp.rdbwhte(pd.DataFrame({"y": y, "x": x, "z": z}), y="y", x="x", z="z")
>>> h > 0
True

References

calonico2025treatment, calonico2020optimal

rdhte

rdhte(data: DataFrame, y: str, x: str, z: Union[str, List[str]], c: float = 0, p: int = 1, h: Optional[Union[float, Tuple[float, float]]] = None, b: Optional[float] = None, kernel: str = 'triangular', bwselect: str = 'mserd', cluster: Optional[str] = None, alpha: float = 0.05, eval_points: Optional[ndarray] = None, n_eval: int = 20, q: Optional[int] = None, vce: Optional[str] = None) -> CausalResult

Conditional average treatment effects (CATE) in sharp RD designs.

Fits y on [1, T] x [1, x-c, ..., (x-c)^p] x [1, z] by kernel weighted least squares on |x - c| <= h and reports CATE(z) = beta_T + z' beta_{T:z}, with robust bias-corrected inference from the order-q fit on the same window.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome and running variable.	required
x	str	Outcome and running variable.	required
z	str or list of str	Moderator(s). A single 0/1, boolean or categorical column gives one effect per level (subgroup analysis, a bandwidth per level); numeric columns give a linear CATE function.	required
c	float	Cutoff.	0
p	int	Polynomial order of the point estimate (0 = local constant).	1
h	float or (float, float)	Bandwidth, or (h_left, h_right); default MSE-optimal (bwselect), per level for subgroups.	None
b	float	Window of the bias-correction fit (default h).	None
kernel	(triangular, uniform, epanechnikov)		'triangular'
bwselect	str	One of the CCT selectors (mserd, msetwo, cerrd, ...).	'mserd'
cluster	str	Cluster column (CR1 variance).	None
alpha	float		0.05
eval_points	array	Moderator values at which to report the CATE (continuous mode).	None
n_eval	int	Number of default evaluation points (quantiles of z).	20
q	int	Order of the bias-correction fit (default p + 1).	None
vce	(hc0, hc1, hc2, hc3, cr1)	Default 'hc3', or 'cr1' with cluster.	'hc0'

Returns:

Type	Description
CausalResult	detail: one row per evaluation point (z_value, cate, cate_bc, se, ci_lower, ci_upper, pvalue); the headline is their equally weighted average. model_info carries the coefficients (coef, coef_bc, se_rb, vcov), bandwidths and a Wald test of homogeneity.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(-1, 1, 2000); z = rng.integers(0, 2, 2000)
>>> y = x + (2 + 3 * z) * (x >= 0) + rng.normal(0, 0.5, 2000)
>>> df = pd.DataFrame({"y": y, "x": x, "z": z})
>>> r = sp.rdhte(df, y="y", x="x", z="z")
>>> r.detail["cate"].round(0).tolist()
[2.0, 5.0]

References

calonico2025treatment, calonico2025rdhte

rdhte_lincom

rdhte_lincom(result: CausalResult, weights: Optional[ndarray] = None, alpha: float = 0.05, linfct: Optional[ndarray] = None) -> Dict[str, Any]

Linear combinations of :func:rdhte effects.

Parameters:

Name	Type	Description	Default
result	CausalResult	Output of :func:rdhte.	required
weights	array	Weights over the rows of result.detail (evaluation points or subgroups).	None
alpha	float		0.05
linfct	array	Rows of weights over the reported coefficients (model_info['coef_names']); takes precedence over weights.	None

Returns:

Type	Description
dict	estimate (from the conventional fit), estimate_bc, se, z, ci, pvalue (robust bias-corrected), and joint (Wald statistic, df, p-value). Scalars for one combination, arrays otherwise.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(-1, 1, 2000); z = rng.integers(0, 2, 2000)
>>> y = x + (2 + 3 * z) * (x >= 0) + rng.normal(0, 0.5, 2000)
>>> r = sp.rdhte(pd.DataFrame({"y": y, "x": x, "z": z}), y="y", x="x", z="z")
>>> out = sp.rdhte_lincom(r, linfct=[-1, 1])
>>> round(out["estimate"])
3

References

calonico2025rdhte

rd_interference

rd_interference(data: DataFrame, y: str, running: str, neighbour_running: str, cutoff: float = 0.0, bandwidth: Optional[float] = None, kernel: str = 'triangular', alpha: float = 0.05) -> RDInterferenceResult

Sharp RDD with network interference (Cabrelli-Marconi 2024).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome.	required
running	str	Own running variable.	required
neighbour_running	str	Average running variable across neighbours (precomputed).	required
cutoff	float		0.0
bandwidth	float	Defaults to IQR of own running variable.	None
kernel	str		'triangular'
alpha	float		0.05

Returns:

Type	Description
RDInterferenceResult

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng(1)
>>> n = 600
>>> r = rng.uniform(-1, 1, size=n)    # own running variable
>>> rn = rng.uniform(-1, 1, size=n)   # avg neighbour running variable
>>> y = (0.5 * (r >= 0) + 0.3 * (rn >= 0)
...      + 0.2 * r + 0.1 * rn + rng.normal(scale=0.3, size=n))
>>> df = pd.DataFrame({"y": y, "run": r, "nbr_run": rn})
>>> res = sp.rd_interference(df, y="y", running="run",
...                          neighbour_running="nbr_run")
>>> type(res).__name__
'RDInterferenceResult'
>>> res.n_obs
600
>>> bool(np.isfinite(res.direct_effect))
True

rdrandinf

rdrandinf(data: DataFrame, y: str, x: str, c: float = 0, wl: Optional[float] = None, wr: Optional[float] = None, statistic: str = 'diffmeans', p: int = 0, covs: Optional[List[str]] = None, kernel: str = 'uniform', n_perms: int = 1000, fuzzy: Optional[str] = None, alpha: float = 0.05, seed: int = 42, *, fuzzy_stat: str = 'itt', nulltau: float = 0.0, bernoulli: Optional[Union[str, float]] = None, ci: Union[bool, Sequence[float], None] = None, d: Optional[float] = None, dscale: float = 0.5, evall: Optional[float] = None, evalr: Optional[float] = None, interfci: Optional[float] = None, vce: str = 'hc3') -> CausalResult

Randomization inference for regression discontinuity designs.

Under the local randomization assumption, units within a small window around the cutoff are treated as if randomly assigned. Inference is based on Fisher's randomization test, with a large-sample approximation reported alongside.

.. versionchanged:: 1.39.0 ⚠️ Results changed in three cases. Re-run anything that used them. See MIGRATION.md.

• wl / wr are now the window's endpoints on the scale of the running variable, as in rdlocrand. They used to be offsets from the cutoff, which is the same thing only when c = 0.
• p > 0 used to residualize the outcome on a polynomial pooled across both sides, which absorbs the effect itself (the score is collinear with treatment inside the window). It now fits a polynomial on each side and reports the difference in intercepts.
• fuzzy= used to report a permutation p-value for the Wald ratio, whose randomization distribution has no finite scale when the permuted first stage passes through zero. The default is now the Anderson-Rubin / intention-to-treat randomization test.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable name.	required
x	str	Running variable name.	required
c	float	RD cutoff value.	0
wl	float	Left and right endpoints of the window, on the scale of the running variable (so wl <= c <= wr). Use :func:rdwinselect to choose a window from covariate balance.	None
wr	float	Left and right endpoints of the window, on the scale of the running variable (so wl <= c <= wr). Use :func:rdwinselect to choose a window from covariate balance.	None
statistic	str	Test statistic: 'diffmeans', 'ksmirnov', 'ranksum', or 'all'. 'ttest' is accepted as an alias for 'diffmeans', matching rdlocrand, where both names select the same statistic. The rank and Kolmogorov-Smirnov statistics need p=0 and the uniform kernel.	'diffmeans'
p	int	Order of the polynomial in the score fitted on each side. The statistic is then the difference in fitted values at the cutoff (or at evall / evalr).	0
covs	list of str	Covariate names to partial out of the outcome before testing.	None
kernel	str	'uniform', 'triangular' or 'epanechnikov' weights. The bandwidth on each side is the distance from the cutoff to that side's window edge.	'uniform'
n_perms	int	Number of randomization draws.	1000
fuzzy	str	Column holding the treatment actually received (fuzzy RD). The estimate is then the Wald ratio inside the window, with the heteroskedasticity-robust 2SLS standard error.	None
alpha	float	Significance level for the confidence interval and the power calculation.	0.05
seed	int	Random seed for reproducibility.	42
fuzzy_stat	(itt, tsls)	Which test backs pvalue in a fuzzy design. 'itt' (alias 'ar') is the randomization test on the reduced form, valid in finite samples: the effect on compliers is zero exactly when the intention-to-treat effect is. 'tsls' is the large-sample z-test on the Wald ratio; no randomization p-value is reported for it.	'itt'
nulltau	float	Treatment effect under the null hypothesis.	0
bernoulli	str or float	Assignment probabilities when the mechanism is independent Bernoulli trials rather than a fixed number of treated units: a column name, or one probability for every unit. Needs p=0 and the uniform kernel.	None
ci	bool or sequence of float	Confidence interval by inverting the randomization test of the difference in means. None / True uses a grid of 201 points over the estimate plus or minus five standard errors; a sequence is used as the grid of effects to test; False skips it. The interval is the range of grid values not rejected at alpha, so its resolution is the grid's.	None
d	float	Effect size for the large-sample power calculation. Defaults to dscale times the control group's outcome standard deviation.	None
dscale	float		0.5
evall	float	Points on the running-variable scale at which the fitted polynomials are evaluated. Both default to the cutoff.	None
evalr	float	Points on the running-variable scale at which the fitted polynomials are evaluated. Both default to the cutoff.	None
interfci	float	Level (such as 0.05) of Rosenbaum's confidence interval under arbitrary interference between units, reported in model_info['interf_ci']. It is the observed difference in means minus the upper and lower quantiles of its randomization distribution, and covers the difference between the statistic and what it would have been had no unit been treated. Needs p=0 and a sharp design.	None
vce	(hc3, hc2, hc1)	Variance of the large-sample test when p > 0: the heteroskedasticity-robust standard error of the regression of the outcome on treatment, the polynomial and their interaction. The default is the one of rdlocrand 3.0. Ignored when p = 0, where the standard error is the HC2 one (Welch's with a uniform kernel). .. versionchanged:: 1.39.0 ⚠️ The large-sample p-value, standard error and power for p > 0 used HC2, as rdlocrand up to 2.0 did. The reference moved to HC3 in 3.0 and so does the default here; vce='hc2' gives the earlier numbers.	'hc3'

Returns:

Type	Description
CausalResult	estimate is the difference in means (or in fitted values at the cutoff when p > 0); pvalue is the randomization p-value of the first requested statistic; ci is the test-inversion interval. model_info carries the large-sample p-value, the power against d, per-side means and standard deviations, and, for fuzzy designs, the intention-to-treat and first-stage estimates.

Notes

What is randomized when p > 0. Treatment is a function of the score, so the randomization that local randomization posits is of scores to units. For a difference in means, weighted or not, that is the same as permuting treatment labels. With a polynomial it is not: the statistic extrapolates each side's fit to the cutoff, and its randomization distribution has to be built from the same extrapolation. For p > 0 this function therefore permutes outcomes against (score, assignment) pairs. rdlocrand up to 2.0 produced much smaller p-values for p > 0; on a design with no effect, 40 observations and p = 1, re-randomizing labels with the scores held fixed rejects a true null 37% of the time at the 5% level, against 4.7% here (tests/reference_parity/test_rdlocrand_extensions_parity.py). The observed statistic and the large-sample p-value agree with rdlocrand to 1e-13.

The large-sample p-value uses the normal distribution and the variance vce selects (HC2, Welch's standard error, when p = 0). With p > 0 rdlocrand 3.0 reports this p-value only and no longer computes a randomization p-value; the one reported here is the outcome-permutation p-value described above.

References

Cattaneo, M.D., Titiunik, R. and Vazquez-Bare, G. (2016). "Inference in Regression Discontinuity Designs under Local Randomization." The Stata Journal, 16(2), 331-367. [@cattaneo2016inference],

[@cattaneo2015randomization], [@cattaneo2024extensions], [@rosenbaum2007interference]

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 500
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.8 * (x >= 0) + 0.5 * x + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({"x": x, "y": y})
>>> res = sp.rdrandinf(df, y="y", x="x", c=0.0, wl=-0.3,
...                    wr=0.3, n_perms=500, seed=42)
>>> round(float(res.estimate), 3)
0.926
>>> res.pvalue is not None
True

rdrbounds

rdrbounds(data: DataFrame, y: str, x: str, c: float = 0, wl: Optional[float] = None, wr: Optional[float] = None, gamma_list: Optional[List[float]] = None, statistic: str = 'ranksum', n_perms: int = 1000, seed: int = 42) -> DataFrame

Rosenbaum sensitivity bounds for RD under local randomization.

Assesses how much hidden bias (departure from random assignment) would be needed to explain away the estimated treatment effect, as R rdlocrand::rdrbounds does. Under Rosenbaum's model with odds ratio Gamma each unit in the window is treated independently with probability Gamma/(1+Gamma) or 1/(1+Gamma). As in R, the bounds are taken over the monotone patterns in which the u units with the largest outcomes get the high probability, u = 1..n: the upper bound is the largest of those p-values and the lower bound the smallest. Gamma = 1 is randomization with the treated share as the common probability.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable name.	required
x	str	Running variable name.	required
c	float	RD cutoff value.	0
wl	float	Window endpoints on the running-variable scale (wl <= c <= wr).	None
wr	float	Window endpoints on the running-variable scale (wl <= c <= wr).	None
gamma_list	list of float	Odds ratios. Defaults to [1, 1.5, 2, 2.5, 3, 4, 5].	None
statistic	(ranksum, diffmeans)	'ranksum' is R's standardised rank sum (control-rank sum minus its mean over sqrt(n0 n1 var(ranks) / n)); 'diffmeans' the difference in means.	'ranksum'
n_perms	int	Bernoulli assignment draws per p-value. The same draws are reused across thresholds and Gamma values (common random numbers, as R reseeds before each p-value).	1000
seed	int	Random seed.	42

Returns:

Type	Description
DataFrame	Columns: gamma, pvalue_upper, pvalue_lower.

Notes

Randomisation p-values: agreement with R is within Monte-Carlo error only. Through 1.28.0 this function split units at the median outcome only -- one threshold instead of the extremum over all of them -- drew a fixed number of treated units rather than Bernoulli assignments, and gave the high probability to the smallest outcomes for the lower bound; R's uminus on its increasing sort selects the largest outcomes, the same patterns as the upper bound, which is what this port follows (with R's draw-to-unit attachment for each bound).

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 500
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.8 * (x >= 0) + 0.5 * x + rng.normal(0, 0.3, n)
>>> df = pd.DataFrame({"x": x, "y": y})
>>> tab = sp.rdrbounds(df, y="y", x="x", c=0.0, wl=-0.3,
...                    wr=0.3, gamma_list=[1.0, 1.5, 2.0],
...                    n_perms=500, seed=42)
>>> tab.shape
(3, 3)
>>> list(tab.columns)
['gamma', 'pvalue_upper', 'pvalue_lower']
>>> bool((tab["pvalue_upper"] >= tab["pvalue_lower"]).all())
True

rdsensitivity

rdsensitivity(data: DataFrame, y: str, x: str, c: float = 0, wlist: Optional[List[float]] = None, nwindows: int = 20, statistic: str = 'diffmeans', p: int = 0, n_perms: int = 500, seed: int = 42, alpha: float = 0.05, plot: bool = False, vce: str = 'hc3') -> DataFrame

Sensitivity of RD estimates across different window widths.

For each window, runs rdrandinf and records the estimate, standard error, and p-value. Optionally produces a plot if matplotlib is available.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable name.	required
x	str	Running variable name.	required
c	float	RD cutoff value.	0
wlist	list of float	Symmetric half-window widths to evaluate. If None, an evenly-spaced grid is generated automatically.	None
nwindows	int	Number of windows when wlist is None.	20
statistic	str	Test statistic for inference.	'diffmeans'
p	int	Polynomial order for adjustment.	0
n_perms	int	Number of permutations per window.	500
seed	int	Random seed.	42
alpha	float	Significance level.	0.05
vce	(hc3, hc2, hc1)	Variance of the large-sample test when p > 0; see :func:rdrandinf.	'hc3'

Returns:

Type	Description
DataFrame	Columns: window, estimate, se, pvalue, ci_lower, ci_upper, significant.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> x = rng.uniform(-1, 1, 500)
>>> y = 0.5 * (x >= 0) + 0.8 * x + rng.normal(0, 0.3, 500)
>>> df = pd.DataFrame({"x": x, "y": y})
>>> sens = sp.rdsensitivity(
...     df, y="y", x="x", wlist=[0.25, 0.5, 0.75],
...     n_perms=200, seed=42,
... )
>>> sens["estimate"].round(3).tolist()
[0.66, 0.834, 1.078]
>>> bool(sens["significant"].all())
True

rdwinselect

rdwinselect(data: DataFrame, x: str, c: float = 0, covs: Optional[List[str]] = None, wmin: Optional[float] = None, wstep: Optional[float] = None, nwindows: int = 10, statistic: str = 'diffmeans', p: int = 0, seed: int = 42, alpha: float = 0.15, *, obsmin: int = 10, wobs: Optional[int] = None, wasymmetric: bool = False, approx: bool = False, n_perms: int = 1000, kernel: str = 'uniform', dropmissing: bool = False, wmasspoints: bool = False, vce: str = 'hc3') -> DataFrame

Data-driven window selection for local randomization RD.

Builds a sequence of nested windows around the cutoff and, in each, tests the balance of every covariate and runs a binomial test on the number of observations either side. The recommended window is the largest one such that it and every window inside it have a minimum balance p-value of at least alpha.

.. versionchanged:: 1.39.0 ⚠️ The window sequence and the p-values changed. Re-run anything that used this function. The default sequence used to start at a fraction of the running variable's range and step to its full extent, so on the U.S. Senate data it tested windows of 5, 15.6, ..., 100 points where the procedure is about windows holding ten, twelve, fourteen observations a side. It now starts at the smallest window with obsmin observations on each side and adds wobs per side per step, as described by Cattaneo, Idrobo and Titiunik (2024). Without covariates the function used to test made-up quantile dummies of the score; it now reports the binomial test only. See MIGRATION.md.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
x	str	Running variable name.	required
c	float	RD cutoff value.	0
covs	list of str	Predetermined covariates to test balance for. Without them only the binomial test is reported and no window is recommended.	None
wmin	float	Half-width of the smallest window. Defaults to the smallest symmetric window with obsmin observations on each side.	None
wstep	float	Increment in half-width. If omitted, each window adds at least wobs observations on each side.	None
nwindows	int	Number of windows to evaluate. Fewer come back if the data run out.	10
statistic	str	Balance statistic: 'diffmeans' (alias 'ttest'), 'ksmirnov', 'ranksum', or 'hotelling'. The first three test each covariate and report the smallest p-value; 'hotelling' is one joint test of all covariates (Hotelling's T-squared), so variable is empty.	'diffmeans'
p	int	Polynomial order for the covariate adjustment model.	0
seed	int	Random seed.	42
alpha	float	Smallest acceptable minimum p-value (level in rdlocrand). It is deliberately larger than 0.05: the concern is failing to detect imbalance, not falsely detecting it.	0.15
obsmin	int	Minimum observations on each side in the smallest window.	10
wobs	int	Observations added on each side per step (5 if neither wobs nor wstep is given).	None
wasymmetric	bool	Let the two sides of the window grow separately.	False
approx	bool	Use the large-sample p-values instead of randomization inference. Deterministic and much faster.	False
n_perms	int	Randomization draws per covariate and window.	1000
kernel	str		'uniform'
dropmissing	bool	Drop rows with a missing covariate before the windows are built.	False
vce	(hc3, hc2, hc1)	Variance of the large-sample balance tests when p > 0; see :func:rdrandinf. Ignored when p = 0. By default the windows are built on every row with an observed score and incomplete rows are dropped inside each window.	'hc3'
wmasspoints	bool	For a discrete running variable: window k runs from the k-th mass point below the cutoff to the k-th at or above it, so each step adds one support point on each side. Cannot be combined with obsmin, wmin, wobs or wstep.	False

Returns:

Type	Description
DataFrame	One row per window: window_left, window_right (endpoints on the running-variable scale), n_left, n_right, p_value (smallest balance p-value across covariates), variable (the covariate attaining it), binom_pvalue and balanced. attrs['recommended_window'] holds (left, right) or None; attrs['recommended_n'] the counts inside it.

Notes

The window sequences, counts, binomial p-values and large-sample balance p-values agree with rdlocrand 1.0 and 3.0 to 1e-9, mass-point windows included (tests/reference_parity/test_rdlocrand_v1_parity.py). Releases 1.1 and 2.0 start the default sequence one observation short below the cutoff (9 with obsmin = 10 on the Senate data) and shift the left edge of mass-point windows by one support point; the help page and Cattaneo, Idrobo and Titiunik (2024, Snippet 2.5) describe what 1.0 and 3.0 do. With wmin= and wstep= all releases agree.

References

[@cattaneo2016inference], [@cattaneo2015randomization], [@cattaneo2024extensions]

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 500
>>> x = rng.uniform(-1, 1, n)
>>> z1 = 0.2 * x + rng.normal(0, 1, n)
>>> df = pd.DataFrame({"x": x, "z1": z1})
>>> tab = sp.rdwinselect(df, x="x", c=0.0, covs=["z1"],
...                      nwindows=5, approx=True)
>>> tab.shape
(5, 8)
>>> list(tab.columns)[:5]
['window_left', 'window_right', 'n_left', 'n_right', 'p_value']
>>> bool(tab["n_left"].iloc[0] >= 10 and tab["n_right"].iloc[0] >= 10)
True

rd_multi_score

rd_multi_score(data: DataFrame, y: str, running_vars: List[str], cutoffs: List[float], bandwidth: Optional[float] = None, kernel: str = 'triangular', alpha: float = 0.05) -> MultiScoreRDResult

Multi-score RDD: treatment if all running variables exceed cutoffs.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome.	required
running_vars	list of str	Multiple running variables.	required
cutoffs	list of float	One cutoff per running variable (same length).	required
bandwidth	float	Defaults to median IQR across running vars.	None
kernel	str		'triangular'
alpha	float		0.05

Returns:

Type	Description
MultiScoreRDResult

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 2000
>>> sat_math = rng.uniform(400, 800, n)
>>> sat_verbal = rng.uniform(400, 800, n)
>>> treat = ((sat_math >= 600) & (sat_verbal >= 600)).astype(int)
>>> y = (10.0 + 0.01 * sat_math + 0.01 * sat_verbal
...      + 5.0 * treat + rng.normal(0, 2, n))
>>> df = pd.DataFrame({"y": y, "math": sat_math, "verbal": sat_verbal})
>>> res = sp.rd_multi_score(
...     df, y="y", running_vars=["math", "verbal"],
...     cutoffs=[600.0, 600.0])
>>> res.n_obs
2000
>>> isinstance(res.boundary_effect, float)
True
>>> print(res.summary())

rd_optimized

rd_optimized(data: DataFrame, y: str, x: str, c: float = 0.0, M: Optional[float] = None, *, fuzzy: Optional[str] = None, M_fuzzy: Optional[float] = None, criterion: str = 'mse', h: Optional[float] = None, sigma2: Optional[float] = None, num_bins: int = 120, n_cells: int = 120, alpha: float = 0.05) -> CausalResult

Sharp regression discontinuity with minimax-optimal linear weights.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome column.	required
x	str	Running variable; units with x >= c are treated. It may be continuous or discrete.	required
c	float	Cutoff.	0.0
M	float	Bound on the absolute second derivative of the conditional mean on each side of the cutoff. If omitted, the Armstrong-Kolesar rule of thumb is used (a global quartic fit on each side), as in :func:rd_honest. The interval is only as credible as this bound: report it and vary it.	None
fuzzy	str	Column with the treatment actually received, for a fuzzy design in which crossing the cutoff changes the probability of treatment. The estimate is then the ratio of the jump in the outcome to the jump in this column.	None
M_fuzzy	float	Bound on the second derivative of the conditional mean of the fuzzy column. Defaults to the same rule of thumb as M.	None
criterion	('mse', 'flci')	"mse" minimises the worst-case mean squared error; "flci" the length of the confidence interval.	"mse"
h	float	Use only observations with |x - c| <= h. By default the window is 2.5 times the optimal local linear bandwidth, which is wide enough for the optimal weights to reach zero inside it.	None
sigma2	float	Outcome variance used to choose the weights (the trade-off depends on M**2 / sigma2). By default it is the average of nearest-neighbour variance estimates on each side of the cutoff. The standard error always uses the unit-level estimates.	None
num_bins	int	If a side has more distinct values of x than this, they are grouped into that many equal-width bins and units in a bin share a weight. The reported bias and standard error are computed from the unit-level data in either case.	120
n_cells	int	Approximate number of pieces, per side, of the piecewise constant second derivative of the least favourable function.	120
alpha	float		0.05

Returns:

Type	Description
CausalResult	model_info holds the bound M, the worst-case bias (max_bias), the unit-level weights (weights, aligned with model_info['index']), the largest distance from the cutoff that carries weight (effective_bandwidth) and the same quantities for local linear regression with a triangular kernel at its optimal bandwidth (local_linear).

Notes

The interval is estimate +/- cv * se where cv is the 1 - alpha quantile of |N(max_bias / se, 1)|, as in :func:rd_honest. It has the stated coverage for every conditional mean function with curvature at most M, whether the running variable is continuous or not. The standard error uses nearest-neighbour variance estimates.

With fuzzy= the weights are still those chosen for the outcome. The interval collects every ratio t for which the weighted sum of Y - t * D is within its own bias-aware critical value of zero, where the worst-case bias uses the bound M + |t| * M_fuzzy. This is an Anderson-Rubin construction: it stays valid when the jump in treatment is small, in which case it can be unbounded (reported as infinite endpoints). se is the delta-method standard error of the ratio and is given for reference only.

Compared with :func:rd_honest, which uses local linear weights under the same bound, the gain in interval length is a few percent for a continuous running variable and can be larger for a discrete one. The two functions share M and should tell the same story.

References

[@imbens2019optimized], [@armstrong2018optimal], [@noack2024biasaware]

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.uniform(-1, 1, 800)
>>> yv = 0.5 * (x >= 0) + np.sin(x) + rng.normal(scale=0.3, size=800)
>>> r = sp.rd_optimized(pd.DataFrame({"y": yv, "x": x}), "y", "x", M=1.0)
>>> bool(r.ci[0] < 0.5 < r.ci[1])
True

rd2d_bw

rd2d_bw(data: DataFrame, y: str, x1: str, x2: str, treatment: str, boundary: Optional[Callable] = None, approach: str = 'location', p: int = 1, kernel: str = 'triangular', *, eval_points: Optional[ndarray] = None, n_eval: int = 1, deriv: Tuple[int, int] = (0, 0), tangvec: Optional[ndarray] = None, kernel_type: str = 'prod', bwselect: str = 'mserd', method: str = 'dpi', vce: str = 'hc1', cluster: Optional[str] = None, fuzzy: Optional[str] = None, bwparam: str = 'main', fitmethod: str = 'joint', masspoints: str = 'check', bwcheck: Union[int, str, None] = 'auto', scaleregul: float = 1.0, scalebiascrct: float = 1.0, stdvars: bool = True, kink_unknown: Union[bool, Tuple[bool, bool]] = False, kink_position: Optional[Sequence] = None, cqt: float = 0.5, distance: Optional[Union[ndarray, List[str]]] = None) -> DataFrame

Bandwidth selection for boundary discontinuity designs.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome, the two running variables and the 0/1 side indicator.	required
x1	str	Outcome, the two running variables and the 0/1 side indicator.	required
x2	str	Outcome, the two running variables and the 0/1 side indicator.	required
treatment	str	Outcome, the two running variables and the 0/1 side indicator.	required
boundary	callable	Boundary x2 = f(x1) for default evaluation points.	None
approach	('location', 'distance')		'location'
p	int		1
kernel	('triangular', 'epanechnikov', 'uniform', 'gaussian')		'triangular'
eval_points	array(k, 2)		None
n_eval	int		1
deriv	(int, int)		(0, 0)
tangvec	array(k, 2)		None
kernel_type	('prod', 'rad')		'prod'
bwselect	str	One of mserd, cerrd, imserd, icerrd, msetwo, certwo, imsetwo, icertwo.	'mserd'
method	('dpi', 'rot')	Bias-bandwidth rule (location approach).	'dpi'
vce	('hc0', 'hc1', 'hc2', 'hc3')		'hc0'
cluster	str		None
fuzzy	str		None
bwparam	('main', 'itt')		'main'
fitmethod	('joint', 'separate')		'joint'
masspoints	('check', 'adjust', 'off')		'check'
bwcheck	(int, 'auto' or None)	'auto' is 20 (location) or 20 + p + 1 (distance).	'auto'
scaleregul	float		1.0
scalebiascrct	float		1.0
stdvars	bool		True
kink_unknown	bool or (bool, bool)		False
kink_position	sequence		None
cqt	float		0.5
distance	array (n, k) or list of str		None

Returns:

Type	Description
DataFrame	One row per evaluation point: b1, b2 and h01, h02, h11, h12 (location) or h0, h1 (distance).

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x1, x2 = rng.uniform(-1, 1, 2000), rng.uniform(-1, 1, 2000)
>>> t = (x1 >= 0).astype(int)
>>> y = x1 + 2 * t + rng.normal(0, 0.3, 2000)
>>> df = pd.DataFrame({"y": y, "x1": x1, "x2": x2, "t": t})
>>> list(sp.rd2d_bw(df, y="y", x1="x1", x2="x2", treatment="t").columns)
['b1', 'b2', 'h01', 'h02', 'h11', 'h12']

References

cattaneo2025boundary

rd2d_plot

rd2d_plot(data: DataFrame, y: str, x1: str, x2: str, treatment: str, boundary: Optional[Callable[..., Any]] = None, result: Optional[CausalResult] = None, plot_type: str = 'scatter', ax: Optional[Any] = None, figsize: Tuple[float, float] = (10, 8)) -> Tuple[Any, Any]

2D boundary RD visualization.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome variable name.	required
x1	str	Running variable names.	required
x2	str	Running variable names.	required
treatment	str	Binary assignment indicator.	required
boundary	callable	Boundary function f(x1) -> x2; None draws the line x1 = 0.	None
result	CausalResult	Result of :func:sp.rd2d; its evaluation points are marked on the scatter, and plot_type='boundary_effects' plots its pointwise bias-corrected effects with robust intervals.	None
plot_type	('scatter', 'heatmap', 'boundary_effects')		'scatter'
ax	matplotlib Axes		None
figsize	tuple		(10, 8)

Returns:

Type	Description
(fig, ax)

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 2000
>>> x1 = rng.uniform(-1, 1, n)
>>> x2 = rng.uniform(-1, 1, n)
>>> treat = (x1 >= 0).astype(int)
>>> y = 1.5 * treat + 0.5 * x1 + 0.3 * x2 + rng.normal(0, 0.5, n)
>>> df = pd.DataFrame({"y": y, "x1": x1, "x2": x2, "treat": treat})
>>> res = sp.rd2d(df, y="y", x1="x1", x2="x2", treatment="treat")
>>> fig, ax = sp.rd2d_plot(
...     df, y="y", x1="x1", x2="x2", treatment="treat", result=res
... )

rd_boost

rd_boost(data: DataFrame, y: str, x: str, c: float = 0, covs: Optional[List[str]] = None, h: Optional[float] = None, n_estimators: int = 200, max_depth: int = 3, learning_rate: float = 0.1, alpha: float = 0.05, seed: int = 42) -> CausalResult

Gradient Boosting for RD — flexible CATE estimation.

Removes a side-specific linear trend in the running variable, so the outcome is carried to the cutoff, then fits separate GBM models on each side and estimates individual-level CATEs as mu_1(z) - mu_0(z). The headline estimate / se / ci are the average effect at the cutoff from a local-linear RD (sp.rdrobust, robust bias-corrected) on the outcome residualised on the covariates with the same GBM, 5-fold cross-fitted, as in sp.rd_flex. CATE standard errors are obtained via bootstrap.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable.	required
x	str	Running variable.	required
c	float	RD cutoff.	0
covs	list of str	Covariate names for heterogeneity.	None
h	float	Bandwidth (auto-selected if None).	None
n_estimators	int	Number of boosting rounds.	200
max_depth	int	Maximum tree depth per round.	3
learning_rate	float	Shrinkage factor.	0.1
alpha	float	Significance level.	0.05
seed	int	Random seed.	42

Returns:

Type	Description
CausalResult	estimate = average effect at the cutoff (see above); detail = per-obs CATE and bootstrap SE; model_info['mean_cate'] is the mean of the fitted CATEs.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 600
>>> x = rng.uniform(-1, 1, n)
>>> z1 = rng.normal(size=n)
>>> z2 = rng.normal(size=n)
>>> treat = (x >= 0).astype(float)
>>> y = (1.0 + 0.8 * z1) * treat + 0.5 * x + rng.normal(0, 0.5, n)
>>> df = pd.DataFrame({"y": y, "x": x, "z1": z1, "z2": z2})
>>> res = sp.rd_boost(
...     df, y="y", x="x", c=0, covs=["z1", "z2"],
...     n_estimators=50, seed=42,
... )
>>> bool(res.ci[0] < 1.0 < res.ci[1])  # the average effect is 1.0
True

rd_cate_summary

rd_cate_summary(data: DataFrame, y: str, x: str, c: float = 0, covs: Optional[List[str]] = None, h: Optional[float] = None, methods: Optional[List[str]] = None, alpha: float = 0.05, seed: int = 42) -> dict

Run multiple ML-RD methods and compare CATE estimates.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable.	required
x	str	Running variable.	required
c	float	RD cutoff.	0
covs	list of str	Covariates for heterogeneity / selection.	None
h	float	Bandwidth (shared across methods).	None
methods	list of str	Subset of ['forest', 'boost', 'lasso']. Default: all three.	None
alpha	float	Significance level.	0.05
seed	int	Random seed.	42

Returns:

Type	Description
dict	Keys = method names, values = CausalResult. Additional key 'comparison' = pd.DataFrame summarising ATEs. Additional key 'heterogeneity_drivers' = top variable importances from the forest (if run).

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 600
>>> x = rng.uniform(-1, 1, n)
>>> z1 = rng.normal(size=n)
>>> z2 = rng.normal(size=n)
>>> treat = (x >= 0).astype(float)
>>> y = (1.0 + 0.8 * z1) * treat + 0.5 * x + rng.normal(0, 0.5, n)
>>> df = pd.DataFrame({"y": y, "x": x, "z1": z1, "z2": z2})
>>> out = sp.rd_cate_summary(
...     df, y="y", x="x", c=0, covs=["z1", "z2"],
...     methods=["forest", "lasso"], seed=42,
... )
>>> sorted(k for k in out if k in ("forest", "lasso"))
['forest', 'lasso']
>>> bool(out["forest"].ci[0] < 1.0 < out["forest"].ci[1])
True

rd_forest

rd_forest(data: DataFrame, y: str, x: str, c: float = 0, covs: Optional[List[str]] = None, h: Optional[float] = None, n_trees: int = 500, min_leaf: int = 20, honesty: bool = True, alpha: float = 0.05, seed: int = 42) -> CausalResult

Causal Forest for RD — heterogeneous treatment effect estimation.

Adapts the Athey-Wager (2019) generalized random forests framework to an RD context. Within bandwidth h, treatment is D = 1(X >= c). The outcome is first carried to the cutoff by removing a side-specific linear trend in the running variable; two random forests then estimate E[Y|Z, D=0] and E[Y|Z, D=1] at the cutoff, and the CATE is their difference.

The headline estimate / se / ci are the average effect at the cutoff from a local-linear RD (sp.rdrobust, robust bias-corrected) on the outcome residualised on the covariates with the same forest, 5-fold cross-fitted, as in sp.rd_flex. The mean of the fitted CATEs is reported separately as model_info['mean_cate'].

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable.	required
x	str	Running variable.	required
c	float	RD cutoff.	0
covs	list of str	Covariate names used as features for heterogeneity detection. Must not include the running variable x.	None
h	float	Bandwidth (uses IK-style automatic selection if None).	None
n_trees	int	Number of trees in each forest.	500
min_leaf	int	Minimum leaf size (larger → more regularisation).	20
honesty	bool	Split-sample (honest) estimation: half the data for tree construction, the other half for leaf predictions.	True
alpha	float	Significance level for confidence intervals.	0.05
seed	int	Random seed.	42

Returns:

Type	Description
CausalResult	estimate = average effect at the cutoff (see above); detail = DataFrame with per-obs CATE and its across-tree SE; model_info includes mean_cate, running_slopes and variable importance.

Notes

Until the 2026-10 fix the forests were fitted to the raw outcome. They see covariates only, so the result was the difference of within-bandwidth means, off by the slope in the running variable times the distance from the cutoff, with a standard error that ignored it.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 800
>>> x = rng.uniform(-1, 1, n)
>>> z1 = rng.normal(size=n)
>>> z2 = rng.normal(size=n)
>>> treat = (x >= 0).astype(float)
>>> y = (1.0 + 0.8 * z1) * treat + 0.5 * x + rng.normal(0, 0.5, n)
>>> df = pd.DataFrame({"y": y, "x": x, "z1": z1, "z2": z2})
>>> res = sp.rd_forest(
...     df, y="y", x="x", c=0, covs=["z1", "z2"],
...     n_trees=50, seed=42,
... )
>>> bool(res.ci[0] < 1.0 < res.ci[1])  # the average effect is 1.0
True
>>> res.estimand
'ATE at cutoff'
>>> list(res.model_info["variable_importance"].keys())[0]
'z1'

rd_lasso

rd_lasso(data: DataFrame, y: str, x: str, c: float = 0, covs: Optional[List[str]] = None, h: Optional[float] = None, kernel: str = 'triangular', cv_folds: int = 5, alpha: float = 0.05) -> CausalResult

LASSO-assisted RD via post-double-selection.

Uses LASSO to select relevant covariates from a potentially large set, then runs a local linear RD with the selected covariates. Follows Belloni, Chernozhukov, and Hansen (2014).

Parameters:

Name	Type	Description	Default
data	DataFrame	Input dataset.	required
y	str	Outcome variable.	required
x	str	Running variable.	required
c	float	RD cutoff.	0
covs	list of str	Candidate covariates (can be large set).	None
h	float	Bandwidth (auto-selected if None).	None
kernel	str	Kernel for local linear regression ('triangular', 'uniform', 'epanechnikov').	'triangular'
cv_folds	int	Cross-validation folds for LASSO penalty selection.	5
alpha	float	Significance level.	0.05

Returns:

Type	Description
CausalResult	estimate = RD treatment effect with LASSO-selected covariates; model_info includes selected_covariates and LASSO paths.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 800
>>> x = rng.uniform(-1, 1, n)
>>> treat = (x >= 0).astype(float)
>>> Z = rng.normal(size=(n, 8))
>>> y = (1.0 * treat + 0.5 * x + 0.7 * Z[:, 0]
...      + rng.normal(0, 0.5, n))
>>> covs = [f"z{i}" for i in range(8)]
>>> df = pd.DataFrame({"y": y, "x": x})
>>> for i, name in enumerate(covs):
...     df[name] = Z[:, i]
>>> res = sp.rd_lasso(df, y="y", x="x", c=0, covs=covs)
>>> round(float(res.estimate), 2)
1.16
>>> res.model_info["n_selected"]
1

rdmc

rdmc(data: DataFrame, y: str, x: str, cutoffs: Optional[List[float]] = None, bandwidth: Optional[float] = None, kernel: str = 'triangular', pooling: str = 'ivw', alpha: float = 0.05, cutoff_var: Optional[str] = None, **rdrobust_kwargs: object) -> RDMultiResult

Multi-cutoff RD design.

Estimates treatment effects at multiple cutoffs and pools them.

Two designs, and they are not the same estimator:

• Unit-specific cutoffs (cutoff_var=): each unit faces its own cutoff, given by a column. This is the design of Cattaneo, Titiunik, Vazquez-Bare & Keele (2016) and what rdmulti::rdmc() implements -- its C argument. Each cutoff is estimated on only the units assigned to it, by delegating to :func:statspai.rdrobust, so it inherits the CCT bandwidth cascade and robust bias correction.
• Shared running variable (cutoffs= alone, the default): one running variable crossed by several thresholds, every unit entering every threshold's local regression under a common rule-of-thumb bandwidth. For that design with the CCT bandwidth and robust inference, use sp.rdms(data, y, x1, cutoff1=[...]).

Before 1.21 only the second was available, while the docstring claimed equivalence to rdmulti::rdmc(). On a design with unit-specific cutoffs and effects of 2.0 / 5.0 / -3.0, the shared-running-variable path returned 0.22 / 0.51 / 0.66 -- because units belonging to one cutoff were pooled into every other cutoff's window, averaging the effects away. Pass cutoff_var= for the R-equivalent estimator.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome variable.	required
x	str	Running variable.	required
cutoffs	list of float	Cutoff values.	None
bandwidth	float	Bandwidth for local polynomial. If None, uses Silverman rule.	None
kernel	str		'triangular'
pooling	str	Pooling method: 'ivw' (inverse-variance weighted) or 'equal'.	'ivw'
alpha	float		0.05

Returns:

Type	Description
RDMultiResult

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 900
>>> running_var = rng.uniform(0, 100, n)
>>> cutoffs = [30, 60]
>>> nearest = np.array([min(cutoffs, key=lambda c: abs(rv - c))
...                     for rv in running_var])
>>> treated = (running_var >= nearest).astype(float)
>>> score = 1.0 + 0.02 * running_var + 2.0 * treated + rng.normal(0, 1, n)
>>> df = pd.DataFrame({"score": score, "running_var": running_var})
>>> result = sp.rdmc(df, y="score", x="running_var", cutoffs=[30, 60])
>>> result.n_cutoffs
2
>>> summary_text = result.summary()
>>> ax = result.plot()

rdms

rdms(data: DataFrame, y: str, x1: str, x2: Optional[str] = None, cutoff1: Union[float, Sequence[float]] = 0, cutoff2: Union[float, Sequence[float]] = 0, treat: Optional[str] = None, bandwidth: Optional[float] = None, kernel: str = 'triangular', alpha: float = 0.05, xnorm: Optional[str] = None, ranges: Optional[Sequence[Tuple[float, float]]] = None, **rdrobust_kwargs: object) -> Union[CausalResult, RDMultiResult]

Multi-score / Geographic RD design at one or several boundary points.

With sequences for cutoff1 and cutoff2 every boundary point is estimated in turn and an :class:RDMultiResult comes back, one row per point, as rdmulti::rdms prints it. xnorm= names a column holding each unit's signed distance to the boundary as a whole (the perpendicular distance, positive for treated units); the RD on that score is reported as the pooled estimate. With scalar cutoffs and no xnorm the single-point :class:CausalResult is returned, as before.

Without x2 the design is one score with cumulative cutoffs: every unit is exposed to each threshold it has passed (a dose that steps up at cutoff1[0], cutoff1[1], ...). Each cutoff is then a sharp RD on x1, on the whole sample or, with ranges, on the units whose score lies in that cutoff's (lower, upper) interval.

Collapses the two-score design onto the score the reference uses -- Euclidean distance to the boundary point, signed by treatment status,

``d = sqrt((x1 - cutoff1)^2 + (x2 - cutoff2)^2) * (2 * treat - 1)``

-- and estimates a sharp RD at d = 0 by delegating to :func:statspai.rdrobust. That is what rdmulti::rdms() does, so the CCT bandwidth cascade, the robust bias correction and the heteroskedasticity-robust variance all come with it.

Following the reference, the reported estimate is the bias-corrected point estimate with the robust interval (Estimate[2] with ci[3, ] in R's return), not the conventional pair; the conventional numbers are in model_info["conventional"].

.. versionchanged:: 1.27.0 ⚠️ Numerical output changed. Re-run anything that used this function. It previously fitted a 2D local linear in (x1, x2) inside a Euclidean window whose width came from Silverman's kernel density rule of thumb, with a homoskedastic variance and no bias correction. On a 6,000-row design with known effects of 1.5 / 2.0 / 2.5 at three boundary points it used 13 / 8 / 27 observations where rdmulti::rdms used 519 / 725 / 743, and returned 1.322 / 4.804 / 2.326 against the reference's 1.396 / 1.778 / 2.463 -- the middle point 170% out, with a standard error of 4.21 against 0.13. See MIGRATION.md.

.. versionadded:: 1.27.0 treat=, the treatment indicator (zvar in R). Treatment on a two-dimensional boundary is not implied by the coordinates, which is why the reference requires it. Omitting it falls back to x1 >= cutoff1 and warns, because that is an assumption about the design rather than a fact about the data.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome variable.	required
x1	str	First running variable (e.g., latitude distance to boundary).	required
x2	str	Second running variable (e.g., longitude distance to boundary). Omit it for one score with cumulative cutoffs.	None
cutoff1	float	Cutoff for x1.	0
cutoff2	float	Cutoff for x2.	0
treat	str	Column holding the 0/1 treatment indicator (R's zvar). Strongly recommended: without it treatment is assumed to be x1 >= cutoff1 and a warning is issued.	None
bandwidth	float	Fixed bandwidth on the signed-distance score. When omitted the CCT MSE-optimal cascade selects it, as the reference does.	None
kernel	str		'triangular'
alpha	float		0.05
xnorm	str	Column with the signed distance to the boundary (positive on the treated side). Adds the pooled estimate: a sharp RD at zero on that score.	None
ranges	sequence of (lower, upper)	Cumulative cutoffs only: for each cutoff, the interval of the score whose units enter its estimate. Use it to keep a neighbouring cutoff out of a cutoff's estimation window. Default: all units.	None
**rdrobust_kwargs	object	Forwarded to :func:statspai.rdrobust (p, bwselect, vce, cluster, covs, ...).	{}

Returns:

Type

Description

CausalResult or RDMultiResult

A :class:RDMultiResult when cutoff1 / cutoff2 are sequences or xnorm is given: cutoff_results holds, per boundary point, the conventional estimate with the robust standard error, interval and p-value; pooled_estimate is the xnorm RD (NaN without xnorm). Otherwise a :class:CausalResult: estimate / se / ci are the bias-corrected point estimate with robust inference. model_info carries the selected bandwidths, effective sample sizes either side, the conventional estimates, and treat_assumed recording whether the treatment convention was inferred.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n = 800
>>> dist_lat = rng.uniform(-5, 5, n)
>>> dist_lon = rng.uniform(-5, 5, n)
>>> treated = (dist_lat >= 0).astype(float)
>>> outcome = (1.0 + 1.5 * treated + 0.1 * dist_lat
...            + 0.1 * dist_lon + rng.normal(0, 1, n))
>>> df = pd.DataFrame({"outcome": outcome,
...                    "dist_lat": dist_lat, "dist_lon": dist_lon})
>>> result = sp.rdms(df, y="outcome", x1="dist_lat", x2="dist_lon",
...                  bandwidth=3.0)
>>> summary_text = result.summary()

rdsampsi

rdsampsi(tau: float, var_left: float = 1.0, var_right: float = 1.0, h_left: float = 0.5, h_right: float = 0.5, alpha: float = 0.05, target_power: float = 0.8, ratio: float = 1.0, data: Optional['pd.DataFrame'] = None, y: Optional[str] = None, x: Optional[str] = None, c: float = 0.0, **rdrobust_kwargs: object) -> RDSampSiResult

Minimum sample size for a given power in an RD design.

Parameters:

Name	Type	Description	Default
tau	float	Treatment effect under the alternative (non-zero).	required
var_left	float	Design mode: outcome variance on each side.	1.0
var_right	float	Design mode: outcome variance on each side.	1.0
h_left	float	Design mode: bandwidth on each side.	0.5
h_right	float	Design mode: bandwidth on each side.	0.5
alpha	float	Significance level of the two-sided test.	0.05
target_power	float	Power to reach.	0.80
ratio	float	Design mode: n_right / n_left.	1.0
data	DataFrame	Data mode: a pilot sample. Requires y and x.	None
y	str	Data mode: outcome and running-variable columns.	None
x	str	Data mode: outcome and running-variable columns.	None
c	float	Data mode: cutoff.	0.0
**rdrobust_kwargs	object	Data mode: options passed to sp.rdrobust.	{}

Returns:

Type	Description
RDSampSiResult	n_left, n_right, n_total (data mode: observations inside the bandwidth window on each side), tau, target_power, alpha.

Notes

Data mode (Cattaneo, Titiunik and Vazquez-Bare 2019, section 2.2): with the robust bias-corrected side variances V-, V+ scaled by n h^(1+2 deriv), the total sample n* is the smallest integer at least the root of power(tau; sqrt(sum V/h^(1+2 deriv) / n)) = target_power. The treated share of the window sample is rho = sqrt(V+)/(sqrt(V-) + sqrt(V+)); the window sample is M = n* / (rho N+/Nh+ + (1-rho) N-/Nh-) and each side is ceil(rho M) / ceil((1-rho) M). Design mode solves the same power equation with n_right = ratio n_left and the design-mode standard error of :func:rdpower; each side is rounded up.

Examples:

>>> import statspai as sp
>>> r = sp.rdsampsi(tau=0.15, target_power=0.80)
>>> (r.n_left, r.n_right, r.n_total)
(489, 489, 978)

References

cattaneo2019power

rdplot

rdplot(data: DataFrame, y: str, x: str, c: float = 0, nbins: Optional[int] = None, binselect: str = 'esmv', p: int = 4, kernel: str = 'uniform', ci_level: float = 0.95, shade_ci: bool = True, donut: float = 0, show_bw: bool = False, h: Optional[float] = None, covs: Optional[List[str]] = None, weights: Optional[str] = None, cluster: Optional[str] = None, hide_ci: bool = False, scatter: bool = True, ax: Optional[Any] = None, figsize: Tuple[float, float] = (10, 7), title: Optional[str] = None, x_label: Optional[str] = None, y_label: Optional[str] = None) -> Tuple[Any, Any]

RD plot: binned scatter with polynomial fit on each side of the cutoff.

The numbers behind the picture are R rdrobust::rdplot's: the number of bins, bin edges and membership, per-bin means, standard errors and t-based intervals, and the kernel-weighted global polynomial on each side (see statspai.rd._rdplot_core). They are returned on the figure as fig.rdplot_data, a dict with R's vars_bins, vars_poly, coef, J, J_IMSE and J_MV.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome and running variable names.	required
x	str	Outcome and running variable names.	required
c	float	Cutoff.	0
nbins	int or (int, int)	Bins per side. If None, chosen by binselect.	None
binselect	str	Bin selection rule, as in R rdplot: 'es' : IMSE-optimal evenly spaced, spacings estimators 'espr' : IMSE-optimal evenly spaced, polynomial regression 'esmv' : mimicking-variance evenly spaced, spacings (default) 'esmvpr' : mimicking-variance evenly spaced, polynomial regression 'qs', 'qspr', 'qsmv', 'qsmvpr' : the quantile-spaced analogues With 20% or more mass points on either side the spacings variants switch to their polynomial-regression versions (R's masspoints='adjust').	'esmv'
p	int	Order of the global polynomial fitted on each side.	4
kernel	(uniform, triangular, epanechnikov)	Kernel weighting the global polynomial over [c - h, c + h]. R's default. (Through 1.28.0 the default read 'triangular' but the argument was ignored and the fit was unweighted.)	'uniform'
ci_level	float	Level of the per-bin t intervals and the polynomial CI band.	0.95
shade_ci	bool	Shade a pointwise CI band around the polynomial (a StatsPAI addition; not shown when covs is given).	True
donut	float	If > 0, shades the donut region |x - c| <= donut.	0
show_bw	bool	If True, shades the sp.rdrobust bandwidth window.	False
h	float	Support of the global polynomial fit on each side (R's h); defaults to the full range of x on that side.	None
covs	list of str	Covariates, adjusted for as R rdplot(covs=, covs_eval='mean') does.	None
weights	str	Observation weights for the polynomial fit.	None
cluster	str	Cluster column: the per-bin intervals become cluster-robust (Stata mean y, vce(cluster c) within each bin, t(G - 1)), and so does the shaded band around the polynomial (CR1). R rdplot has no cluster option; without it the intervals assume independent observations.	None
hide_ci	bool	If True, suppress all interval displays.	False
scatter	bool	Show binned means.	True
ax	matplotlib Axes		None
figsize	tuple		(10, 7)
title	str		None
x_label	str		None
y_label	str		None

Returns:

Type	Description
(fig, ax)	fig.rdplot_data carries the numbers.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 1000
>>> x = rng.uniform(-1, 1, n)
>>> y = 0.5 * x + 2.0 * (x >= 0) + rng.normal(0, 0.4, n)
>>> df = pd.DataFrame({'y': y, 'x': x})
>>> fig, ax = sp.rdplot(df, y='y', x='x', c=0, title='RD plot')
>>> fig.rdplot_data["J"]  # bins left / right of the cutoff
(23, 23)
>>> fig.savefig('rd_plot.png')

Typical flow: visualise first, then estimate with :func:rdrobust:

>>> result = sp.rdrobust(df, y='y', x='x', c=0)
>>> abs(result.estimate - 2.0) < 0.3
True

rdplotdensity

rdplotdensity(data: DataFrame, x: str, c: float = 0, p: int = 2, n_grid: int = 50, h: Optional[Union[float, Tuple[float, float]]] = None, ci_level: float = 0.95, hist: bool = True, nbins: int = 30, ax: Optional[Any] = None, figsize: Tuple[float, float] = (10, 7), title: Optional[str] = None) -> Tuple[Any, Any]

Density of the running variable on each side of the cutoff.

The Python counterpart of R rddensity::rdplotdensity(rddensity(X), X): on each side, the local-polynomial density estimator of Cattaneo, Jansson & Ma (R lpdensity) is evaluated on an evenly spaced grid over [c - 3 h_l, c + 3 h_r] (clipped to the data) with the rddensity bandwidths, order p for the curve and p + 1 for the robust bias-corrected band, the triangular kernel, the mass-point-adjusted empirical CDF and R's side scaling n_side / (n - 1). The numbers are returned on the figure as fig.rdplotdensity_data (Estl / Estr with R's columns grid, bw, nh, f_p, f_q, se_p, se_q).

Through 1.28.0 this used a rule-of-thumb bandwidth, a per-side empirical CDF that rescaled each curve by n / n_side (roughly doubling both densities), and a heuristic standard error.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
x	str	Running variable.	required
c	float	Cutoff.	0
p	int	Local polynomial order (rddensity's p).	2
n_grid	int	Grid points per side (R's plotN, whose default is 10).	50
h	float or (float, float)	Bandwidths (h_left, h_right); by default those of :func:sp.rddensity at the same c and p.	None
ci_level	float	Level of the robust bias-corrected band f_q ± z se_q.	0.95
hist	bool	Overlay a density-scaled histogram.	True
nbins	int	Histogram bins per side.	30
ax	matplotlib Axes		None
figsize	tuple		(10, 7)
title	str		None

Returns:

Type	Description
(fig, ax)	fig.rdplotdensity_data carries the numbers.

References

cattaneo2020simple

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> df = pd.DataFrame({"x": rng.uniform(-1, 1, 500)})
>>> fig, ax = sp.rdplotdensity(df, x="x", c=0)
>>> est = fig.rdplotdensity_data["Estl"]
>>> bool(abs(est["f_p"][-1] - 0.5) < 0.2)  # U(-1, 1) has density 1/2
True

fit

fit(data: Any = None, y: Any = None, x: Any = None, c: Any = 0, *, method: str = 'rdrobust', **kwargs: Any) -> Any

Alias for :func:_rd_dispatch. See sp.rd.__doc__ for usage.