install.packages(c("rdrobust",
"rddensity",
"lmtest",
"sandwich"))Sharp Regression Discontinuity Design —————————————————–
1 Install and load packages
library(rdrobust)
library(rddensity)
library(lmtest)
library(sandwich)
packageVersion("rdrobust")
#> [1] '4.0.0'
packageVersion("rddensity")
#> [1] '3.0'
packageVersion("lmtest")
#> [1] '0.9.40'
packageVersion("sandwich")
#> [1] '3.1.1'2 Load and check the data
sharp_data <- read.csv("C:/Users/.../sharp_data.csv")head(sharp_data)
#> AGNO AGNO_c D YKS
#> 1 44.64775 -25.35225 0 318.1068
#> 2 47.31709 -22.68291 0 323.1706
#> 3 51.43038 -18.56962 0 325.9956
#> 4 53.78935 -16.21065 0 345.2440
#> 5 57.87032 -12.12968 0 339.0280
#> 6 59.29873 -10.70127 0 342.5197
tail(sharp_data)
#> AGNO AGNO_c D YKS
#> 95 85.29534 15.29534 1 428.4752
#> 96 86.53061 16.53061 1 422.0859
#> 97 86.91692 16.91692 1 434.5603
#> 98 87.67141 17.67141 1 429.4159
#> 99 90.51579 20.51579 1 447.7176
#> 100 91.61646 21.61646 1 427.7910YKS: Öğrencinin Yüksek Öğretim Sınavı Puanı
AGNO: Öğrencinin Ağırlıklı Genel Not Ortalaması
AGNO_c: Öğrencinin eşik değerine (70) göre ortalanmış AGNO’su
D: Öğrencinin vakıftan burs alıp almadığı
3 Check manipulation of the running variable
Code
density_test <- rddensity(X = sharp_data$AGNO, # running variable
c = 70) # cutoff
summary(density_test)
#>
#> Manipulation testing using local polynomial density estimation.
#>
#> Number of obs = 100
#> Model = unrestricted
#> Kernel = triangular
#> BW method = estimated
#> VCE method = jackknife
#>
#> c = 70 Left of c Right of c
#> Number of obs 51 49
#> Eff. Number of obs 39 23
#> Order est. (p) 2 2
#> Order bias (q) 3 3
#> BW est. (h) 7.921 6.287
#>
#> Method T P > |T|
#> Robust 0.3411 0.733
#>
#>
#> P-values of binomial tests (H0: p=0.5).
#>
#> Window Length <c >=c P>|T|
#> 6.087 + 6.087 33 20 0.0984
#> 6.290 + 6.109 34 20 0.0759
#> 6.494 + 6.131 35 20 0.0581
#> 6.698 + 6.153 37 22 0.0674
#> 6.902 + 6.176 37 22 0.0674
#> 7.106 + 6.198 37 23 0.0925
#> 7.309 + 6.220 38 23 0.0722
#> 7.513 + 6.242 38 23 0.0722
#> 7.717 + 6.264 38 23 0.0722
#> 7.921 + 6.287 39 23 0.0559
rdplotdensity(rdd = density_test, # rddenisty() result
X = sharp_data$AGNO,
noPlot = TRUE) # running variable
#> $Estl
#> Call: lpdensity
#>
#> Sample size 51
#> Polynomial order for point estimation (p=) 2
#> Order of derivative estimated (v=) 1
#> Polynomial order for confidence interval (q=) 3
#> Kernel function triangular
#> Scaling factor 0.515151515151515
#> Bandwidth method user provided
#>
#> Use summary(...) to show estimates.
#>
#> $Estr
#> Call: lpdensity
#>
#> Sample size 49
#> Polynomial order for point estimation (p=) 2
#> Order of derivative estimated (v=) 1
#> Polynomial order for confidence interval (q=) 3
#> Kernel function triangular
#> Scaling factor 0.494949494949495
#> Bandwidth method user provided
#>
#> Use summary(...) to show estimates.
#>
#> $Estplot
4 Simple mean difference
4.1 Regression equation
\[ YKS_i = \beta_0 + \beta_1 D_i + \epsilon_i \]
4.2 Estimation
rd_mean <- lm(YKS ~ D, # model
data = sharp_data) # data
coef(summary(rd_mean))
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 358.64499 1.747540 205.22847 6.715047e-131
#> D 62.90158 2.496486 25.19605 1.319258e-444.3 Plot
5 Linear functional form - same regression line on both sides
5.1 Regression equation
\[ YKS_i = \beta_0 + \beta_1 D_i + \beta_2 AGNO_i + \epsilon_i \]
5.2 Estimation
rd_parallel <- lm(YKS ~ D + AGNO,
data = sharp_data)
coef(summary(rd_parallel))
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 254.407587 10.9705711 23.190004 2.365823e-41
#> D 40.134864 2.9842957 13.448689 6.771265e-24
#> AGNO 1.624911 0.1698836 9.564851 1.167380e-155.3 Plot
6 Linear functional form - different regression lines on both sides
6.1 Regression equation
\[ YKS_i = \beta_0 + \beta_1 D_i + \beta_2 AGNO_{c,i} + \beta_3 D_i * AGNO_{c,i} + \epsilon_i \]
6.2 Estimation
rd_lin <- lm(YKS ~ D + AGNO_c + I(D*AGNO_c),
data = sharp_data)
coef(summary(rd_lin))
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 371.524724 1.7707948 209.806762 1.315444e-129
#> D 41.603947 2.8390879 14.653983 3.176609e-26
#> AGNO_c 2.201514 0.2246551 9.799531 4.009082e-16
#> I(D * AGNO_c) -1.169992 0.3200142 -3.656064 4.179514e-046.3 Plot
7 Quadratic functional form - different regression lines on both sides
7.1 Regression equation
\[ YKS_i = \beta_0 + \beta_1 D_i + \beta_2 AGNO_{c,i} + \beta_3 AGNO_{c,i}^2 + \beta_4 D_i * AGNO_{c,i} + \beta_5 D_i * AGNO_{c,i}^2 + \epsilon_i \]
7.2 Estimation
rd_quad <- lm(YKS ~ D + AGNO_c + I(AGNO_c^2) + I(D*AGNO_c) + I(D*AGNO_c^2),
data = sharp_data)
coef(summary(rd_quad))
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 372.81671614 2.58779043 144.0675844 4.515711e-112
#> D 39.86339303 4.17595749 9.5459288 1.697315e-15
#> AGNO_c 2.63252821 0.66488251 3.9593886 1.460167e-04
#> I(AGNO_c^2) 0.01979083 0.02870461 0.6894651 4.922292e-01
#> I(D * AGNO_c) -1.46493163 1.01144536 -1.4483547 1.508454e-01
#> I(D * AGNO_c^2) -0.02677161 0.04705384 -0.5689570 5.707423e-017.3 Plot
8 {rdrobust} bandwidth, estimation and plotting
8.1 Linear functional form
8.1.1 Bandwidth
bw_lin <- rdbwselect(y = sharp_data$YKS, # outcome
x = sharp_data$AGNO, # running variables
c = 70, # cutoff
vce = "hc0", # robust standard error estimator
p = 1, # polynomial form (1 = linear)
kernel = "triangular", # triangular kernel weighting
bwselect = "mserd")$bws[1,1]
# mserd: one common MSE-optimal bandwidth on two sides
# msetwo: two separate MSE-optimal bandwidth on two sides
bw_lin
#> [1] 5.3723798.1.2 Estimation
rdrobust_lin <- rdrobust(y = sharp_data$YKS, # outcome
x = sharp_data$AGNO, # running variable
c = 70, # cutoff
p = 1, # polynomial form (1 = linear)
vce = "hc0", # robust standard error estimator
h = bw_lin) # bandwidth
summary(rdrobust_lin)
#> Call: rdrobust
#>
#> Sharp RD estimates using local polynomial regression.
#>
#> Number of Obs. 100
#> BW type Manual
#> Kernel Triangular
#> VCE method HC0
#>
#> Left Right
#> Number of Obs. 51 49
#> Eff. Number of Obs. 27 17
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 5.372 5.372
#> BW bias (b) 5.372 5.372
#> rho (h/b) 1.000 1.000
#> Unique Obs. 51 49
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 39.431 5.718 0.000 [26.783 , 54.720]
#> =====================================================================
# estimate and standard error
c("tau" = rdrobust_lin$coef[1],
"se" = rdrobust_lin$se[1])
#> tau se
#> 39.431198 5.0555298.1.3 Plot
8.2 Quadratic functional form
8.2.1 Bandwidth
bw_quad <- rdbwselect(y = sharp_data$YKS, # outcome
x = sharp_data$AGNO, # running variables
c = 70, # cutoff
vce = "hc0", # robust standard error estimator
p = 2, # polynomial form (2 = quadratic)
kernel = "triangular", # triangular kernel weighting
bwselect = "mserd")$bws[1,1]
# mserd: one common MSE-optimal bandwidth on two sides
# msetwo: two separate MSE-optimal bandwidth on two sides
bw_quad
#> [1] 7.4027378.2.2 Estimation
rdrobust_quad <- rdrobust(y = sharp_data$YKS, # outcome
x = sharp_data$AGNO, # running variable
c = 70, # cutoff
p = 2, # polynomial form (1 = quadratic)
vce = "hc0", # robust standard error estimator
h = bw_quad) # bandwidth
summary(rdrobust_quad)
#> Call: rdrobust
#>
#> Sharp RD estimates using local polynomial regression.
#>
#> Number of Obs. 100
#> BW type Manual
#> Kernel Triangular
#> VCE method HC0
#>
#> Left Right
#> Number of Obs. 51 49
#> Eff. Number of Obs. 38 27
#> Order est. (p) 2 2
#> Order bias (q) 3 3
#> BW est. (h) 7.403 7.403
#> BW bias (b) 7.403 7.403
#> rho (h/b) 1.000 1.000
#> Unique Obs. 51 49
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 41.205 5.055 0.000 [24.996 , 56.654]
#> =====================================================================
# estimate and standard error
c("tau" = rdrobust_quad$coef[1],
"se" = rdrobust_quad$se[1])
#> tau se
#> 41.204934 6.1034588.2.3 Plot
9 Back to lm() but add bandwidth
9.1 Bandwidth and kernel weights
# function to calculate kernel weights
triangular_kernel_weights <- function(x, cutoff = 0, bandwidth) {
u <- (x - cutoff) / bandwidth # standardized distance u
K <- ifelse(abs(u) <= 1, 1 - abs(u), 0) # triangular kernel evaluation K(u)
weights <- K / sum(K) # normalized weights (sum to 1)
return(weights)
} # triangular_kernel_weightsCalculate the optimal bandwidth.
sharp_data_bw_lin <- subset(x = sharp_data,
subset = abs(AGNO_c) <= bw_lin)
# calculate triangular kernel weights
tkweights_bw_lin <- triangular_kernel_weights(x = sharp_data_bw_lin$AGNO_c,
cutoff = 0,
bandwidth = bw_lin)9.2 Plot kernel weights
9.3 Estimation
rd_lin <- lm(YKS ~ D + AGNO_c + I(D*AGNO_c),
weights = tkweights_bw_lin,
data = sharp_data_bw_lin)
print(coeftest(rd_lin, vcov = vcovHC(rd_lin, type = "HC0")))
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 374.81275 1.87841 199.5376 < 2.2e-16 ***
#> D 39.43120 5.05553 7.7996 1.481e-09 ***
#> AGNO_c 3.66923 0.84634 4.3354 9.557e-05 ***
#> I(D * AGNO_c) -3.90765 1.86415 -2.0962 0.04244 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# compare to {rdronbust}
c("tau" = rdrobust_lin$coef[1],
"se" = rdrobust_lin$se[1])
#> tau se
#> 39.431198 5.0555299.4 Plot
10 Bandwidth sensitivity analysis
10.1 Linear functional form with interaction and half optimal bandwidth
10.1.1 Re-estimate kernel weights
bw_lin_half <- bw_lin / 2
bw_lin_half
#> [1] 2.686189
sharp_data_bw_lin_half <- subset(x = sharp_data,
subset = abs(AGNO_c) <= bw_lin_half)
# calculate triangular kernel weights
tkweights_bw_lin_half <- triangular_kernel_weights(x = sharp_data_bw_lin_half$AGNO_c,
cutoff = 0,
bandwidth = bw_lin_half)10.1.2 Estimation
rd_lin_bw_half <- lm(YKS ~ D + AGNO_c + I(D*AGNO_c),
weights = tkweights_bw_lin_half,
data = sharp_data_bw_lin_half)
print(coeftest(rd_lin_bw_half, vcov = vcovHC(rd_lin_bw_half, type = "HC0")))
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 375.1250 1.9502 192.3488 < 2.2e-16 ***
#> D 41.6082 7.6175 5.4622 2.862e-05 ***
#> AGNO_c 3.6316 1.5014 2.4187 0.02578 *
#> I(D * AGNO_c) -7.2880 6.7937 -1.0728 0.29682
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# compare
rdrobust_lin_bw_half <- rdrobust(y = sharp_data$YKS,
x = sharp_data$AGNO,
c = 70,
p = 1,
vce = "hc0",
h = bw_lin_half)
c("tau" = rdrobust_lin_bw_half$coef[1],
"se" = rdrobust_lin_bw_half$se[1])
#> tau se
#> 41.608209 7.61750410.1.3 Plot
10.2 Linear functional form with interaction and twice optimal bandwidth
10.2.1 Re-estimate kernel weights
bw_lin_twice <- 2 * bw_lin
bw_lin_twice
#> [1] 10.74476
sharp_data_bw_lin_twice <- subset(x = sharp_data,
subset = abs(AGNO_c) <= bw_lin_twice)
# calculate triangular kernel weights
tkweights_bw_lin_twice <- triangular_kernel_weights(x = sharp_data_bw_lin_twice$AGNO_c,
cutoff = 0,
bandwidth = bw_lin_twice)10.2.2 Estimation
rd_lin_bw_twice <- lm(YKS ~ D + AGNO_c + I(D*AGNO_c),
weights = tkweights_bw_lin_twice,
data = sharp_data_bw_lin_twice)
print(coeftest(rd_lin_bw_twice, vcov = vcovHC(rd_lin_bw_twice, type = "HC0")))
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 373.0194 1.6875 221.0465 < 2.2e-16 ***
#> D 38.7804 3.9550 9.8055 4.439e-15 ***
#> AGNO_c 2.6038 0.4849 5.3699 8.498e-07 ***
#> I(D * AGNO_c) -1.2346 0.8324 -1.4831 0.1422
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# compare
rdrobust_lin_bw_twice <- rdrobust(y = sharp_data$YKS,
x = sharp_data$AGNO,
c = 70,
p = 1,
vce = "hc0",
h = bw_lin_twice)
c("tau" = rdrobust_lin_bw_twice$coef[1],
"se" = rdrobust_lin_bw_twice$se[1])
#> tau se
#> 38.780363 3.95496410.2.3 Plot
11 Quadratic functional form with interaction and optimal bandwidth
11.1 Bandwidth and kernel weights
sharp_data_bw_quad <- subset(x = sharp_data,
subset = abs(AGNO_c) <= bw_quad)
# calculate triangular kernel weights
tkweights_bw_quad <- triangular_kernel_weights(x = sharp_data_bw_quad$AGNO_c,
cutoff = 0,
bandwidth = bw_quad)11.2 Estimation
rd_quad <- lm(YKS ~ D + AGNO_c + I(AGNO_c^2) + I(D*AGNO_c) + I(D*AGNO_c^2),
weights = tkweights_bw_quad,
data = sharp_data_bw_quad)
print(coeftest(rd_quad, vcov = vcovHC(rd_quad, type = "HC0")))
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 375.19104 2.10566 178.1823 < 2.2e-16 ***
#> D 41.20493 6.10346 6.7511 7.128e-09 ***
#> AGNO_c 4.55342 1.78487 2.5511 0.01335 *
#> I(AGNO_c^2) 0.28904 0.30626 0.9438 0.34913
#> I(D * AGNO_c) -8.51092 4.29929 -1.9796 0.05242 .
#> I(D * AGNO_c^2) 0.55130 0.61159 0.9014 0.37103
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# compare to {rdronbust}
c("tau" = rdrobust_quad$coef[1],
"se" = rdrobust_quad$se[1])
#> tau se
#> 41.204934 6.10345811.3 Plot
