Biraz Rahatlayalım: Regresyon Discontinuity Tasarımları

Dünyamız, insanların birbirlerine uyguladığı kurallar ve kanunlar üzerine kurulu. Bu kurallar ve kanunlar bazen keyfi şekillerde oluşturulabiliyor. Ancak bu keyfi kararlar bize deney yapma imkanı sağlayabiliyor.

RD (regresyon discontinuity) tasarımı, bu kurallardan faydalanan bir yöntemdir. Fuzzy ve sharp olmak üzere iki farklı şekli vardır. Sharp tasarımı, gözlemlediğimiz bir olayda kuralların tedavi ve kontrol gruplarını kesin olarak belirlemek üzerine kuruludur. Fuzzy tasarım ise bu grupları ayırabilmek için araç değişkenler gibi bir yapıya ihtiyaç duyar.

Kitabın 6.1.1 grafiğinin tekrarı (replication)

6.1.1 üç ayrı grafikten oluşuyor. Doğrusal, doğrusal olmayan, süreksiz olmakla karıştırılan doğrusal olmayan RD tasarımları. İlk olarak bu üç tasarım için veri oluşturuyoruz.

1. adım: Her biri veri için 100 gözlem belirleyin.

n=100

2. adım: Rastgele x verisi oluşturun.

set.seed(1149)
x <- runif(n)

head(x)

## [1] 0.901608043 0.508481772 0.003583934 0.644593220 0.051421000 0.273405757

3. adım: koşullara göre 3 tane y oluşturun. Eşik değeri 0.5 olarak belirlenmiştir.

y_dogrusal <- x + (x > 0.5) * 0.25 + rnorm(n, mean = 0, sd = 0.1)
y_dogrusal_olmayan <- 0.5 * sin(6 * (x - 0.5)) + 0.5 + (x > 0.5) * 0.25 + rnorm(n, mean = 0, sd = 0.1)
y_hatalı <- 1 / (1 + exp(-25 * (x - 0.5))) + rnorm(n, mean = 0, sd = 0.1)

x ve Y’yi doğrusal bir grafik çizmek için bir data frame içinde birleştirin.

df <- data.frame(x = x, y_dogrusal = y_dogrusal)

Bu dataframe’i eşik değerine göre ikiye bölelim.

df_alt <- df[df$x < 0.5, ]
df_ust <- df[df$x > 0.5, ]

1. grafik, doğrusal sharp tasarım

library(plotrix)

A <- plot(df_alt$x, df_alt$y_dogrusal, xlim = c(0, 1), ylim = c(0, 2), pch = 16, cex = 0.5, col = "black", main = "A. Doğrusal E[Y|X]", xlab = "x", ylab = "Y")
points(df_ust$x, df_ust$y_dogrusal, pch = 16, cex = 0.5, col = "black")
lm_alt <- lm(y_dogrusal ~ x, data = df_alt)
lm_ust <- lm(y_dogrusal ~ x, data = df_ust)
ablineclip(lm_alt, x1 = 0,x2 = .5, col = "black", lwd = 2)
ablineclip(lm_ust, x1 = .5,x2 = 1, col = "black", lwd = 2)
abline(v = 0.5, lty = 2)

2. Grafik için doğrusal olmayan tasarım.

df_dogrusal_olmayan <- data.frame(x = x, y_dogrusal_olmayan = y_dogrusal_olmayan)

Dogrusal olmayan veri setini ikiye ayıralım

df_dogrusal_olmayan_alt <- df_dogrusal_olmayan[df_dogrusal_olmayan$x < 0.5, ]
df_dogrusal_olmayan_ust <- df_dogrusal_olmayan[df_dogrusal_olmayan$x > 0.5, ]

Bu sefer plot fonksiyonu yerine ggplot2 paketini kullanabiliriz. (Not: Bu ders notunun amacı kullanılan paketleri öğretmek değil. Paketlerin detayları için ayrı bir ders hazırlanabilir.)

library(ggplot2)

B<- ggplot() +
  geom_smooth(data = df_dogrusal_olmayan_alt, aes(x = x, y = y_dogrusal_olmayan), method = "lm", formula = y ~ poly(x, 2), se = FALSE, color = "black")+
  geom_point(data = df_dogrusal_olmayan_alt, aes(x = x, y = y_dogrusal_olmayan), size = 1, shape = 16, color = "black") +
  geom_smooth(data = df_dogrusal_olmayan_ust, aes(x = x, y = y_dogrusal_olmayan), method = "lm", formula = y ~ poly(x, 2), se = FALSE, color = "black") +
  geom_point(data = df_dogrusal_olmayan_ust, aes(x = x, y = y_dogrusal_olmayan), size = 1, shape = 16, color = "black") +
  labs(title = "B. Doğrusal Olmayan E[Y|X]", y = "Y", x = "x") +
  theme_bw() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  xlim(0, 1)

3. grafik RD ile karıştırılma ihtimali olan tasarım.

df_hatalı <- data.frame(x = x, y_hatalı = y_hatalı)

df_hatalı_alt <- df_hatalı[df_hatalı$x < 0.5, ]
df_hatalı_ust <- df_hatalı[df_hatalı$x >= 0.5, ]

Quadratic fonksiyon yerine, bir süreksizlik fonksiyonu belirlenmiş,

discontinuity <- function(x) {
  1 / (1 + exp(-25 * (x - 0.5)))
}

C<- ggplot() +
  geom_smooth(data = df_hatalı_alt, aes(x = x, y = y_hatalı), method = "lm", se = FALSE, color = "black") +
  geom_smooth(data = df_hatalı_ust, aes(x = x, y = y_hatalı), method = "lm", se = FALSE, color = "black") +
  geom_function(fun = discontinuity, linetype = "dashed") +
  geom_point(data = df_hatalı, aes(x = x, y = y_hatalı), size = 1, shape = 16, color = "black") +
  labs(title = "C. Doğrusal olmama durumunun süreksiz olmakla karıştırılması", y = "Y", x = "x") +
  theme_bw() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  xlim(0, 1)

Bu gösterim yerine, 2 grafiği birleştirerek göstermek isterseniz cowplot paketini kullanabilirsiniz.

library(cowplot)

cowplot::plot_grid(B, C, ncol = 1, align = "h")

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

Farklı Trend Fonksiyonları

\(E[Y_{0i} | x_i]\) ve \(E[Y_{1i} | x_i]\) için farklı trend fonksiyonları yazılabilir.

\[ E[Y_{0i} | x_i] = f_0(x_i) = \alpha + \beta_{01} \tilde{x_i} + \beta_{02} \tilde{x_i}^2 + .... + \beta_{0p} \tilde{x_i}^p \]

\[ E[Y_{1i} | x_i] = f_1(x_i) = \alpha + \rho + \beta_{11} \tilde{x_i} + \beta_{12} \tilde{x_i}^2 + .... + \beta_{1p} \tilde{x_i}^p \]

burada \(\tilde{x_i} = x_i - x_0\), yani \(x_i\) \(x_0\)’nun merkezinde normalleştirilmiştir böylece \(D_i\)’ın katsayısı tedavi (treatment effect) etkisini verir.

Regresyon modelini türetirsek,

\[ E[Y_{i} | x_i] = E[Y_{0i} | x_i] + (E[Y_{1i} | x_i] - E[Y_{0i} | x_i]) \cdot D_i \] \[ E[Y_{i} | x_i] = \alpha + \beta_{01} \tilde{x_i} + \beta_{02} \tilde{x_i}^2 + .... + \beta_{0p} \tilde{x_i}^p + (\alpha + \rho + \beta_{11} \tilde{x_i} + \beta_{12} \tilde{x_i}^2 + .... + \beta_{1p} \tilde{x_i}^p - \alpha + \beta_{01} \tilde{x_i} + \beta_{02} \tilde{x_i}^2 + .... + \beta_{0p} \tilde{x_i}^p) \cdot D_i \]

\[ E[Y_{i} | x_i] = \alpha + \beta_{01} \tilde{x_i} + \beta_{02} \tilde{x_i}^2 + .... + \beta_{0p} \tilde{x_i}^p + \rho + \beta_{11} \tilde{x_i}\cdot D_i + \beta_{12} \tilde{x_i}^2\cdot D_i + .... + \beta_{1p} \tilde{x_i}^p\cdot D_i - \beta_{01} \tilde{x_i}\cdot D_i - \beta_{02} \tilde{x_i}^2\cdot D_i - .... - \beta_{0p} \tilde{x_i}^p\cdot D_i \]

\[ E[Y_{i} | x_i] = \alpha + \beta_{01} \tilde{x_i} + \beta_{02} \tilde{x_i}^2 + .... + \beta_{0p} \tilde{x_i}^p + \rho + (\beta_{11}- \beta_{01}) \tilde{x_i}\cdot D_i + (\beta_{12}- \beta_{02}) \tilde{x_i}^2\cdot D_i + .... + (\beta_{1p} - \beta_{0p}) \tilde{x_i}^p\cdot D_i \] {}

\[ E[Y_{i} | x_i] = \alpha + \beta_{01} \tilde{x_i} + \beta_{02} \tilde{x_i}^2 + .... + \beta_{0p} \tilde{x_i}^p + \rho + \beta_1^{\ast} \tilde{x_i}\cdot D_i + \beta_2^{\ast} \tilde{x_i}^2\cdot D_i + .... + \beta_p^{\ast} \tilde{x_i}^p\cdot D_i + \eta_i \]

Şekil 6.1.2 tekrarı (replication)

Mostly Harmless Econometrics kitabındaki Şekil 6.1.2 örneği, Lee’nin (2008) parti görevinin yeniden seçilme olasılıkları üzerindeki etkisini incelediği keskin RD tasarımını göstermektedir. Lee, ABD Temsilciler Meclisi’nde bir koltuk için Demokrat adayın partisinin geçen sefer sandalyeyi kazanması durumunda bir sonraki seçim için avantajının olup olmadığıyla ilgilenmektedir, Temsilciler Meclisi üyelerinin dikkat çeken başarısı, temsilcilerin ofislerini ve kaynaklarını kendileri veya partileri için avantaj elde etmek için kullanıp kullanmadıkları sorusunu gündeme getirmektedir. Bu varsayım makul görünmektedir, ancak görevde olanların başarısının gerçek bir seçim avantajını yansıtması gerekmemektedir. Görevde olanlar, adaylar ve partiler olarak kazanma potansiyeli gösteren kişilerdir ve seçmenleri memnun etme veya oyları geri alma konusunda daha başarılı olabilirler.

Lee, görevde olmanın nedensel etkisini yakalamak için bir önceki seçimdeki göreli oy paylarını kazanma olasılığına bağlı olarak inceliyor. Özellikle, bir seçimi kazananın \(D_i = 1 (x_i \geq 0)\) şeklinde belirlendiği gerçeğinden yararlanıyor; burada \(x_i\), kazananın oy payı marjını temsil etmektedir. \(D_i\), \(x_i\)’nin deterministik bir fonksiyonu olduğundan, \(x_i\) dışında başka bir karıştırıcı değişkenin olmadığına dikkat edilmelidir. Bu, RD tasarımının sinyal özelliğidir.

1- Gerekli paketler

library(haven)
library(data.table)
library(ggplot2)
library(gridExtra)

Verisetini https://economics.mit.edu/people/faculty/josh-angrist/mhe-data-archive adresinden indirebilirsiniz. Bu linke gittiğinizde Lee (2008)’in altında Lee2008.zip dosyası indirilebilir.

download.file('https://economics.mit.edu/people/faculty/josh-angrist/mhe-data-archive', 'Lee2008.zip')

lee <- data.table(read_dta('individ_final.dta'))

library(dplyr)
library(DT)
datatable(lee)

2- Verisetini sadece 2 değişkene indirin. Oy farkı (Running variable, x, ve sonuç değişkeni (Y).

df <- lee %>%
  select(myoutcomenext, difshare) %>%
  na.omit()

head(df)

##    myoutcomenext    difshare
## 1:             0  0.06148815
## 2:             0 -0.06148815
## 3:             1  0.10486948
## 4:             0 -0.10486948
## 5:             0 -0.53572100
## 6:             0  0.16446409

3- D değişkenini oluşturun.

df <- df %>%
  mutate(D = as.numeric(difshare >= 0))

head(df)

##    myoutcomenext    difshare D
## 1:             0  0.06148815 1
## 2:             0 -0.06148815 0
## 3:             1  0.10486948 1
## 4:             0 -0.10486948 0
## 5:             0 -0.53572100 0
## 6:             0  0.16446409 1

D değişkeni, Oy farkı (difshare) pozitifse 1 değilse 0 şeklinde oluşturuldu.

4- 4. dereceden logit regresyonu tahmin edin.

logit <- glm(formula = myoutcomenext ~ poly(difshare, degree = 4) +
               poly(difshare, degree = 4) * D,
             family = binomial(link = "logit"),
             data = df)
df <- df %>%
  mutate(pmyoutcomenext = predict(logit, type = "response"))

head(df)

##    myoutcomenext    difshare D pmyoutcomenext
## 1:             0  0.06148815 1    0.695785959
## 2:             0 -0.06148815 0    0.057708113
## 3:             1  0.10486948 1    0.738774903
## 4:             0 -0.10486948 0    0.029672513
## 5:             0 -0.53572100 0    0.001221896
## 6:             0  0.16446409 1    0.778284984

5- 0.005 aralıklarla yerel ortalamalar oluşturun.

breaks <- round(seq(-1, 1, by = 0.005), 3)
df <- df  %>%
  mutate(i005 = cut(difshare, breaks = breaks, labels = head(breaks, -1), right = TRUE),
         i005 = as.numeric(as.character(i005))) %>%
  group_by(i005) %>%
  mutate(m_next = mean(myoutcomenext),
         mp_next = mean(pmyoutcomenext)) %>%
  ungroup() %>%
  filter(i005 > -0.251 & i005 < 0.251)

head(df)

## # A tibble: 6 × 7
##   myoutcomenext difshare     D pmyoutcomenext   i005 m_next mp_next
##           <dbl>    <dbl> <dbl>          <dbl>  <dbl>  <dbl>   <dbl>
## 1             0   0.0615     1         0.696   0.06  0.72    0.697 
## 2             0  -0.0615     0         0.0577 -0.065 0.0594  0.0570
## 3             1   0.105      1         0.739   0.1   0.678   0.737 
## 4             0  -0.105      0         0.0297 -0.105 0.0556  0.0308
## 5             0   0.164      1         0.778   0.16  0.792   0.777 
## 6             0  -0.164      0         0.0130 -0.165 0       0.0133

6- Grafikler

KO <- ggplot(data = df, aes(x = i005)) +
  geom_point(aes(y = m_next)) +
  geom_line(aes(y = mp_next, group = i005 >= 0)) +
  geom_vline(xintercept = 0, linetype = 'longdash') +
  xlab('Demokratların kazanma oy oranları marjini, t seçimi') +
  ylab('Kazanma Olasılığı, t+1 Seçimi') +
  ggtitle('A')

KO

2. grafik için iki ayrı grup oluşturun

df2 <- lee %>%
  mutate(i005 = cut(difshare, breaks = breaks, labels = head(breaks, -1), right = TRUE),
         i005 = as.numeric(as.character(i005))) %>%
  group_by(i005) %>%
  summarize(m_vic = mean(mofficeexp, na.rm = TRUE),
            mp_vic = mean(mpofficeexp, na.rm = TRUE)) %>%
  ungroup() %>%
  filter(i005 > -0.251 & i005 < 0.251)

head(df2)

## # A tibble: 6 × 3
##     i005  m_vic mp_vic
##    <dbl>  <dbl>  <dbl>
## 1 -0.25  0.0513 0.0995
## 2 -0.245 0.315  0.106 
## 3 -0.24  0      0.114 
## 4 -0.235 0      0.121 
## 5 -0.23  0.0698 0.130 
## 6 -0.225 0      0.138

OM <- ggplot(data = df2, aes(x = i005)) +
  geom_point(aes(y = m_vic)) +
  geom_line(aes(y = mp_vic, group = i005 >= 0)) +
  geom_vline(xintercept = 0, linetype = 'longdash') +
  xlab('Demokratların kazanma oy oranları marjini, t seçimi') +
  ylab('t seçimiyle birlikte geçmiş kazanma sayısı') +
  ggtitle('B')

OM

lee.p <- grid.arrange(KO, OM, ncol = 1)

Bulanık Kesikli Rassal Deney (Fuzzy RD)

Bulanık Kesikli Rassal Deney (Fuzzy RD), tedavinin olasılığı veya beklenen değerindeki süreksizliklerden faydalanarak bir ortak değişkene bağlı olarak tedavinin etkisini değerlendirir. Bu araştırma tasarımı, tedavi durumunu deterministik olarak açıp kapatmak yerine, süreksizliğin tedavi durumu için araçsal bir değişkene dönüştüğü bir yaklaşımı benimser. Bunun nasıl çalıştığını anlamak için, \(D_i\)’nin tedavi durumunu önceki gibi temsil etmesine izin verelim, ancak burada \(D_i\) artık deterministik olarak eşik geçiş kuralıyla ilişkili değildir. Bunun yerine, tedavi olasılığı eşik değerinde sıçrama gösterir.

STAR deneyi

Eğitim alanında öncü bir rastgele atama çalışması olan Tennessee STAR Deneyi, ilkokuldaki daha küçük sınıfların etkilerini tahmin etmek amacıyla tasarlanmıştır. Eğitim ekonomistleri, sınıf ortamının özellikleriyle çocukların öğrenme arasında nedensel bağlantılar kurmaya yönelik uzun bir geleneğe sahiptir. Bu, “eğitim üretimi” adını verdiğimiz bir araştırma alanıdır. Bu terminoloji, okul ortamının özelliklerini paraya mal olan girdiler olarak düşündüğümüz, okulların ürettiği çıktının ise öğrencilerin öğrenmesi olduğu gerçeğini yansıtır. Eğitim üretimi araştırmalarında temel soru, maliyetler göz önüne alındığında hangi girdilerin en fazla öğrenmeyi ürettiğidir. Sınıf mevcudu, en pahalı girdilerden biridir, çünkü daha küçük sınıflar ancak daha fazla öğretmen istihdam edilerek elde edilebilir. Bu nedenle, daha küçük sınıfların maliyetinin daha yüksek öğrenci başarısıyla uyumlu olup olmadığını bilmek önemlidir. STAR Deneyi, bu soruya cevap aramayı amaçlamaktadır.

Deneysel olmayan verileri kullanan birçok eğitim üretimi çalışması, sınıf mevcudu ile öğrenci öğrenmesi arasında çok az veya hiç bağlantı olmadığını göstermektedir. Bu nedenle, okul sistemlerinin, öğretmen maliyetinden tasarruf ederek başarıda herhangi bir azalma olmadan daha az öğretmen istihdam etmeleri mümkün olabilir. Ancak, sınıf mevcudu ile öğrenci başarısı arasındaki gözlemlenen ilişki, göründüğü gibi yorumlanmamalıdır, çünkü daha zayıf öğrenciler genellikle kasıtlı olarak daha küçük sınıflara atanır. Rastgele bir deney, elmaları elmalarla karşılaştırmamızı sağlayarak (farklı büyüklükte sınıflara atanmış öğrencileri karşılaştırarak) bu sorunun üstesinden gelir. Tennessee STAR Deneyi’nin sonuçları, daha küçük sınıflar için güçlü ve kalıcı bir getiri olduğunu göstermektedir (Finn ve Achilles, 1990 ve Krueger, 1999’da STAR verilerinin ekonometrik analizi için orijinal çalışmaya bakınız).

Bu deney aklaşık 12 milyon dolara mal oldu ve 1985-1986 döneminde bir anaokulu grubunda uygulandı. Çalışma, öğrencilerin üçüncü sınıfa gelene kadar dört yıl süren bir süreyi kapsadı ve yaklaşık 11.600 çocuğu içeriyordu. 1985-1986 döneminde Tennessee’deki normal sınıflardaki ortalama sınıf mevcudu yaklaşık 22.3 idi. Deney, öğrencileri üç farklı tedavi grubuna atadı: 13-17 çocuklu küçük sınıflar, 22-25 çocuklu normal sınıflar ve tam zamanlı veya yarı zamanlı öğretmen yardımcısı olan normal sınıflar. Deneye katılmayı seçen okullarda her grupta en az üç sınıf bulunuyordu.

Rastgele bir deney hakkında sorulan ilk soru, randomizasyonun deneklerin özelliklerini farklı tedavi grupları arasında başarılı bir şekilde dengeliyip dengelemediğidir. Bu değerlendirmeyi yapmak için, tedavi öncesi sonuçları veya diğer ortak değişkenleri gruplar arasında karşılaştırmak yaygındır. Ne yazık ki, STAR verileri tedavi öncesi test puanlarını içermese de, çocukların ırk ve yaş gibi özelliklerine bakmamıza izin vermektedir. Krueger’ın (1999) çalışmasından alınan Tablo 2.2.1, bu değişkenlerin ortalamalarını karşılaştırır. Tablodaki öğrenci özellikleri arasında ücretsiz öğle yemeği alıp almama, öğrenci ırkı ve öğrenci yaşının yer aldığını görebiliriz. Üç sınıf türü arasında bu özelliklerdeki farklılıklar küçüktür ve son sütundaki p değerlerine göre hiçbiri istatistiksel olarak önemli düzeyde farklı değildir. Bu, rastgele atamanın amaçlandığı gibi çalıştığını gösterir.

Tablo 2.2.1 Replikasyonu

Veri seti Kruger (1999) https://economics.mit.edu/people/faculty/josh-angrist/mhe-data-archive ’den indirilebir.

Stata veri setini R’a yüklemek için haven paketini indirmeyi unutmayın.

library(haven)
webstar <- read_dta("C:/Users/hutku/Downloads/webstar.dta")

datatable(webstar)

## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.html

srace değişkeninin açıklamaları

value label 1 white 2 black 3 asian 4 hispanic 5 am. indian 6 other 9 missing

srace değişkeni NA’leri çıkar.

webstar <- data.frame(webstar)

webstar <- webstar[!is.na(webstar$srace), ]

Gerekli değişiklikleri yapın.

library(dplyr)
webstar <- as.data.frame(webstar)
webstar <- webstar %>% mutate(white_asian = case_when(srace==1 | srace==3 ~ 1 , TRUE ~ 0))

datatable(webstar)

## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.html

cltypek değişkeninin açıklamaları.

value label 1 small class 2 regular class 3 regular + aide class

webstar <- webstar[!is.na(webstar$cltypek), ]

library(car)
res.ftest <- fligner.test(white_asian ~ cltypek, data = webstar)
res.ftest

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  white_asian by cltypek
## Fligner-Killeen:med chi-squared = 2.7615, df = 2, p-value = 0.2514

sesk değişkenin açıklamarı.

value label 1 free lunch 2 non-free lunch

webstar <- webstar[!is.na(webstar$sesk), ]

webstar <- webstar %>% mutate(free_lunch = case_when(sesk==1 ~ 1 , sesk==2 ~ 0))

webstar %>%
  group_by(cltypek) %>%
  summarise_at(c("white_asian", "free_lunch", "treadssk", "tmathssk"), mean, na.rm = TRUE)

## # A tibble: 3 × 5
##   cltypek                  white_asian free_lunch treadssk tmathssk
##   <dbl+lbl>                      <dbl>      <dbl>    <dbl>    <dbl>
## 1 1 [small class]                0.682      0.471     441.     491.
## 2 2 [regular class]              0.675      0.477     435.     483.
## 3 3 [regular + aide class]       0.660      0.503     435.     483.

res.ftest <- fligner.test(free_lunch ~ cltypek, data = webstar)
res.ftest

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  free_lunch by cltypek
## Fligner-Killeen:med chi-squared = 3.2153, df = 2, p-value = 0.2004

Star deneyinin yorumları bu bölümle ilgili değil. Bu deneyin açıklamaları Bölüm 2’de bulunabilir.

The Angrist and Lavy (1999) replikasyonu

Verisetlerini indirin

Şu internet sitesinden data4 ve data5 dta (stata) uzantılı veri setlerini indirebilirsiniz.

https://dataverse.harvard.edu/file.xhtml?persistentId=doi:10.7910/DVN/XRSUJU/0FCG1X&version=1.0

library(haven)
final4 <- read_dta("C:/Users/hutku/Downloads/final4.dta")

datatable(final4)

final5 <- read_dta("C:/Users/hutku/Downloads/final5.dta")

datatable(final5)

Angrist ve Lavy (1999) makalesinde, İbn Meymun’un kuralının İsrail devlet okullarında sınıf mevcudu belirlemek için kullanıldığından bahsediyorlar. Bu kurala göre, en fazla 40 öğrenci bir sınıfa kaydedilebilir. Makalede, bu kuralın varyasyonlarına dayanarak sınıf mevcudunun öğrenci başarısı üzerindeki etkisini tahmin etmek mümkün olduğunu belirtiyorlar.

İbn Meymun’un kuralına göre, sınıf mevcudu 40 öğrenciye kadar kayıtla bire bir artar. Ancak 41 öğrenci kaydolduğunda, ortalama sınıf mevcudu 20.5’e düşer. Benzer şekilde, 80 öğrenci kaydolduğunda ortalama sınıf mevcudu 40 olurken, 81 öğrenci kaydolduğunda ortalama sınıf mevcudu 27’ye düşer.

Bu varyasyonları kullanarak, araştırmacılar İsrailli öğrencilerin başarısı üzerinde sınıf mevcudunun etkisini tahmin etmek için bu kuralı kullanmışlardır. Bu, sınıf mevcudu değişkeninin doğal bir deneyim gibi işlev görmesini sağlar, çünkü kayıt sayısının belirli bir eşik değeri aşıldığında sınıf mevcudu beklenmedik şekilde değişir.

Angrist ve Lavy, bu tahminleri kullanarak sınıf mevcudunun öğrenci başarısı üzerindeki etkisini değerlendirmişler ve bu etkinin pozitif olduğunu bulmuşlardır. Daha küçük sınıflarda öğrenim gören öğrencilerin genellikle daha iyi performans gösterdiği sonucuna varmışlardır. Bu bulgular, daha küçük sınıfların öğrenci başarısını artırabileceği fikrini desteklemektedir.

Mostly Harmless’da bulunmayan ama makale içinde verilen tablolara da yer vermeye çalışacağız. Öncelikle özet tablosunu hazırlayın.

library(tidyr)

final4 %>% 
  select(classize, c_size, tip_a, verbsize, mathsize, avgverb, avgmath) %>%
  rename("Enrollment" = "c_size", "Percent disadvantaged" = "tip_a", "Class size" = "classize", "Reading size" = "verbsize", "Math size" = "mathsize", "Average verbal" = "avgverb", "Average math" = "avgmath") %>%
  mutate_at(c('Enrollment', 'Class size','Math size', 'Reading size', 'Average math','Average verbal', 'Percent disadvantaged'), as.numeric) %>%
  summarise(across(where(is.numeric), .fns = 
                     list(Mean = ~mean(.,na.rm=TRUE),
                          S.D. = ~sd(.,na.rm=TRUE),
                          Q10 = ~quantile(., 0.10, na.rm=TRUE),
                          Q25 = ~quantile(., 0.25, na.rm=TRUE),
                          Q50 = ~quantile(., 0.25, na.rm=TRUE),
                          Q75 = ~quantile(., 0.75, na.rm=TRUE),
                          Q90 = ~quantile(., 0.90, na.rm=TRUE)))) %>%
  pivot_longer(everything(), names_sep='_', names_to=c('variable', '.value')) %>%
  mutate(across(where(is.numeric), round, 2))

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(where(is.numeric), round, 2)`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
## 
##   # Previously
##   across(a:b, mean, na.rm = TRUE)
## 
##   # Now
##   across(a:b, \(x) mean(x, na.rm = TRUE))

## # A tibble: 7 × 8
##   variable               Mean  S.D.   Q10   Q25   Q50   Q75   Q90
##   <chr>                 <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Class size             30.3  6.41  22    26    26    35    38  
## 2 Enrollment             78.4 37.9   30    51    51   102   128  
## 3 Percent disadvantaged  13.9 13.4    2     4     4    19    35  
## 4 Reading size           27.6  6.55  19    24    24    32    36  
## 5 Math size              28.1  6.56  19    24    24    33    36  
## 6 Average verbal         72.5  7.99  62.2  67.7  67.7  78.2  82  
## 7 Average math           68.9  8.77  57.5  63.6  63.6  75.0  79.4

final5 %>% 
  select(classize, c_size, tip_a, verbsize, mathsize, avgverb, avgmath) %>%
  rename("Enrollment" = "c_size", "Percent disadvantaged" = "tip_a", "Class size" = "classize", "Reading size" = "verbsize", "Math size" = "mathsize", "Average verbal" = "avgverb", "Average math" = "avgmath") %>%
  mutate_at(c('Enrollment', 'Class size','Math size', 'Reading size', 'Average math','Average verbal', 'Percent disadvantaged'), as.numeric) %>%
  summarise(across(where(is.numeric), .fns = 
                     list(Mean = ~mean(.,na.rm=TRUE),
                          S.D. = ~sd(.,na.rm=TRUE),
                          Q10 = ~quantile(., 0.10, na.rm=TRUE),
                          Q25 = ~quantile(., 0.25, na.rm=TRUE),
                          Q50 = ~quantile(., 0.25, na.rm=TRUE),
                          Q75 = ~quantile(., 0.75, na.rm=TRUE),
                          Q90 = ~quantile(., 0.90, na.rm=TRUE)))) %>%
  pivot_longer(everything(), names_sep='_', names_to=c('variable', '.value')) %>%
  mutate(across(where(is.numeric), round, 2))

## # A tibble: 7 × 8
##   variable               Mean  S.D.   Q10   Q25   Q50   Q75   Q90
##   <chr>                 <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Class size             30.0  6.6   21    26    26    35    38  
## 2 Enrollment             77.9 39.1   31    50    50   100   128. 
## 3 Percent disadvantaged  14.1 13.6    2     4     4    20    35  
## 4 Reading size           27.3  6.62  19    23    23    32    36  
## 5 Math size              27.7  6.68  19    23    23    33    36  
## 6 Average verbal         74.4  8.08  64.2  69.9  69.9  79.8  83.3
## 7 Average math           67.3 10.0   54.8  61.1  61.1  74.1  79.4

Buradaki sonuçlar, makale de bulunan sonuçlardan biraz farklı. Yazarların notlarından anlaşılacağı üzere veriseti üzerinde bir subset yapılmış.

final4ek<- final4 %>% 
  filter(classize > 1, classize < 45, c_size>5)

final4ek %>% 
  select(classize, c_size, tip_a, verbsize, mathsize, avgverb, avgmath) %>%
  rename("Enrollment" = "c_size", "Percent disadvantaged" = "tip_a", "Class size" = "classize", "Reading size" = "verbsize", "Math size" = "mathsize", "Average verbal" = "avgverb", "Average math" = "avgmath") %>%
  mutate_at(c('Enrollment', 'Class size','Math size', 'Reading size', 'Average math','Average verbal', 'Percent disadvantaged'), as.numeric) %>%
  summarise(across(where(is.numeric), .fns = 
                     list(Mean = ~mean(.,na.rm=TRUE),
                          S.D. = ~sd(.,na.rm=TRUE),
                          Q10 = ~quantile(., 0.10, na.rm=TRUE),
                          Q25 = ~quantile(., 0.25, na.rm=TRUE),
                          Q50 = ~quantile(., 0.25, na.rm=TRUE),
                          Q75 = ~quantile(., 0.75, na.rm=TRUE),
                          Q90 = ~quantile(., 0.90, na.rm=TRUE)))) %>%
  pivot_longer(everything(), names_sep='_', names_to=c('variable', '.value')) %>%
  mutate(across(where(is.numeric), round, 2))

## # A tibble: 7 × 8
##   variable               Mean  S.D.   Q10   Q25   Q50   Q75   Q90
##   <chr>                 <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Class size             30.3  6.35  22    26    26    35    38  
## 2 Enrollment             78.3 37.9   30    51    51   102   128. 
## 3 Percent disadvantaged  13.9 13.4    2     4     4    19    35  
## 4 Reading size           27.7  6.54  19    24    24    32    36  
## 5 Math size              28.1  6.55  19    24    24    33    36  
## 6 Average verbal         72.5  7.99  62.2  67.7  67.7  78.2  82  
## 7 Average math           68.9  8.77  57.5  63.6  63.6  75.0  79.4

final5ek<- final5 %>% 
  filter(classize > 1, classize < 45, c_size>5)

final5ek %>% 
  select(classize, c_size, tip_a, verbsize, mathsize, avgverb, avgmath) %>%
  rename("Enrollment" = "c_size", "Percent disadvantaged" = "tip_a", "Class size" = "classize", "Reading size" = "verbsize", "Math size" = "mathsize", "Average verbal" = "avgverb", "Average math" = "avgmath") %>%
  mutate_at(c('Enrollment', 'Class size','Math size', 'Reading size', 'Average math','Average verbal', 'Percent disadvantaged'), as.numeric) %>%
  summarise(across(where(is.numeric), .fns = 
                     list(Mean = ~mean(.,na.rm=TRUE),
                          S.D. = ~sd(.,na.rm=TRUE),
                          Q10 = ~quantile(., 0.10, na.rm=TRUE),
                          Q25 = ~quantile(., 0.25, na.rm=TRUE),
                          Q50 = ~quantile(., 0.25, na.rm=TRUE),
                          Q75 = ~quantile(., 0.75, na.rm=TRUE),
                          Q90 = ~quantile(., 0.90, na.rm=TRUE)))) %>%
  pivot_longer(everything(), names_sep='_', names_to=c('variable', '.value')) %>%
  mutate(across(where(is.numeric), round, 2))

## # A tibble: 7 × 8
##   variable               Mean  S.D.   Q10   Q25   Q50   Q75   Q90
##   <chr>                 <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Class size             29.9  6.54  21    26    26    35    38  
## 2 Enrollment             77.9 39.1   31    50    50   100   129. 
## 3 Percent disadvantaged  14.1 13.6    2     4     4    20    35  
## 4 Reading size           27.3  6.58  19    23    23    32    36  
## 5 Math size              27.7  6.64  19    23    23    33    36  
## 6 Average verbal         74.4  8.08  64.2  69.8  69.8  79.8  83.3
## 7 Average math           67.3 10.0   54.8  61.1  61.1  74.1  79.4

Tabloda Discontinuity (Süreksizlik) örnekleminin sonuçları da verilmiş. Enrollment’ın 36-45 aralığı

final4ekdis<- final4ek %>% 
  filter(c_size>35 & c_size<46 | c_size>75 & c_size<86 | c_size>115 & c_size<125)

final4ekdis %>% 
  select(classize, c_size, tip_a, verbsize, mathsize, avgverb, avgmath) %>%
  rename("Enrollment" = "c_size", "Percent disadvantaged" = "tip_a", "Class size" = "classize", "Reading size" = "verbsize", "Math size" = "mathsize", "Average verbal" = "avgverb", "Average math" = "avgmath") %>%
  mutate_at(c('Enrollment', 'Class size','Math size', 'Reading size', 'Average math','Average verbal', 'Percent disadvantaged'), as.numeric) %>%
  summarise(across(where(is.numeric), .fns = 
                     list(Mean = ~mean(.,na.rm=TRUE),
                          S.D. = ~sd(.,na.rm=TRUE),
                          Q10 = ~quantile(., 0.10, na.rm=TRUE),
                          Q25 = ~quantile(., 0.25, na.rm=TRUE),
                          Q50 = ~quantile(., 0.25, na.rm=TRUE),
                          Q75 = ~quantile(., 0.75, na.rm=TRUE),
                          Q90 = ~quantile(., 0.90, na.rm=TRUE)))) %>%
  pivot_longer(everything(), names_sep='_', names_to=c('variable', '.value')) %>%
  mutate(across(where(is.numeric), round, 2))

## # A tibble: 7 × 8
##   variable               Mean  S.D.   Q10   Q25   Q50   Q75   Q90
##   <chr>                 <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Class size             31.0  7.24  21    25    25    38    40  
## 2 Enrollment             77.9 29.6   41    43    43   116   119  
## 3 Percent disadvantaged  13.0 12.4    1     4     4    18    32  
## 4 Reading size           28.2  7.71  18    22    22    35    37.2
## 5 Math size              28.6  7.69  18    23    23    35    38  
## 6 Average verbal         72.4  7.88  61.9  67.0  67.0  78.4  81.7
## 7 Average math           68.7  9.14  56.8  62.6  62.6  75.4  79.7

final5ekdis<- final5ek %>% 
  filter(c_size>35 & c_size<46 | c_size>75 & c_size<86 | c_size>115 & c_size<125)

final5ekdis %>% 
  select(classize, c_size, tip_a, verbsize, mathsize, avgverb, avgmath) %>%
  rename("Enrollment" = "c_size", "Percent disadvantaged" = "tip_a", "Class size" = "classize", "Reading size" = "verbsize", "Math size" = "mathsize", "Average verbal" = "avgverb", "Average math" = "avgmath") %>%
  mutate_at(c('Enrollment', 'Class size','Math size', 'Reading size', 'Average math','Average verbal', 'Percent disadvantaged'), as.numeric) %>%
  summarise(across(where(is.numeric), .fns = 
                     list(Mean = ~mean(.,na.rm=TRUE),
                          S.D. = ~sd(.,na.rm=TRUE),
                          Q10 = ~quantile(., 0.10, na.rm=TRUE),
                          Q25 = ~quantile(., 0.25, na.rm=TRUE),
                          Q50 = ~quantile(., 0.25, na.rm=TRUE),
                          Q75 = ~quantile(., 0.75, na.rm=TRUE),
                          Q90 = ~quantile(., 0.90, na.rm=TRUE)))) %>%
  pivot_longer(everything(), names_sep='_', names_to=c('variable', '.value')) %>%
  mutate(across(where(is.numeric), round, 2))

## # A tibble: 7 × 8
##   variable               Mean  S.D.   Q10   Q25   Q50   Q75   Q90
##   <chr>                 <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Class size             30.7  7.43  21    24    24    38    40  
## 2 Enrollment             75.8 29.2   40.4  43    43    85   119  
## 3 Percent disadvantaged  13.7 13.2    2     4     4    17    36  
## 4 Reading size           28.1  7.35  18    21    21    35    38  
## 5 Math size              28.5  7.46  18    22    22    35    38  
## 6 Average verbal         74.5  8.22  63.8  69.6  69.6  80.5  83.6
## 7 Average math           67.0 10.3   54.5  60.8  60.8  74.1  80.0

Bu subset daha sonra, süreksizlik regresyon tablolarında da kullanılacaktır.

Meymun_rule <- function(x) {x / (floor((x - 1)/40) + 1)}

final4Enrollmentmeans <- final4ek %>%
  group_by(c_size) %>%
  summarise(mean_class_size = mean(classize, na.rm = TRUE))

ggplot(data = final4Enrollmentmeans, aes(x = c_size)) + geom_line(aes(y = mean_class_size)) +
           stat_function(fun      = Meymun_rule,
                         linetype = "dashed")           +
           expand_limits(y = 0)                         +
           scale_x_continuous(breaks = seq(0, 240, 40)) +
           ylab("Class size")                           +
           xlab("Enrollment count")                     +
           ggtitle("B. Fourth grade")

final5Enrollmentmeans <- final5ek %>%
  group_by(c_size) %>%
  summarise(mean_class_size = mean(classize, na.rm = TRUE))

ggplot(data = final5Enrollmentmeans, aes(x = c_size)) + geom_line(aes(y = mean_class_size)) +
           stat_function(fun      = Meymun_rule,
                         linetype = "dashed")           +
           expand_limits(y = 0)                         +
           scale_x_continuous(breaks = seq(0, 240, 40)) +
           ylab("Class size")                           +
           xlab("Enrollment count")                     +
           ggtitle("A. Fifth grade")

library(sandwich)
library(lmtest)
library(AER)

library(estimatr)

Tüm Veri Regresyon sonuçları

Makele içinde veri setinde yapılan bazı değişikliklerle başlayalın.

final5ek = final5ek %>% 
       mutate(avgmath = case_when(avgmath >100 ~ avgmath - 100,
                             TRUE ~ avgmath))

final5ek <- filter(final5ek,c_leom==1 & c_pik<3)

final5ek <- filter(final5ek,classize>1 & classize<45 & c_size>5)

final5ek <- filter(final5ek,mathsize>0)

ols1 <- lm_robust(avgmath ~ classize, data = final5ek, clusters = final5ek$schlcode,se_type = "stata")

ols2 <- lm_robust(avgmath ~ classize + tipuach , data = final5ek, clusters = final5ek$schlcode,se_type = "stata")

ols3 <- lm_robust(avgmath ~ classize + tipuach + c_size , data = final5ek, clusters = final5ek$schlcode, se_type = "stata")

Süreksizlik Örneklemi

olsdis1 <- lm_robust(avgmath ~ classize, data = final5ekdis, clusters = final5ekdis$schlcode,se_type = "stata")

olsdis2 <- lm_robust(avgmath ~ classize + tipuach , data = final5ekdis, clusters = final5ekdis$schlcode,se_type = "stata")

olsdis3 <- lm_robust(avgmath ~ classize + tipuach + c_size , data = final5ekdis, clusters = final5ekdis$schlcode, se_type = "stata")

library(huxtable)

huxreg(ols1, ols2, ols3,
       error_format = "[{statistic}]", 
       note         = "{stars}. T statistics in brackets.", number_format = 2)

huxreg(olsdis1, olsdis2, olsdis3,
       error_format = "[{statistic}]", 
       note         = "{stars}. T statistics in brackets.", number_format = 2)

Araç Değişken Regresyonları

Araç değişkeni oluşturun.

final5ek$f <- final5ek$c_size / (floor((final5ek$c_size - 1) / 40) + 1)

final5ekdis$f <- final5ekdis$c_size / (floor((final5ekdis$c_size - 1) / 40) + 1)

iv1 <- iv_robust(avgmath ~ classize  | f , data = final5ek, clusters = final5ek$schlcode,se_type = "stata")

iv2 <- iv_robust(avgmath ~ classize + tipuach  | f + tipuach , data = final5ek, clusters = final5ek$schlcode,se_type = "stata")

iv3 <- iv_robust(avgmath ~ classize + tipuach + c_size  | f + tipuach + c_size , data = final5ek, clusters = final5ek$schlcode,se_type = "stata")

huxreg(iv1, iv2, iv3,
       error_format = "[{statistic}]", 
       note         = "{stars}. T statistics in brackets.")

Süreksiz örneklem araç regresyonları

ivdis1 <- iv_robust(avgmath ~ classize  | f , data = final5ekdis, clusters = final5ekdis$schlcode,se_type = "stata")

ivdis2 <- iv_robust(avgmath ~ classize + tipuach  | f + tipuach , data = final5ekdis, clusters = final5ekdis$schlcode,se_type = "stata")

ivdis3 <- iv_robust(avgmath ~ classize + tipuach + c_size  | f + tipuach + c_size , data = final5ekdis, clusters = final5ekdis$schlcode,se_type = "stata")

huxreg(ivdis1, ivdis2, ivdis3,
       error_format = "[{statistic}]", 
       note         = "{stars}. T statistics in brackets.")

Son olarak 3. modellerin hepsini bir tabloda gösterip karşılaştıralım.

huxreg(ols3, olsdis3, iv3, ivdis3,
       error_format = "[{statistic}]", 
       note         = "{stars}. T statistics in brackets.")

Angrist ve Lavy’nin bahsettiği “Mostly Harmless” kitabındaki 6.2.1 tablosunda, beşinci sınıf matematik puanları için tahminler kontrolsüz ve kontrollü durumda karşılaştırılmıştır. Kontrolsüz tahminlerde, sınıf mevcudu ile test puanları arasında güçlü bir pozitif ilişki olduğu görülmektedir.

Ancak, kontrol olarak dezavantajlı öğrencilerin yüzdesi dahil edildiğinde, pozitif ilişki anlamlılığını kaybetmektedir. Yani, daha büyük sınıflardaki öğrencilerin yüksek puanlarının, dezavantajlı öğrenci oranındaki farklılıklardan kaynaklanabileceği sonucuna varılmaktadır.

Ancak, 2SLS (iki aşamalı en küçük kareler) tahminleri, OLS (en küçük kareler) tahminlerinden farklı bir sonuç vermektedir. 2SLS tahminleri, daha küçük sınıfların test puanlarını artırdığını göstermektedir. Bu sonuç, Tennessee STAR randomize çalışmasının sonuçlarıyla uyumlu olduğunu ifade etmektedir.

Sonuç olarak, bahsedilen çalışmalardan elde edilen bulgular, sınıf mevcudunun öğrenci başarısı üzerindeki etkisinin karmaşık olduğunu göstermektedir. Daha küçük sınıfların avantajları ve dezavantajları dikkate alınmalı ve bu konuda daha fazla araştırma yapılması gerekmektedir.