Regression Discontinuity Design

Causality at the Cutoff

Bogdan G. Popescu
bogdan.popescu@johncabot.edu

John Cabot University

Intro

How can we identify causal effects using observational data?

Some of the ways to conduct causal analysis with observational data, we need to run

  • differences in differences
  • regression discontinuity design

RDD

RDD are about arbitrary rules and thresholds determining access to policy programs

  • subjects are in the program if they have scores above a specific threshold
  • subject are not in the program if they have scores below a specific threshold

Or the other way around

Running or forcing variable - index that measures eligibility for the program (e.g. scores)
Cutoff / cutpoint / threshold - the number that the determines access to the program

RDD

RDD

Student ID Exam Score Scholarship Awarded? First-Year GPA
101 83 No 3.2
102 84 No 3.3
103 85 Yes 3.6
104 86 Yes 3.5
105 87 Yes 3.7
106 84.5 No 3.4
107 85.1 Yes 3.6

RDD

Hypothetical tutoring program

Students take an entrance exam

Those who score 70 or lower get a free tutor for the year

Students then take an exit exam at the end of the year

RDD

Hypothetical tutoring program

Show the code 
# Fake program data
set.seed(1234)
num_students <- 1000
tutoring <- tibble(
  id = 1:num_students,
  entrance_exam = rbeta(num_students, shape1 = 7, shape2 = 2),
  exit_exam = rbeta(num_students, shape1 = 5, shape2 = 3)
) |> 
  mutate(entrance_exam = entrance_exam * 100,
         tutoring = entrance_exam <= 70) |> 
  mutate(exit_exam = exit_exam * 40 + 10 * tutoring + entrance_exam / 2) |> 
  mutate(tutoring_fuzzy = ifelse(entrance_exam > 60 & entrance_exam < 80,
                                 sample(c(TRUE, FALSE), n(), replace = TRUE),
                                 tutoring)) |> 
  mutate(tutoring_text = factor(tutoring, levels = c(FALSE, TRUE),
                                labels = c("No tutor", "Tutor")),
         tutoring_fuzzy_text = factor(tutoring_fuzzy, levels = c(FALSE, TRUE),
                                      labels = c("No tutor", "Tutor"))) |> 
  mutate(entrance_centered = entrance_exam - 70)
Show the code 
ggplot(tutoring, aes(x = entrance_exam, y = tutoring_text, fill = tutoring_text)) +
  geom_vline(xintercept = 70, linewidth = 2, color = "#FFDC00") + 
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4,
             position = position_jitter(width = 0, height = 0.15, seed = 1234)) + 
  labs(x = "Entrance exam score", y = NULL) + 
  guides(fill = "none") +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  theme_bw(base_size = 28)

Intuition

Hypothetical tutoring program

The people right before and right after the threshold are very similar

The idea is that we can measure the effect of the program by examining students just around the threshold

Show the code 
ggplot(tutoring, aes(x = entrance_exam, y = tutoring_text, fill = tutoring_text)) +
  geom_vline(xintercept = 70, linewidth = 2, color = "#FFDC00") + 
    annotate(geom = "rect", fill = "grey50", alpha = 0.25, ymin = -Inf, ymax = Inf,
           xmin = 70 - 5,  xmax = 70 + 5) +
    annotate(geom = "rect", fill = "grey50", alpha = 0.5, ymin = -Inf, ymax = Inf,
           xmin = 70 - 2,  xmax = 70 + 2) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4,
             position = position_jitter(width = 0, height = 0.15, seed = 1234)) + 
  labs(x = "Entrance exam score", y = NULL) + 
  guides(fill = "none") +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  theme_bw(base_size = 28)

Intuition

Hypothetical tutoring program

The people right before and right after the threshold are very similar

The idea is that we can measure the effect of the program by examining students just around the threshold.

Show the code 
ggplot(tutoring, aes(x = entrance_exam, y = tutoring_text, fill = tutoring_text)) +
  geom_vline(xintercept = 70, linewidth = 2, color = "#FFDC00") + 
    annotate(geom = "rect", fill = "grey50", alpha = 0.25, ymin = -Inf, ymax = Inf,
           xmin = 70 - 5,  xmax = 70 + 5) +
    annotate(geom = "rect", fill = "grey50", alpha = 0.5, ymin = -Inf, ymax = Inf,
           xmin = 70 - 2,  xmax = 70 + 2) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4,
             position = position_jitter(width = 0, height = 0.15, seed = 1234)) + 
  labs(x = "Entrance exam score", y = NULL) + 
  guides(fill = "none") +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  theme_bw(base_size = 28)+
    coord_cartesian(xlim = c(70 - 5, 70 + 5))

Intuition

Hypothetical tutoring program

Show the code 
ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4)+
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  #guides(fill = "none") +
    guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  theme_bw(base_size = 28)+
    coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))

Intuition

Hypothetical tutoring program

Show the code 
ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4)+
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  #guides(fill = "none") +
    guides(fill = guide_legend(reverse = TRUE)) +
  geom_smooth(data = filter(tutoring, entrance_exam <= 70),
              method = "lm",
              color="red",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  geom_smooth(data = filter(tutoring, entrance_exam > 70),
              method = "lm", 
              color="black",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
    scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  theme_bw(base_size = 28)+
    coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))

Intuition

Hypothetical tutoring program

Show the code 
library("broom")
tutoring <- tutoring |>
  mutate(tutoring = relevel(factor(tutoring), ref = "TRUE"))  # "TRUE" = treated is reference

rdd_tutoring <- lm(exit_exam ~ entrance_centered + tutoring, data = tutoring) |> 
  tidy()

effect_control <- filter(rdd_tutoring, term == "(Intercept)")$estimate
late <- filter(rdd_tutoring, term == "tutoringFALSE")$estimate
effect_treatment <- effect_control + late


ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  #guides(fill = "none") +
    guides(fill = guide_legend(reverse = TRUE)) +
  geom_smooth(data = filter(tutoring, entrance_exam <= 70),
              method = "lm",
              color="red",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  geom_smooth(data = filter(tutoring, entrance_exam > 70),
              method = "lm", 
              color="black",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  theme_bw(base_size = 28)+
  annotate(geom = "segment", x = 70, xend = 70, 
           y = effect_control, yend = effect_treatment,
           size = 4, color = "red") +
  annotate(geom = "label", x = 75, y = effect_treatment - (late / 2),
           label = "δ\n(causal\n effect)",
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
    coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))

Intuition

Hypothetical tutoring program

The \(\delta\) measures the gap in outcome for people on both sides of the cutpoint

\(\delta\) is also called LATE (local average treatment effects)

The size of the gap depends on how you draw the lines on each side of the cutoff

The type of lines you choose can change the estimate of \(\delta\)

  • parametric vs. non-parametric
  • bandwidths
  • kernels

Intuition

Hypothetical tutoring program

Show the code 
ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  #guides(fill = "none") +
    guides(fill = guide_legend(reverse = TRUE)) +
  geom_smooth(data = filter(tutoring, entrance_exam <= 70),
              method = "lm",
              color="red",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  geom_smooth(data = filter(tutoring, entrance_exam > 70),
              method = "lm", 
              color="black",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  theme_bw(base_size = 28)+
  annotate(geom = "segment", x = 70, xend = 70, 
           y = effect_control, yend = effect_treatment,
           size = 4, color = "red") +
  annotate(geom = "label", x = 75, y = effect_treatment - (late / 2),
           label = "δ\n(causal\n effect)",
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
    coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))

Intuition

Hypothetical tutoring program

Show the code 
program_data <- tutoring %>%
  mutate(entrance_centered = 
           entrance_exam - 70)
model1 <- lm(exit_exam ~ 
               entrance_centered + tutoring,
             data = program_data)

x<-tidy(model1)
late<-round(x$estimate[x$term=="tutoringFALSE"],2)
const<-round(x$estimate[x$term=="(Intercept)"],2)


ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  #guides(fill = "none") +
    guides(fill = guide_legend(reverse = TRUE)) +
  geom_smooth(data = filter(tutoring, entrance_exam <= 70),
              method = "lm",
              color="red",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  geom_smooth(data = filter(tutoring, entrance_exam > 70),
              method = "lm", 
              color="black",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  theme_bw(base_size = 28)+
  annotate(geom = "segment", x = 70, xend = 70, 
           y = effect_control, yend = effect_treatment,
           size = 4, color = "red") +
  annotate(geom = "label", x = 75, y = const + (late / 2),
           label = paste("δ=", late),
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
    coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))

Intuition

Hypothetical tutoring program

Show the code 
# 1. Fit cubic models
fit_control <- lm(exit_exam ~ poly(entrance_exam, 3, raw = TRUE),
                  data = filter(tutoring, entrance_exam < 70))

fit_treat <- lm(exit_exam ~ poly(entrance_exam, 3, raw = TRUE),
                data = filter(tutoring, entrance_exam >= 70))

# 2. Get tidy versions
tidy_control <- tidy(fit_control)
tidy_treat <- tidy(fit_treat)

cutoff <- 70
effect_control <- predict(fit_control, newdata = data.frame(entrance_exam = cutoff))
effect_treatment <- predict(fit_treat, newdata = data.frame(entrance_exam = cutoff))
delta <- round(effect_treatment-effect_control, 2)


# Create prediction data for plotting
x_vals <- tibble(entrance_exam = seq(30, 100, length.out = 300))

preds <- x_vals %>%
  mutate(
    exit_exam = ifelse(
      entrance_exam < 70,
      predict(fit_control, newdata = x_vals),
      predict(fit_treat, newdata = x_vals)
    )
  )

ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4) +
  geom_line(data = preds, aes(x = entrance_exam, y = exit_exam),
            inherit.aes = FALSE, color = "black", size = 1.5) +
  geom_vline(xintercept = 70, color = "#FFDC00", size = 2) +
  theme_bw(base_size = 28)+
scale_fill_manual(values = c("Tutor" = "red", "No tutor" = "black"), name = NULL,
                  guide = guide_legend(reverse = TRUE))+
  labs(x = "Entrance exam score", y = "Exit exam score")+
    annotate(geom = "label", x = 75, y = const + (late / 2),
           label = paste("δ=", delta),
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
      coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))+
      guides(fill = guide_legend(reverse = TRUE))

Intuition

Hypothetical tutoring program

Show the code 
library(dplyr)

# Add sine-transformed variable
tutoring <- tutoring |>
  mutate(sin_term = sin(0.2 * entrance_exam))

# Control group: entrance_exam < 70
fit_trig_control <- nls(
  exit_exam ~ a + b * entrance_exam + c * sin(d * entrance_exam),
  data = filter(tutoring, entrance_exam < 70),
  start = list(a = 55, b = 0.4, c = 15, d = 0.4)
)

# Treatment group: entrance_exam >= 70
fit_trig_treat <- nls(
  exit_exam ~ a + b * entrance_exam + c * sin(d * entrance_exam),
  data = filter(tutoring, entrance_exam >= 70),
  start = list(a = 55, b = 0.4, c = 15, d = 0.4)
)

x_vals <- tibble(entrance_exam = seq(30, 100, length.out = 300))

preds_trig <- x_vals |>
  mutate(
    exit_exam = ifelse(
      entrance_exam < 70,
      predict(fit_trig_control, newdata = x_vals),
      predict(fit_trig_treat, newdata = x_vals)
    )
  )

cutoff <- 70

# Create a new data frame for x = 70
cutoff_df <- data.frame(entrance_exam = cutoff)

# Predicted outcome from each model
effect_control <- predict(fit_trig_control, newdata = cutoff_df)
effect_treatment <- predict(fit_trig_treat, newdata = cutoff_df)

# Compute the treatment effect (jump at cutoff)
delta <- round(effect_treatment-effect_control,2)

ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4) +
  geom_line(data = preds_trig, aes(x = entrance_exam, y = exit_exam),
            inherit.aes = FALSE, color = "black", size = 1.5) +
  geom_vline(xintercept = 70, color = "#FFDC00", size = 2) +
  theme_bw(base_size = 28)+
scale_fill_manual(values = c("Tutor" = "red", "No tutor" = "black"), name = NULL,
                  guide = guide_legend(reverse = TRUE))+
  labs(x = "Entrance exam score", y = "Exit exam score")+
    annotate(geom = "label", x = 75, y = const + (late / 2),
           label = paste("δ=", delta),
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
      coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))+
      guides(fill = guide_legend(reverse = TRUE))

Lines and Fit

Parametric line should fit the data pretty well

Non-parametric lines use the data to find the best line, often with windows and moving averages.

For example, one common method is LOESS (Locally Estimated Scatterplot Smoothing)

Measuring the Gap Parametrically

The gap can be easily measured by centering the running variable around the threshold

Show the code 
head(tutoring, 5)

\[ y = \beta_0 + \beta_1 \text{Running Variable (centered)} + \beta_2 \text{Treatment} \]

Measuring the Gap Parametrically

And this is how we can display the results

Show the code 
library(modelsummary)
program_data <- tutoring |> 
  mutate(entrance_centered = 
           entrance_exam - 70)

#Running model
model1 <- lm(exit_exam ~  entrance_centered + tutoring, data = program_data)

#Make table
modelsummary(model1, stars = TRUE, gof_map = c("nobs", "r.squared"))
(1)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 70.306***
(0.501)
entrance_centered 0.514***
(0.027)
tutoringFALSE -10.968***
(0.802)
Num.Obs. 1000
R2 0.273

Measuring the Gap Parametrically

And this is how we can display the results on the graph

Show the code 
x<-tidy(model1)
gap<-round(x$estimate[x$term=="tutoringFALSE"],2)
const<-round(x$estimate[x$term=="(Intercept)"],2)

ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  #guides(fill = "none") +
    guides(fill = guide_legend(reverse = TRUE)) +
  geom_smooth(data = filter(tutoring, entrance_exam <= 70),
              method = "lm",
              color="red",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  geom_smooth(data = filter(tutoring, entrance_exam > 70),
              method = "lm", 
              color="black",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  theme_bw(base_size = 28)+
  annotate(geom = "segment", x = 70, xend = 70, 
           y = effect_control, yend = effect_treatment,
           size = 4, color = "red") +
  annotate(geom = "label", x = 75, y = const + (late / 2),
           label = paste("δ=", gap),
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
    coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))

Measuring the Gap Nonparametrically

We need to calculate the gap nonparametrically using rdrobust

library(rdrobust)
rd_out <- rdrobust(y = tutoring$exit_exam,
                   x = tutoring$entrance_exam,
                   c = 70)

#summary(rd_out)
#Extracting conventional coefficient
gap_rdrobust<-rd_out$coef["Conventional", "Coeff"]

Measuring the Gap Nonparametrically

And this is how we can display the results on the graph

Show the code 
model_control <- lm(exit_exam ~ entrance_exam,
                    data = filter(tutoring, entrance_exam <= 70))

model_treat <- lm(exit_exam ~ entrance_exam,
                  data = filter(tutoring, entrance_exam > 70))

effect_control <- predict(model_control, newdata = data.frame(entrance_exam = 70))
effect_treatment <- predict(model_treat, newdata = data.frame(entrance_exam = 70))

late <- round(gap_rdrobust,2)

ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  #guides(fill = "none") +
    guides(fill = guide_legend(reverse = TRUE)) +
  geom_smooth(data = filter(tutoring, entrance_exam <= 70),
              method = "loess",
              color="red",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  geom_smooth(data = filter(tutoring, entrance_exam > 70),
              method = "loess", 
              color="black",
              #aes(color = tutoring_text),
              size = 2, show.legend = FALSE) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  theme_bw(base_size = 28)+
  annotate(geom = "segment", x = 70, xend = 70, 
           y = effect_control, yend = effect_treatment,
           size = 4, color = "red") +
    annotate(geom = "label", x = 75, y = const + (late / 2),
           label = paste("δ=", late),
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
    coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))

Bandwidths

What we care about is the area right around the cutoff.

The observations far away do not matter because they are too diffent and thus, they are not comparable

That is why we need to focus on bandwidth or the window around cutoff

Bandwidths

These are the observations within a 5-point bandwidth

Show the code 
#Marking obs. within a 5-point bandwidth
tutoring <- tutoring %>%
  mutate(
    in_bw_5 = abs(entrance_exam - 70) <= 5,
    in_bw_10 = abs(entrance_exam - 70) <= 10
  )


#Running model
program_data_bw_5<-subset(tutoring, in_bw_5==TRUE)
model2 <- lm(exit_exam ~  entrance_centered + tutoring, data = program_data_bw_5)

x<-tidy(model2)
gap<-round(x$estimate[x$term=="tutoringFALSE"],2)
const<-round(x$estimate[x$term=="(Intercept)"],2)


ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.2) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  #guides(fill = "none") +
    guides(fill = guide_legend(reverse = TRUE)) +
  theme_bw(base_size = 28)+
    geom_point(data = filter(tutoring, in_bw_5),
           aes(x = entrance_exam, y = exit_exam),
           size = 3, pch = 21, color = "white", alpha = 1)+
geom_smooth(data = filter(tutoring, in_bw_5 & entrance_exam <= 70),
            method = "lm",
            color = "red",
            size = 2,
            show.legend = FALSE) +
geom_smooth(data = filter(tutoring, in_bw_5 & entrance_exam > 70),
            method = "lm",
            color = "black",
            size = 2,
            show.legend = FALSE)+
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  annotate(geom = "segment", x = 70, xend = 70, 
           y = const, yend = const + gap,
           size = 4, color = "red") +
    annotate(geom = "label", x = 75, y = const + (gap / 2),
           label = paste("δ=", gap),
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
    coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))

Bandwidths

These are the observations within a 10-point bandwidth

Show the code 
#Running model
program_data_bw_10<-subset(tutoring, in_bw_10==TRUE)
model3 <- lm(exit_exam ~  entrance_centered + tutoring, data = program_data_bw_10)

x<-tidy(model3)
gap<-round(x$estimate[x$term=="tutoringFALSE"],2)
const<-round(x$estimate[x$term=="(Intercept)"],2)


ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.2) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  #guides(fill = "none") +
    guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  theme_bw(base_size = 28)+
    geom_point(data = filter(tutoring, in_bw_10),
           aes(x = entrance_exam, y = exit_exam),
           size = 3, pch = 21, color = "white", alpha = 1)+
  geom_smooth(data = filter(tutoring, in_bw_10 & entrance_exam <= 70),
            method = "lm",
            color = "red",
            size = 2,
            show.legend = FALSE) +
geom_smooth(data = filter(tutoring, in_bw_10 & entrance_exam > 70),
            method = "lm",
            color = "black",
            size = 2,
            show.legend = FALSE)+
  annotate(geom = "segment", x = 70, xend = 70, 
           y = const, yend = const + gap,
           size = 4, color = "red") +
    annotate(geom = "label", x = 75, y = const + (gap / 2),
           label = paste("δ=", gap),
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
    coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))

Bandwidths

Optimal Bandwidth

But which ones do we choose?

The rdrobust package allows us to decide on the bandwidth:

Show the code 
rd_out <- rdrobust(y = tutoring$exit_exam,
                   x = tutoring$entrance_exam,
                   c = 70)

bw_opt <- rd_out$bws[1]  # optimal bandwidth for conventional estimate
bw_opt
[1] 7.615756

Bandwidths

Optimal Bandwidth

These are the results using the 7.62 bandwidth.

Show the code 
#Running model
program_data_opt<-subset(tutoring, abs(entrance_centered)<bw_opt)
model4 <- lm(exit_exam ~  entrance_centered + tutoring, data = program_data_opt)

x<-tidy(model4)
gap<-round(x$estimate[x$term=="tutoringFALSE"],2)
const<-round(x$estimate[x$term=="(Intercept)"],2)


ggplot(tutoring, aes(x = entrance_exam, y = exit_exam, fill = tutoring_text)) +
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.2) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  #guides(fill = "none") +
    guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  theme_bw(base_size = 28)+
    geom_point(data = program_data_opt,
           aes(x = entrance_exam, y = exit_exam),
           size = 3, pch = 21, color = "white", alpha = 1)+
    geom_smooth(data = filter(program_data_opt, entrance_exam <= 70),
              method = "lm",
              color = "red",
              size = 2,
              show.legend = FALSE) +
  geom_smooth(data = filter(program_data_opt, entrance_exam > 70),
              method = "lm",
              color = "black",
              size = 2,
              show.legend = FALSE) +
  annotate(geom = "segment", x = 70, xend = 70, 
           y = const, yend = const + gap,
           size = 4, color = "red") +
    annotate(geom = "label", x = 75, y = const + (gap / 2),
           label = paste("δ=", gap),
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
    coord_sf(xlim = c(30, 100), 
           ylim = c(40, 90))

Bandwidths

It’s a good idea to have multiple bandwidths

  • some that make sense
  • some that are decided by the optimal bandwidth algorithm

Bandwidths

Add either table or graph

Kernel

The other important part of RD is choosing the kernel

Kernel = method for assigning importance to observations based on distance to the cutoff

Bandwidths

Kernel: Uniform

Show the code 
library(ggnewscale)
# Define the boxcar (uniform) kernel function
rec <- function(x) (abs(x) <= 1) * 0.5

# Set cutoff and bandwidth
cutoff <- 70
bandwidth <- 40

# For rectangular kernel overlay
kernel_df <- data.frame(
  entrance_exam = seq(cutoff - bandwidth - 10, cutoff + bandwidth + 10, length.out = 200)
)
kernel_df$u <- (kernel_df$entrance_exam - cutoff) / bandwidth
kernel_df$boxcar_weight <- rec(kernel_df$u)
kernel_df$boxcar_weight_rescaled <- kernel_df$boxcar_weight * 50 + 40  # visual rescaling

# Apply kernel weights to main data
tutoring$boxcar_weight <- with(tutoring, {
  u <- (entrance_exam - cutoff) / bandwidth
  rec(u)
})

# Run weighted regression
model4 <- lm(exit_exam ~ entrance_centered + tutoring, 
             data = tutoring, 
             weights = boxcar_weight)

x <- tidy(model4)
gap <- round(x$estimate[x$term == "tutoringFALSE"], 2)
const <- round(x$estimate[x$term == "(Intercept)"], 2)

# Plot with rectangular kernel overlay
ggplot() + 
  geom_point(data = tutoring, 
             aes(x = entrance_exam, 
                 y = exit_exam, 
                 fill = tutoring_text,
                 size = boxcar_weight),
             shape = 21, 
             alpha = 0.2) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
  geom_smooth(
    data = filter(tutoring, entrance_exam <= 70),
    aes(x = entrance_exam, y = exit_exam),
    method = "lm",
    color = alpha("red", 0.2),
    fill = alpha("red", 0.2),
    size = 2,
    show.legend = FALSE
  ) +
  geom_smooth(
    data = filter(tutoring, entrance_exam > 70),
    aes(x = entrance_exam, y = exit_exam),
    method = "lm",
    color = alpha("black", 0.2),
    fill = alpha("black", 0.2),
    alpha = 0.2,
    size = 2,
    show.legend = FALSE
  ) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  scale_size_continuous(name = "Uniform\nWeight", range = c(2, 6)) +
  theme_bw(base_size = 28) +
  coord_cartesian(xlim = c(20, 120), ylim = c(40, 90)) +
  new_scale_color() +
  
  # Rectangle-shaped kernel overlay using geom_segment()
  geom_segment(
    aes(x = cutoff - bandwidth, xend = cutoff + bandwidth,
        y = 0.5 * 50 + 40, yend = 0.5 * 50 + 40,
        color = "Boxcar\nKernel"),
    size = 2
  ) +
  geom_segment(aes(x = cutoff - bandwidth, xend = cutoff - bandwidth,
                 y = 40, yend = 0.5 * 50 + 40),
             color = "blue", size = 1) +
  geom_segment(aes(x = cutoff + bandwidth, xend = cutoff + bandwidth,
                 y = 40, yend = 0.5 * 50 + 40),
             color = "blue", size = 1)+
  
  scale_color_manual(name = NULL, values = c("Boxcar\nKernel" = "blue"))+
    annotate(geom = "segment", x = 70, xend = 70, 
           y = const, yend = const + gap,
           size = 4, color = "red") +
  annotate(geom = "label", x = 75, y = const + (gap / 2),
           label = paste("δ=", gap),
           size = 5,
           color = "red", fill = alpha("white", 0.8))+
    guides(
  size = guide_legend(order = 1),
  fill = guide_legend(order = 2),
  color = guide_legend(order = 3)
)

Bandwidths

Kernel: Triangular

Show the code 
# Set cutoff and bandwidth
cutoff <- 70
bandwidth <- 40  # adjust as needed

# Create a data frame to draw the triangular kernel
kernel_df <- data.frame(
  entrance_exam = seq(cutoff - bandwidth, cutoff + bandwidth, length.out = 200)
)

# Triangular kernel weights for visualization
kernel_df$tri_weight <- with(kernel_df, {
  u <- (entrance_exam - cutoff) / bandwidth
  ifelse(abs(u) <= 1, 1 - abs(u), 0)
})

# Rescale for visual overlay
kernel_df$tri_weight_rescaled <- kernel_df$tri_weight * 50 + 40  # adjust scaling as needed

# Compute triangular weights in main dataset
tutoring$tri_weight <- with(tutoring, {
  u <- (entrance_exam - cutoff) / bandwidth
  ifelse(abs(u) <= 1, 1 - abs(u), 0)
})

# Run weighted regression using triangular kernel
model4 <- lm(exit_exam ~ entrance_centered + tutoring, 
             data = tutoring, 
             weights = tri_weight)

# Extract estimates for annotation
x <- tidy(model4)
gap <- round(x$estimate[x$term == "tutoringFALSE"], 2)
const <- round(x$estimate[x$term == "(Intercept)"], 2)


ggplot()+ 
  geom_point(data = tutoring, 
             aes(x = entrance_exam, 
                 y = exit_exam, 
                 fill = tutoring_text,
                 size = tri_weight),
             shape = 21, 
             alpha = 0.2) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
geom_smooth(
  data = filter(tutoring, entrance_exam <= 70),
  aes(x = entrance_exam, y = exit_exam),
  method = "lm",
  color = alpha("red", 0.2),        # line color
  fill = alpha("red", 0.2),          # confidence interval fill
  size = 2,
  show.legend = FALSE
) +
geom_smooth(
  data = filter(tutoring, entrance_exam > 70),
  aes(x = entrance_exam, y = exit_exam),
  method = "lm",
  color = alpha("black", 0.2),        # line color
  fill = alpha("black", 0.2),          # confidence interval fill
  alpha = 0.2,
  size = 2,
  show.legend = FALSE
) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  scale_size_continuous(name = "Triangular\nWeight", range = c(2, 6)) +  # legend control
  theme_bw(base_size = 28)+
  annotate(geom = "segment", x = 70, xend = 70, 
           y = const, yend = const + gap,
           size = 4, color = "red") +
    annotate(geom = "label", x = 75, y = const + (gap / 2),
           label = paste("δ=", gap),
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
  coord_cartesian(xlim = c(20, 120), ylim = c(40, 90)) +
  new_scale_color() +
  geom_line(
    data = kernel_df,
    aes(x = entrance_exam, 
        y = tri_weight_rescaled,
        color = "Triangular\n Kernel"), # <-- also mapped to show dashed line in legend
    inherit.aes = FALSE,
    size = 2
  )+
  scale_color_manual(name = NULL, values = c("Triangular\n Kernel" = "blue"))+
    guides(
  size = guide_legend(order = 1),
  fill = guide_legend(order = 2),
  color = guide_legend(order = 3)
)

Bandwidths

Kernel: Epanechnikov

Show the code 
kernel_df <- data.frame(
  entrance_exam = seq(cutoff - bandwidth, cutoff + bandwidth, length.out = 200)
)

u <- (kernel_df$entrance_exam - cutoff) / bandwidth
kernel_df$epa_weight <- ifelse(abs(u) <= 1, 0.75 * (1 - u^2), 0)
kernel_df$epa_weight_rescaled <- kernel_df$epa_weight * 50 + 40  # adjust as needed


# Compute Epanechnikov weights
tutoring$epa_weight <- with(tutoring, {
  u <- (entrance_exam - cutoff) / bandwidth
  ifelse(abs(u) <= 1, 0.75 * (1 - u^2), 0)
})

model4 <- lm(exit_exam ~ entrance_centered + tutoring, 
             data = tutoring, 
             weights = epa_weight)

x<-tidy(model4)
gap<-round(x$estimate[x$term=="tutoringFALSE"],2)
const<-round(x$estimate[x$term=="(Intercept)"],2)


ggplot()+ 
  geom_point(data = tutoring, 
             aes(x = entrance_exam, 
                 y = exit_exam, 
                 fill = tutoring_text,
                 size = epa_weight),
             shape = 21, 
             alpha = 0.2) +
  geom_vline(xintercept = 70, size = 2, color = "#FFDC00") + 
  labs(x = "Entrance exam score", y = "Exit exam score") + 
geom_smooth(
  data = filter(tutoring, entrance_exam <= 70),
  aes(x = entrance_exam, y = exit_exam),
  method = "lm",
  color = alpha("red", 0.2),        # line color
  fill = alpha("red", 0.2),          # confidence interval fill
  size = 2,
  show.legend = FALSE
) +
geom_smooth(
  data = filter(tutoring, entrance_exam > 70),
  aes(x = entrance_exam, y = exit_exam),
  method = "lm",
  color = alpha("black", 0.2),        # line color
  fill = alpha("black", 0.2),          # confidence interval fill
  alpha = 0.2,
  size = 2,
  show.legend = FALSE
) +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  scale_color_manual(values = c("black", "red")) +
  scale_size_continuous(name = "Epanechnikov\nWeight", range = c(2, 6)) +  # legend control
  theme_bw(base_size = 28)+
  annotate(geom = "segment", x = 70, xend = 70, 
           y = const, yend = const + gap,
           size = 4, color = "red") +
    annotate(geom = "label", x = 75, y = const + (gap / 2),
           label = paste("δ=", gap),
           size=5,
           color = "red", fill = alpha(c("white"),0.8))+
  coord_cartesian(xlim = c(20, 120), ylim = c(40, 90)) +
  new_scale_color() +
  geom_line(
    data = kernel_df,
    aes(x = entrance_exam, 
        y = epa_weight_rescaled,
        color = "Epanechnikov\n Kernel"), # <-- also mapped to show dashed line in legend
    inherit.aes = FALSE,
    size = 2
  )+
  scale_color_manual(name = NULL, values = c("Epanechnikov\n Kernel" = "blue"))+
  guides(
  size = guide_legend(order = 1),
  fill = guide_legend(order = 2),
  color = guide_legend(order = 3)
)

Practical Advice

In the end there are many different moving parameter parameters.

The best thing to do is to use rdrobust and present all the the results

New

Show the code 
# Uniform
# Run rdrobust with all results
rd1 <- rdrobust(y = tutoring$exit_exam, x = tutoring$entrance_exam, c = 70,
               kernel = "Uniform")
# Create a named list of models (each method with its estimate and se)
models <- list(
  "Conventional" = rd1$Estimate[1],
  "Bias-Corrected" = rd1$Estimate[2],
  "Robust Bias-Corrected" = rd1$Estimate[3]
)

ses <- list(
  "Conventional" = rd1$se[1],
  "Bias-Corrected" = rd1$se[2],
  "Robust Bias-Corrected" = rd1$se[3]
)

# Build a modelsummary-style data frame
library(tibble)


results1 <- tibble::tibble(
  Statistic = c("Estimate", "Std. Error", "Kernel", "Obs.", "Bandwidth"),
  `Conventional` = c(
    sprintf("%.3f", models$Conventional),
    sprintf("%.3f", ses$Conventional),
    rd1$kernel,  # keep as-is (a string)
    sum(rd1$N),
    sprintf("%.3f",rd1$bws[1])
  ),
  `Bias-Corrected` = c(
    sprintf("%.3f", models$`Bias-Corrected`),
    sprintf("%.3f", ses$`Bias-Corrected`),
    rd1$kernel,
    sum(rd1$N),
    sprintf("%.3f",rd1$bws[2])
  ),
  `Robust Bias-Corrected` = c(
    sprintf("%.3f", models$`Robust Bias-Corrected`),
    sprintf("%.3f", ses$`Robust Bias-Corrected`),
    rd1$kernel,
    sum(rd1$N),
    sprintf("%.3f",rd1$bws[2])
  )
)

#Triangular

# Run rdrobust with all results
rd2 <- rdrobust(y = tutoring$exit_exam, x = tutoring$entrance_exam, c = 70,
               kernel = "Triangular")
# Create a named list of models (each method with its estimate and se)
models <- list(
  "Conventional" = rd2$Estimate[1],
  "Bias-Corrected" = rd2$Estimate[2],
  "Robust Bias-Corrected" = rd2$Estimate[3]
)

ses <- list(
  "Conventional" = rd2$se[1],
  "Bias-Corrected" = rd2$se[2],
  "Robust Bias-Corrected" = rd2$se[3]
)

# Build a modelsummary-style data frame
library(tibble)


results2 <- tibble::tibble(
  Statistic = c("Estimate", "Std. Error", "Kernel", "Obs.", "Bandwidth"),
  `Conventional` = c(
    sprintf("%.3f", models$Conventional),
    sprintf("%.3f", ses$Conventional),
    rd2$kernel,  # keep as-is (a string)
    sum(rd2$N),
    sprintf("%.3f",rd2$bws[1])
  ),
  `Bias-Corrected` = c(
    sprintf("%.3f", models$`Bias-Corrected`),
    sprintf("%.3f", ses$`Bias-Corrected`),
    rd2$kernel,
    sum(rd2$N),
    sprintf("%.3f",rd2$bws[2])
  ),
  `Robust Bias-Corrected` = c(
    sprintf("%.3f", models$`Robust Bias-Corrected`),
    sprintf("%.3f", ses$`Robust Bias-Corrected`),
    rd2$kernel,
    sum(rd2$N),
    sprintf("%.3f",rd2$bws[2])
  )
)

#Epanechnikov

# Run rdrobust with all results
rd3 <- rdrobust(y = tutoring$exit_exam, x = tutoring$entrance_exam, c = 70,
               kernel = "Epanechnikov")
# Create a named list of models (each method with its estimate and se)
models <- list(
  "Conventional" = rd3$Estimate[1],
  "Bias-Corrected" = rd3$Estimate[2],
  "Robust Bias-Corrected" = rd3$Estimate[3]
)

ses <- list(
  "Conventional" = rd3$se[1],
  "Bias-Corrected" = rd3$se[2],
  "Robust Bias-Corrected" = rd3$se[3]
)

# Build a modelsummary-style data frame
library(tibble)


results3 <- tibble::tibble(
  Statistic = c("Estimate", "Std. Error", "Kernel", "Obs.", "Bandwidth"),
  `Conventional` = c(
    sprintf("%.3f", models$Conventional),
    sprintf("%.3f", ses$Conventional),
    rd3$kernel,  # keep as-is (a string)
    sum(rd1$N),
    sprintf("%.3f",rd3$bws[1])
  ),
  `Bias-Corrected` = c(
    sprintf("%.3f", models$`Bias-Corrected`),
    sprintf("%.3f", ses$`Bias-Corrected`),
    rd3$kernel,
    sum(rd1$N),
    sprintf("%.3f",rd3$bws[2])
  ),
  `Robust Bias-Corrected` = c(
    sprintf("%.3f", models$`Robust Bias-Corrected`),
    sprintf("%.3f", ses$`Robust Bias-Corrected`),
    rd3$kernel,
    sum(rd3$N),
    sprintf("%.3f",rd3$bws[2])
  )
)

# Create a separator row
separator_uniform <- tibble::tibble(
  Statistic = "<b>Uniform Kernel</b>",
  Conventional = "", `Bias-Corrected` = "", `Robust Bias-Corrected` = ""
)

separator_triangular <- tibble::tibble(
  Statistic = "<b>Triangular Kernel</b>",
  Conventional = "", `Bias-Corrected` = "", `Robust Bias-Corrected` = ""
)

separator_epanechnikov <- tibble::tibble(
  Statistic = "<b>Epanechnikov Kernel</b>",
  Conventional = "", `Bias-Corrected` = "", `Robust Bias-Corrected` = ""
)

results <- rbind(
  separator_uniform,
  results1,
  separator_triangular,
  results2,
  separator_epanechnikov,
  results3
)
# Show using modelsummary custom table
modelsummary::datasummary_df(results, output = "markdown")
Statistic Conventional Bias-Corrected Robust Bias-Corrected
Uniform Kernel
Estimate -9.901 -10.296 1.758
Std. Error 1.758 1.758 2.028
Kernel Uniform Uniform Uniform
Obs. 1000 1000 1000
Bandwidth 5.648 10.607 10.607
Triangular Kernel
Estimate -9.992 -10.228 1.708
Std. Error 1.708 1.708 2.049
Kernel Triangular Triangular Triangular
Obs. 1000 1000 1000
Bandwidth 7.616 11.669 11.669
Epanechnikov Kernel
Estimate -9.897 -10.097 1.686
Std. Error 1.686 1.686 2.018
Kernel Epanechnikov Epanechnikov Epanechnikov
Obs. 1000 1000 1000
Bandwidth 7.066 11.389 11.389

Limitations and Caveats

The gap at the bandwidth is the LATE (Local Average Treatment Effect).

Thus, the findings have high internal validity and low external validity:

  • the findings are limited in scope
  • these are local empirical results

Placebo Test

Falsification via Fake Cutoffs

A placebo test checks for discontinuities at points where no treatment should occur.

It helps ensure that:

  • The observed treatment effect at the true cutoff is not an artifact of underlying trends or model misspecification.
  • There are no jumps in the outcome at arbitrary thresholds.

Placebo Test

Example Strategy

  • Re-run the RDD at other values of the running variable (e.g. 60, 65, 75, 80).
  • If \(\delta\) is non-zero at those points, the original result may be spurious.
Show the code 
library(purrr)
# Fake placebo cutoffs
placebo_cutoffs <- c(60, 65, 70, 75, 80)

# Store estimates
placebo_results <- map_dfr(placebo_cutoffs, function(cut) {
  fit <- rdrobust(y = tutoring$exit_exam, x = tutoring$entrance_exam, c = cut)
  
  tibble(
    cutoff = cut,
    estimate = fit$coef[1, 1],  # [row 1 = Conventional, column 1 = estimate]
    se = fit$se[1, 1]           # [row 1 = Conventional, column 1 = std error]
  )
})

ggplot(placebo_results, aes(x = cutoff, y = estimate)) +
  geom_point(size = 4, color = "red") +
  geom_errorbar(aes(ymin = estimate - 1.96 * se,
                    ymax = estimate + 1.96 * se),
                width = 1.5, color = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
  labs(title = "Placebo Tests at Fake Cutoffs",
       x = "Fake Cutoff",
       y = "Estimated Treatment Effect (δ)") +
  theme_bw(base_size = 18)

Covariate Balance

Are Subjects Comparable Around the Cutoff?

One key assumption of RDD is that units just above and just below the cutoff are similar in all aspects except treatment assignment.

This means pre-treatment covariates should be balanced across the threshold.

Why It Matters

  • Helps validate the continuity assumption
  • Suggests that selection bias is unlikely
  • Strengthens causal interpretation of \(\delta\)

Covariate Balance

How to Test It

Run the same RDD logic, but use pre-treatment covariates as outcomes.

If there’s no jump, the assumption holds.

Caveat 1

Manipulation at the threshold

People might be aware of the cutoff so they might modify their behavior to get in/out of program

For example, students could try to get a slightly lower grade to benefit from the tutoring program.

This means that people at the cutoff are not representative of the control and treatment group.

Caveat 1

Manipulation Test

The way to test for manipulation at the cutoff is the McCrary density test.

This is an example of Non-manipulation

Show the code 
library(dplyr)

set.seed(123)
num_students <- 1000
tutoring <- tibble(
  id = 1:num_students,
  entrance_exam = rbeta(num_students, shape1 = 7, shape2 = 2),
  exit_exam_base = rbeta(num_students, shape1 = 5, shape2 = 3)
) |> 
  mutate(
    entrance_exam = entrance_exam * 100,

    # Assignment to tutoring
    tutoring_assignment = entrance_exam <= 70,

    # Fuzzy compliance behavior
    tutoring_fuzzy = case_when(
      tutoring_assignment & runif(n()) < 0.85 ~ TRUE,
      !tutoring_assignment & runif(n()) < 0.15 ~ TRUE,
      TRUE ~ FALSE
    ),
    # Make tutoring have a **positive effect**
    exit_exam = exit_exam_base * 40 + 10 * tutoring_assignment + entrance_exam / 2,
    exit_exam_fuzzy = exit_exam_base * 40 + 10 * tutoring_fuzzy + entrance_exam / 2,
    # Labels
    tutoring_text = factor(tutoring_assignment, levels = c(FALSE, TRUE), labels = c("No tutor", "Tutor")),
    tutoring_fuzzy_text = factor(tutoring_fuzzy, levels = c(FALSE, TRUE), labels = c("No tutor", "Tutor")),
    entrance_centered = entrance_exam - 70
  )
[1] 1

Caveat 1

Manipulation Test

The way to test for manipulation at the cutoff is the McCrary density test.

This is an example of Manipulation

Show the code 
set.seed(123)  # for reproducibility
num_students <- 1500  # a bit more to have room for thinning

tutoring2 <- tibble(
  id = 1:num_students,
  entrance_exam = rbeta(num_students, shape1 = 7, shape2 = 2),
  exit_exam = rbeta(num_students, shape1 = 5, shape2 = 3)
) |> 
  mutate(entrance_exam = entrance_exam * 100) |>

  # 🔥 Add manipulation: fewer people just above the cutoff
  filter(!(entrance_exam > 70 & entrance_exam < 75 & runif(n()) < 0.5)) |>  # randomly drop 50% just above cutoff
  
  mutate(
    # Assignment
    tutoring_assignment = entrance_exam <= 70,
    
    # Fuzzy compliance
    tutoring_fuzzy = case_when(
      tutoring_assignment & runif(n()) < 0.85 ~ TRUE,
      !tutoring_assignment & runif(n()) < 0.15 ~ TRUE,
      TRUE ~ FALSE
    ),
    
    # Outcome
    exit_exam = exit_exam * 40 + 10 * tutoring_fuzzy + entrance_exam / 2,
    
    # Labels
    tutoring_text = factor(tutoring_assignment, levels = c(FALSE, TRUE), labels = c("No tutor", "Tutor")),
    tutoring_fuzzy_text = factor(tutoring_fuzzy, levels = c(FALSE, TRUE), labels = c("No tutor", "Tutor")),
    entrance_centered = entrance_exam - 70
  )


library(rddensity)
mtest <- rddensity(X = tutoring2$entrance_exam, c=70)
man<-rdplotdensity(mtest,X = tutoring2$entrance_exam)$Estplot+
  geom_vline(xintercept=70, linetype="dashed") +
  labs(x="CD4 Count", y="Density")

Noncompliance

The logic for noncompliance means that people on the margin of the cutoff might end up in/out of the program.

Show the code 
ggplot(tutoring, aes(x = entrance_exam, y = tutoring_fuzzy_text, fill = entrance_exam <= 70)) +
  geom_vline(xintercept = 70, linewidth = 2, color = "#FFDC00") + 
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4,
             position = position_jitter(width = 0, height = 0.15, seed = 1234)) + 
  labs(x = "Entrance exam score", y = NULL) + 
  guides(fill = "none") +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  theme_bw(base_size = 28)

Compliance

The logic for noncompliance means that people on the margin of the cutoff might end up in/out of the program.

Show the code 
ggplot(tutoring, aes(x = entrance_exam, y = tutoring_text, fill = tutoring_text)) +
  geom_vline(xintercept = 70, linewidth = 2, color = "#FFDC00") + 
  geom_point(size = 3, pch = 21, color = "white", alpha = 0.4,
             position = position_jitter(width = 0, height = 0.15, seed = 1234)) + 
  labs(x = "Entrance exam score", y = NULL) + 
  guides(fill = "none") +
  scale_fill_manual(values = c("black", "red"), name = NULL) +
  theme_bw(base_size = 28)

Noncompliance

Addressing it

The way to address noncompliance is to use fuzzy discontinuities:

  • Use an instrument for which side of the cutoff people should be on

Thus, the effect will only be for compliers near the cutoff

Compliance

Person Type If Assigned to Treatment If Assigned to Control
Complier Takes treatment Does not take treatment
Always-Taker Takes treatment Takes treatment
Never-Taker Does not take treatment Does not take treatment
Defier Does not take treatment Takes treatment

Example

Noncompliance

Addressing it

So, we can use an instrument = above/below cutoff to address non-compliance

The instrument has to meet three important criteria

  • Relevance

    • The cutoff changes the likelihood of receiving the treatment.
      → People below the cutoff are more likely to get tutoring.

  • Exclusion

    • The cutoff affects the outcome only through the treatment.
      → Crossing the threshold doesn’t directly change the outcome except by changing tutoring status.

  • Exogeneity

    • Unobserved confounders between the tutoring program and exit exam scores are unrelated to the cutoff.

Noncompliance

Addressing it

Noncompliance

Addressing it

The way to deal with noncompliance is to use a Fuzzy RD

We can do it in the following way:

rdrobust(y = tutoring$exit_exam, x = tutoring$entrance_exam, 
         c = 70, fuzzy = tutoring$tutoring_fuzzy)
Fuzzy RD estimates using local polynomial regression.

Number of Obs.                 1000
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                  246          754
Eff. Number of Obs.             187          350
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                  13.056       13.056
BW bias (b)                  19.608       19.608
rho (h/b)                     0.666        0.666
Unique Obs.                     246          754

First-stage estimates.

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional    -0.680     0.072    -9.382     0.000    [-0.822 , -0.538]    
        Robust         -         -    -8.029     0.000    [-0.869 , -0.528]    
=============================================================================

Treatment effect estimates.

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional    11.955     2.030     5.889     0.000     [7.976 , 15.934]    
        Robust         -         -     4.940     0.000     [7.450 , 17.250]    
=============================================================================

Noncompliance

Addressing it

Show the code 
# Uniform
# Run rdrobust with all results
rd1 <- rdrobust(y = tutoring$exit_exam, x = tutoring$entrance_exam, c = 70,
                fuzzy = tutoring$tutoring_fuzzy,
               kernel = "Uniform")
# Create a named list of models (each method with its estimate and se)
models <- list(
  "Conventional" = rd1$Estimate[1],
  "Bias-Corrected" = rd1$Estimate[2],
  "Robust Bias-Corrected" = rd1$Estimate[3]
)

ses <- list(
  "Conventional" = rd1$se[1],
  "Bias-Corrected" = rd1$se[2],
  "Robust Bias-Corrected" = rd1$se[3]
)

# Build a modelsummary-style data frame
library(tibble)


results1 <- tibble::tibble(
  Statistic = c("Estimate", "Std. Error", "Kernel", "Obs.", "Bandwidth"),
  `Conventional` = c(
    sprintf("%.3f", models$Conventional),
    sprintf("%.3f", ses$Conventional),
    rd1$kernel,  # keep as-is (a string)
    sum(rd1$N),
    sprintf("%.3f",rd1$bws[1])
  ),
  `Bias-Corrected` = c(
    sprintf("%.3f", models$`Bias-Corrected`),
    sprintf("%.3f", ses$`Bias-Corrected`),
    rd1$kernel,
    sum(rd1$N),
    sprintf("%.3f",rd1$bws[2])
  ),
  `Robust Bias-Corrected` = c(
    sprintf("%.3f", models$`Robust Bias-Corrected`),
    sprintf("%.3f", ses$`Robust Bias-Corrected`),
    rd1$kernel,
    sum(rd1$N),
    sprintf("%.3f",rd1$bws[2])
  )
)

#Triangular

# Run rdrobust with all results
rd2 <- rdrobust(y = tutoring$exit_exam, x = tutoring$entrance_exam, c = 70,
                                fuzzy = tutoring$tutoring_fuzzy,
               kernel = "Triangular")
# Create a named list of models (each method with its estimate and se)
models <- list(
  "Conventional" = rd2$Estimate[1],
  "Bias-Corrected" = rd2$Estimate[2],
  "Robust Bias-Corrected" = rd2$Estimate[3]
)

ses <- list(
  "Conventional" = rd2$se[1],
  "Bias-Corrected" = rd2$se[2],
  "Robust Bias-Corrected" = rd2$se[3]
)

# Build a modelsummary-style data frame
library(tibble)


results2 <- tibble::tibble(
  Statistic = c("Estimate", "Std. Error", "Kernel", "Obs.", "Bandwidth"),
  `Conventional` = c(
    sprintf("%.3f", models$Conventional),
    sprintf("%.3f", ses$Conventional),
    rd2$kernel,  # keep as-is (a string)
    sum(rd2$N),
    sprintf("%.3f",rd2$bws[1])
  ),
  `Bias-Corrected` = c(
    sprintf("%.3f", models$`Bias-Corrected`),
    sprintf("%.3f", ses$`Bias-Corrected`),
    rd2$kernel,
    sum(rd2$N),
    sprintf("%.3f",rd2$bws[2])
  ),
  `Robust Bias-Corrected` = c(
    sprintf("%.3f", models$`Robust Bias-Corrected`),
    sprintf("%.3f", ses$`Robust Bias-Corrected`),
    rd2$kernel,
    sum(rd2$N),
    sprintf("%.3f",rd2$bws[2])
  )
)

#Epanechnikov

# Run rdrobust with all results
rd3 <- rdrobust(y = tutoring$exit_exam, x = tutoring$entrance_exam, c = 70,
                                fuzzy = tutoring$tutoring_fuzzy,
               kernel = "Epanechnikov")
# Create a named list of models (each method with its estimate and se)
models <- list(
  "Conventional" = rd3$Estimate[1],
  "Bias-Corrected" = rd3$Estimate[2],
  "Robust Bias-Corrected" = rd3$Estimate[3]
)

ses <- list(
  "Conventional" = rd3$se[1],
  "Bias-Corrected" = rd3$se[2],
  "Robust Bias-Corrected" = rd3$se[3]
)

# Build a modelsummary-style data frame
library(tibble)


results3 <- tibble::tibble(
  Statistic = c("Estimate", "Std. Error", "Kernel", "Obs.", "Bandwidth"),
  `Conventional` = c(
    sprintf("%.3f", models$Conventional),
    sprintf("%.3f", ses$Conventional),
    rd3$kernel,  # keep as-is (a string)
    sum(rd1$N),
    sprintf("%.3f",rd3$bws[1])
  ),
  `Bias-Corrected` = c(
    sprintf("%.3f", models$`Bias-Corrected`),
    sprintf("%.3f", ses$`Bias-Corrected`),
    rd3$kernel,
    sum(rd1$N),
    sprintf("%.3f",rd3$bws[2])
  ),
  `Robust Bias-Corrected` = c(
    sprintf("%.3f", models$`Robust Bias-Corrected`),
    sprintf("%.3f", ses$`Robust Bias-Corrected`),
    rd3$kernel,
    sum(rd3$N),
    sprintf("%.3f",rd3$bws[2])
  )
)

# Create a separator row

separator_fs <- tibble::tibble(
  Statistic = "<b>First stage</b>",
  Conventional = "", `Bias-Corrected` = "", `Robust Bias-Corrected` = ""
)

separator_ss <- tibble::tibble(
  Statistic = "<b>Treatment effect estimates</b>",
  Conventional = "", `Bias-Corrected` = "", `Robust Bias-Corrected` = ""
)

separator_uniform <- tibble::tibble(
  Statistic = "<b>Uniform Kernel</b>",
  Conventional = "", `Bias-Corrected` = "", `Robust Bias-Corrected` = ""
)

separator_triangular <- tibble::tibble(
  Statistic = "<b>Triangular Kernel</b>",
  Conventional = "", `Bias-Corrected` = "", `Robust Bias-Corrected` = ""
)

separator_epanechnikov <- tibble::tibble(
  Statistic = "<b>Epanechnikov Kernel</b>",
  Conventional = "", `Bias-Corrected` = "", `Robust Bias-Corrected` = ""
)



# First-stage estimates
first_stage_rd1 <- tibble::tibble(
  Statistic = c("Estimate", "Std. Error"),
  `Conventional` = c(
    sprintf("%.3f", rd1$tau_T[1]),
    sprintf("%.3f", rd1$se_T[1])),
  
  `Bias-Corrected` = c(
    sprintf("%.3f", rd1$tau_T[2]),
    sprintf("%.3f", rd1$se_T[2])),
  
  `Robust Bias-Corrected` = c(
    sprintf("%.3f", rd1$tau_T[3]),
    sprintf("%.3f", rd1$se_T[3]))
)


# First-stage estimates
first_stage_rd2 <- tibble::tibble(
  Statistic = c("Estimate", "Std. Error"),
  `Conventional` = c(
    sprintf("%.3f", rd2$tau_T[1]),
    sprintf("%.3f", rd2$se_T[1])),
  
  `Bias-Corrected` = c(
    sprintf("%.3f", rd2$tau_T[2]),
    sprintf("%.3f", rd2$se_T[2])),
  
  `Robust Bias-Corrected` = c(
    sprintf("%.3f", rd2$tau_T[3]),
    sprintf("%.3f", rd2$se_T[3]))
)

# First-stage estimates
first_stage_rd3 <- tibble::tibble(
  Statistic = c("Estimate", "Std. Error"),
  `Conventional` = c(
    sprintf("%.3f", rd3$tau_T[1]),
    sprintf("%.3f", rd3$se_T[1])),
  
  `Bias-Corrected` = c(
    sprintf("%.3f", rd3$tau_T[2]),
    sprintf("%.3f", rd3$se_T[2])),
  
  `Robust Bias-Corrected` = c(
    sprintf("%.3f", rd3$tau_T[3]),
    sprintf("%.3f", rd3$se_T[3]))
)

results_final1 <- rbind(
  separator_uniform,
  separator_fs,
  first_stage_rd1,
  separator_ss,
  results1)

results_final2 <- rbind(separator_triangular,
  separator_fs,
  first_stage_rd2,
  separator_ss,
  results2)
  
results_final3 <- rbind(separator_epanechnikov,
  separator_fs,
  first_stage_rd3,
  separator_ss,
  results3)

# Show using modelsummary custom table
modelsummary::datasummary_df(results_final1, output = "markdown")
Statistic Conventional Bias-Corrected Robust Bias-Corrected
Uniform Kernel
First stage
Estimate -0.633 -0.625 -0.625
Std. Error 0.092 0.092 0.110
Treatment effect estimates
Estimate 18.112 18.512 4.144
Std. Error 4.144 4.144 4.984
Kernel Uniform Uniform Uniform
Obs. 1000 1000 1000
Bandwidth 6.123 10.533 10.533

Noncompliance

Addressing it

Show the code 
# Show using modelsummary custom table
modelsummary::datasummary_df(results_final2, output = "markdown")
Statistic Conventional Bias-Corrected Robust Bias-Corrected
Triangular Kernel
First stage
Estimate -0.680 -0.699 -0.699
Std. Error 0.072 0.072 0.087
Treatment effect estimates
Estimate 16.644 16.644 2.899
Std. Error 2.899 2.899 3.546
Kernel Triangular Triangular Triangular
Obs. 1000 1000 1000
Bandwidth 13.090 19.608 19.608

Noncompliance

Addressing it

Show the code 
# Show using modelsummary custom table
modelsummary::datasummary_df(results_final3, output = "markdown")
Statistic Conventional Bias-Corrected Robust Bias-Corrected
Epanechnikov Kernel
First stage
Estimate -0.660 -0.672 -0.672
Std. Error 0.077 0.077 0.092
Treatment effect estimates
Estimate 17.159 17.440 3.214
Std. Error 3.214 3.214 3.902
Kernel Epanechnikov Epanechnikov Epanechnikov
Obs. 1000 1000 1000
Bandwidth 10.779 16.937 16.937

Summary: Regression Discontinuity Design

Goal:
Estimate local causal effect of treatment at a known cutoff using a continuous running variable.

Key Assumptions:
- Continuity: Potential outcomes \(Y(0)\) and \(Y(1)\) change smoothly at the cutoff.
- No Precise Manipulation: Units cannot perfectly sort just above or below the cutoff.

Diagnostics & Robustness:
- McCrary Test: Check for manipulation in the running variable’s density.
- Covariate Balance: No jumps in pre-treatment covariates at cutoff.
- Placebo Tests: No effects at fake cutoffs.

Interpretation:
- Estimate is Local Average Treatment Effect (LATE) at the cutoff.
- RDD is valid only near the cutoff—not population-wide.