Regression Discontinuity Design

# Fake program data
set.seed(1234)   # makes the random numbers repeatable (so results are the same each time)
num_students <- 1000

tutoring <- tibble(
  id = 1:num_students,  # create 1000 students, each with a unique id
  
  # entrance exam scores: drawn from a beta distribution (skewed toward higher scores)
  entrance_exam = rbeta(num_students, shape1 = 7, shape2 = 2),
  
  # exit exam base scores: drawn from another beta distribution
  exit_exam = rbeta(num_students, shape1 = 5, shape2 = 3)
) |> 
  mutate(
    # scale entrance exam from 0–1 to 0–100 to make it look like a real exam score
    entrance_exam = entrance_exam * 100,
    
    # assign tutoring: students with entrance exam ≤ 70 get a tutor
    tutoring = entrance_exam <= 70
  ) |> 
  mutate(
    # create final exit exam scores:
    # start with a base (exit_exam * 40),
    # add +10 points if the student had tutoring,
    # and add a benefit from prior knowledge (half of entrance score)
    exit_exam = exit_exam * 40 + 10 * tutoring + entrance_exam / 2
  ) |> 
  mutate(
    # introduce some “fuzziness” (real life isn’t perfect!)
    # for students with entrance exam between 60 and 80,
    # randomly decide whether they got tutoring or not
    tutoring_fuzzy = ifelse(
      entrance_exam > 60 & entrance_exam < 80,
      sample(c(TRUE, FALSE), n(), replace = TRUE),
      tutoring
    )
  ) |> 
  mutate(
    # create nice text labels for charts and tables
    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(
    # center the entrance exam around 70 (the cutoff)
    # so 0 = exactly at the tutoring rule
    entrance_centered = entrance_exam - 70
  )

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)

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)

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

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_cartesian(xlim = c(30, 100),
           ylim = c(40, 90))

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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))

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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))

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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))

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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))

# 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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))+
      guides(fill = guide_legend(reverse = TRUE))

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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))+
      guides(fill = guide_legend(reverse = TRUE))

head(tutoring, 5)

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

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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))

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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))

#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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))

#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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))

#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_cartesian(xlim = c(30, 100), 
           ylim = c(40, 90))

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 = "Uniform\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("Uniform\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)
)

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

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

library(dplyr)
set.seed(123)  # Set random seed for reproducibility
num_students <- 1000  # Number of simulated students

# Create a tibble (data frame) with simulated student data
tutoring <- tibble(
  id = 1:num_students,  # Unique ID for each student

  # Simulate entrance exam scores using a Beta distribution (scaled 0–100)
  entrance_exam = rbeta(num_students, shape1 = 7, shape2 = 2),

  # Simulate base ability for exit exam (before adding treatment effects)
  exit_exam_base = rbeta(num_students, shape1 = 5, shape2 = 3)
) |> 
  mutate(
    entrance_exam = entrance_exam * 100,  # Convert entrance exam scores to 0–100 scale

    # Deterministic assignment rule:
    # Students with entrance_exam ≤ 70 are assigned to receive tutoring
    tutoring_assignment = entrance_exam <= 70,

    # Fuzzy compliance: not everyone follows their assignment perfectly
    tutoring_fuzzy = case_when(
      # 85% of assigned students actually take tutoring
      tutoring_assignment & runif(n()) < 0.85 ~ TRUE,
      # 15% of unassigned students take tutoring anyway
      !tutoring_assignment & runif(n()) < 0.15 ~ TRUE,
      # Everyone else does not take tutoring
      TRUE ~ FALSE
    ),

    # Construct the true exit exam score assuming tutoring has a positive effect
    # - Base score scaled to 0–40
    # - +10 points for being assigned tutoring
    # - +0.5 points per entrance exam point
    exit_exam = exit_exam_base * 40 + 10 * tutoring_assignment + entrance_exam / 2,

    # Same as above but using actual (fuzzy) tutoring status
    exit_exam_fuzzy = exit_exam_base * 40 + 10 * tutoring_fuzzy + entrance_exam / 2,

    # Label variables for clarity (for use in plots/tables)
    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")),

    # Center entrance exam scores around 70 for regression interpretation
    entrance_centered = entrance_exam - 70
  )

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="Entrance exam score", y="Density")

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)

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)