---
title: "LRErdd: Regression Discontinuity Designs as Local Randomized Experiments"
author: "Ibon Tamayo, Alessandra Mattei, Fabrizia Mealli and Marie-Abele Bind"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{LRErdd: Regression Discontinuity Designs as Local Randomized Experiments}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  warning = FALSE,
  message = FALSE
)
set.seed(27)
```

## 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:

* **Design phase.** A subset of functions focuses on the selection of suitable
  subpopulations around the threshold for which valid causal inference can be
  drawn. These functions provide summary statistics of pre- and post-treatment
  variables by treatment status, and select subpopulations where pre-treatment
  variables are well balanced between treatment groups, using
  randomization-based tests with adjustment for multiplicities. Functions for
  a visual inspection of the results are also provided.

* **Analysis phase.** A second set of functions draws inference on causal
  effects for the selected subpopulations using randomization-based modes of
  inference. Specifically, the Fisher Exact p-value (FEP) and Neyman
  approaches are implemented for the analysis of both *sharp* and *fuzzy* RD
  designs.

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

```{r install, eval = FALSE}
install.packages("LRErdd")
```

```{r load}
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`).

```{r data}
data(grants)
dim(grants)
head(grants[, c("S", "Z", "A", "W", "dropout", "sex", "hsgrade")])
```

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.

```{r create}
# 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:

```{r fullprint}
rdo$full_print()
```

## 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:

```{r summarystats}
rdo$summary_statistics_mean()
```

### 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:

* `"percentage"`: considering the range of the forcing variable as 100%, the
  bandwidth is defined as a percentage of that range.
* `"balunit"`: a symmetric bandwidth is defined by its total width in the
  units of the forcing variable.
* `"unbalunit"`: an asymmetric bandwidth, defined by the lower and upper
  bounds provided in `range_value`.

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):

```{r bandwidth}
rdo$bandwidth_selection(
  typerange   = "percentage",
  range_value = 25,
  num_it      = 100,
  plot        = TRUE
)
```

### 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"`):

```{r distplot}
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:

```{r summarybw}
rdo$summary_bandwidth(buffers = c(500, 1000, 1500), num_it = 200)
```

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:

* `'Sharp FEP'`: Fisher's exact p-value approach for sharp RD designs within
  the selected subpopulations.
* `'Sharp Neyman'`: Neyman approach for sharp RD designs.
* `'Fuzzy FEP'`: Fisher's exact p-value approach for fuzzy RD designs.
* `'Fuzzy Neyman'`: Neyman approach for fuzzy RD designs.

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):

```{r sharpfep}
rdo$causal_effect(
  method      = "Sharp FEP",
  buffers     = c(500, 1000, 1500),
  num_it      = 100,
  outcome     = "dropout",
  typeoutcome = "binary",
  plot        = TRUE
)
```

### 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):

```{r sharpneyman}
rdo$causal_effect(
  method      = "Sharp Neyman",
  buffers     = c(500, 1000, 1500),
  num_it      = 100,
  outcome     = "dropout",
  typeoutcome = "binary",
  cin         = 0.05,
  plot        = TRUE
)
```

### 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):

```{r fuzzyneyman}
rdo$causal_effect(
  method      = "Fuzzy Neyman",
  buffers     = c(500, 1000, 1500),
  num_it      = 100,
  outcome     = "dropout",
  typeoutcome = "binary",
  treatm_cov  = "W",
  plot        = TRUE
)
```

### 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:

```{r fuzzyfep}
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
)
```

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:

```{r shiny, eval = FALSE}
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](Squeme.jpg){width=100%}

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](https://doi.org/10.1214/15-AOAS809)

## Session info

```{r session}
sessionInfo()
```