install.packages(c("rdrobust",
"rddensity",
"ivreg",
"lmtest",
"sandwich"))Fuzzy Regression Discontinuity Design —————————————————–
1 Install and load packages
library(rdrobust)
library(rddensity)
library(ivreg)
library(lmtest)
library(sandwich)
packageVersion("rdrobust")
#> [1] '4.0.0'
packageVersion("rddensity")
#> [1] '3.0'
packageVersion("ivreg")
#> [1] '0.6.8'
packageVersion("lmtest")
#> [1] '0.9.40'
packageVersion("sandwich")
#> [1] '3.1.1'2 Check data
head(fuzzy_data)
#> ogrenci_no okul_no TKS TKS_c Z p D SM AGNO
#> 1 1 1 28 -2 0 0.165 0 28 2.325117
#> 2 2 1 28 -2 0 0.165 0 28 2.264310
#> 3 3 1 28 -2 0 0.165 0 28 2.364069
#> 4 4 1 28 -2 0 0.165 0 28 2.411930
#> 5 5 1 28 -2 0 0.165 0 28 2.297855
#> 6 6 1 28 -2 0 0.165 0 28 2.329940
tail(fuzzy_data)
#> ogrenci_no okul_no TKS TKS_c Z p D SM AGNO
#> 3027 3027 100 31 1 1 0.755 1 16 3.034225
#> 3028 3028 100 31 1 1 0.755 1 16 3.081669
#> 3029 3029 100 31 1 1 0.755 1 16 2.873357
#> 3030 3030 100 31 1 1 0.755 1 16 3.071537
#> 3031 3031 100 31 1 1 0.755 1 16 3.135065
#> 3032 3032 100 31 1 1 0.755 1 16 2.9816773 Aggregate data for alternative analyses
# aggregate variables by okul_no
fuzzy_school_data <- aggregate(
cbind(AGNO, D, TKS, TKS_c, Z) ~ okul_no,
data = fuzzy_data,
FUN = mean)
# derive the weight variable (row count per okul_no)
weights <- aggregate(
AGNO ~ okul_no,
data = fuzzy_data,
FUN = length)$AGNO
fuzzy_school_data$weights <- weights4 Check manipulation of the running variable
Code
density_test <- rddensity(X = fuzzy_school_data$TKS,
c = 30)
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 = 30 Left of c Right of c
#> Number of obs 42 58
#> Eff. Number of obs 42 58
#> Order est. (p) 2 2
#> Order bias (q) 3 3
#> BW est. (h) 15 15
#>
#> Method T P > |T|
#> Robust 1.3767 0.1686
#> Warning in summary.CJMrddensity(density_test): There are repeated observations.
#> Point estimates and standard errors have been adjusted. Use option
#> massPoints=FALSE to suppress this feature.
#>
#> P-values of binomial tests (H0: p=0.5).
#>
#> Window Length / 2 <c >=c P>|T|
#> 9.000 21 37 0.0479
#> 9.667 21 37 0.0479
#> 10.333 26 40 0.1089
#> 11.000 29 44 0.1006
#> 11.667 29 44 0.1006
#> 12.333 35 50 0.1284
#> 13.000 37 54 0.0929
#> 13.667 37 54 0.0929
#> 14.333 41 57 0.1293
#> 15.000 42 58 0.1332
rdplotdensity(rdd = density_test,
X = fuzzy_school_data$TKS,
noPlot = TRUE)
#> $Estl
#> Call: lpdensity
#>
#> Sample size 42
#> 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.424242424242424
#> Bandwidth method user provided
#>
#> Use summary(...) to show estimates.
#>
#> $Estr
#> Call: lpdensity
#>
#> Sample size 58
#> 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.585858585858586
#> Bandwidth method user provided
#>
#> Use summary(...) to show estimates.
#>
#> $Estplot
5 Bandwidth and data filtering
# bandwidth
bw <- rdbwselect(y = fuzzy_data$AGNO, # outcome
x = fuzzy_data$TKS, # running variables
c = 30, # cutoff
vce = "cr2", # cluster-robust estimator
p = 1, # polynomial form (1 = linear)
kernel = "uniform", # uniform to match results, also could be triangular
bwselect = "mserd",
cluster = fuzzy_data$okul_no)$bws
#> Warning in rdbwselect(y = fuzzy_data$AGNO, x = fuzzy_data$TKS, c = 30, vce =
#> "cr2", : Mass points detected in the running variable.
h_optimal <- bw[1,1]
# data filtering
fuzzy_data_subset <- subset(fuzzy_data,
abs(TKS_c) <= h_optimal)
fuzzy_school_data_subset <- subset(fuzzy_school_data,
abs(TKS_c) <= h_optimal)6 Manual 2SLS with lm()
6.1 Regression equation
6.1.1 Crude way - no standard errors
First stage:
\[ D_j = \alpha_0 + {\color{red}{\alpha_1}} Z_j + \alpha_2 TKS_{c,j} + \alpha_3 Z_j*TKS_{c,j} + v_j \]
Second stage:
\[ AGNO_{ij} = \beta_0 + {\color{red}{\beta_1}} Z_j + \beta_2 TKS_{c,j} + \beta_3 Z_j*TKS_{c,j} + \epsilon_{ij} \] Causal effect:
\[ \tau = {\color{red}{\beta_1}} / {\color{red}{\alpha_1}} \]
6.1.2 Proper way with standard errors
First stage:
\[ D_j = \alpha_0 + \alpha_1 Z_j + \alpha_2 TKS_{c,j} + \alpha_3 Z_j*TKS_{c,j} + v_j \]
Second stage:
\[ AGNO_{ij} = \beta_0 + {\color{red}{\beta_1}} \hat D_j + \beta_2 TKS_{c,j} + \beta_3 Z_j*TKS_{c,j} + \epsilon_{ij} \] Causal effect:
\[ \tau = {\color{red}{\beta_1}} \]
6.2 Estimation
6.2.1 Using student-level data
stage_1 <- lm(D ~ Z + TKS_c + Z:TKS_c,
data = fuzzy_data_subset)
fuzzy_data_subset$D_hat <- predict(stage_1)
stage_2 <- lm(AGNO ~ D_hat + TKS_c + Z:TKS_c,
data = fuzzy_data_subset)
# heteroskedasticity- and cluster-robust standard errors
coeftest(stage_2,
vcov = vcovCL(stage_2,
type = "HC0",
cluster = ~okul_no))
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.334733 0.248922 9.3794 < 2e-16 ***
#> D_hat 0.546485 0.326475 1.6739 0.09452 .
#> TKS_c -0.044874 0.071774 -0.6252 0.53200
#> TKS_c:Z 0.072218 0.070717 1.0212 0.30744
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.2.2 Using school-level data
stage_1 <- lm(D ~ Z + TKS_c + Z:TKS_c,
data = fuzzy_school_data_subset,
weights = weights)
fuzzy_school_data_subset$D_hat <- predict(stage_1)
stage_2 <- lm(AGNO ~ D_hat + TKS_c + Z:TKS_c,
data = fuzzy_school_data_subset,
weights = weights)
# heteroskedasticity-robust standard errors
coeftest(stage_2,
vcov = vcovCL(stage_2,
type = "HC0"))
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.334733 0.248922 9.3794 1.692e-09 ***
#> D_hat 0.546485 0.326475 1.6739 0.1071
#> TKS_c -0.044874 0.071774 -0.6252 0.5377
#> TKS_c:Z 0.072218 0.070717 1.0212 0.3173
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 17 2SLS with {ivreg}
7.1 Estimation
7.1.1 Using student-level data
fit_ivreg <- ivreg(AGNO ~ D + TKS_c + Z:TKS_c | Z + TKS_c + Z:TKS_c,
data = fuzzy_data_subset)
# heteroskedasticity and cluster-robust standard errors
coeftest(fit_ivreg,
vcov = vcovCL(
fit_ivreg,
cluster = ~okul_no,
type = "HC0"))
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.334733 0.147635 15.8142 < 2.2e-16 ***
#> D 0.546485 0.184816 2.9569 0.003193 **
#> TKS_c -0.044874 0.044321 -1.0125 0.311606
#> TKS_c:Z 0.072218 0.042770 1.6885 0.091681 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 17.1.2 Using school-level data
fit_ivreg <- ivreg(AGNO ~ D + TKS_c + Z:TKS_c | Z + TKS_c + Z:TKS_c,
data = fuzzy_school_data_subset,
weights = weights)
# heteroskedasticity-robust standard errors
coeftest(fit_ivreg,
vcov = vcovCL(
fit_ivreg,
type = "HC0"))
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.334733 0.147635 15.8142 3.415e-14 ***
#> D 0.546485 0.184816 2.9569 0.006873 **
#> TKS_c -0.044874 0.044321 -1.0125 0.321418
#> TKS_c:Z 0.072218 0.042770 1.6885 0.104268
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 18 2SLS with {rdrobust}
8.1 Estimation
8.1.1 Using student-level data
fit_rdrobust <- rdrobust(y = fuzzy_data$AGNO, # outcome
x = fuzzy_data$TKS_c, # running variables
c = 0, # cutoff
fuzzy = fuzzy_data$D, # assignment status
p = 1, # 1: linear functional form
h = h_optimal, # bandwidth
kernel = "uniform", # uniform kernel weights
cluster = fuzzy_data$okul_no,
vce = "cr2") # standard error adjustment
#> Warning in rdrobust(y = fuzzy_data$AGNO, x = fuzzy_data$TKS_c, c = 0, fuzzy =
#> fuzzy_data$D, : Mass points detected in the running variable.
summary(fit_rdrobust)
#> Call: rdrobust
#>
#> Fuzzy RD estimates using local polynomial regression. Std. errors are clustered (100 clusters).
#>
#> Number of Obs. 3032
#> BW type Manual
#> Kernel Uniform
#> VCE method CR2
#>
#> Left Right
#> Number of Obs. 895 2137
#> Eff. Number of Obs. 221 635
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 4.506 4.506
#> BW bias (b) 4.506 4.506
#> rho (h/b) 1.000 1.000
#> Clusters (g) 42 58
#> Unique Obs. 15 15
#>
#> First-stage estimates.
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> =====================================================================
#> Rd Effect 0.433 0.152 0.879 [-1.390 , 1.624]
#> =====================================================================
#>
#> Treatment effect estimates.
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 0.546 1.599 0.110 [-0.177 , 1.739]
#> =====================================================================
c("tau" = fit_rdrobust$coef[1],
"se" = fit_rdrobust$se[1])
#> tau se
#> 0.5464854 0.21175478.1.2 Using school-level data
fit_rdrobust <- rdrobust(y = fuzzy_school_data$AGNO,
x = fuzzy_school_data$TKS_c,
c = 0,
fuzzy = fuzzy_school_data$D,
p = 1,
h = h_optimal,
kernel = "uniform",
weight = fuzzy_school_data$weights,
vce = "hc0")
#> Warning in rdrobust(y = fuzzy_school_data$AGNO, x = fuzzy_school_data$TKS_c, :
#> Mass points detected in the running variable.
summary(fit_rdrobust)
#> Call: rdrobust
#>
#> Fuzzy RD estimates using local polynomial regression.
#>
#> Number of Obs. 100
#> BW type Manual
#> Kernel Uniform
#> VCE method HC0
#>
#> Left Right
#> Number of Obs. 42 58
#> Eff. Number of Obs. 8 20
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 4.506 4.506
#> BW bias (b) 4.506 4.506
#> rho (h/b) 1.000 1.000
#> Unique Obs. 15 15
#>
#> First-stage estimates.
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> =====================================================================
#> Rd Effect 0.433 0.184 0.854 [-1.129 , 1.363]
#> =====================================================================
#>
#> Treatment effect estimates.
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 0.546 1.962 0.050 [0.001 , 1.562]
#> =====================================================================
c("tau" = fit_rdrobust$coef[1],
"se" = fit_rdrobust$se[1])
#> tau se
#> 0.5464854 0.18148619 {rdrobust} with triangular kernel weights
9.1 Estimation
bw <- rdbwselect(y = fuzzy_data$AGNO, # outcome
x = fuzzy_data$TKS, # running variables
c = 30, # cutoff
vce = "cr2", # cluster-robust estimator
p = 1, # polynomial form (1 = linear)
kernel = "triangular",
bwselect = "mserd",
cluster = fuzzy_data$okul_no)$bws
#> Warning in rdbwselect(y = fuzzy_data$AGNO, x = fuzzy_data$TKS, c = 30, vce =
#> "cr2", : Mass points detected in the running variable.
h_optimal <- bw[1,1]
fit_rdrobust <- rdrobust(y = fuzzy_school_data$AGNO,
x = fuzzy_school_data$TKS_c,
c = 0,
fuzzy = fuzzy_school_data$D,
p = 1,
h = h_optimal,
kernel = "triangular",
weight = fuzzy_school_data$weights,
vce = "hc0")
#> Warning in rdrobust(y = fuzzy_school_data$AGNO, x = fuzzy_school_data$TKS_c, :
#> Mass points detected in the running variable.
summary(fit_rdrobust)
#> Call: rdrobust
#>
#> Fuzzy RD estimates using local polynomial regression.
#>
#> Number of Obs. 100
#> BW type Manual
#> Kernel Triangular
#> VCE method HC0
#>
#> Left Right
#> Number of Obs. 42 58
#> Eff. Number of Obs. 10 26
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 5.948 5.948
#> BW bias (b) 5.948 5.948
#> rho (h/b) 1.000 1.000
#> Unique Obs. 15 15
#>
#> First-stage estimates.
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> =====================================================================
#> Rd Effect 0.520 -0.058 0.954 [-1.197 , 1.129]
#> =====================================================================
#>
#> Treatment effect estimates.
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 0.523 2.982 0.003 [0.259 , 1.252]
#> =====================================================================
c("tau" = fit_rdrobust$coef[1],
"se" = fit_rdrobust$se[1])
#> tau se
#> 0.5232477 0.13566209.2 Plot {rdrobust}
9.3 Plot student-level
9.4 Plot school-level
