The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

LRErdd: Regression Discontinuity Designs as Local Randomized Experiments

Ibon Tamayo, Alessandra Mattei, Fabrizia Mealli and Marie-Abele Bind

2026-07-19

Introduction

The LRErdd package provides a set of functions for the design and analysis of Regression Discontinuity (RD) Designs as local randomized experiments within the potential outcome approach, as formalized in Li, Mattei and Mealli (2015).

The key insight underlying this approach is that, under a local randomization assumption, valid causal inference can be drawn for subpopulations of treated and control units whose realized values of the forcing variable fall in a neighborhood of the threshold and who have similar background characteristics.

The package supports the two phases of an RD study:

There are two ways of working with the package:

  1. Programmatically in R, through the RegressionDiscontinuityClass R6 class, which encapsulates the whole workflow. This option is illustrated in this vignette and is the recommended one for larger datasets and reproducible analyses.
  2. Interactively, through a user-friendly Shiny application shipped with the package (see the last section of this vignette).

Installation

install.packages("LRErdd")
library(LRErdd)

The case study: Italian university grants and student dropout

We illustrate the workflow with data from an Italian study concerning the effects of university grants on student dropout (Li, Mattei and Mealli, 2015). Eligibility for the grant is determined by a measure of the economic situation of the student’s family (S, the forcing variable) falling below a threshold of 15,000 euros, which makes it a natural regression discontinuity design. Because not all eligible students applied for and received the grant, the design is fuzzy: the treatment received (W) does not coincide with the eligibility status (Z).

data(grants)
dim(grants)
#> [1] 15984    21
head(grants[, c("S", "Z", "A", "W", "dropout", "sex", "hsgrade")])
#>          S Z A W dropout sex hsgrade
#> 1 16833.09 0 1 0       0   1      81
#> 2  4002.53 1 1 1       0   1      70
#> 3 16905.84 0 1 0       0   1      70
#> 4 11836.94 1 1 1       0   1      60
#> 5 11935.00 1 1 1       0   0      70
#> 6  2015.45 1 1 1       1   0      70

The dataset contains information on 15,984 students. The main variables are the forcing variable S, the eligibility indicator Z, the grant application indicator A, the treatment received W, and the outcome dropout (1 if the student dropped out). The remaining columns are pre-treatment covariates (sex, high school type and grade, cohort, university, and field of study).

Creating a RegressionDiscontinuityClass object

The whole analysis is driven by an object of class RegressionDiscontinuityClass. To create it, we need the dataset (a data frame), the name of the forcing variable, the threshold, whether the treatment group is below ("lower") or above ("upper") the threshold, and the vector of covariate names used to assess balance.

# Covariates of interest
cov <- c("HSTech", "hsgrade", "Y2005")

rdo <- RegressionDiscontinuityClass$new(
  data       = grants,
  forcing    = "S",
  threshold  = 15000,
  Z1S1       = "lower",
  covariates = cov
)

The full_print() method summarizes the information stored in the object:

rdo$full_print()
#> [1] "------------------------------------"
#>  [1] "Names of data frame -> S"             
#>  [2] "Names of data frame -> Z"             
#>  [3] "Names of data frame -> A"             
#>  [4] "Names of data frame -> W"             
#>  [5] "Names of data frame -> dropout"       
#>  [6] "Names of data frame -> sex"           
#>  [7] "Names of data frame -> HSHumanity"    
#>  [8] "Names of data frame -> HSScience"     
#>  [9] "Names of data frame -> HSTech"        
#> [10] "Names of data frame -> HSOther"       
#> [11] "Names of data frame -> hsgrade"       
#> [12] "Names of data frame -> Y2004"         
#> [13] "Names of data frame -> Y2005"         
#> [14] "Names of data frame -> Y2006"         
#> [15] "Names of data frame -> University"    
#> [16] "Names of data frame -> Humanity"      
#> [17] "Names of data frame -> Science"       
#> [18] "Names of data frame -> Social.Science"
#> [19] "Names of data frame -> BioMed"        
#> [20] "Names of data frame -> Tech"          
#> [21] "Names of data frame -> Other"         
#> [1] "Forcing variable -> S"
#> [1] "Threshold -> 15000"
#> [1] "Z1S1 -> lower"
#> [1] "Covariate -> HSTech"  "Covariate -> hsgrade" "Covariate -> Y2005"  
#> [1] "------------------------------------"

Design phase

Summary statistics

The summary_statistics_mean() method reports the mean of each selected covariate, overall and by treatment assignment status, together with the group sizes:

rdo$summary_statistics_mean()
#>      Name      All     Z=0      Z=1
#> 1  HSTech     0.39    0.36     0.40
#> 2 hsgrade    81.25   81.94    81.00
#> 3   Y2005     0.34    0.36     0.34
#> 4       N 15984.00 4281.00 11703.00

Bandwidth selection

The bandwidth_selection() method evaluates the influence of a candidate bandwidth around the threshold on the balance of the covariates. Three ways of defining the range are available through the typerange argument:

The method returns the number of selected records with the limits of the selected range, a table with initial and bandwidth-specific p-values (raw and adjusted for multiplicity) from randomization-based balance tests, and, optionally, a histogram of the forcing variable with the threshold (black) and the selected limits (red):

rdo$bandwidth_selection(
  typerange   = "percentage",
  range_value = 25,
  num_it      = 100,
  plot        = TRUE
)
#> $RangeInfo
#>   Number of records Down limit Up limit Difference
#> 1              7584   11250.25 18749.75     7499.5
#> 
#> $pvalues
#>   Variable Initial Pvalue Initial Adj pvalue Bandwidth specific Pvalue
#> 1   HSTech           0.00               0.03                      0.27
#> 2  hsgrade           0.00               0.03                      0.96
#> 3    Y2005           0.03               0.09                      0.05
#>   Bandwidth specific Adj pvalue
#> 1                          0.64
#> 2                          1.00
#> 3                          0.15
#> 
#> $RangePlot

Distribution of a covariate before and after the selection

The distribution_plot() method compares the distribution of a covariate between treatment groups in the complete dataset and within the selected bandwidth. The covariate can be binary (typecov = "binary") or continuous (typecov = "conti"):

rdo$distribution_plot(
  typerange   = "percentage",
  range_value = 25,
  num_it      = 100,
  covariate   = "sex",
  typecov     = "binary"
)

Comparing several candidate bandwidths

The summary_bandwidth() method implements randomization-based tests with multiplicity adjustment for assessing balance in the background variables between treatment groups, for several candidate bandwidths at once. It returns a table of adjusted p-values by covariate and bandwidth; large p-values indicate subpopulations where the covariates are well balanced:

rdo$summary_bandwidth(buffers = c(500, 1000, 1500), num_it = 200)
#>     Names p-value[Bandwidth=500] p-value[Bandwidth=1000]
#> 1  HSTech                   1.00                   0.900
#> 2 hsgrade                   0.58                   1.000
#> 3   Y2005                   1.00                   0.385
#>   p-value[Bandwidth=1500]
#> 1                   1.000
#> 2                   0.955
#> 3                   0.830

In practice we recommend a larger number of iterations (e.g. num_it = 1000 or more) for the final analysis; a smaller value is used here to keep the vignette light.

Analysis phase: inference on causal effects

The causal_effect() method is the main function of the package. Four modes of inference are implemented through the method argument:

In all cases, the user selects the bandwidths (buffers), the outcome and its type ("binary" or "numeric"), and the number of iterations used to approximate the randomization distribution.

Sharp FEP

Under the sharp interpretation, we test the sharp null hypothesis of no treatment effect for the eligibility status. The method reports, for each bandwidth, the observed difference in average outcomes by assignment status and the Fisher exact p-value, together with a histogram of the randomization distribution of the test statistic (the dashed line marks the observed value):

rdo$causal_effect(
  method      = "Sharp FEP",
  buffers     = c(500, 1000, 1500),
  num_it      = 100,
  outcome     = "dropout",
  typeoutcome = "binary",
  plot        = TRUE
)
#> $table
#>   Bandwidth    N Difference in average outcomes by treatment status\nStatistic
#> 1       500  528                                                       -0.0595
#> 2      1000 1042                                                       -0.0260
#> 3      1500 1577                                                       -0.0656
#>   Absolute value of difference in average outcomes p-value
#> 1                                           0.0595   0.161
#> 2                                           0.0260   0.386
#> 3                                           0.0656   0.009
#> 
#> $plot

Sharp Neyman

The Neyman approach estimates the average causal effect of eligibility with its standard error and a confidence interval at the level defined by cin (the alpha level, 0.05 by default):

rdo$causal_effect(
  method      = "Sharp Neyman",
  buffers     = c(500, 1000, 1500),
  num_it      = 100,
  outcome     = "dropout",
  typeoutcome = "binary",
  cin         = 0.05,
  plot        = TRUE
)
#> $table
#>   Bandwidth    N     ACE Neyman SE       L       U
#> 1       500  528 -0.0595    0.0423 -0.1425  0.0235
#> 2      1000 1042 -0.0260    0.0302 -0.0853  0.0332
#> 3      1500 1577 -0.0656    0.0246 -0.1139 -0.0173
#> 
#> $plot

Fuzzy Neyman

In the fuzzy design, the treatment received (W, provided through the treatm_cov argument) does not coincide with the assignment. The Neyman approach reports both the intention-to-treat effect and the complier average causal effect (CACE):

rdo$causal_effect(
  method      = "Fuzzy Neyman",
  buffers     = c(500, 1000, 1500),
  num_it      = 100,
  outcome     = "dropout",
  typeoutcome = "binary",
  treatm_cov  = "W",
  plot        = TRUE
)
#> $table
#>    Bandwidth Estimand     ACE Neyman SE       L       U
#> 2        500    ITT.W  0.3794    0.0289  0.3228  0.4361
#> 3        500    ITT.Y -0.0595    0.0423 -0.1425  0.0235
#> 4        500     CACE -0.1568    0.1119 -0.3761  0.0625
#> 5       1000    ITT.W  0.3932    0.0206  0.3529  0.4336
#> 6       1000    ITT.Y -0.0260    0.0302 -0.0853  0.0332
#> 7       1000     CACE -0.0662    0.0769 -0.2170  0.0846
#> 8       1500    ITT.W  0.3778    0.0165  0.3454  0.4101
#> 9       1500    ITT.Y -0.0656    0.0246 -0.1139 -0.0173
#> 10      1500     CACE -0.1736    0.0657 -0.3025 -0.0448
#> 
#> $plot

Fuzzy FEP

The Fuzzy FEP approach tests the null hypothesis of no effect for compliers. Noncompliance can be one-sided (sided = "onesided") or two-sided (sided = "twosided"). Three test statistics are used: the IV estimate, the maximum likelihood estimate (MLE), and the posterior median of the CACE. The num_it_fuzzy argument controls the number of iterations used to compute the MLE/posterior mode of the complier average causal effects:

rdo$causal_effect(
  method       = "Fuzzy FEP",
  buffers      = c(500, 1000, 1500),
  num_it       = 100,
  outcome      = "dropout",
  typeoutcome  = "binary",
  sided        = "twosided",
  treatm_cov   = "W",
  num_it_fuzzy = 5,
  plot         = TRUE
)
#> $table
#> $table[[1]]
#>   Bandwidth Statistic: IV estimate of CACE Statistic: MLE of CACE
#> 2       500                            0.0                      0
#> 3      1000                            0.2                      0
#> 4      1500                            0.0                      0
#>   Statistic: Posterior median of CACE
#> 2                                 0.0
#> 3                                 0.2
#> 4                                 0.0
#> 
#> $table[[2]]
#>         CACE.IV     CACE.MLE       CACE.PM
#> 2   0.121401861  0.120962946  0.1085443334
#> 3  -0.076923077 -0.077347660 -0.0509405694
#> 4  -0.020833333 -0.021235691 -0.0184509592
#> 5  -0.058316374 -0.057989495 -0.0746075755
#> 6   0.040497137  0.040750630  0.0373806480
#> 7   0.083178864  0.082580819  0.0729336064
#> 8   0.024202655  0.024636647  0.0199358256
#> 9   0.001076548  0.001474292 -0.0052912478
#> 10  0.026738593  0.027119947  0.0169827540
#> 11 -0.025550824 -0.025180402 -0.0230312374
#> 12  0.039389529  0.039074604  0.0363220111
#> 13  0.021868573  0.021503205  0.0290408361
#> 14 -0.068539020 -0.068833109 -0.0694434905
#> 15  0.004873657  0.005239995 -0.0005615341
#> 16 -0.019087523 -0.018808781 -0.0175841706
#> 
#> $table[[3]]
#>         [,1]       [,2]        [,3]
#>  -0.15682699 -0.1564848 -0.13780520
#>  -0.06619532 -0.0658559 -0.09382708
#>  -0.17363510 -0.1733321 -0.17730450
#> 
#> 
#> $plot

Note that the Fuzzy FEP approach is computationally more demanding than the other three; for final analyses with a large number of iterations, expect longer running times.

The Shiny application

All the functionality described above is also available through a user-friendly Shiny application. It can be launched with:

open_LRErdd_framework()

The dashboard presents two sections: the input options on the left panel, managed through a drop-down list covering all the steps of the workflow, and the results on the right, in the result window.

Workflow of the LRErdd Shiny application
Workflow of the LRErdd Shiny application

The steps mirror the programmatic workflow:

  1. Load data. The user can choose the grants dataset included in the package, or upload a dataset from an Excel (.xlsx) or R (.RData) file.
  2. Define settings. The forcing variable, the threshold, and whether the treatment group is below or above the threshold.
  3. Summary statistics. Means of the selected variables by treatment assignment status.
  4. Bandwidth selection. Interactive selection of the bandwidth (as a percentage of the range, symmetric, or asymmetric), with balance tables and plots updated dynamically.
  5. Summary bandwidth selection. Randomization-based balance tests with multiplicity adjustment for several candidate bandwidths.
  6. Inference on causal effects. The four modes of inference (Sharp FEP, Sharp Neyman, Fuzzy FEP, Fuzzy Neyman) with their specific settings.

References

Li, F., Mattei, A. and Mealli, F. (2015). Bayesian inference for regression discontinuity designs with application to the evaluation of Italian university grants. The Annals of Applied Statistics, 9(4), 1906–1931. doi:10.1214/15-AOAS809

Session info

sessionInfo()
#> R version 4.3.3 (2024-02-29)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
#> 
#> locale:
#> [1] C
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] LRErdd_0.1.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.4       jsonlite_1.8.8     dplyr_1.1.4        compiler_4.3.3    
#>  [5] highr_0.10         gtools_3.9.5       tidyselect_1.2.0   jquerylib_0.1.4   
#>  [9] scales_1.3.0       yaml_2.3.8         fastmap_1.1.1      ggplot2_3.4.4     
#> [13] R6_2.5.1           labeling_0.4.3     generics_0.1.3     knitr_1.45        
#> [17] tibble_3.2.1       munsell_0.5.0      bslib_0.6.1        pillar_1.9.0      
#> [21] RColorBrewer_1.1-3 rlang_1.1.3        utf8_1.2.4         cachem_1.0.8      
#> [25] xfun_0.41          sass_0.4.8         cli_3.6.2          withr_2.5.0       
#> [29] magrittr_2.0.3     digest_0.6.34      grid_4.3.3         cowplot_1.1.3     
#> [33] lifecycle_1.0.4    vctrs_0.6.5        evaluate_0.23      glue_1.7.0        
#> [37] farver_2.1.1       fansi_1.0.5        colorspace_2.1-0   rmarkdown_2.25    
#> [41] tools_4.3.3        pkgconfig_2.0.3    htmltools_0.5.7

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.