install.packages(c("rdrobust", "rddensity", "lmtest", "sandwich"))Sharp Regression Discontinuity Design
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'1 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
2 Simple mean difference
2.1 Regression equation
\[ YKS = \beta_0 + \beta_1 D + \epsilon \]
2.2 Estimation
m1 <- lm(YKS ~ D, # model
data = sharp_data) # data
coef(summary(m1))
#> 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-442.3 Plotting
Code
# plot the model
plot(x = sharp_data$AGNO, # running variable
y = sharp_data$YKS, # outcome
col = rgb(148/255, 163/255, 184/255, alpha = 0.75),
main = "YKS ~ D",
xlab = "AGNO (atama değişkeni)",
ylab = "YKS_puanı", pch = 19)
abline(v = 70, lty = 2, col = "red")
x_left <- seq(from = min(sharp_data$AGNO),
to = 70,
length.out = 100)
y_left <- predict(object = m1,
newdata = data.frame(AGNO = x_left, D = 0))
lines(x = x_left,
y = y_left,
col = "blue",
lwd = 2)
x_right <- seq(from = 70,
to = max(sharp_data$AGNO),
length.out = 100)
y_right <- predict(object = m1,
newdata = data.frame(AGNO = x_right, D = 1))
lines(x = x_right,
y = y_right,
col = "blue",
lwd = 2)
3 Linear functional form - same regression line on both sides
3.1 Regression equation
\[ YKS = \beta_0 + \beta_1 D + \beta_2 AGNO + \epsilon \]
3.2 Estimation
m2 <- lm(YKS ~ D + AGNO,
data = sharp_data)
coef(summary(m2))
#> 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-153.3 Plotting
Code
plot(x = sharp_data$AGNO,
y = sharp_data$YKS,
col = rgb(148/255, 163/255, 184/255, alpha = 0.75),
main = "YKS ~ D + AGNO",
xlab = "AGNO (atama değişkeni)",
ylab = "YKS_puanı", pch = 19)
abline(v = 70, lty = 2, col = "red")
x_left <- seq(from = min(sharp_data$AGNO),
to = 70,
length.out = 100)
y_left <- predict(object = m2,
newdata = data.frame(AGNO = x_left, D = 0))
lines(x = x_left,
y = y_left,
col = "blue",
lwd = 2)
x_right <- seq(from = 70,
to = max(sharp_data$AGNO),
length.out = 100)
y_right <- predict(object = m2,
newdata = data.frame(AGNO = x_right, D = 1))
lines(x = x_right,
y = y_right,
col = "blue",
lwd = 2)
4 Linear functional form - different regression lines on both sides
4.1 Regression equation
\[ YKS = \beta_0 + \beta_1 D + \beta_2 AGNO_c + \beta_3 D * AGNO_c + \epsilon \]
4.2 Estimation
m3 <- lm(YKS ~ D + AGNO_c + D*AGNO_c,
data = sharp_data)
coef(summary(m3))
#> 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
#> D:AGNO_c -1.169992 0.3200142 -3.656064 4.179514e-044.3 Plotting
Code
plot(x = sharp_data$AGNO_c,
y = sharp_data$YKS,
col = rgb(148/255, 163/255, 184/255, alpha = 0.75),
main = "YKS ~ D + AGNO_c + D * AGNO_c",
xlab = "AGNO (atama değişkeni - ortalanmış)",
ylab = "YKS_puanı", pch = 19)
abline(v = 0, lty = 2, col = "red")
x_left <- seq(from = min(sharp_data$AGNO_c),
to = 0,
length.out = 100)
y_left <- predict(m3,
newdata = data.frame(AGNO_c = x_left, D = 0))
lines(x = x_left,
y = y_left,
col = "blue",
lwd = 2)
x_right <- seq(from = 0,
max(sharp_data$AGNO_c), length.out = 100)
y_right <- predict(object = m3,
newdata = data.frame(AGNO_c = x_right, D = 1))
lines(x_right, y_right, col = "blue", lwd = 2)
5 Quadratic functional form - different regression lines on both sides
5.1 Regression equation
\[ YKS = \beta_0 + \beta_1 D + \beta_2 AGNO_c + \beta_3 AGNO_c^2 + \beta_4 D * AGNO_c + \beta_5 D * AGNO_c^2 + \epsilon \]
5.2 Estimation
m4 <- lm(YKS ~ D + AGNO_c + I(AGNO_c^2) + D*AGNO_c + D*I(AGNO_c^2),
data = sharp_data)
coef(summary(m4))
#> 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
#> D:AGNO_c -1.46493163 1.01144536 -1.4483547 1.508454e-01
#> D:I(AGNO_c^2) -0.02677161 0.04705384 -0.5689570 5.707423e-015.3 Plotting
Code
plot(x = sharp_data$AGNO_c,
y = sharp_data$YKS,
col = rgb(148/255, 163/255, 184/255, alpha = 0.75),
main = "YKS ~ D + AGNO_c + AGNO_c^2 + \n D*AGNO_c + D*AGNO_c^2",
xlab = "AGNO (atama değişkeni - ortalanmış)",
ylab = "YKS_puanı", pch = 19)
abline(v = 0,
lty = 2,
col = "red")
x_left <- seq(from = min(sharp_data$AGNO_c),
to = 0,
length.out = 100)
y_left <- predict(object = m4,
newdata = data.frame(AGNO_c = x_left, D = 0))
lines(x = x_left,
y = y_left,
col = "blue",
lwd = 2)
x_right <- seq(from = 0,
to = max(sharp_data$AGNO_c),
length.out = 100)
y_right <- predict(object = m4,
newdata = data.frame(AGNO_c = x_right, D = 1))
lines(x = x_right,
y = y_right,
col = "blue",
lwd = 2)
6 {rdrobust} bandwidth, estimation and plotting
6.1 Linear functional form
6.1.1 Bandwidth
bw_m3 <- 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
# mserd: one common MSE-optimal bandwidth on two sides
# msetwo: two separate MSE-optimal bandwidth on two sides6.1.2 Estimation
rd_m3 <- 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_m3[1,1]) # bandwidth
summary(rd_m3)
#> 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" = rd_m3$coef[1], "se" = rd_m3$se[1])
#> tau se
#> 39.431198 5.0555296.1.3 Plotting
Code
plot(rd_m3,
y = sharp_data$YKS,
x_run = sharp_data$AGNO)
6.2 Quadratic functional form
6.2.1 Bandwidth
bw_m4 <- 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
# mserd: one common MSE-optimal bandwidth on two sides
# msetwo: two separate MSE-optimal bandwidth on two sides6.2.2 Estimation
rd_m4 <- 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_m4[1,1]) # bandwidth
summary(rd_m4)
#> 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" = rd_m4$coef[1], "se" = rd_m4$se[1])
#> tau se
#> 41.204934 6.1034586.2.3 Plotting
Code
plot(rd_m4,
y = sharp_data$YKS,
x_run = sharp_data$AGNO)
7 Linear functional form with interaction and optimal bandwidth
7.1 Bandwidth and kernel weights
# function to calculate kernel weights
triangular_kernel_weights <- function(x, cutoff = 0, bandwidth) {
# standardized distance u
u <- (x - cutoff) / bandwidth
# triangular kernel evaluation K(u)
K <- ifelse(abs(u) <= 1, 1 - abs(u), 0)
# normalized weights (sum to 1)
weights <- K / sum(K)
return(weights)
} # triangular_kernel_weights
# calculate Imbens-Kalyanaraman optimal bandwidth
IKbandwidth(Y = sharp_data$YKS,
X = sharp_data$AGNO,
cutpoint = 70,
kernel = "triangular")
#> [1] 7.239205
# or get the bandwidth from the rdrobust package
h_opt_m3 <- bw_m3[1,1]
cat("Optimal Bant Genişliği (h):", h_opt_m3, "\n")
#> Optimal Bant Genişliği (h): 5.372379
sharp_data_m3_bw <- subset(x = sharp_data,
subset = AGNO >= (70 - h_opt_m3) & AGNO <= (70 + h_opt_m3))
# calculate triangular kernel weights
tkweights <- triangular_kernel_weights(x = sharp_data_m3_bw$AGNO_c,
cutoff = 0,
bandwidth = h_opt_m3)7.2 Plot kernel weights
Code
plot(y = tkweights,
x = sharp_data_m3_bw$AGNO_c,
ylab = "Triangular Kernel Weights",
xlab = "AGNO_c")
7.3 Estimation
Code
# estimate the model
m3_bw <- lm(YKS ~ D + AGNO_c + D*AGNO_c,
weights = tkweights,
data = sharp_data_m3_bw)
print(coeftest(m3_bw, vcov = vcovHC(m3_bw, 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 ***
#> D:AGNO_c -3.90765 1.86415 -2.0962 0.04244 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
c("tau" = rd_m3$coef[1], "se" = rd_m3$se[1]) # compare to rdrobust
#> tau se
#> 39.431198 5.0555297.4 Plotting
Code
# plot the model
# weight the observations on the plot as they were weighted in the model
abs_agno <- abs(sharp_data$AGNO_c)
tkweights_full <- ifelse(abs_agno <= h_opt_m3, 1 - (abs_agno / h_opt_m3), 0)
norm_weights <- (tkweights_full - min(tkweights_full)) / (max(tkweights_full) - min(tkweights_full))
point_sizes <- 0.8 + (norm_weights * 1.6)
plot(x = sharp_data$AGNO_c,
y = sharp_data$YKS,
col = ifelse(abs(sharp_data$AGNO_c) <= h_opt_m3,
rgb(148/255, 163/255, 184/255, alpha = 0.75),
"pink"),
pch = 19,
cex = point_sizes,
main = "YKS ~ D + AGNO_c + D * AGNO_c",
xlab = "AGNO (atama değişkeni - ortalanmış)",
ylab = "YKS_puanı")
abline(v = 0,
lty = 2,
col = "red")
abline(v = -h_opt_m3,
lty = 3,
col = "gray40")
abline(v = h_opt_m3,
lty = 3,
col = "gray40")
grid_left <- seq(from = -h_opt_m3,
to = 0,
length.out = 100)
pred_left <- predict(object = m3_bw,
newdata = data.frame(AGNO_c = grid_left, D = 0))
lines(x = grid_left,
y = pred_left,
col = "blue",
lwd = 2)
grid_right <- seq(from = 0,
to = h_opt_m3,
length.out = 100)
pred_right <- predict(object = m3_bw,
newdata = data.frame(AGNO_c = grid_right, D = 1))
lines(x = grid_right,
y = pred_right,
col = "blue",
lwd = 2)
8 Bandwidth sensitivity analysis
8.1 Linear functional form with interaction and half optimal bandwidth
8.1.1 Re-estimate kernel weights
h_opt_m3 <- bw_m3[1,1] / 2
cat("Optimal Bant Genişliği (h):", h_opt_m3, "\n")
#> Optimal Bant Genişliği (h): 2.686189
sharp_data_m3_bw <- subset(x = sharp_data,
subset = AGNO >= (70 - h_opt_m3) & AGNO <= (70 + h_opt_m3))
# calculate triangular kernel weights
tkweights <- triangular_kernel_weights(x = sharp_data_m3_bw$AGNO_c,
cutoff = 0,
bandwidth = h_opt_m3)8.1.2 Estimation
m3_bw <- lm(YKS ~ D + AGNO_c + D*AGNO_c,
weights = tkweights,
data = sharp_data_m3_bw)
print(coeftest(m3_bw, vcov = vcovHC(m3_bw, 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 *
#> D:AGNO_c -7.2880 6.7937 -1.0728 0.29682
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 18.1.3 Plotting
Code
# weight the observations on the plot as they were weighted in the model
abs_agno <- abs(sharp_data$AGNO_c)
tkweights_full <- ifelse(abs_agno <= h_opt_m3, 1 - (abs_agno / h_opt_m3), 0)
norm_weights <- (tkweights_full - min(tkweights_full)) / (max(tkweights_full) - min(tkweights_full))
point_sizes <- 0.8 + (norm_weights * 1.6)
plot(x = sharp_data$AGNO_c,
y = sharp_data$YKS,
col = ifelse(abs(sharp_data$AGNO_c) <= h_opt_m3,
rgb(148/255, 163/255, 184/255, alpha = 0.75),
"pink"),
pch = 19,
cex = point_sizes,
main = "YKS ~ D + AGNO_c + D * AGNO_c",
xlab = "AGNO (atama değişkeni - ortalanmış)",
ylab = "YKS_puanı")
abline(v = 0,
lty = 2,
col = "red")
abline(v = -h_opt_m3,
lty = 3,
col = "gray40")
abline(v = h_opt_m3,
lty = 3,
col = "gray40")
grid_left <- seq(from = -h_opt_m3,
to = 0,
length.out = 100)
pred_left <- predict(object = m3_bw,
newdata = data.frame(AGNO_c = grid_left, D = 0))
lines(x = grid_left,
y = pred_left,
col = "blue",
lwd = 2)
grid_right <- seq(from = 0,
to = h_opt_m3,
length.out = 100)
pred_right <- predict(object = m3_bw,
newdata = data.frame(AGNO_c = grid_right, D = 1))
lines(x = grid_right,
y = pred_right,
col = "blue",
lwd = 2)
8.2 Linear functional form with interaction and twice optimal bandwidth
8.2.1 Re-estimate kernel weights
h_opt_m3 <- 2 * bw_m3[1,1]
cat("Optimal Bant Genişliği (h):", h_opt_m3, "\n")
#> Optimal Bant Genişliği (h): 10.74476
sharp_data_m3_bw <- subset(x = sharp_data,
subset = AGNO >= (70 - h_opt_m3) & AGNO <= (70 + h_opt_m3))
# calculate triangular kernel weights
tkweights <- triangular_kernel_weights(x = sharp_data_m3_bw$AGNO_c,
cutoff = 0,
bandwidth = h_opt_m3)8.2.2 Estimation
m3_bw <- lm(YKS ~ D + AGNO_c + D*AGNO_c,
weights = tkweights,
data = sharp_data_m3_bw)
print(coeftest(m3_bw, vcov = vcovHC(m3_bw, 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 ***
#> D:AGNO_c -1.2346 0.8324 -1.4831 0.1422
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 18.2.3 Plotting
Code
# plot the model
# weight the observations on the plot as they were weighted in the model
abs_agno <- abs(sharp_data$AGNO_c)
tkweights_full <- ifelse(abs_agno <= h_opt_m3, 1 - (abs_agno / h_opt_m3), 0)
norm_weights <- (tkweights_full - min(tkweights_full)) / (max(tkweights_full) - min(tkweights_full))
point_sizes <- 0.8 + (norm_weights * 1.6)
plot(x = sharp_data$AGNO_c,
y = sharp_data$YKS,
col = ifelse(abs(sharp_data$AGNO_c) <= h_opt_m3,
rgb(148/255, 163/255, 184/255, alpha = 0.75),
"pink"),
pch = 19,
cex = point_sizes,
main = "YKS ~ D + AGNO_c + D * AGNO_c",
xlab = "AGNO (atama değişkeni - ortalanmış)",
ylab = "YKS_puanı")
abline(v = 0,
lty = 2,
col = "red")
abline(v = -h_opt_m3,
lty = 3,
col = "gray40")
abline(v = h_opt_m3,
lty = 3,
col = "gray40")
grid_left <- seq(from = -h_opt_m3,
to = 0,
length.out = 100)
pred_left <- predict(object = m3_bw,
newdata = data.frame(AGNO_c = grid_left, D = 0))
lines(x = grid_left,
y = pred_left,
col = "blue",
lwd = 2)
grid_right <- seq(from = 0,
to = h_opt_m3,
length.out = 100)
pred_right <- predict(object = m3_bw,
newdata = data.frame(AGNO_c = grid_right, D = 1))
lines(x = grid_right,
y = pred_right,
col = "blue",
lwd = 2)
9 Quadratic functional form with interaction and optimal bandwidth
9.1 Bandwidth and kernel weights
# get the bandwidth from the rdrobust package
h_opt_m4 <- bw_m4[1,1]
cat("Optimal Bant Genişliği (h):", h_opt_m4, "\n")
#> Optimal Bant Genişliği (h): 7.402737
sharp_data_m4_bw <- subset(x = sharp_data,
subset = AGNO >= (70 - h_opt_m4) & AGNO <= (70 + h_opt_m4))
# calculate triangular kernel weights
tkweights <- triangular_kernel_weights(x = sharp_data_m4_bw$AGNO_c,
cutoff = 0,
bandwidth = h_opt_m4)9.2 Estimation
m4_bw <- lm(YKS ~ D + AGNO_c + I(AGNO_c^2) + D*AGNO_c + D*I(AGNO_c^2),
weights = tkweights,
data = sharp_data_m4_bw)
print(coeftest(m4_bw, vcov = vcovHC(m4_bw, 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
#> D:AGNO_c -8.51092 4.29929 -1.9796 0.05242 .
#> D:I(AGNO_c^2) 0.55130 0.61159 0.9014 0.37103
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
c("tau" = rd_m4$coef[1], "se" = rd_m4$se[1]) # compare to rdrobust
#> tau se
#> 41.204934 6.1034589.3 Plotting
Code
# plot the model
abs_agno <- abs(sharp_data$AGNO_c)
tkweights_full <- ifelse(abs_agno <= h_opt_m4, 1 - (abs_agno / h_opt_m4), 0)
norm_weights <- (tkweights_full - min(tkweights_full)) / (max(tkweights_full) - min(tkweights_full))
point_sizes <- 0.8 + (norm_weights * 1.6)
plot(x = sharp_data$AGNO_c,
y = sharp_data$YKS,
col = ifelse(abs(sharp_data$AGNO_c) <= h_opt_m4,
rgb(148/255, 163/255, 184/255, alpha = 0.75),
"pink"),
pch = 19,
cex = point_sizes,
main = "YKS ~ D + AGNO_c + AGNO_c^2 + \n D*AGNO_c + D*AGNO_c^2",
xlab = "AGNO (atama değişkeni - ortalanmış)",
ylab = "YKS_puanı")
abline(v = 0,
lty = 2,
col = "red")
abline(v = -h_opt_m4,
lty = 3,
col = "gray40")
abline(v = h_opt_m4,
lty = 3,
col = "gray40")
grid_left <- seq(from = -h_opt_m4,
to = 0,
length.out = 100)
pred_left <- predict(object = m4_bw,
newdata = data.frame(AGNO_c = grid_left, D = 0))
lines(x = grid_left,
y = pred_left,
col = "blue",
lwd = 2)
grid_right <- seq(from = 0,
to = h_opt_m4,
length.out = 100)
pred_right <- predict(object = m4_bw,
newdata = data.frame(AGNO_c = grid_right, D = 1))
lines(x = grid_right,
y = pred_right,
col = "blue",
lwd = 2)