day15-NaturalExperimentsRDD.Rmd

---
title: |
 | Natural Experiments and Discontinuities
 | ICPSR 2021 Session 1
date: '`r format(Sys.Date(), "%B %d, %Y")`'
author: |
  | Jake Bowers, Ben Hansen, Tom Leavitt
author-meta: Jake Bowers, Ben Hansen, Tom Leavitt
bibliography:
 - 'BIB/MasterBibliography.bib'
fontsize: 10pt
geometry: margin=1in
graphics: yes
biblio-style: authoryear-comp
colorlinks: true
biblatexoptions:
  - natbib=true
output:
  beamer_presentation:
    slide_level: 2
    keep_tex: true
    latex_engine: xelatex
    citation_package: biblatex
    template: icpsr.beamer
    incremental: true
    includes:
        in_header:
           - defs-all.sty
    md_extensions: +raw_attribute-tex_math_single_backslash+autolink_bare_uris+ascii_identifiers+tex_math_dollars
    pandoc_args: [ "--csl", "chicago-author-date.csl" ]
---


<!-- To show notes  -->
<!-- https://stackoverflow.com/questions/44906264/add-speaker-notes-to-beamer-presentations-using-rmarkdown -->

```{r setup1_env, echo=FALSE, include=FALSE}
library(here)
source(here::here("rmd_setup.R"))
opts_chunk$set(echo = TRUE, digits = 4)
```

```{r setup2_loadlibs, echo=FALSE, include=FALSE}
## Load all of the libraries that we will use when we compile this file
## We are using the renv system. So these will all be loaded from a local library directory
library(dplyr)
library(ggplot2)
library(coin)
library(RItools)
library(optmatch)
library(estimatr)
library(dplyr)
library(MASS)
library(magrittr)
library(tidyr)
library(haven)
library(parallel)
library(rdd)
```

## Today

 1. Agenda: Natural experiments and discontinuities (with and without regression).
 3. Questions arising from the reading or assignments or life?

## So far:

- Q: How can we argue that the relationship between $Z$ and $Y$ is $Z
  \rightarrow Y$? 
- A: By excluding alternative explanations for this relationship.
- Q: How do we exclude alternative explanations?
- A: Many approaches
  1. Randomize $Z$;
  2. Randomize an instrument for $Z$;
  3. Remove specific $X \rightarrow Z$ relationships via stratification and/or weighting.
  4. (Today:) use a discontinuity as a kind of instrument for $Z$.

# Discontinuity Designs

## General Introduction and Setup

Today we are going to use the data on close US House of Representatives races
1942--2008 used in  @caugheysekhon2011.^[The full replication data is
available for download
\href{http://sekhon.berkeley.edu/rep/RDReplication.zip}{here}. But we read it
directly below.] @caugheysekhon2011 engage in a debate whose
participants seek to identify the causal effect of the so-called "incumbency
advantage." That is, what effect does a candidate's status as an incumbent
have on whether or not that candidate wins an election? Obviously, whether or
not a candidate is an incumbent is *not* randomly assigned.

Let's first load the data:

```{r}
rdd_data <- read_dta("http://jakebowers.org/Matching/RDReplication.dta") %>%
  filter(Use == 1) ## Use is indicator for whether unit is included in RD incumbency advantage sample
```

## The Running Variable

The "running variable" is called `DifDPct`, which is defined as the Democratic margin of victory or defeat in the election; in other words, DifDPct is the difference between the percentage of all votes that were cast for the leading Democrat in the race and the percentage cast for the leading non-Democrat. Races in which no Democrat ran or in which the top two vote-getters were both Democrats are coded as missing.

```{r}
running_var <- matrix(c("DifDPct", "Democrat Margin of Victory"),
  ncol = 2,
  byrow = TRUE
)

dimnames(running_var) <- list(1, c("Running Variable", "Description"))

kable(running_var)
```

## The Treatment Variable

The treatment variable is whether or not the Democratic candidate wins the election or not. If the candidate wins the election, then that candidate is assigned to "treatment." If the candidate loses the election, then he or she is assigned to "control.''

```{r}
treatment <- matrix(c("DemWin", "Democrat Wins Election"),
  ncol = 2,
  byrow = TRUE
)
dimnames(treatment) <- list(1, c("Treatment", "Description"))
kable(treatment)
```

Now let's quickly look at the empirical distribution of the treatment variable:

```{r}
table(rdd_data$DemWin)
```

## Outcome Variables

In @caugheysekhon2011, the primary outcome variables of interest are as
follows: whether a democrat wins the next election, the proportion voting for
a democrat in the next election, and the democratic vote margin in the next
election.

```{r echo=FALSE}
dvs <- matrix(c(
  "DWinNxt", "Dem Win t + 1",
  "DPctNxt", "Dem % t + 1",
  "DifDPNxt", "Dem % Margin t + 1"
),
ncol = 2,
byrow = TRUE
)

dimnames(dvs) <- list(
  seq(
    from = 1,
    to = 3,
    by = 1
  ),
  c("Outcome", "Description")
)

kable(dvs)
```

## Outcomes and Running Variables

Effect of a democrat winning on democrat votes in next election? Alternative explanations? Reasonable comparisons?

```{r echo=FALSE, out.width=".7\\textwidth", warnings=FALSE}
g <- ggplot(rdd_data,aes(x=DifDPct, y=DPctNxt,color=DemWin)) +
    geom_point()
print(g)
```


## Baseline Covariates

The relevant baseline covariates (all measured prior to the realization of the running variable) are:

```{r echo=FALSE}
covs <- matrix(c(
  "DWinPrv", "Dem Win t - 1",
  "DPctPrv", "Dem % t - 1",
  "DifDPPrv", "Dem % Margin t - 1",
  "IncDWNOM1", "Inc's D1 NOMINATE",
  "DemInc", "Dem Inc in Race",
  "NonDInc", "Rep Inc in Race",
  "PrvTrmsD", "Dem's # Prev Terms",
  "PrvTrmsO", "Rep's # Prev Terms",
  "RExpAdv", "Rep Experience Adv",
  "DExpAdv", "Dem Experience Adv",
  "ElcSwing", "Partisan Swing",
  "CQRating3", "CQ Rating {-1, 0, 1}",
  "DSpndPct", "Dem Spending %",
  "DDonaPct", "Dem Donation %",
  "SoSDem", "Dem Sec of State",
  "GovDem", "Dem Governor",
  "DifPVDec", "Dem Pres % Margin", ## average over decade
  "DemOpen", "Dem-held Open Seat",
  "NonDOpen", "Rep-held Open Seat",
  "OpenSeat", "Open Seat",
  "VtTotPct", "Voter Turnout %",
  "GovWkPct", "Pct Gov't Worker",
  "UrbanPct", "Pct Urban",
  "BlackPct", "Pct Black",
  "ForgnPct", "Pct Foreign Born"
),
ncol = 2,
byrow = TRUE
)

dimnames(covs) <- list(
  seq(
    from = 1,
    to = 25,
    by = 1
  ),
  c("Covariate", "Description")
)

kable(covs)
```

## Local Randomization Framework

Let the index $i \in \left\{1, \dots , n\right\}$ run over the $n$ experimental units. In the context of a regression discontinuity design, let $R_i$ denote the random "score variable" (also known as "running variable''). The random assignment variable, $Z_i \in \left\{0, 1\right\}$, which indicates whether subject $i$ is assigned to treatment, $Z_i = 1$, or to control, $Z_i = 0$, is a deterministic function of $R_i$. In the context of the incumbency advantage study, we can define $Z_i$ as follows:
\begin{align*}
Z_i \equiv \mathbbm{I}\left[R_i > 0 \right],
\end{align*}
where $\mathbbm{I}\left[\cdot\right]$ is an indicator function that is $1$ if the argument to the function is true and $0$ if false. In this particular study, if the margin of victory is greater than $0$, then candidate $i$ is treated ($Z_i = 1$), and if not, then candidate $i$ is assigned to control ($Z_i = 0$).

\medskip

Let $W_0 = [\underline{r}, \overline{r}]$, where $\underline{r} < r_0 < \overline{r}$, denote the window around the cutpoint (or threshold value), $r_0$, that sorts units into treatment or control.

\medskip

How to choose the window?

## Optimal Bandwidth Selection {.allowframebreaks}

We know from @hansensales2015, @rosenbaum2008, @berger1988 that "if a researcher
pre-specifies a sequence of hypotheses and corresponding level-$\alpha$ tests,
tests those hypotheses in order, and stops testing after the first non-rejected
hypothesis, then the probability of incorrectly rejecting at least one correct
hypothesis is at most $\alpha$" \citep[p. 185]{hansensales2015}.
\smallskip

As applied to bandwidth selection in the RDD context, the SIUP implies that one should start with a set of candidate bandwidths and sequentially test for covariate balance (beginning from either the largest candidate bandwidth or the smallest candidate bandwidth).
\smallskip

Let's specify a set of candidate bandwidths and then sequentially test covariate balance. Before actually testing, though, we want to specify a balance criterion and then maximize effective sample size subject to that criterion.

```{r}

bal_fmla <- reformulate(covs[1:25], response = "DemWin")

candidate_bands <- seq(
  from = -5,
  to = 5,
  by = .1
)
```

Now let's first filter our dataset and check for balance in the largest candidate bandwidth spanning from $-5$ to $5$.

```{r}
lower_bound <- seq(from = -5, to = -0.1, by = 0.1)

upper_bound <- seq(from = 0.1, to = 5, by = 0.1) %>%
  sort(decreasing = TRUE)

rdd_dataA <- rdd_data
rdd_dataA %<>% filter(DifDPct > lower_bound[1] & DifDPct < upper_bound[1])

rdd_dataA %$% summary(DifDPct)

xBalance(
  fmla = bal_fmla,
  data = rdd_dataA,
  report = "chisquare.test"
)

xb1 <- xBalance(
  fmla = bal_fmla,
  data = rdd_dataA,
  report = "all"
)
xb1$results
```

Now let's write a function to perform this same procedure over all candidate
bandwidth sizes beginning with the largerst candidate bandwidth and
subsequently testing smaller and smaller bandwidths in order.

```{r}

chi_squared_balance <- function(i,
                                running_var,
                                bal_fmla,
                                data) {

  # Preliminaries
  suppressMessages(stopifnot(require(dplyr, quietly = TRUE)))
  suppressMessages(stopifnot(require(parallel, quietly = TRUE)))
  suppressMessages(stopifnot(require(magrittr, quietly = TRUE)))
  suppressMessages(stopifnot(require(RItools, quietly = TRUE)))

  lower_bound <- seq(from = -5, to = -0.1, by = 0.1)

  upper_bound <- seq(from = 0.1, to = 5, by = 0.1) %>%
    sort(decreasing = TRUE)

  data %<>% filter(running_var > lower_bound[i] & running_var < upper_bound[i])

  # Effective Sample Size
  ess <- nrow(data)

  p_value <- xBalance(
    fmla = bal_fmla,
    data = data,
    report = "chisquare.test"
  )$overall[[3]]

  bands <- cbind(ess, p_value, lower_bound[i], upper_bound[i])

  return(bands)
}

is <- seq(
  from = 1,
  to = length(seq(
    from = -5,
    to = -0.1,
    by = 0.1
  )),
  by = 1
)
```

Now use the function:

```{r}
cl <- makeCluster(parallel::detectCores())

band_df <- data.frame(t(parSapply(cl, is,
  chi_squared_balance,
  running_var = rdd_data$DifDPct,
  bal_fmla = bal_fmla,
  data = rdd_data
))) %>%
  rename(
    ess = X1,
    p_value = X2
  )

parallel::stopCluster(cl)

kable(band_df)
```

And look at the results

```{r}

band_df$limits <- paste(round(band_df$X3,1),round(band_df$X4,1),sep=",")
g2 <- ggplot(band_df,aes(x=p_value,y=ess,label=limits))+
    geom_point() +
    geom_label()
print(g2)
```

```{r}
band_df[47,c("ess","p_value","X3","X4")]
```


```{r}

g + xlim(c(-.4,.4))

```


## Outcome Analysis and the Exclusion Restriction {.allowframebreaks}

Using the discontinuity as an instrument brings us back to instrumental variables.

\medskip

Let's assume that within a certain bandwidth around the cutpoint, $W_0$,
the assumption of a local randomized experiment obtains. In other words,
probabilties are uniformly distributed on $R_i \in W_0$ for all $i$. Yet even
if this assumption obtains, the running variable, $R$, might still relate to
potential outcomes through a mechanism other than $Z \given R \in W_0$.

\medskip

The exclusion restriction is a strong assumption in the RDD context since $Z$
is a deterministic function of the running variable, $R$, which implies that
treatment and control groups, by construction, will be imbalanced on the
running variable. Thus, if $R$ relates to $(y_{c}, y_{t})$, through a mechanism
other than $Z \given R \in W_0$, then the exlcusion restriction is violated.

\medskip

If we think the exclusion restriction holds, then we can simply perform outcome
analysis within a window around the cutpoint where local randomization
presumably holds:

```{r}

rdd_data47 <-  rdd_data %>% filter(DifDPct > lower_bound[47] & DifDPct < upper_bound[47])

rdd_data47 %>% nrow()

lm_robust(formula = DPctNxt ~ DemWin, data = rdd_data47)
```

We coud also perform permutation inference within the window around the cutpoint under the assumption that local randomization holds (useful especially if the window contains relatively few observations, skewed outcomes, etc..):

```{r}

set.seed(1:5)
sharp_null_dist <- replicate(
  n = 10^3,
  expr = coef(lm(formula = DPctNxt ~ sample(DemWin), data = rdd_data47))[[2]]
)

obsstat <- coef(lm(formula = DPctNxt ~ DemWin, data = rdd_data47))[["DemWin"]]

upper_p_val <- mean(sharp_null_dist >= obsstat)
lower_p_val <- mean(sharp_null_dist <= obsstat)
two_sided_p_val <- min(1, 2 * min(upper_p_val, lower_p_val))

two_sided_p_val

rdd_data47$DemWinF <- factor(rdd_data47$DemWin)

## Alternatively
wctest <- wilcox_test(DPctNxt~DemWinF,data=rdd_data47,distribution=approximate(nresample=1000))
ttest <- oneway_test(DPctNxt~DemWinF,data=rdd_data47,distribution=approximate(nresample=1000))
pvalue(wctest)
pvalue(ttest)
```

## Outcome Analysis and the Exclusion Restriction {.allowframebreaks}

What can we do if the exclusion restriction is violated? Or worse, if the
running variable, $R$, confounds the relationship between $Z$ and $Y$? One
approach is to model potential outcomes as a function of the running variable,
and then to "de-trend" (or "transform'') the outcome variable and to
subsequently make the claim of "residual ignorability."

For example, we could do the following:
```{r}
tmp_lm <- lm(DPctNxt ~ DifDPct, data = rdd_data47)

rdd_data47 %<>% mutate(resid_DPctNxt = resid(tmp_lm))
```

```{r}
## An influential point robust loess smoother.
ggplot(
  data = rdd_data47,
  mapping = aes(x = DifDPct, y = DPctNxt)
) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE) +
  geom_smooth(method = rlm, se = FALSE, colour = "black") +
  geom_smooth(
    method = loess, se = FALSE, colour = "red",
    method.args = list(deg = 1, span = 1 / 3, family = "symmetric")
  ) +
  xlab("Democratic Margin of Victory (Running Variable)") +
  ylab("Democratic Vote Percentage in Next Election (Outcome Variable)") +
  ggtitle("Relationship between Running Variable and Outcome") +
  geom_vline(xintercept = 0, linetype = "dashed")
```

## Outcome Analysis

We can now perform outcome analysis on the detrended outcome variable:

```{r}
lm_robust(formula = resid_DPctNxt ~ DemWin, data = rdd_data47)
```

@saleshansen2020, however, propose robust regression, which is less sensitive to violations of the regression model's assumptions --- still removing linear trend here.

```{r}
tmp_rlm <- rlm(DPctNxt ~ DifDPct, data = rdd_data47)
rdd_data47 %<>% mutate(rlm_resid_DPctNxt = resid(tmp_rlm))
lm_robust(formula = rlm_resid_DPctNxt ~ DemWin, data = rdd_data47)

wctest2 <- wilcox_test(rlm_resid_DPctNxt~DemWinF,data=rdd_data47,distribution=approximate(nresample=1000))
ttest2<- oneway_test(rlm_resid_DPctNxt~DemWinF,data=rdd_data47,distribution=approximate(nresample=1000))
pvalue(wctest2)
pvalue(ttest2)

```

## Matching to create variable "bandwidths"

```{r}
mdist <- match_on(bal_fmla,data=rdd_data)
pm1 <- pairmatch(mdist,data=rdd_data,remove.unmatchables = TRUE)
summary(pm1)
```

## Continuity in Potential Outcomes Framework {.allowframebreaks}

An alternative approach to regression discontinuity is common in the econometrics literature \citep[see, e.g.,][]{imbenslemieux2008,lee2008,calonicoetal2014,hahnetal2001,leelemieux2010,mccrary2008}. As @hahnetal2001 and @imbenslemieux2008 state, the estimand (i.e., the quantity we seek to estimate) is:
\begin{align*}
\lim \limits_{r \downarrow 0}\E\left[Y \given R = r > 0\right] - \lim \limits_{r \uparrow 0}\E\left[Y \given R = r < 0 \right] &  \\
& = \E\left[Y_t \given R = 0\right] - \E\left[Y_c \given R = 0 \right]  \\
& = \E\left[Y_t - Y_c \given R = 0 \right]
\end{align*}

This estimand assumes that the (1) left and right limits of $\E\left[Y \given R \right]$ exist as $R$ approaches the cutpoint $r = 0$ and (2) conditional expectation of the outcome given the running variable is left-continuous and right-continuous at the cutpoint $r = 0$, where in general a function $f$ is continuous from the right at $a$ if $\lim \limits_{x \downarrow a} f(x) = f(a)$ and is continuous from the left if $\lim \limits_{x \uparrow a} f(x) = f(a)$.


## Continuity in Potential Outcomes Framework {.allowframebreaks}

We don't actually need to know the true form of $\E[\left[Y \given R =
r\right]$. Instead, we can use methods, such as local regression, to
approximate the function of $\E[\left[Y \given R = r\right]$ at values of $R$
below and above the cutpoint.


```{r}
rdd_data %<>% arrange(DifDPct)

rdd_data_cp <- dplyr::select(.data = rdd_data, DifDPct, DPctNxt, DemWin) %>%
  mutate(cp = ifelse(test = DifDPct < 0, yes = "Below Cutpoint",
    no = "Above Cutpoint"
  ))
```

```{r echo=FALSE}
ggplot(
  data = rdd_data_cp,
  mapping = aes(
    x = DifDPct,
    y = DPctNxt
  )
) +
  geom_point() +
  geom_smooth(method = lm, se = FALSE) +
  geom_smooth(method = rlm, se = FALSE, colour = "black") +
  geom_smooth(method = loess, se = FALSE, colour = "red", method.args = list(deg = 1, span = 1 / 2, family = "symmetric")) +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 4)
  ) +
  xlab("Democratic Margin of Victory (Running Variable)") +
  ylab("Democratic Vote Percentage in Next Election (Outcome Variable)") +
  ggtitle("Relationship between Running Variable and Outcome") +
  facet_wrap(
    facets = ~cp,
    nrow = 1,
    ncol = 2
  )
```

## Continuity in Potential Outcomes Framework {.allowframebreaks}


Let's first use a simple linear model to estimate the relationship on both sides of the cutpoint. Here we are not using a window.

```{r}

lm_predict_below <- predict(
  object = lm(formula = DPctNxt ~ DifDPct, data = rdd_data, subset = DifDPct < 0),
  newdata = data.frame(DifDPct = 0)
)

lm_predict_above <- predict(
  object = lm(formula = DPctNxt ~ DifDPct, data = rdd_data, subset = DifDPct > 0),
  newdata = data.frame(DifDPct = 0)
)

lm_predict_above - lm_predict_below
```

## Continuity in Potential Outcomes Framework {.allowframebreaks}

We could also do the same thing with LOESS regression:

```{r}

loess_predict_below <- predict(
  object = loess(
    formula = DPctNxt ~ DifDPct,
    data = rdd_data,
    subset = DifDPct < 0,
    surface = "direct"
  ),
  newdata = data.frame(DifDPct = 0)
)

loess_predict_above <- predict(
  object = loess(
    formula = DPctNxt ~ DifDPct,
    data = rdd_data,
    subset = DifDPct > 0,
    surface = "direct"
  ),
  newdata = data.frame(DifDPct = 0)
)

loess_predict_above - loess_predict_below
```

We could also use a $p$th order polynomial. For example:

```{r}

lm_poly_predict_below <- predict(
  object = lm(
    formula = DPctNxt ~ I(DifDPct^3) + I(DifDPct^2) + DifDPct,
    data = rdd_data, subset = DifDPct < 0
  ),
  newdata = data.frame(DifDPct = 0)
)

lm_poly_predict_above <- predict(
  object = lm(
    formula = DPctNxt ~ I(DifDPct^3) + I(DifDPct^2) + DifDPct,
    data = rdd_data, subset = DifDPct > 0
  ),
  newdata = data.frame(DifDPct = 0)
)

lm_poly_predict_above - lm_poly_predict_below
```

This next calculates a bandwidth and then does what we just did above (all in one step):

```{r}

rdest1 <- RDestimate(DPctNxt~DifDPct,data=rdd_data, se.type="HC3")
summary(rdest1)

plot(rdest1)
```


```{r}
library(rdrobust)

rdrest1 <- rdrobust(y=rdd_data$DPctNxt,x=rdd_data$DifDPct)
summary(rdrest1)


rdplot(y=rdd_data$DPctNxt,x=rdd_data$DifDPct)

```


\citet[647]{cattaneotitiunikvazquez-bare2017} view the potential outcomes as a random sample from an infinite superpopulation. How does this fit with this course's emphasis on finite population, design-based causal inference?


## Summary

 - A discontinuity can be an instrument --- creating a "natural experiment" (We
   recommend @dunning2012 for more on natural experiments.)
 - Notice that we could have matched observations pre-vs-post within a window,
   to strengthen the argument for the "as randomized" assumption.
 - One can use a model to help the argument that the pre- and post-observations
   are comparable (for example, removing trends).  
 - "Regression discontinuity" tends to imply use of a discontinity plus a
   continuity assumption for estimation --- estimation of a difference in
   curves at the cut-point.
 - So, there are two major approaches to the use of discontinuities for causal
   inference: 
   1. It is a  "natural experiment" that allows a kind of instrumental
      variables estimation plus or minus stratification or trend adjustment;
   2. It is like an instrument in the middle of a continuous and smooth process
      that we aim to approximate via flexible data models.

## References