Fuzzy regression discontinuity

Fuzzy parametric estimation

Let’s make an instrument! We’ll also center the running variable just like we did with sharp regression discontinuity:

tutoring_centered <- tutoring %>% 
  mutate(entrance_centered = entrance_exam - 70,
         below_cutoff = entrance_exam <= 70)
tutoring_centered

## # A tibble: 1,000 x 7
##       id entrance_exam tutoring tutoring_text exit_exam entrance_centered below_cutoff
##    <dbl>         <dbl> <lgl>    <chr>             <dbl>             <dbl> <lgl>       
##  1     1          92.4 FALSE    No tutor           78.1            22.4   FALSE       
##  2     2          72.8 FALSE    No tutor           58.2             2.77  FALSE       
##  3     3          53.7 TRUE     Tutor              62.0           -16.3   TRUE        
##  4     4          98.3 FALSE    No tutor           67.5            28.3   FALSE       
##  5     5          69.7 TRUE     Tutor              54.1            -0.288 TRUE        
##  6     6          68.1 TRUE     Tutor              60.1            -1.93  TRUE        
##  7     7          86.0 FALSE    No tutor           73.1            16.0   FALSE       
##  8     8          85.7 TRUE     Tutor              76.7            15.7   FALSE       
##  9     9          85.9 FALSE    No tutor           57.8            15.9   FALSE       
## 10    10          89.5 FALSE    No tutor           79.9            19.5   FALSE       
## # … with 990 more rows

Now we have a new column named below_cutoff that we’ll use as an instrument. Most of the time this will be the same as the tutoring column, since most people are compliers. But some people didn’t comply, like person 8 here who was not below the cutoff but still used the tutoring program.

Before using the instrument, let’s first run a model that assumes the cutoff is sharp. As we did with the sharp parametric analysis, we’ll include two explanatory variables:

\[ \text{Exit exam} = \beta_0 + \beta_1 \text{Entrance exam score}_\text{centered} + \beta_2 \text{Tutoring program} + \epsilon \]

We’ll use a bandwidth of 10:

# Bandwidth ±10
model_sans_instrument <- lm(exit_exam ~ entrance_centered + tutoring,
                            data = filter(tutoring_centered,
                                          entrance_centered >= -10 & 
                                            entrance_centered <= 10))
tidy(model_sans_instrument)

## # A tibble: 3 x 5
##   term              estimate std.error statistic   p.value
##   <chr>                <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)         59.3      0.503     118.   9.75e-313
## 2 entrance_centered    0.511    0.0665      7.69 1.17e- 13
## 3 tutoringTRUE        11.5      0.744      15.4  1.77e- 42

Here, the coefficient for tutoringTRUE shows the size of the jump, which is 11.5. This means that participating in the tutoring program causes an increase of 11.5 points on the final exam for people in the bandwidth.

BUT THIS IS WRONG. This is not a sharp discontinuity, so we can’t actually do this. Instead, we need to run a 2SLS model that includes our instrument in the first stage, which will then remove the endogeneity built into participation in the program. We’ll estimate this set of models:

\[ \begin{aligned} \widehat{\text{Tutoring program}} &= \gamma_0 + \gamma_1 \text{Entrance exam score}_\text{centered} + \gamma_2 \text{Below cutoff} + \omega \\ \text{Exit exam} &= \beta_0 + \beta_1 \text{Entrance exam score}_\text{centered} + \beta_2 \widehat{\text{Tutoring program}} + \epsilon \end{aligned} \]

We could manually run the first stage model, generate predicted tutoring and then use those predicted values in the second stage model like we did in the instrumental variables example, but that’s tedious and nobody wants to do all that work. We’ll use iv_robust() from the estimatr package instead.

model_fuzzy <- iv_robust(
  exit_exam ~ entrance_centered + tutoring | entrance_centered + below_cutoff,
  data = filter(tutoring_centered, entrance_centered >= -10 & entrance_centered <= 10)
)
tidy(model_fuzzy)

##                term estimate std.error statistic  p.value conf.low conf.high  df   outcome
## 1       (Intercept)    60.14     1.018      59.1 9.7e-200    58.14     62.14 400 exit_exam
## 2 entrance_centered     0.44     0.099       4.4  1.4e-05     0.24      0.63 400 exit_exam
## 3      tutoringTRUE     9.74     1.912       5.1  5.4e-07     5.98     13.50 400 exit_exam

Based on this model, using below_cutoff as an instrument, we can see that the coefficient for tutoringTRUE is different now! It’s 9.74, which means that the tutoring program causes an average increase of 9.74 points on the final exam for compliers in the bandwidth.

Notice that last caveat. Because we’re working with regression discontinuity, we’re estimating a local average treatment effect (LATE) for people in the bandwidth. Because we’re working with instrumental variables, we’re estimating the LATE for compliers only. That means our fuzzy regression discontinuity result here is doubly robust.

If we compare this fuzzy result to the sharp result, we can see a sizable difference:

# gof_omit here will omit goodness-of-fit rows that match any of the text. This
# means 'contains "IC" OR contains "Low" OR contains "Adj" OR contains "p.value"
# OR contains "statistic" OR contains "se_type"'. Basically we're getting rid of
# all the extra diagnostic information at the bottom
modelsummary(list("No instrument (wrong)" = model_sans_instrument,
                  "Fuzzy RD (bw = 10)" = model_fuzzy),
             gof_omit = 'IC|Log|Adj|p\\.value|statistic|se_type',
             stars = TRUE) %>% 
  # Add a background color to row 5
  row_spec(5, background = "#F5ABEA")

We can (and should!) do all the other things that we talked about in the regression discontinuity example, like modifying the bandwidth, adding polynomial terms, and so forth to see how robust the finding is. But we won’t do any of that here.

Fuzzy nonparametric estimation

We can also use nonparametric methods to measure the size of the fuzzy gap at the cutoff. We’ll use rdrobust() just like we did in the sharp example. The only difference is that we have to add one extra argument. That’s it!

To do fuzzy estimation with rdrobust(), use the fuzzy argument to specify the treatment column (or tutoring in our case). Importantly (and confusingly! this took me waaaaay too long to figure out!), you do not need to specify an instrument (or even create one!). All you need to specify is the column that indicates treatment status—rdrobust() will do all the above/below-the-cutoff instrument stuff behind the scenes for you.

rdrobust(y = tutoring$exit_exam, x = tutoring$entrance_exam, 
         c = 70, fuzzy = tutoring$tutoring) %>% 
  summary()

## Call: rdrobust
## 
## Number of Obs.                 1000
## BW type                       mserd
## Kernel                   Triangular
## VCE method                       NN
## 
## Number of Obs.                 238         762
## Eff. Number of Obs.            170         347
## Order est. (p)                   1           1
## Order bias  (q)                  2           2
## BW est. (h)                 12.985      12.985
## BW bias (b)                 19.733      19.733
## rho (h/b)                    0.658       0.658
## Unique Obs.                    238         762
## 
## =============================================================================
##         Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
## =============================================================================
##   Conventional     9.683     1.893     5.116     0.000     [5.973 , 13.393]    
##         Robust         -         -     4.258     0.000     [5.210 , 14.095]    
## =============================================================================

That’s all! Using nonparametric methods, with a triangular kernel and a bandwidth of ±12.96, the causal effect of the tutoring program for compliers in the bandwidth is 9.683.

We can (and should!) do all the other nonparametric robustness checks that we talked about in the regression discontinuity example, like modifying the bandwidth (ideal, half, double) and messing with the kernel (uniform, triangular, Epanechnikov) to see how robust the finding is. But again, we won’t do any of that here.