Sharp RDD

library(tidyverse)
library(targets)
library(modelsummary)
library(tinytable)
library(scales)
library(ggdag)
library(dagitty)
library(parameters)
library(rdrobust)
library(rddensity)

tar_config_set(
  store = here::here('_targets'),
  script = here::here('_targets.R')
)

rdd <- tar_read(data_rdd_sharp)

invisible(list2env(tar_read(graphic_functions), .GlobalEnv))
set_annotation_fonts()

gof_map <- tribble(
  ~raw,            ~clean,               ~fmt, ~omit,
  "nobs",          "N",                  0,    FALSE,
  "r.squared",     "$R^2$",          2,    FALSE,
  "adj.r.squared", "\\(R^2\\) adjusted", 2,    FALSE
)

rdd |> 
  ggplot(aes(
    x = entrance_score, y = as.logical(experiential_learning), 
    color = factor(experiential_learning)
  )) +
  geom_point(
    size = 0.5, alpha = 0.5,
    position = position_jitter(width = 0, height = 0.25, seed = 1234)
  ) +
  geom_vline(xintercept = 75) +
  labs(
    x = "Entrance test score", 
    y = "Participated in\nexperiential learning program"
  ) +
  scale_color_manual(values = c(clrs$Prism[2], clrs$Prism[8]), guide = "none") +
  theme_np()

ggplot(rdd, aes(x = entrance_score, fill = as.logical(experiential_learning))) +
  geom_histogram(binwidth = 1, color = "white", boundary = 0) +
  geom_vline(xintercept = 75) +
  labs(x = "Entrance exam score", y = "Count", fill = "In program") +
  scale_fill_manual(values = c(clrs$Prism[2], clrs$Prism[8])) +
  theme_np()

test_density <- rddensity(rdd$entrance_score, c = 75)
# summary(test_density)

plot_density_test <- rdplotdensity(
  rdd = test_density,
  X = rdd$entrance_score,
  type = "both"
)

ggplot(rdd, aes(x = entrance_score, y = prosocial_index, color = as.logical(experiential_learning))) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_smooth(data = filter(rdd, entrance_score <= 75, entrance_score >= 70),
              method = "lm", se = FALSE, linewidth = 2) +
  geom_smooth(data = filter(rdd, entrance_score > 75, entrance_score <= 80),
              method = "lm", se = FALSE, linewidth = 2) +
  geom_vline(xintercept = 75.5) +
  coord_cartesian(xlim = c(65, 85), ylim = c(60, 100)) +
  labs(x = "Entrance test score", y = "Prosocial index", color = "Participated") +
  scale_color_manual(values = c(clrs$Prism[2], clrs$Prism[8])) +
  theme_np(base_size = 16)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

blah <- rdplot(y = rdd$prosocial_index, x = rdd$entrance_score, c = 75, hide = TRUE)
## [1] "Mass points detected in the running variable."
blah$rdplot +
  theme_np() +
  labs(title = NULL, x = "Entrance test score", y = "Prosocial index")

model_rdrobust <- rdrobust(y = rdd$prosocial_index, x = rdd$entrance_score, c = 75)

model_rdrobust |>
  summary()
## Sharp RD estimates using local polynomial regression.
## 
## Number of Obs.                  967
## BW type                       mserd
## Kernel                   Triangular
## VCE method                       NN
## 
## Number of Obs.                  440          527
## Eff. Number of Obs.             254          307
## Order est. (p)                    1            1
## Order bias  (q)                   2            2
## BW est. (h)                   7.904        7.904
## BW bias (b)                  12.811       12.811
## rho (h/b)                     0.617        0.617
## Unique Obs.                      25           26
## 
## =============================================================================
##         Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
## =============================================================================
##   Conventional     6.278     1.728     3.634     0.000     [2.892 , 9.664]     
##         Robust         -         -     2.869     0.004     [1.921 , 10.206]    
## =============================================================================

tidy(model_rdrobust)
## # A tibble: 3 × 7
##   term           estimate std.error statistic  p.value conf.low conf.high
##   <chr>             <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 Conventional       6.28      1.73      3.63 0.000279     2.89      9.66
## 2 Bias-Corrected     6.06      1.73      3.51 0.000449     2.68      9.45
## 3 Robust             6.06      2.11      2.87 0.00412      1.92     10.2

model_bw5 <- lm(
  prosocial_index ~ entrance_score + experiential_learning,
  data = filter(rdd, entrance_score >= (75 - 5) & entrance_score <= (75 + 5))
)

model_bw10 <- lm(
  prosocial_index ~ entrance_score + experiential_learning,
  data = filter(rdd, entrance_score >= (75 - 10) & entrance_score <= (75 + 10))
)

models <- list(
  "BW = 5" = model_bw5,
  "BW = 10" = model_bw10,
  "rdrobust" = model_rdrobust
)

modelsummary(
  models,
  statistic = "[{conf.low}, {conf.high}]",
  fmt = 2,
  coef_map = c(
    `experiential_learning` = "Experiential learning",
    `Conventional` = "Experiential learning"
  ),
  gof_map = gof_map,
  add_columns = tibble("True effect" = 8.5) |> magrittr::set_attr("position", 2)
)

	True effect	BW = 5	BW = 10	rdrobust
Experiential learning	8.50	9.07	8.82	6.28
		[5.80, 12.35]	[6.41, 11.22]	[2.89, 9.66]
N		431	720
$R^2$		0.27	0.35
\(R^2\) adjusted		0.26	0.35

treatment_effects <- enframe(models) |> 
  mutate(params = map(value, \(x) tidy(x, conf.int = TRUE))) |> 
  unnest(params) |> 
  filter(term %in% c("experiential_learning", "Conventional")) |> 
  mutate(name = factor(name, levels = names(models)))

ggplot(treatment_effects, aes(x = estimate, y = fct_rev(name))) +
  geom_vline(xintercept = 8.5, color = clrs$Prism[2]) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high)) +
  labs(x = "Treatment effect", y = NULL) +
  theme_np(base_size = 16) +
  theme(panel.grid.major.y = element_blank())