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 Load and check the data
fuzzy_data <- read.csv("C:/Users/.../fuzzy_data.csv")head(fuzzy_data)
#> okul_no AGNO TKS TKS_c Z D
#> 1 1 2.325117 28 -2 0 0
#> 2 1 2.264310 28 -2 0 0
#> 3 1 2.364069 28 -2 0 0
#> 4 1 2.411930 28 -2 0 0
#> 5 1 2.297855 28 -2 0 0
#> 6 1 2.329940 28 -2 0 0
tail(fuzzy_data)
#> okul_no AGNO TKS TKS_c Z D
#> 3027 100 3.034225 31 1 1 1
#> 3028 100 3.081669 31 1 1 1
#> 3029 100 2.873357 31 1 1 1
#> 3030 100 3.071537 31 1 1 1
#> 3031 100 3.135065 31 1 1 1
#> 3032 100 2.981677 31 1 1 1
aggregate(TKS ~ okul_no,
data = fuzzy_data,
FUN = mean)
#> okul_no TKS
#> 1 1 28
#> 2 2 30
#> 3 3 20
#> 4 4 23
#> 5 5 21
#> 6 6 41
#> 7 7 16
#> 8 8 21
#> 9 9 25
#> 10 10 26
#> 11 11 19
#> 12 12 35
#> 13 13 16
#> 14 14 35
#> 15 15 45
#> 16 16 35
#> 17 17 43
#> 18 18 29
#> 19 19 32
#> 20 20 29
#> 21 21 20
#> 22 22 17
#> 23 23 42
#> 24 24 23
#> 25 25 43
#> 26 26 31
#> 27 27 36
#> 28 28 40
#> 29 29 15
#> 30 30 44
#> 31 31 35
#> 32 32 18
#> 33 33 44
#> 34 34 18
#> 35 35 41
#> 36 36 27
#> 37 37 42
#> 38 38 37
#> 39 39 41
#> 40 40 32
#> 41 41 42
#> 42 42 21
#> 43 43 36
#> 44 44 27
#> 45 45 43
#> 46 46 32
#> 47 47 34
#> 48 48 41
#> 49 49 18
#> 50 50 43
#> 51 51 31
#> 52 52 38
#> 53 53 30
#> 54 54 18
#> 55 55 31
#> 56 56 33
#> 57 57 21
#> 58 58 24
#> 59 59 18
#> 60 60 25
#> 61 61 24
#> 62 62 42
#> 63 63 26
#> 64 64 32
#> 65 65 16
#> 66 66 32
#> 67 67 19
#> 68 68 24
#> 69 69 22
#> 70 70 34
#> 71 71 30
#> 72 72 34
#> 73 73 31
#> 74 74 42
#> 75 75 36
#> 76 76 17
#> 77 77 38
#> 78 78 20
#> 79 79 19
#> 80 80 18
#> 81 81 22
#> 82 82 42
#> 83 83 38
#> 84 84 32
#> 85 85 20
#> 86 86 44
#> 87 87 35
#> 88 88 37
#> 89 89 40
#> 90 90 16
#> 91 91 30
#> 92 92 36
#> 93 93 20
#> 94 94 40
#> 95 95 38
#> 96 96 33
#> 97 97 37
#> 98 98 29
#> 99 99 35
#> 100 100 31AGNO: Öğrencinin Ağırlıklı Genel Not Ortalaması
TKS: Toplam kayıt sayısı
TKS_c: Okulun eşik değerine (30) göre ortalanmış TKS’sı
Z: Okul TKS’sının eşik değerini (30) aşıp aşmadığı
D: Okulda sınıfların bölünüp bölünmediği
3 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(
TKS ~ okul_no,
data = fuzzy_data,
FUN = length)$TKS
weights
#> [1] 28 30 20 23 21 41 16 21 25 26 19 35 16 35 45 35 43 29 32 29 20 17 42 23 43
#> [26] 31 36 40 15 44 35 18 44 18 41 27 42 37 41 32 42 21 36 27 43 32 34 41 18 43
#> [51] 31 38 30 18 31 33 21 24 18 25 24 42 26 32 16 32 19 24 22 34 30 34 31 42 36
#> [76] 17 38 20 19 18 22 42 38 32 20 44 35 37 40 16 30 36 20 40 38 33 37 29 35 31
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
h_optimal <- rdbwselect(y = fuzzy_data$AGNO, # outcome
x = fuzzy_data$TKS, # running variables
c = 30, # cutoff
vce = "cr1", # cluster-robust estimator
p = 1, # polynomial form (1 = linear)
kernel = "uniform", # kernel weights
bwselect = "mserd",
cluster = fuzzy_data$okul_no)$bws[1,1]
#> Warning in rdbwselect(y = fuzzy_data$AGNO, x = fuzzy_data$TKS, c = 30, vce =
#> "cr1", : Mass points detected in the running variable.
# cr1 for CR1 cluster-robust variance estimator with degrees-of-freedom correction at the cluster level
# cr2 for CR2 cluster-robust variance estimator with leverage adjustments (Bell McCaffrey).
# 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 Method 1
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:
\[ {\color{red}{\tau}} = {\color{red}{\beta_1}} / {\color{red}{\alpha_1}} \]
6.1.2 Method 2
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:
\[ {\color{red}{\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", # kernel weights
cluster = fuzzy_data$okul_no,
vce = "cr1") # 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 CR1
#>
#> 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.560 4.560
#> BW bias (b) 4.560 4.560
#> 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.172 0.864 [-1.216 , 1.449]
#> =====================================================================
#>
#> Treatment effect estimates.
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 0.546 1.829 0.067 [-0.056 , 1.619]
#> =====================================================================
c("tau" = fit_rdrobust$coef[1],
"se" = fit_rdrobust$se[1])
#> tau se
#> 0.5464854 0.19354308.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.560 4.560
#> BW bias (b) 4.560 4.560
#> 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
h_optimal <- rdbwselect(y = fuzzy_data$AGNO, # outcome
x = fuzzy_data$TKS, # running variables
c = 30, # cutoff
vce = "cr1", # cluster-robust estimator
p = 1, # polynomial form (1 = linear)
kernel = "triangular",
bwselect = "mserd",
cluster = fuzzy_data$okul_no)$bws[1,1]
#> Warning in rdbwselect(y = fuzzy_data$AGNO, x = fuzzy_data$TKS, c = 30, vce =
#> "cr1", : Mass points detected in the running variable.
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. 13 30
#> Order est. (p) 1 1
#> Order bias (q) 2 2
#> BW est. (h) 6.021 6.021
#> BW bias (b) 6.021 6.021
#> 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.526 -0.023 0.982 [-1.155 , 1.129]
#> =====================================================================
#>
#> Treatment effect estimates.
#>
#> =====================================================================
#> Point Robust Inference
#> Estimate z P>|z| [ 95% C.I. ]
#> ---------------------------------------------------------------------
#> RD Effect 0.518 3.159 0.002 [0.285 , 1.218]
#> =====================================================================
c("tau" = fit_rdrobust$coef[1],
"se" = fit_rdrobust$se[1])
#> tau se
#> 0.5180220 0.12999369.2 Plot {rdrobust}
9.3 Plot student-level
9.4 Plot school-level
