R Dersi 14

R Markdown

İstatistik kısmını sonlandırdık. Şu ana kadar veriyi bulma veya oluşturma, temizleme ve dönüştürme, analiz etme, görseleştirme ve yorumlama aşamalarının hepsini öğrenmiş olduk. Geriye sadece raporlama kısmı kaldı. R Markdown da bu noktada devreye giriyor. R Markdown; yazı, resim, tablo, kod gibi farklı elementleri tek bir dökümanda birleştirmemize yarar. Böylece web sitesine koymak için blog, paylaşmak için yazı, sunum yazababilirsiniz. Ayrıca kullanarak da formüller, tablolar koyup istediğimiz formatta yapabiliriz. Mesela makaleler falan hep bu şekilde yazılır.

R Markdown dökümanı oluşturduğunuz zaman tepede YAML Header dediğimiz bir kısım çıkar. Yazar ismi, çıktı formatı ve gün gibi bilgileri burada yazarız. Kullanılacak template gibi şeyler de bu kısımda belirtilir.

Eğer düz yazı girmek istiyorsak normal şekilde yazabiliriz. R scriptin aksine comment gibi # koymaya gerek yoktur. Düz yazı yazar gibi girebiliriz. Yazıyı farklı formatlarda yazabiliriz. Mesela heading için farklı seviyeler var.

# Heading 1

## Heading 2

### Heading 3

Listelemek için satır başlarına - veya *

Liste 1
Liste 2
Liste 3

Texti bold yapmak için yazıyı ** ** arasına almak lazım

Italic yapmak içinse * * arasına alırız

Link koymak içinse Link ismi burada yapılır ancak boşluk bırakılmadan.

Ancak burada kod girmek için farklı bir sistem kullanmamız gerekir. Kod girmek için ise {r} arasına kodumuzu yazarız. Bunu oluşturmanın basit bir yolu Command Option I tuşlarına basmak. Böylece doğrudan bir code chunk oluşturabiliriz.

# Kod buraya yazılacak

Kod oluştururken curly bracket içinde r diye belirtmiştik. Bunun yanına farklı optionlar ekleyebiliriz. Bunlar

include: Eğer bunu FALSE seçersek ne kodu ne de outputu dökümana yansır. Ancak kod yine de çalıştırılır.
eval: Bunu FALSE seçersek yazdığımız kod otomatik olarak çalıştırılmaz.
echo: Sadece kodu göstermemek ama kodun çıktısını göstermek istiyorsak bunu FALSE yaparız.
message: Eğer kod bir mesaj gösteriyorsa saklamak için FALSE seçeriz.
warning: Eğer kod bir mesaj warning saklamak için FALSE seçeriz.

Döküman içine resim koymak istiyorsak yine paranteler arasında boşluk bırakmadan.

![Yazı buraya] (image_url veya image_path)

Tablo oluşturmak için bu şekilde bir yöntem kullanırız ancak bu bizim için gerekli değil. Genellikle tabloyu paketler ile oluşturup direkt alacağız.

Header 1	Header 2
Cell 1	Cell 2

t-test tablosu ile

data(mtcars)
library(broom)
library(knitr)
library(kableExtra)
t_test_result <- t.test(mpg ~ am, data = mtcars)
tidy_result <- tidy(t_test_result)

kable(tidy_result, format = "html", 
      caption = "T-Test Results") %>%
  kable_styling(bootstrap_options = 
                  c("striped", "hover", "condensed", "responsive"))

T-Test Results
estimate	estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high	method	alternative
-7.244939	17.14737	24.39231	-3.767123	0.0013736	18.33225	-11.28019	-3.209684	Welch Two Sample t-test	two.sided

Assumption testler için de bir örneğe bakalım

shapiro_test <- shapiro.test(mtcars$mpg)
tidy_shapiro <- tidy(shapiro_test)
kable(tidy_shapiro, format = "markdown", 
      caption = "Shapiro-Wilk Normality Test Results")

Shapiro-Wilk Normality Test Results
statistic	p.value	method
0.9475647	0.1228814	Shapiro-Wilk normality test

Daha geniş customize seçenekleri için de farklı paketler kullanabiliriz.

library(kableExtra)
kable(tidy_result, format = "markdown", 
      caption = "Enhanced T-Test Results") %>%
  kable_styling(full_width = FALSE)

Enhanced T-Test Results
estimate	estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high	method	alternative
-7.244939	17.14737	24.39231	-3.767123	0.0013736	18.33225	-11.28019	-3.209684	Welch Two Sample t-test	two.sided

Böylece tabloyu da oluşturmuş oluyoruz. Ancak regresyon tabloları için farklı yöntemler de var. Bunlar daha güçlü.

library(stargazer)
model1 <- lm(mpg ~ wt + hp, data = mtcars)
model2 <- lm(mpg ~ wt + qsec, data = mtcars)

stargazer(model1, model2, 
          type = "latex", 
          title = "Regression Results", 
          float = FALSE, 
          omit.table.layout = "n",
          notes = "")

% Table created by stargazer v.5.2.3 by Marek Hlavac, Social Policy Institute. E-mail: marek.hlavac at gmail.com % Date and time: Fri, Oct 04, 2024 - 11:47:55

library(modelsummary)
modelsummary(list("Model 1" = model1, 
                  "Model 2" = model2), 
             output = "kableExtra", fmt = "%.3f")

	Model 1	Model 2
(Intercept)	37.227	19.746
	(1.599)	(5.252)
wt	−3.878	−5.048
	(0.633)	(0.484)
hp	−0.032
	(0.009)
qsec		0.929
		(0.265)
Num.Obs.	32	32
R2	0.827	0.826
R2 Adj.	0.815	0.814
AIC	156.7	156.7
BIC	162.5	162.6
Log.Lik.	−74.326	−74.360
F	69.211	69.033
RMSE	2.47	2.47

model <- lm(mpg ~ wt + hp, data = mtcars)
tidy(model)

## # A tibble: 3 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  37.2      1.60        23.3  2.57e-20
## 2 wt           -3.88     0.633       -6.13 1.12e- 6
## 3 hp           -0.0318   0.00903     -3.52 1.45e- 3

Report paketi farklı bir alternatif olarak karşımıza çıkıyor

library(report)
report(t.test(mtcars$mpg ~ mtcars$am))

## Effect sizes were labelled following Cohen's (1988) recommendations.
## 
## The Welch Two Sample t-test testing the difference of mtcars$mpg by mtcars$am
## (mean in group 0 = 17.15, mean in group 1 = 24.39) suggests that the effect is
## negative, statistically significant, and large (difference = -7.24, 95% CI
## [-11.28, -3.21], t(18.33) = -3.77, p = 0.001; Cohen's d = -1.41, 95% CI [-2.26,
## -0.53])

Döküman içinde tablolara, resimlere, figürlere ve sectionlara da referans vermemiz mümkündür. Bunları bookdown adını verdiğimiz bir sistemi ve paketi kullanarak yapıyoruz. Burada bahsedeceğim kısım çok kısıtlı olacak. Eğer merak ediyorsanız daha detaylı bakabilirsiniz.

Resimlere referans vermek için:

Logistic Regression Plot Daha sonra bu referansı kullanmak için:

Refer to Figure @ref(fig:logistic) to see the logistic regression plot.

Figürlere referans vermek için:

plot(mtcars$wt, mtcars$mpg)

Scatter Plot of MPG vs Weight

Daha sonra bu referansı kullanmak için:

See Figure @ref(fig:fig1) for the scatter plot.

Referans

Sectionlara referans verme

Daha sonra bu referansı kullanmak için:

Refer to Section @ref(sec:introduction) for more details.

Tablolara referans verme

t_test_result <- t.test(mpg ~ am, data = mtcars)
tidy_result <- broom::tidy(t_test_result)
knitr::kable(tidy_result, format = "markdown", 
             caption = "T-Test Results")

T-Test Results
estimate	estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high	method	alternative
-7.244939	17.14737	24.39231	-3.767123	0.0013736	18.33225	-11.28019	-3.209684	Welch Two Sample t-test	two.sided

Daha sonra bu referansı kullanmak için:

See Table @ref(tab:summary) for the summary statistics.

Daha önceden özellikle formüller için Latex kullanabileceğimizi söylemiştim. Formülleri line içinde göstermek için tek dolar isareti içine alırız. Mesela $6^3*8$ gibi yapabiliriz. Ancak eğer formülü ayrı bir kısım olarak göstermek istiyorsak cift dolar isareti içine alırız.

\[6^3*8\]

Yazı yazarken genelde modellerimizin formülün koyarız. Böylece okuyanların bütün metod kısmını okumadan formüle bakıp ne yaptığımızı anlamasını sağlarız. Regresyonun formülünü kolay bir şekilde yazmak için paket var.

Equatiomatic ile de formülleri çıkarabiliyoruz. Bunları Latex ile çıkardığı için doğrudan kullanabiliyoruz.

library(equatiomatic)
linear_model <- lm(mpg ~ wt + hp, data = mtcars)
logistic_model <- glm(am ~ wt + hp, 
                      data = mtcars, family = binomial)

extract_eq(linear_model)

\[ \operatorname{mpg} = \alpha + \beta_{1}(\operatorname{wt}) + \beta_{2}(\operatorname{hp}) + \epsilon \]

extract_eq(logistic_model)

\[ \log\left[ \frac { P( \operatorname{am} = \operatorname{1} ) }{ 1 - P( \operatorname{am} = \operatorname{1} ) } \right] = \alpha + \beta_{1}(\operatorname{wt}) + \beta_{2}(\operatorname{hp}) \]

Ayrı bir code chunk girmeden, kodumuzu doğrudan text içinde de çalıştırmamız mümkün. Bu mantığı Latex ile benzer düşünebilirsiniz.

The mean of the mpg column in the mtcars dataset is 20.090625.

Bunu daha advanced şekilde yapıyorsak

The t-statistic for comparing MPG between automatic and manual cars is -3.77.

Yazı referansı vermek içinse Latex’ten yararlanırız. Bunun için bib dediğimiz bir döküman gerekir. Bunu farklı programla oluşturabilirsiniz. Zotera falan buna yarar. Referansın latex kodunu da Google Scholar üstünden falan rahatça bulabilirsiniz.

Referans vermek için

Markdownda bu şekilde referans verebilirsiniz (@johnson_introduction_2018).

Causal Inference

Regresyon konusunda correlation ve causation farkına değinmiştik. Mesela dondurma satışlarının artması ile boğulma vakaları arasında bir pozitif ilişkli olduğunu varsayalım. Burada akılımıza ilk gelen sonuç dondurma tüketimi insanların fiziksel durumunu olumsuz etkiliyor ve boğulmalara sebep oluyor gibi bir çıkarım mı olur yoksa başka bir bağlantı olduğunu mu düşünürüz?

Z değişkeni veya confounder burada devreye giriyor. Z değişkeni X ve Y değişkenleri arasındaki ilişkiyi düzenleyen değişkene denir. Mesela bu örnekte o değişken yaz. Mantıklı bir yaklaşım yaparak X ve Y arasında sebep-sonuç ilişkisinin olmadığını ama başka bir değişkenin bu ikisini de artırabileceğini söylüyoruz.

Regresyon bu korelasyon meselesinin önüne geçiyor ama bu sonuçtan illa nedensellik çıkaramıyoruz. Nedenselliğe daha yakın bir yorum yapmamızı sağlıyor ancak yine de bağımsız ve bağımlı değişkenler arasındaki ilişkyi ölçüyor. Tam olarak bir nedensellik estabilsh etmeyebiliyor.

Burada farkında olmamız gereken farklı noktalar var. Burada 3 tanesinden bahseceğim ancak elbette bunlarla sınırlı değil. Bir örnek üstünden gidelim. Aklımızda bir soru olsun. Partilerin kampanya harcamaları ve seçime katılım oranları arasındaki ilişkiyi ölçmek isteyelim. Bunun için 60 ülkeden son 5 seçimi aldığımız bir veri kullanalım. Regresyon ile ölçtüğümüz zaman sonuç istatistiksel olarak anlamlı ve pozitif çıkabilir. Burada yorum olarak şunu yapmamız lazım: seçim harcamalarının yüksek olduğu seçimlerde seçime katılım artmış. Seçim harcaması seçime katılımı etkilemez diyen null hipotezi reddediyoruz ve bu yorumu yapıyoruz.

Şimdi de bir partinin kampanya başı olduğunuzu düşünün. Bu yazıyı okuyorsunuz. Bu durumda kesin olarak harcamaları artırırsak seçime katılım artar diyebilir miyiz? Yani parti politikası olarak bunu yapmayı önerir misiniz?

Cevap evet olabilir ama causal inference burada devreye giriyor ve bize bir dur bekle diyor. 3 düşünme mantığından bahsetmiştim. Burada inceleyelim:

İlk olarak sahip olmamız gereken düşünce tipi coutnerfactual thinking. Eğer X olmasaydı ne olurdu sorusunu sormamız lazım. Mesela partiler seçim harcamasını hiç yapmasaydı katılım ne olurdu? Bu durumu establish edebiliyor muyuz?
İkinci olaraksa burada tek faktör seçim harcaması mı? Yoksa farkında olmadığımız bir Z değişkeni mi var? Mesela iki adayın oylarının yakın olduğu bir seçimde partiler öne çıkmak için daha çok harcama yapar ve seçim yakın olduğu için de seçmenler daha kritik bir seçim diye katılımı artırabilir. Ancak bu seçim kampanyasına harcanan paradan bağımsız olarak gerçekleşir.
Son düşünce bu durum için geçerli değil ama cause and effect ilişkisi kurmak istiyorsak effect causedan sonra gelmeli. Yani ilişkinin yönünü ve zamanlamasını doğru bilmeliyiz.

Bu 3 düşünce geçerli olmadan nedensellik ilişkisi kurmak hatalı olacaktır. İki değişken arasındaki ilişkiyi regresyonla falan ölçeriz ama nedensellik kurmak için bundan daha fazlası gerekir. Dersin bu kısmında nedensellik kurmamıza yardımcı olabilen bazı yöntemlerden bahsedeceğim. Bu yöntemler

Randomized Controlled Trials
Natural Experiments
Matching
Difference in Differences
Instumental Variables
Regression Discontinuity Design
Panel Data
Synthetic Control

Bu konulara geçmeden önce nedensellik ilişkilerinin görselleştirilmesi ile ilgili bir konuya değinmek lazım. Directed Acyclic Graphs veya DAG nedenselliğin yönünün ve ne yapıda olduğunun görsel olarak gösterilmesidir. Daha basitçe ve kolay anlaşılan şekilde nedensellik ilişkisini görebiliriz. Bu tarz grafiklerde noktalara node, çizgilere de edge denir. Network analizinde de aynı terimleri kullanırız.

Bu grafikte mediator değişken dediğimiz kavramın bir etkisini görüyoruz. X değişkeninin Y değişkeni üzerine etkisi var ancak bu sadece Z değişkeni aracılığıyla oluyor. Mesela eğitim ve oy verme davranışı arasında bu tarz bir ilişki olabilir. Eğitim politik bilgiyi artırıyor ve bu da oy verme davranışını etkiliyor.

library(dagitty)
dag_mediator <- dagitty("dag {
  X -> Z -> Y
}")
plot(dag_mediator)

## Plot coordinates for graph not supplied! Generating coordinates, see ?coordinates for how to set your own.

Burada ise confounder örneğini görüyoruz. U değişkeni hem X hem de Y değişkenini etkiliyor.

dag_confounder <- dagitty("dag {
  X -> Y
  U -> X
  U -> Y
}")
plot(dag_confounder)

## Plot coordinates for graph not supplied! Generating coordinates, see ?coordinates for how to set your own.

Collider değişken ise en sorunlu görebileceğimiz örnek. Bir Z değişkeni hem X hem de W tarafından etklieniyor ancak Y değişkenini etkilemiyor. Ancak biz Z değişkenini kontrol olarak analizimize eklersek W ve Y arasında bir yol açmış oluruz ve aslında işe yaramayan bir değişken etkili gibi görünebilir.

dag_collider <- dagitty("dag {
  X -> Y
  X -> Z
  W -> Z
}")
plot(dag_collider)

## Plot coordinates for graph not supplied! Generating coordinates, see ?coordinates for how to set your own.

Randomized Control Trials

Bunları tamamladıktan sonra ilk konumuz olan RCT’lere odaklanalım. RCT deneysel bir araştırma yöntemidir ve bir müdalalenin etkisini grupları treatment ve control olmak üzere ikiye bölüp inceler. Klasik deneysel yöntem budur. Bu noktada çok önemli bir nokta bu ayrımı rassal olarak yapmamız gereklilğidir. Eğer belli bir değişken üstünden ayrım yaparsak sadece intervention değişkenini değil diğer değişkeni de ölçmüş oluruz.

Bu yöntem gold standard diyebileciğimiz bir araştırma metodu. Hem causal ilişki kurulmasını sağlıyor, hem selection biası minimuma indiriyor hem de iki grup arasındaki ayrımı net bir şekilde görebiliyoruz. Zaten bütün bilim dallarında kullanılmasının sebebi de etkili olması.

Mesela ABD’de çok popüler olan bir yöntem. Belli bir bölgede yaşayan kişiler rastgele seçilir. Bazılarına “Bu seçimde oy kullanın” yazan mesajlar gönderilir. Diğer kişilere ise gönderilmez. Böylece rassal olarak treatment ve control grupları oluşturulmuş olur. Seçim sonrasındaysa kişilerin oy verip vermediğine dair bilgiler public olarak paylaşıldığı için bu kişilerin oy kullanıp kullanmadığına bakılır. Böylece deneyimiin ne kadar etkili olduğunu anlarız.

Burada counterfactual thinking dediğimiz yaklaşımın etkisini görüyoruz. Ya bunu vermeseydik ne olurdu sorusunun cevabı control group. Randomization ile de internval validity sağlamış oluyoruz. Böylece outcome değişkenini etkileyecek etmenlerin sayısını minimuma indiriyoruz. Hem observe ettiğimiz hem de etmediğimiz değişkenlerin etkisi azalıyor. ANCAK external validity yani genelleştirme konusunda sıkıntı yaşayabiliriz. Mesela bu yukarıdaki deney ABD’de işe yararken başka ülkelerde kültürel farklardan dolayı ters tepebilir.

set.seed(353)
n <- 1000
rct_data <- data.frame(
  subject_id = 1:n,
  treatment = sample(c(0, 1), n, replace = TRUE),
  outcome = rnorm(n, mean = 50, sd = 10)
)
rct_data$outcome[rct_data$treatment == 1] <- 
  rct_data$outcome[rct_data$treatment == 1] + 5
t.test(outcome ~ treatment, data = rct_data)

## 
##  Welch Two Sample t-test
## 
## data:  outcome by treatment
## t = -8.6646, df = 977.96, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -6.679756 -4.212771
## sample estimates:
## mean in group 0 mean in group 1 
##        49.77776        55.22402

Güzel bir yöntem olmasına rağmen eksikleri de var. Hem pahalı hem de pratik olmayabiliyorlar. Birçok konu için yapması zor oluyor. Aynı zamanda etik sorunlar da var. Kişilere seçim öncesi yalan bilgi vermek etik açıdan sorunlu. Kişilerin de onayı gerekebileceği analizler olduğu için belli durumlar haricinde çok tercih edemiyoruz.

Natural Experiment

Natural Experiment ise kişiler veya subjectler deneyi tasarlayan kişinin müdahalesi olmadan treatment ve control gruplarına ayrıldığı zaman gerçekleşir. Bu doğal olarak da (coğrafi) lab müdahalesinden bağımsız ama insan yapımı (devletlerin bölünmesi) olarak da gerçekleşebilir. Deneyi randomize edebileceğimiz bir ortam yoksa veya geçmişte olan bu tarz bir bölünmeyi inceliyorsak bu yönteme başvurabiliriz. Kısaca randomness real-world tarafından oluyor, bizim tarafımızdan değil.

Doğu Almanya broadcast öğrneği.

Acemoğlu

Burada beklentimiz yine bu iki grubu bölen şeyin araştırdığımız konudan bağımsız olması. Yani incelediğimiz treatment ve outcome değişkenlerinin dışında bir dağılım yapılmış olmalı. Bu dağılımın ise as-if random olduğunu varsayıyoruz.

Dağılım yapılması elimizde olmadığı için sanki random gibi dağılım yapılmış olduğunu varsayıyoruz. Yine RCT’ler gibi treatment ve control grubunun ne olduğu ve bunlara nasıl etki edeceği belli olmalı.

Comparability de iki grubun arasında oluşan farkın treatment dışında bir sebepten kaynaklanmaması gerekliliği ve eğer bu treatment olmasaydı iki grubun birbirine benzer olacağı assumptionı.

Natural experimentların en büyük eksiklerinden birisi bulmanın zorluğudur. Grupların birbirinden ayrılma sebebi ile araştırdığımız outcome’ın etkisinin bağımsız olduğu konuları bulmak zor olabiliyor. Aynı zamanda iki grubun bu ayrım dışında birbirine benzer olduğu konusunda kesinlik sağlaması zor olabiliyor.

Bunların genelleştirmesi de zor olabiliyor. Başka contextlerde dağılımın nasıl olacağını bilemediğimiz için sadece olay üstünde yorum yapabiliyoruz.

Bazense elimizde direkt araştırma verisi olabilir. Mesela bir anket veya bölgeler, insanlar hakkında bilgi veren başka bir veri. Bu gibi durumlarda quasi-experiment dediğimiz yöntemleri uygularız. Bu yöntemlerde ne elimizde önceden bildiğimiz bir ayrım ne de kendi elimizde kişileri böldüğümüz bir ayrım vardır. Veri üstünden ayrımı kendimiz yapmaya çalışırız ancak bu kişilerin sonuçlarının ne olduğu bellidir. Bu quasi-exprimental yöntemlerde öğreneceğimiz ilk konu da matching.

Matching

Matching observational çalışmalarda randomized experimentmış gibi bir analiz yapmamızı sağlar. Bunu da deneysel çalışmaların özelliklerini mimick ederek yapar. Elde ettiğimiz veriden özellikleri yakın olan kişileri bularak bunları treatment ve control gruplarına atamak. Bu birbirine benzer kişiler bizim treatment değişkenizimde ayrışacakalar ancak bunun dışında birbirlerine benzer olacak. Aynı RCT’lerde olduğu gibi. Burada önemli bir assumption cofounder olabilecek değişkenlerin unobserve edilmiş olmaması. Eğer önemli değişkenleri kaçırırsak random assignmentın dışına çıkabiliriz.

Mesela bir ülkede zorunlu oy kullanma kanunu olmasının o ülkedeki kişilerin oy kullanma ihtimalini artırıp artırmadığına bakalım. Bu durumd bizim treatment değerimiz oy kullanmanın zorunlu olduğu bir ülkede yaşamak. Ancak biz doğrudan bu ülkeleri ikiye bölüp karşılaştırma yapamayız çünkü ülkelerin kendi karakteristikleri bu konuda etkili olabilir. Ondan dolayı öncelikle insanları eğitim, gelir, yaş gibi değerlerini alıyoruz ve matching sistemi bunlarda en yakın olan kişileri de iki farklı grupta alıyor ve treatment ile control gruplarına atıyor. Böylece hem diğer anlamda birbirine yakın kişileri seçmiş oluyoruz hem de treatment grubumuz belli bir ayrım oluyor.

Bunu yaparken kullandığı yöntem genel olarak Propensity Score olur. Burada unitlerin treatment grubunda olma ihtimallerini belli sosyo-ekonomik değişkenlere göre inceleyeceğiz. Böylece birbirlerine benzer skorlara sahip kişileri karşılaştırmış olacağız.

set.seed(123)
data <- data.frame(
  compulsory_voting = rbinom(1000, 1, 0.5),
  turnout = rbinom(1000, 1, 0.7),            
  income = runif(1000, 20000, 100000),      
  education = sample(c("High School", "College", 
                       "Postgraduate"), 1000, 
                     replace = TRUE),
  age = runif(1000, 18, 80),                
  urban = rbinom(1000, 1, 0.6)              
)

Burada kişinin sosyolojik değerlerine bağlı olarak compulsory votingde olma ihtimalini hesapladık. Bundan sonra benzer skorlara sahip kişileri colpulsory voting olan ve olmayan ülkelerde eşleştirecez. Böylece bu kişilerin hepsi eğitimli zengin olacak ancak bazıları compulsory voting olan ülkede yaşarken diğerleri öbür ülkeden olacak.

ps_model <- glm(compulsory_voting ~ 
                  income + education + age + urban, 
                data = data, family = binomial)
data$propensity_score <- predict(ps_model, 
                                 type = "response")

library(MatchIt)
matchit_model <- matchit(compulsory_voting ~ 
                           income + education 
                         + age + urban, 
                         data = data, method = 
                           "nearest", distance 
                         = "logit")
summary(matchit_model)

## 
## Call:
## matchit(formula = compulsory_voting ~ income + education + age + 
##     urban, data = data, method = "nearest", distance = "logit")
## 
## Summary of Balance for All Data:
##                       Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance                     0.4977        0.4884          0.1902     1.0610
## income                   61281.0371    58534.0020          0.1163     1.0483
## educationCollege             0.3245        0.3846         -0.1283          .
## educationHigh School         0.3367        0.2781          0.1240          .
## educationPostgraduate        0.3387        0.3373          0.0031          .
## age                         48.8371       48.4923          0.0193     1.0493
## urban                        0.6227        0.6095          0.0273          .
##                       eCDF Mean eCDF Max
## distance                 0.0514   0.0791
## income                   0.0359   0.0713
## educationCollege         0.0601   0.0601
## educationHigh School     0.0586   0.0586
## educationPostgraduate    0.0015   0.0015
## age                      0.0097   0.0327
## urban                    0.0133   0.0133
## 
## Summary of Balance for Matched Data:
##                       Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance                     0.4977        0.4906          0.1454     1.1165
## income                   61281.0371    59372.1803          0.0808     1.0734
## educationCollege             0.3245        0.3671         -0.0910          .
## educationHigh School         0.3367        0.2860          0.1073          .
## educationPostgraduate        0.3387        0.3469         -0.0171          .
## age                         48.8371       48.4934          0.0193     1.0437
## urban                        0.6227        0.6085          0.0293          .
##                       eCDF Mean eCDF Max Std. Pair Dist.
## distance                 0.0391   0.0730          0.1458
## income                   0.0266   0.0649          0.8584
## educationCollege         0.0426   0.0426          0.5502
## educationHigh School     0.0507   0.0507          0.4850
## educationPostgraduate    0.0081   0.0081          0.8829
## age                      0.0100   0.0325          1.0471
## urban                    0.0142   0.0142          0.9751
## 
## Sample Sizes:
##           Control Treated
## All           507     493
## Matched       493     493
## Unmatched      14       0
## Discarded       0       0

matched_data <- match.data(matchit_model)

plot(matchit_model, type = "qq")

Bu kişileri eşleştirdikten sonra elimizde birbirine benzer iki grup var. Bu kişiler sosyolojik değerleri benzeyen ancak treatment olarak ayrılan iki grup oluşturdu. Aynı RCT’lerde istediğimiz gibi grup oluşturmuş olduk. Bundan sonra ise regresyon analiz ile grupları karşılaştırabiliriz.

model_outcome <- glm(turnout ~ compulsory_voting, 
                     data = matched_data, 
                     family = binomial)
summary(model_outcome)

## 
## Call:
## glm(formula = turnout ~ compulsory_voting, family = binomial, 
##     data = matched_data)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.72373    0.09604   7.536 4.85e-14 ***
## compulsory_voting  0.28235    0.13989   2.018   0.0436 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1199.8  on 985  degrees of freedom
## Residual deviance: 1195.8  on 984  degrees of freedom
## AIC: 1199.8
## 
## Number of Fisher Scoring iterations: 4

Burada bir sonuç bulamadık ama bu veriyi kendiğimiz kurduğumuz için geldi. Gerçek dünya verisinde gelme ihtimali yüksekti.

Bu yöntem güzel gözükse bile eksikleri de var. Mesela uygun eşler bulunmazsa gereksiz yere veri kaybına neden olabilir. Başka bir eksik yanı ise incelediğimiz değişkenlere olan bağlılığımız. Eğer incelemedeğimiz önemli bir değişken verse sonuçarı önemli ölçüde etkileyebilir. Aynı zamanda sonuçlar seçtiğimiz metoda da çok bağlı olacağı için burada da dikkat etmemiz gerekir. Son olarak iki grup arasında balans tutmayabilir bu da sonuçları saptırır.

Difference in Differences

Bu yöntem Rubin Causal Model olarak da bilien potential outcomes framework temeline kuruludur. Bu mantığa göre incelediğimiz her unit farklı treatment kondisyonları altında potansiyel sonuçlara sahiptir ancak biz bunların sadece bir tanesini inceleyebiliriz.

DiD modeli ise bunun üstesinden treatmentın veriliş tarihine göre grupları karşılaştırarak verir. Treatment alan ve almayan grubun farkını treatment aldığı süreden itibaren karşılaştırır ve aralarında fark var mı yok mu ona bakar. Buradaki mantık kontrol grubunun, treatment grubu treatment almasaydı neye benzeyeceğini göstermesi.

$Y_{it} = \alpha + \beta_1 D_i + \beta_2 T_t + \delta (D_i \times T_t) + \varepsilon_{it}$

D1 treatment grubun dummy variable’ı, T ise zamanın dummy variable’ı. Delta ise zaman ve grubun ortak etkisi.

$\hat{\delta} = (\bar{Y}_{T,1} - \bar{Y}_{T,0}) - (\bar{Y}_{C,1} - \bar{Y}_{C,0})$

Formül kısaca böyle de verilebilir. Treatment grubunun treatment aldıktan sonraki haliyle önceki halinin farkına bakıyoruz. Aynısını kontrol grubu için de yapıyoruz. Sonra bu ikisinin arasındaki farka bakarız. Zaten modelin de adı buradan geliyor. Farkların arasındaki farka bakıyoruz.

Bu modelin işe yaraması için çok önemli bir assumption da Parallel Trends assumption. Bu eğer treatment olmasaydı iki grubun da paralel bir şekilde hareket edeceğini yani aralarında değişim açısından bir fark olmayacağını söyler.

Aynı zamanda treatment uygulanmadan önce sonucun etkilenmemesi gerekir. Buna no anticipation deriz. Yani beklentiye göre pozisyon alınması bu çalışmanın etkisini azaltır.

Bunlara ek olarak Stable Unit Treatment Value Assumption (SUTVA) dediğimiz bir assumption da var. Bir grubun treatmentı diğer grubu etkilememeli. Yani gruplar birbirinden tamamen bağımsız olmalı. Bu natural experiment ile benzer gelecektir. Zaten DiD natural experiment yaparken kullandığımız analiz yöntemlerinden birisidir.

Son assumption ise common shocks. Buna göre dış şoklar her iki grubu da etkilemeli. Eğer gruplar birbirinden farklı etkileniyorsa bu violate edilir.

Bunun çok meşhur bir örneği Ekonomi alanında Nobel de kazanan bir çalışma. Bu çalışma 1992 yılında asgari ücretin $4.25 olduğu iki eyaleti karşılaştırıyor: New Jersey ve Pennsylvania. O yıl içinde asgari ücret NJ’de 5 dolara yükseltiliyor. Bu çalışmada hem asgari ücret yükseltilmesinden önce hem de sonra iki eyaletin de istidham oranları karşılaştırılıyor. Beklentilerin aksine işe alımlarda ciddi bir düşüş yaşanmıyor ve daha iyi bir değişim bile oluyor.

library(tidyverse)
library(did)
set.seed(123)
did_data_large <- expand_grid(
  year = 2010:2019,
  state = paste0("State_", LETTERS[1:10]) 
) %>%
  mutate(
    id = row_number(),
    treated = as.numeric(state %in% c("State_B", 
                                      "State_D", 
                                      "State_F") 
                         & year >= 2015),
    time = year - 2010,
    group = as.numeric(state %in% c("State_B", 
                                    "State_D", 
                                    "State_F")),
    gname = ifelse(group == 1, 2015, 0),
    outcome = 50 + 2 * time + 5 * group + 10 * 
      treated + rnorm(n(), mean = 0, sd = 5)
  )

Parallel trendsi kontrol etmek için:

ggplot(did_data_large %>% filter(year < 2015) %>% 
         group_by(year, group) %>% 
         summarise(mean_outcome = mean(outcome)), 
       aes(x = year, y = mean_outcome, color = factor(group))) +
  geom_line(size = 1.2) +
  labs(x = "Year", y = "Average Outcome", color = "Group") +
  scale_color_discrete(labels = c("Control", "Treated")) +
  scale_x_continuous(breaks = 2010:2014) +
  theme_minimal() +
  ggtitle("Parallel Trends Pre-Treatment (Averaged by Group)")

Bu veride bu assumption violate edilmiş diyebiliriz. Trendlerin paralel gitmediğini görüyoruz. Ancak bu simüle veri olduğu için sanki burada bir sıkıntı yokmuş gibi bakalım.

did_data_large <- did_data_large %>%
  mutate(post_treatment = ifelse(year >= 2015, 1, 0))
lm_model <- lm(outcome ~ group * post_treatment + 
                 factor(state), data = did_data_large)
summary(lm_model)

## 
## Call:
## lm(formula = outcome ~ group * post_treatment + factor(state), 
##     data = did_data_large)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.2785  -4.0954  -0.0822   3.6607  13.2730 
## 
## Coefficients: (1 not defined because of singularities)
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           53.6229     1.9179  27.959  < 2e-16 ***
## group                  7.8496     2.8231   2.780  0.00664 ** 
## post_treatment        11.2230     1.3562   8.276 1.27e-12 ***
## factor(state)State_B  -3.2278     2.5372  -1.272  0.20665    
## factor(state)State_C   0.2961     2.5372   0.117  0.90736    
## factor(state)State_D  -0.9006     2.5372  -0.355  0.72346    
## factor(state)State_E  -0.1681     2.5372  -0.066  0.94733    
## factor(state)State_F       NA         NA      NA       NA    
## factor(state)State_G   1.6889     2.5372   0.666  0.50738    
## factor(state)State_H  -1.3455     2.5372  -0.530  0.59724    
## factor(state)State_I  -0.2263     2.5372  -0.089  0.92913    
## factor(state)State_J   0.8262     2.5372   0.326  0.74548    
## group:post_treatment   7.7900     2.4760   3.146  0.00226 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.673 on 88 degrees of freedom
## Multiple R-squared:  0.7181, Adjusted R-squared:  0.6829 
## F-statistic: 20.38 on 11 and 88 DF,  p-value: < 2.2e-16

Bu sonuçlara baktığımız zaman treated bizim odaklanmak istediğimiz etki. Etki pozitif görünüyor ve p-değeri 0.05 altında olduğu için null hipotezi reddediyoruz.

Did etkisine bakmak için

ggplot(did_data_large %>% 
         group_by(year, group) %>% 
         summarise(mean_outcome = mean(outcome)), 
       aes(x = year, y = mean_outcome, color = 
             factor(group))) +
  geom_line(size = 1.2) +
  geom_vline(xintercept = 2015, linetype = 
               "dashed", color = "red") +
  labs(x = "Year", y = "Average Outcome", 
       color = "Group") +
  scale_color_discrete(labels = 
                         c("Control", "Treated")) +
  scale_x_continuous(breaks = 
                       c(2012, 2015, 2018)) +
  theme_minimal() +
  ggtitle("DiD Visualization (Averaged by Group)")

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

Son olarak da daha düzenli bir analiz için did paketini kullanabiliriz.

did_data <- expand_grid(
  year = 2010:2019,
  state = c("A", "B")
) %>%
  mutate(
    id = row_number(),
    treated = as.numeric(state %in% c("B") & year >= 2015),
    time = year - 2010,
    group = as.numeric(state %in% c("B")),
    gname = ifelse(group == 1, 2015, 0),
    outcome = 50 + 2 * time + 5 * group + 10 * 
      treated + rnorm(n(), mean = 0, sd = 5)
  )

Burada değişkenleri artık analize göre girebilirsiniz ancak control_group argümanı ve panel argümanı verinin tipine bağlı olarak değişiyor.

did_result_large <- att_gt(
  yname = "outcome",
  gname = "gname",
  idname = "id",
  tname = "year",
  xformla = ~1, 
  data = did_data_large,
  control_group = "nevertreated",
  panel = FALSE 
)

## Warning in pre_process_did(yname = yname, tname = tname, idname = idname, : Be aware that there are some small groups in your dataset.
##   Check groups: 2015.

summary(did_result_large)

## 
## Call:
## att_gt(yname = "outcome", tname = "year", idname = "id", gname = "gname", 
##     xformla = ~1, data = did_data_large, panel = FALSE, control_group = "nevertreated")
## 
## Reference: Callaway, Brantly and Pedro H.C. Sant'Anna.  "Difference-in-Differences with Multiple Time Periods." Journal of Econometrics, Vol. 225, No. 2, pp. 200-230, 2021. <https://doi.org/10.1016/j.jeconom.2020.12.001>, <https://arxiv.org/abs/1803.09015> 
## 
## Group-Time Average Treatment Effects:
##  Group Time ATT(g,t) Std. Error [95% Simult.  Conf. Band]  
##   2015 2011   0.7143     4.1252       -9.2640     10.6927  
##   2015 2012  -7.1224     3.6453      -15.9402      1.6953  
##   2015 2013   3.9654     3.1076       -3.5517     11.4824  
##   2015 2014   1.3298     5.2260      -11.3114     13.9710  
##   2015 2015  13.1605     5.3819        0.1421     26.1788 *
##   2015 2016   4.1650     5.4166       -8.9373     17.2673  
##   2015 2017   5.7902     6.1050       -8.9772     20.5576  
##   2015 2018   8.9478     5.0289       -3.2166     21.1122  
##   2015 2019   3.2021     5.2135       -9.4089     15.8131  
## ---
## Signif. codes: `*' confidence band does not cover 0
## 
## P-value for pre-test of parallel trends assumption:  0.32932
## Control Group:  Never Treated,  Anticipation Periods:  0
## Estimation Method:  Doubly Robust

Bu sonuçlar bize daha açıklaycı bir resim çiziyor. 2015 yılı significant gözüküyor. Bu da 2015 yılında treatmentın etkili olduğunun bir göstergesi. Ancak sonasında bu signigicance kaybolmuş. Bu da ya etkisinin azaldığını ya da önceki kadar önemli bir fark olmadığının göstergesi. Zaten graik de bize benzer bir veri vermişti. Etkiler ise treatment sonrasında pozitif duruyor.

Bu modelin eksikleri arasında outlier verilerin potansiyel olarak ciddi etkisi, external validitiynin düşük olması var. Aynı zamanda anticipation olmamasını da her zaman kabul etmemi zor çünkü sosyal bir alan ve etkilenmesi çok olası. Parallels trend assumtion da ihlal edilmesi çok muhtemel. Çok güçlü ve etkili bir model olmasına rağmen bu tarz güçlü assumptionlara bağlı olması onun bir eksikliği denebilir.

Instrumental Variables

İstatistik analizlerinde endogeneity dediğimiz bir problem vardır. Mesela eğitim ve siyasi katılım arasındaki ilişkiyi inceliyoruz. İki değişkeni regresyona soktuk ve olumlu sonuç bulduk. Bu tahmin ettiğimiz gibi daha eğitimli kişilerin daha yüksek siyasi katılıma sahip olduğunu mu gösterir?

Belki. Bu ikisini birden etkileyen bazı faktörler olabilir. Mesela kişinin motivasyonu onun hem eğitimini hem de siyasi katılımını etkileyebilir. Eğer biz bu değişkeni dışarıda bırakırsak analiz sonuçlarımız sapmış olur.

Endogeneity bağımsız değişkenlerimizin error ile korele olmasıdır. Bunun 3 temel sebebi vardır:

Omitted variable bias: Eğer önemli bir değişkeni verimize eklemezek coefficientların esitmationı biased olur.
Reverse causality veya Simultaneity: Hem X Y’yi hem de Y X’i etkilediği durumlarda olur.
Measurement error: Bağımsız değişkenimizin ölçümünde hata varsa yine bias olarak karşımıza çıkar.

Instrumental variable ise bu noktada devreye girer ve endogeneity sorununu çözmeye çalışır. IV’yi bir Z değişkeni olarak düşünebiliriz. Endogenous değişkenle yani X ile korele ama error term ile korele olmayan değişkenlere denir. Böylece X’in error term ile ilişkili olmayan kısmını izole etmiş oluruz ve böylece X’in Y’ye causal etkisini ölçebiliriz.

Bunu iki aşamalı bir süreçle yaparız. Önce bu enstrümanı kullanarak X’i tahmin ederiz. Sonrasındaysa bu tahminimizi Y’yi tahmin etmek için kullanırız.

Eğitim ve siyasi katılım arasındaki ilişkiye dönelim. En yakın okula olan mesafeyi instrument olarak düşünelim. Bu eğitim seviyesini etkileyebilecek bir faktör ancak doğrudan siyasi katılım üstüne bir etkisi yok. Siyasi katılım üstüne etkisi sadece eğitim aracılığıyla oluyor. Bu da onu iyi bir instrument yapar.

Bir instrumentın başarılı olarak alınması için farklı gereklilikler vardır. Bunlar:

Instrument olarak aldığımız Z değeri, X ile güçlü bir ilişkiye sahip olmalı. Yani Z’yi kullanarak X’i tahmin edebilmeliyiz.
Aynı zamanda error term ile korele olmamalıdır. Yani bu değişken, Y değişkenini doğrudan etkilememli, etkisi sadece X aracılığıyla olmalı. Okula yakınlık > eğitim > siyasi katılım olarak düşünebiliriz.

Daha önceden de söylediğim gibi bunu 2 aşama ile yaparız. Önce bu enstrümanı kullanarak X’i tahmin ederiz. Bu X’in varyasyonunu izole etmemize yarar. Böylece kullandığımız instrumentın etkisiyle değişen X’i öçlmüş oluruz. Sonrasındaysa bu tahminimizi Y’yi tahmin etmek için kullanırız. Böylece de endogeneity problemi olmadan X’i Y’i tahmin etmek için kullanmış oluruz. Buna da Two-Stage Least Squares (2SLS) metodu deriz ve IV altında en çok kullanılan metod budur.

library(AER)
set.seed(123)
n <- 1000
z <- rnorm(n) 
x <- 0.5 * z + rnorm(n)
y <- 2 * x + rnorm(n) 
data <- data.frame(x = x, y = y, z= z)

Bunu iki şekilde yapabiliriz.

first_stage <- lm(x ~ z, data = data)
summary(first_stage)

## 
## Call:
## lm(formula = x ~ z, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0279 -0.6914  0.0043  0.7087  3.2911 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.04105    0.03183    1.29    0.198    
## z            0.58805    0.03211   18.32   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.006 on 998 degrees of freedom
## Multiple R-squared:  0.2516, Adjusted R-squared:  0.2508 
## F-statistic: 335.4 on 1 and 998 DF,  p-value: < 2.2e-16

second_stage <- lm(y ~ fitted(first_stage), data = data)
summary(second_stage)

## 
## Call:
## lm(formula = y ~ fitted(first_stage), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.6982 -1.6365 -0.0032  1.5333  7.6762 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.01847    0.07183  -0.257    0.797    
## fitted(first_stage)  1.96757    0.12277  16.027   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.263 on 998 degrees of freedom
## Multiple R-squared:  0.2047, Adjusted R-squared:  0.2039 
## F-statistic: 256.9 on 1 and 998 DF,  p-value: < 2.2e-16

iv_model <- ivreg(y ~ x | z, data = data)
summary(iv_model)

## 
## Call:
## ivreg(formula = y ~ x | z, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.81289 -0.63188 -0.01815  0.66055  3.46755 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.01847    0.03111  -0.594    0.553    
## x            1.96757    0.05318  37.002   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9801 on 998 degrees of freedom
## Multiple R-Squared: 0.8508,  Adjusted R-squared: 0.8506 
## Wald test:  1369 on 1 and 998 DF,  p-value: < 2.2e-16

Bunu instrument kullanmadan yapacağımız regresyonla karşılaştıralım

ols_model <- lm(y ~ x, data = data)
summary(ols_model)

## 
## Call:
## lm(formula = y ~ x, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8344 -0.6389 -0.0293  0.6647  3.4311 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.02074    0.03098  -0.669    0.503    
## x            2.01243    0.02663  75.563   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9787 on 998 degrees of freedom
## Multiple R-squared:  0.8512, Adjusted R-squared:  0.8511 
## F-statistic:  5710 on 1 and 998 DF,  p-value: < 2.2e-16

Instrumentın işe yarayıp yaramadığını test edebiliriz. Bunun için instrumentın X’i tahmin etmesinde kullandığımız F istatistiğine bakarız. Eğer 10’un üstündeyse güçlü diyebiliriz.

summary(first_stage)$fstatistic

##    value    numdf    dendf 
## 335.4354   1.0000 998.0000

Eğer birden çok instrumentımız olsaydı Hansen J-test kullanarak bunların exogenous olup olmadığına bakardık. Amacımız bağımsız olmaları çünkü.

summary(iv_model, diagnostics = TRUE)$diagnostics

##                  df1 df2   statistic      p-value
## Weak instruments   1 998 335.4354362 8.049720e-65
## Wu-Hausman         1 997   0.9534897 3.290695e-01
## Sargan             0  NA          NA           NA

Diğer yöntemler gibi bunun da eksikleri var elbet. İlk olarak iyi bir instrument bulmak çok zor. Bulan kişi bundan doktora tezi bile çıkartır. Hem Y’yi doğrudan etkilmeyen hem de X ile güçlü ilişkiye sahip olması kriteri bunu zorlaştırıyor. Bir çok instrument mantıklı gelmesine rağmen zayıf ilişkiye sahip. Böylece hem regresyon biased oluyor hem de standard errorler büyüyor. Bu da analizin gücünü zayıflatıyor. Hatta bu durumda düz regresyondan daha kötü sonuç alabiliyoruz.

Başka bir sorun ise overfitting. Eğer açıklayıcı değişkenden daha çok instrument kullanırsak bu analizi o veriye göre dizayn etmemize sebep olabilir ve böylece başka bir veride çalışması ve gücü azalmış olur.

Regression Discontinuity Design

Regression Discontinuity Design (RDD) deney yapmanın mümkün olmadığı veya etik bulunmadığı durumlarda kullanabileceğimiz bir yöntemdir. Quasi experimental bir metoddur ve continious olarak assume edebileceğimiz bir değişkende bir cutoff belirleme üstüne kuruludur. Bu metoda göre cutoff noktasının hemen üstü ve hemen altında bulunan unitler veya kişiler biribrine çok benzerdir. Ancak bir gruba treatment uygulanırken diğerine uygulanmaz. İki grubun ayrışması da bu noktada olur. Treatment edilen ve kontrol grubu arasında bir fark oluşursa da treatmentın etkili olduğunu söyleyebiliriz. Buradaki değişkene running veya forcing variable deriz.

Mesela bir sınavda 85 puanı cutoff point olarak belirleyelim. Bunun üstünde alanlara burs verilsin, diğerlerine verilmesin. Daha sonradan bu öğrencilerin yıl sonu performansı ölçülsün. Eğer burs verilen öğrenciler burs verilmeyen öğrencilere göre daha başarılı olduysa bu bursun eğitimi güçlendirici etkisini anlamamıza yarar. En başta 86 ve 84 alan öğrenciler çok benzerken daha sonrasında fark gözlemlediysek bunun sebebi treatment olur. Yani sadece altında veya üstünde kalanı değil bu thresholda çok yakın olanları karşılaştırırız.

Özetle 4 basic componentı vardır:

Running variable(X): Continious olur ve unitlerimizin treatment alıp almayacağını belirleyen değişkendir.
Cutoff point(C): Bu unitlerin treatment alıp almayacağını belirleyen threshold değeridir.
Treatment(T): Kişilerin treatment alıp almadığını gösterir.
Outcome(Y): Treatmentın etkisini ölçtüğümüz değişkendir.

İki tip RDD vardır diyebiliriz.

Sharp RDD: Treatment ayrımı keskin bir şekilde cutoff point ile uygulanır. Eğer unit cutoff üstündeyse alır ama altındaysa almazgibi.
Fuzzy RDD: Treatment ayrımı önceki kadar keskin değildir. Mesela düşük puan alan ancak geliri de yüksek olmayan öğrenciler burs alabilir. Bu durumda treatmentın uygulanma ihtimali cutoff noktasında çok yükselir ama 0’dan 1’e değil de 0.30 ile 0.6 arasında gibi.

Siyaset biliminden bir örnek verecek olursak İngiltere’de adayların bölgelerine bakabiliriz. Seçim bölgelerinde %51 ile kazanan aday ile %48 ile kaybeden arasında çok ciddi bir fark olmayacaktır ancak bu iki aday yakın olmasına rağmen biri seçimi kazanmış diğeri kaybetmiş oluyor. Bunu adayların gelecekti siyasi pozisyonlarının değişimi üstüne etkisine bakarsak seçimi kazanmanın nedensel etkisini gösterebilir.

Veya gerçek bir çalışma yine İngiltere’de adaylara bakıyor. Bu şekilde kazanan ve kaybeden adayların malvarlığının ne kadar değiştiğini inceliyor. Böylece seçim kazanmanın gelire etkisini incelemiş oluyor.

Buradaki temel assumption treatment uygulanan ve uygulanmayan kişilerin treatment haricine çok benzer oldukları. Cutoff çevresinde değer almaları da bu assumptionı destekler nitelikte.

İkinci bir assumption ise bu cutoffun treatment olmadan outcome’ı etkilemeyecek olması. Mesela sınav sonuçlarında en kötü öğrencilere bile 84 verme değil daha doğru ve continious bir ayrım yapılması lazım.

library(rdrobust)
library(tidyverse)
library(rdd)

## Loading required package: Formula

set.seed(123)
N <- 1000
X <- runif(N, -1, 1) 
Y <- 3 + 2 * X + 5 * (X >= 0) + rnorm(N)

Görsel olarak inceleyebiliriz.

ggplot(data.frame(X, Y), aes(x = X, y = Y)) +
  geom_point(alpha = 0.3) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_smooth(data = subset(data.frame(X, Y), X < 0),
              method = "lm", se = FALSE, color = "blue") +
  geom_smooth(data = subset(data.frame(X, Y), X >= 0), 
              method = "lm", se = FALSE, color = "red") +
  theme_minimal() +
  labs(title = "RDD Plot", x = "Running Variable", 
       y = "Outcome")

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

rdd_model <- RDestimate(Y ~ X, cutpoint = 0)
summary(rdd_model)

## 
## Call:
## RDestimate(formula = Y ~ X, cutpoint = 0)
## 
## Type:
## sharp 
## 
## Estimates:
##            Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|)  
## LATE       0.8810      879          5.074     0.1526      33.25    2.161e-242
## Half-BW    0.4405      448          5.095     0.2221      22.94    1.962e-116
## Double-BW  1.7620     1000          5.120     0.1315      38.94     0.000e+00
##               
## LATE       ***
## Half-BW    ***
## Double-BW  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## F-statistics:
##            F     Num. DoF  Denom. DoF  p
## LATE       2892  3         875         0
## Half-BW    1153  3         444         0
## Double-BW  3879  3         996         0

RDD’nin temel amacı local average treatment effect (LATE) dediğimiz kavramı tahmin etmek yani tam olarak cutoff pointlerde çevresinde olan farkı inceliyor kalanını değil. Bu sonuçlara göre treatment cutoff noktasında etkili olmuş.

Bandwith dediğimiz kavram ise cutoff’un ne kadar çevresini almamız gerektiğine bakar. Küçük bandwith karşılaştırma gücünü artırır ancak sample size azalır. RDestimate otomatik olarak seçer ama başka bandwithler de deneyebiliriz.

rdd_model2 <- RDestimate(Y ~ X, cutpoint = 0, bw = 5)
summary(rdd_model2)

## 
## Call:
## RDestimate(formula = Y ~ X, cutpoint = 0, bw = 5)
## 
## Type:
## sharp 
## 
## Estimates:
##            Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|)
## LATE        5.0       1000          5.129     0.1296      39.56    0       
## Half-BW     2.5       1000          5.124     0.1304      39.31    0       
## Double-BW  10.0       1000          5.131     0.1294      39.64    0       
##               
## LATE       ***
## Half-BW    ***
## Double-BW  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## F-statistics:
##            F     Num. DoF  Denom. DoF  p
## LATE       4080  3         996         0
## Half-BW    3982  3         996         0
## Double-BW  4122  3         996         0

rdd_result <- rdrobust(Y, X)
summary(rdd_result)

## Sharp RD estimates using local polynomial regression.
## 
## Number of Obs.                 1000
## BW type                       mserd
## Kernel                   Triangular
## VCE method                       NN
## 
## Number of Obs.                  507          493
## Eff. Number of Obs.             186          178
## Order est. (p)                    1            1
## Order bias  (q)                   2            2
## BW est. (h)                   0.358        0.358
## BW bias (b)                   0.574        0.574
## rho (h/b)                     0.624        0.624
## Unique Obs.                     507          493
## 
## =============================================================================
##         Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
## =============================================================================
##   Conventional     5.150     0.244    21.083     0.000     [4.671 , 5.628]     
##         Robust         -         -    17.919     0.000     [4.638 , 5.777]     
## =============================================================================

Bu metod ise bizim için en optimal gördüğü bandwithi kendisi seçer. Bu RDestimate metoduna göre daha güçlü olabilir.

Placebo test dediğimiz şey ise belirlediğimiz cutoff dışında bir cutoff belrleyip analizi yapmak oluyor. Eğer yine significant sonuç çıkarsa bizim RDD assumptionları ile sorun olduğunuzu gösterir.

rdd_model3 <- RDestimate(Y ~ X, cutpoint = 0.7)
summary(rdd_model3)

## 
## Call:
## RDestimate(formula = Y ~ X, cutpoint = 0.7)
## 
## Type:
## sharp 
## 
## Estimates:
##            Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|)   
## LATE       0.2781     272           -0.2627   0.2514      -1.0448  0.29614    
## Half-BW    0.1390     145           -0.2111   0.3719      -0.5677  0.57025    
## Double-BW  0.5562     424           -0.3783   0.1945      -1.9452  0.05175   .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## F-statistics:
##            F        Num. DoF  Denom. DoF  p        
## LATE        2.1519  3         268         1.881e-01
## Half-BW     0.6758  3         141         8.635e-01
## Double-BW  18.5743  3         420         4.958e-11

Bu sonuçlar significant çıkmadığı için iyi bir cutoff seçtik diyebiliriz.

Elbet her metod gibi eksikleri var. Sample size özellikle cutoff çevresinde büyük olmalı. İki grubu karşılaştırmak istiyorsak 5-10 gibi küçük sample bizim işimize yaramayacaktır. İkinci olarak da cutoff noktasından uzaklaştıkça sonuçlarımızın genellenebilirliği azalır. Bu durumda unitler birbirlerinden çok farklı olabilir ve sırf bu farklarından dolayı sonucu elde etmiş olabilirler. Bizim amacımız benzer olanları karşılaştırmak.

Panel Data

Ekonometri alanının canı ciğeri olan bir metod ise panel data. Buradaki causal inference metodlarının aksine kullanımı çok daha yaygın bir yöntem. Başka adı ise cross-sectional longitudinal veridir. Cross sectional veri aynı dönem içinde incelenmiş farklı unitler iken longitudinal(time-series) veri bir unitin farklı zamanlarda incelenmesidir. Yani panel data farklı unitlerin farklı zamanlarda incelenmesidir.

Mesela 2019 yılında Türkiye’nin illerindeki işsizlik oranlarını alıyorsak bu cross-sectional veridir. Ankara’nın 2010-2019 arası işsizlik versini alıyorsak bu longitudinal veridir. Ancak Türkiye illerinin 2010-2019 arası işsizlik versini alıyorsak bu panel data olur.

Bu objeleri hem kendi aralarında hem de kendi içlerinde incelememiz causal inference için önemli bir güç oluşturur. Nedensellik çıkarımı yapmak istiyorsak unitlerin time-invariant confounderları olabilir. Mesela bir öğrencinin dönem içindeki performasını takip etmek istiyoruz. Öğrencinin zekası burada time-invariant confounder olarak alınabilir. Ne kadar eğitim veririsek de değiştirsek de bu değişmeyeceği için bunu kontrol etmemiz gerekir ve panel data ile bunu kontrol edebiliriz. Ancak bunlar unitler arasında değiştiği için bizim için önemlidir.

Aynı zamanda temporal etkileri de görebiliriz. Böylece neyin önce geldiğini ve bunun nelere neden olabileceğini anlayabiliriz. Böylece deneylerde yaptığımız gibi X ve Y arasında neden sonuç ilişkisi kurabiliriz.

Mesela demokrasinin ekonomik büyümeye etkisi olup olmadığınab bakmak istiyorsak 50 ülkeyi 30 yıl boyunca inceleyebiliriz. Burada ülkenin o yıl demokrasi olup olmadığı, GDP artışı ve ticari açıklık gibi belli değişkenleri ekleyebiliriz.

Bunu da istersek coğrafya gibi değişmeyen faktörleri kontrol ederek yapabiliriz. Burada panel data ve genel olarak multi-level modeller ile ilgili çok temel bir ayrıma değinmek gerekiyor.

Temel olarak 3 tip panel data vardır:

Pooled OLS: Datanın panel formatını ignore edip cross sectional gibi davranmak için bunu kullanırız. Veriyi temel olarak anlamak için iyi olabilir ama zaman içinde değişmeyen faktörler mevcutsa biased sonuçlar almamıza sebep olur.
Fixed Effects(FE): Time-invariant yani zaman içinde değişmeyen karakteristikleri her unite kendi inteceptini vererek kontol eder. Buradaki temel assumption unitler arasındaki farkın unitlere özel değişkenlerle açıklanabileceğidir. Böylece her unit kendi içinde incelenmiş olur ve unitlerin değişmez ve baz özelliklerinin etkisi kaybolmuş olur. Yani zeki bir öğrenciyi genelle karşılaştırmak yerine kendi içinde değerlendiririz. Bunu da her unit için ayri bir dummy koyarak yapar.
Random Effects(RE): Bu model ise bu etkilerin random olduğunu ve sonucu etkilemeyeceğini varsayar. Yani burada independent değişkenler ile korele ve bizim observe edemediğimiz bir etkinin olduğunu söylüyor. Öyle bir değişken olsa bile bunun bir pattern ile değil de random şekilde dağıldığını ve bunun bize sorun yaratmayacağını assume ediyor. Burada ise random bir intercept ve kişiler ile zaman arasındaki varyasyonu kapsayan bir error term ekliyor.

Bu iki metod arasında seçim yapmak bizim teorimize kalmış. Eğer biz observe edemediğimiz ama önemli olabilecek faktörler olduğunu düşünüyorsak fixed effects kullanmalıyız. Ancak bunlar random ve etkili değil diyorsak da random effects kullanmalıyız. Bu iki metod arasında seçim yapmak için Hausman Test kullanırız. Bu test unique errorlerin yani random effectin regresorler ile korele olup olmadığına bakar. Eğer significant çıkarsa fixed effects çıkmazsa random effects kullanırız.

Fixed effects kullanmanın bir dezavantajı generaliziblity’i azaltmasıdır. Her unitin kendi intercepti olduğu için sadece ele alınan birimlerin içindeki değişimlere odaklanır. Bu diğer birimlere benzer souçlar çıkartmayı zorlaştırır.

Yani, fixed effects model, ele alınan birimlerin içindeki gözlemlenmeyen faktörleri daha iyi kontrol ederken, random effects modeli daha geniş genelleme imkanı sunar ama daha güçlü varsayımlar gerektirir.

library(plm)
set.seed(123)
countries <- rep(c("A", "B", "C", "D", "E"), 
                 each = 10)
years <- rep(2010:2019, 5)
democracy <- c(rep(0, 25), rep(1, 25))
gdp_growth <- rnorm(50, mean = 2, sd = 1) + 
  democracy * 0.5
data <- data.frame(country = countries, year = years, 
                   democracy = democracy, 
                   gdp_growth = gdp_growth)
pdata <- pdata.frame(data, index = c("country", "year"))

veya plm ile gelen hazır bir veri de kullanabiliriz.

data("Grunfeld", package = "plm")

pooled_model <- plm(inv ~ value + capital, 
                    data = Grunfeld, 
                    model = "pooling")
summary(pooled_model)

## Pooling Model
## 
## Call:
## plm(formula = inv ~ value + capital, data = Grunfeld, model = "pooling")
## 
## Balanced Panel: n = 10, T = 20, N = 200
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -291.6757  -30.0137    5.3033   34.8293  369.4464 
## 
## Coefficients:
##                Estimate  Std. Error t-value  Pr(>|t|)    
## (Intercept) -42.7143694   9.5116760 -4.4907 1.207e-05 ***
## value         0.1155622   0.0058357 19.8026 < 2.2e-16 ***
## capital       0.2306785   0.0254758  9.0548 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    9359900
## Residual Sum of Squares: 1755900
## R-Squared:      0.81241
## Adj. R-Squared: 0.8105
## F-statistic: 426.576 on 2 and 197 DF, p-value: < 2.22e-16

Sonuçları yorumlayalım

fe_model <- plm(inv ~ value + capital, 
                data = Grunfeld, 
                model = "within")
summary(fe_model)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = inv ~ value + capital, data = Grunfeld, model = "within")
## 
## Balanced Panel: n = 10, T = 20, N = 200
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -184.00857  -17.64316    0.56337   19.19222  250.70974 
## 
## Coefficients:
##         Estimate Std. Error t-value  Pr(>|t|)    
## value   0.110124   0.011857  9.2879 < 2.2e-16 ***
## capital 0.310065   0.017355 17.8666 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    2244400
## Residual Sum of Squares: 523480
## R-Squared:      0.76676
## Adj. R-Squared: 0.75311
## F-statistic: 309.014 on 2 and 188 DF, p-value: < 2.22e-16

Sonuçları yorumlayalım

re_model <- plm(inv ~ value + capital, 
                data = Grunfeld, 
                model = "random")
summary(re_model)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = inv ~ value + capital, data = Grunfeld, model = "random")
## 
## Balanced Panel: n = 10, T = 20, N = 200
## 
## Effects:
##                   var std.dev share
## idiosyncratic 2784.46   52.77 0.282
## individual    7089.80   84.20 0.718
## theta: 0.8612
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -177.6063  -19.7350    4.6851   19.5105  252.8743 
## 
## Coefficients:
##               Estimate Std. Error z-value Pr(>|z|)    
## (Intercept) -57.834415  28.898935 -2.0013  0.04536 *  
## value         0.109781   0.010493 10.4627  < 2e-16 ***
## capital       0.308113   0.017180 17.9339  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    2381400
## Residual Sum of Squares: 548900
## R-Squared:      0.7695
## Adj. R-Squared: 0.76716
## Chisq: 657.674 on 2 DF, p-value: < 2.22e-16

Sonuçları yorumlayalım

phtest(fe_model, re_model)

## 
##  Hausman Test
## 
## data:  inv ~ value + capital
## chisq = 2.3304, df = 2, p-value = 0.3119
## alternative hypothesis: one model is inconsistent

Bu test bize random effects modelin daha iyi olabileceğini gösteriyor.

Basit regresyon yöntemlerine göre daha güçlü olsa bile yine de endogeity önemli bir sorun olarak kalıyor. Bu da IV yöntemlerini kullanmayı gerektirebilir. Ayrıca measurement error yani değişkenlerin ölçümünde yaşananbielcek bir sorun heö fe hem de re için önemli biaslar oluşturabilir. Ancak en ciddi görebileceğimiz sorun attrition dediğimiz bir mesele. Burada unitler ölçülmediği için sampledan düşebilirler ve bu da önemli etkilere sahip olabilir. Özellikle de unitlerin düşme sebebi random değilse ve bir etkenin sonucuysa bu analimizi olumsuz olarak etkiler.

Synthetic Control

En son olaraksa en kompleks metod olan synthetic control var. Diğerleri klasik regresyon metodlarının benzeri ve 1-2 farklı analiz yapılmış haliydi. Bu ise çok daha ileri ve henüz bile gelişmekte olan bir alan. Amacı bir interventionın etkisini anlamak. Bunun ten bilinen örneklerinden birisi California’nın tütüne karşı yürüttüğü problemin başarılı olup olmadığını görmek. Ancak ne yazzık ki counterfactual olarak ne olabileceğini göremiyoruz. Elimizde sadece gerçeklik var.

Sentetik Kontol ise bu sorunun önüne treat edilen unitin sentetik bir versiyounu oluşturarak geçiyor. Bunu da treat edilmemiş ve benzer özelliklere sahip unitlerin ağırlıklı kombinasyonunu kullanarak yapıyor. Böylece o unite en benzer treat edilmemiş bir uniti oluşturmaya ve intervention olmasaydı ne olurduyu anlamaya çalışıyor.

Bunu yapmak için predictor değişkenlerini iyi seçmemiz lazım. Bunları belirledikten sonra unitiöize benzer treat edilmemiş unitler bulmak lazım. Sonra algoritma ile unite en benzer uniti oluşturuyoruz ve son olarak da bu iki grubu karşılaştırıyoruz. Bu modelin önemli bir assumptionı da inceleyemediğimiz faktörlerin zaman içinde sabir kalması. Sonuç olarak biz treatment öncesi için bir sentetik karşılaştırma oluşturuyoruz. Eğer burada outcome’ı etikleyecek başka değişikler olursa bu bizim sonuçlarımıza yansır.

En popüler örneklerden birisi Calofirnia ve tütün programıydı. Burada:

Treated unit: California
Outcome: Kişi naşına satılan sigara
Intervention: 988 yılında başlanan tütün programı
Donor pool: Tütün programı uygulamayan diğer 38 ABD eyaleti
Predictor değişkenler: Önceki sigara satışları ve demokrafik ile ekonomik faktörler

Bu çalışmada araştırmacılar sentetik California oluşturyor ve programın sigara satışlarını azalttığını yani başarışı olduğunu buluyor

library(Synth)

## ##
## ## Synth Package: Implements Synthetic Control Methods.

## ## See https://web.stanford.edu/~jhain/synthpage.html for additional information.

data("basque")
dataprep.out <- dataprep(
  foo = basque,
  predictors = c("gdpcap","invest"),  
  predictors.op = "mean",  
  dependent = "gdpcap",  
  unit.variable = "regionno", 
  time.variable = "year",
  treatment.identifier = 17,  
  controls.identifier = c(2, 3, 4, 5, 6, 7, 
                          8, 9, 10, 12),
  time.predictors.prior = 1955:1975, 
  time.optimize.ssr = 1955:1975,
  unit.names.variable = "regionname", 
  time.plot = 1955:1997 
)

## 
##  Missing data- treated unit; predictor: invest ; for period: 1955 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1956 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1957 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1958 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1959 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1960 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1961 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1962 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1963 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data - control unit: 2 ; predictor: invest ; for period: 1955 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data - control unit: 2 ; predictor: invest ; for period: 1956 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.

Burada diğer testlerde olduğu gibi p-değeri almıyoruz.

synth.out <- synth(dataprep.out)

## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.2450866 
## 
## solution.v:
##  0.9668934 0.03310656 
## 
## solution.w:
##  0.004177178 0.005942501 0.006044762 0.8816474 0.00425588 0.004660736 0.004400995 0.003692502 0.08197436 0.003203686

path.plot(synth.res = synth.out, 
          dataprep.res = dataprep.out, 
          Ylab = "GDP per capita", 
          Xlab = "Year", 
          Legend = c("Basque Country", 
                     "Synthetic Basque Country"))

abline(v = 1975, col = "blue", lty = 2, lwd = 2)

Burada hangi unitlerin en çok kullanıma sokulduğunu görebiliriz. Bu da sentetik kontrolün nasıl oluştuğunu görmek için işe yarar.

synth.tables <- synth.tab(dataprep.res = 
                            dataprep.out, 
                          synth.res = synth.out)
synth.tables$tab.w

##    w.weights             unit.names unit.numbers
## 2      0.004              Andalucia            2
## 3      0.006                 Aragon            3
## 4      0.006 Principado De Asturias            4
## 5      0.882       Baleares (Islas)            5
## 6      0.004               Canarias            6
## 7      0.005              Cantabria            7
## 8      0.004        Castilla Y Leon            8
## 9      0.004     Castilla-La Mancha            9
## 10     0.082               Cataluna           10
## 12     0.003            Extremadura           12

Klasik p-değeri kullanmadığımız için place testi yaparak sonucun rassal mı yoksa veriyi güzel kurduğumuz için mi oluşutuğunu görebiliriz.

dataprep.placebo <- dataprep(
  foo = basque,
  predictors = c("gdpcap", "invest"),  
  predictors.op = "mean",  
  dependent = "gdpcap",  
  unit.variable = "regionno", 
  time.variable = "year",
  treatment.identifier = 2,  
  controls.identifier = c(3, 4, 5, 6, 7, 8, 9, 10, 12),
  time.predictors.prior = 1955:1975, 
  time.optimize.ssr = 1955:1975,
  unit.names.variable = "regionname", 
  time.plot = 1955:1997 
)

## 
##  Missing data- treated unit; predictor: invest ; for period: 1955 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1956 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1957 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1958 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1959 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1960 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1961 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1962 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data- treated unit; predictor: invest ; for period: 1963 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data - control unit: 3 ; predictor: invest ; for period: 1955 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.
## 
##  Missing data - control unit: 3 ; predictor: invest ; for period: 1956 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.

placebo.out <- synth(dataprep.placebo)

## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 0.0003042486 
## 
## solution.v:
##  1 0 
## 
## solution.w:
##  0.06014675 0.05804176 0.01448733 0.07992888 0.05403698 0.09519529 0.124022 0.01976704 0.494374

path.plot(synth.res = placebo.out, 
          dataprep.res = dataprep.placebo, 
          Ylab = "GDP per capita", 
          Xlab = "Year", 
          Legend = c("Control Region 2", 
                     "Synthetic Control Region 2"))
abline(v = 1975, col = "blue", lty = 2, lwd = 2)

Burada iki grup arasındaki fark çok güçlü görünmüyor. Bu da sentetik kontrolümüzün işe yaradığına bir işaret.

Modelin önemli eksikleri de var elbet. Güzel bir pre-intervention fit olması önemli. Bunu sağlamazsak iyi bir karşılaştırma yapamayız. Intervention kontrol edilen unitleri etkilememeli ve bunların kararları tamamen bağımsız olmalı. Mesela sigara yasağı ile insanlar yan eyalete gidip sigara alıyorsa bu onu etkileyen bir faktördür. Ayrıca treatment dışında bir etki olmamalı ve hepsi sabit kalıyor diyebilmeiyiz. Mesela ülke çapında başlayan bir sigara kampanyası veya ekonomik kriz bunu etkileyebilir.

R Dersleri - Ali Onur Gitmez

R Dersi 14

Ali Onur Gitmez

2024-07-15

R Markdown

Referans

Causal Inference

Randomized Control Trials

Natural Experiment

Matching

Difference in Differences

Instrumental Variables

Regression Discontinuity Design

Panel Data

Synthetic Control

SON