R ve Analiz

Bu yazıyı yazmaktaki amacım R kullanarak bilimsel çalışma yapmak isteyen fakat bir sebepten başlayamayanlara aslında sürecin ne kadar kolay olduğunu göstermek. Metin içinde kullanılan bilgileri farklı kaynaklardan elde ettim:

1. Paketlerin kullanımını, fonksiyonların içeriğini paketlerde yer alan vinyetlerden okuyarak öğrendim.
2. Buradaki bilgilerin bazılarını Datacamp’de aldığım R kurslardan edindim. Bugüne kadar tamamladığım kursların (career track) isimlerini en sonda yazdım. Derslerin detaylarına bağlantılardan ulaşabilirsiniz.
3. Burada bahsi geçen konuların şahsım adına çok zorlaştığı noktalarda Stack Exchange, Stack Overflow ve ücretsiz R sitelerinde yazılanları takip ederek sonuca ulaştım. Hangi sayfa benim için en kolay anlatımı yaptıysa (ben nasıl öğrendiysem), ona atıf yaptım. Kullandığım bütün kaynakları kodlardan önce paylaştım.

Önce hangi kodun ne işe yaradığını yorumlarla yanına yazdım, ancak kopyala-yapıştır ile konsola attığımda çok uğraştırdı. O yüzden daha anlaşılır olması amacıyla kodların altına ekledim.

Buradaki temel anlatımından elde edilecek faydayı büyütmek, farklı boyutlara getirmek tamamiyle sizin çabanızla ilgili. R kullanarak istatistiki analize giriş yapmak isteyenler için faydalı olması dileğiyle.

Haydi başlayalım!

R Paketleri

Burada yapacağımız işlemlere başlamadan önce aşağıdaki R paketlerinin kurulması ve çalıştırılması gerekli. RStudio kullanacağımız için readxl paketi zaten veri aktarımı yapılırken etkinleşeceği için aşağıdaki tabloya eklemedim. Her bir paket için kurulum ve etkinleştirme komutları aşağıdaki tabloda. Kopyala yapıştır yapıp konsola aktarabileceğiniz gibi R alışkanlığınızın artması için tek tek yazmanızı tavsiye ederim.

Veri seti

MODAVICA 2020’de sunduğum “Effect of Life Cycles on Capital Expenditures: Evidence from Borsa Istanbul” bildirisinde kullandığım veri setinin buradaki anlatıma göre düzenlenmiş hali. Setin yapısı aşağıda sunulmuştur. Veri setimizin ilk üç sütunu tanımsal değerler olduğu için bazı fonksiyonlarda aralıkları kullanacağız.

Veri setinin aktarımı

Öncelikle aşağıdaki veri setini bilgisayarınıza indirin.

(yeni sekmede açılacak)

Şahsi tercihim RStudio’da bir proje oluşturmanız. Bunu iki sebeple söylüyorum:

1. .RData dosyasınız sadece burada yaptıklarınızla ilgili olacaktır. Veri setini bozsanız bile başka çalışmalara ait dosyalarınız yerli yerinde duracaktır (deneyim diyelim).
2. En son aşamada kullandığımız not defteri de açtığınız projenin bir parçası olacağı için komutları hatırlamak daha kolay olacak (bu sitenin hiçbir kazanç amacı yoktur).

RStudio’daki “Import Dataset” butonu ile dosyayı indirdiğiniz konumdan içeri vseti ismiyle aktarın.

Çalışma Planımız

1. Deskriptif istatistik hesaplayacağız
   1. Özet istatistikkleri raporlayacağız.
   2. BD, X₁, KD₁,KD₂,KD₃,KD₄ ve KD₅ için korelasyon matrisi hazırlayacağız.
   3. X₁ değişkenini baz alıp gruplar için farklar testi yapacağız
   4. Grupları kullanarak toplam ve yıllara göre grafikler hazırlayacağız.
2. BD = X₁+KD₁+KD₂+KD₃+KD₄+KD₅ modelini yıl ve birim için kukla değişkenler kullanarak sabit etkiler ile test edeceğiz.
   1. Modeli standart hataların düzelterek raporlayacağız.
3. Aynı modeli kantil regresyon ile test edeceğiz.
   1. Kantil regresyon modeline bootstrapping yapacağız.
   2. Kantil regresyon modelinde R² hesaplayacağız.
   3. Kantil Regresyon modelini toplam olarak raporlayacağız.
4. Bütün elde ettiğimiz sonuçları paylaşılabilir hale getireceğiz.

Deskriptif İstatistikler

Özet İstatistiki

Öncelikle veri setimizde ne neymiş bir bakalım. Gözlemlerimin tanımlarının (ilk üç sütun) herhangi bir istatistiki karşılığı olmadığı için istastiki özetimize eklemiyoruz. summary fonksiyonu içindeki vseti[,4:10] tanımının da amacı ilk üçü istemediğimizi belirtmek. Daha kapsamlı özet istatistikler için bu makaleyi kullanabilirsiniz. Bana bu kadar sıfır yetmez daha fazla isterim diyorsanız da üçüncü komutu çalıştırın.

desk<-summary(vseti[,4:10])
desk
options(digits = 7)

Korelasyon Matrisi

Korelasyon matrisini iki şekilde elde edebiliriz R’da. Birinci yöntem R’ın kendi korelasyon matris fonksiyonu olan “cor” ile yapılabilir ancak bu fonksiyona dair esas sorun ilişkilerin istatistiki anlamlılığı raporlanmaz.

kormat<-cor(vseti[,4:10])
kormat

Peki nedir bu sorunun çözümü? Hmisc paketi. Aşağıdaki komutları çalıştırdığınızda karşınıza gelen korelasyon matrisinde iki kısım olacak. Birincisi değişkenler arası korelasyon, ikincisi korelasyonların istatistiki anlamlılığı.

kormat_pv<-rcorr(as.matrix((vseti[,4:10])))
kormat_pv

Farkların Anlamlığı Testi

Elimizdeki veri setindeki X1 değişkeninde 0.00, 0.25, 0.50, 0.75 ve 1.00 değerlerinden 39’ar adet olmak üzere 195 gözlem var. Hiçbir birim için eksik gözlemimiz yok. X1 değişkeni içindeki grupların Y ortalamaları arasında bir fark olup olmadığını inceleyelim şimdi. Değerleri nasıl gruplayacağız.

Grup12 <- vseti %>%
filter(X1 == 0 | X1 == 0.25) %>%
select(X1, BD)

Grup13 <- vseti %>%
filter(X1 == 0 | X1 == 0.5) %>%
select(X1, BD)

Grup14 <- vseti %>%
filter(X1 == 0 | X1 == 0.75) %>%
select(X1, BD)

Grup15 <- vseti %>%
filter(X1 == 0 | X1 == 1) %>%
select(X1, BD)

Grup23 <- vseti %>%
filter(X1 == 0.25 | X1 == 0.5) %>%
select(X1, BD)

Grup24 <- vseti %>%
filter(X1 == 0.25 | X1 == 0.75) %>% 
select(X1, BD)

Grup25 <- vseti %>%
filter(X1 == 0.25 | X1 == 1) %>%
select(X1, BD)

Grup34 <- vseti %>%
filter(X1 == 0.5 | X1 == 0.75) %>%
select(X1, BD)

Grup35 <- vseti %>% 
filter(X1 == 0.5 | X1 == 1) %>%
select(X1, BD)

Grup45 <- vseti %>%
filter(X1 == 0.75 | X1 == 1) %>%
select(X1, BD)

Üstteki komutları kısaca izah etmek isterim. dplyr paketini etkinleştirdiğimiz için %>% operatörü ile bütün komutu ile tek bir satırda yazmak yerine parçalara bölmüş oluyoruz.

Grup12 <- vseti %>%
filter(X1 == 0 | X1 == 0.25) %>%
select(X1, BD)

# Grup12 tablosunu oluştur ve vseti isimli veri setini kullan
# vseti içerisindeki X1 değişkeni için 0 ve 0.25 olan gözlemleri filtrele
# vseti içindeki 0 ve 0.25 olan X1 gözlemlerini için X1 ve BDS değişkenlerini Grup12 tablosuna yerleştir.

Grupları oluşturduğumuza göre anlamlılık testlerimiz yapalım her bir grup için (UC Business Analytics R Programming Guide).

t.test(BD ~ X1, data = Grup12)
t.test(BD ~ X1, data = Grup13)
t.test(BD ~ X1, data = Grup14)
t.test(BD ~ X1, data = Grup15)
t.test(BD ~ X1, data = Grup23)
t.test(BD ~ X1, data = Grup24)
t.test(BD ~ X1, data = Grup25)
t.test(BD ~ X1, data = Grup34)
t.test(BD ~ X1, data = Grup35)
t.test(BD ~ X1, data = Grup45)

Grafiklerin Oluşturulması

Grafikleri oluşturmak için önce her bir grubun ve grupların yıllara göre ortalamasını veri setimize ekleyelim.

vseti<- vseti %>% 
group_by(X1) %>%
mutate(BD_X1_ort=mean(BD))
 

# vseti kullanarak vseti'ne işlemleri yap
# vseti'ni X1 kullanarak grupla 
# vseti'nin sonuna BD'nin X1 gruplarına göre ortalamasını ekle

vseti <- vseti %>%
group_by(X1, FY) %>% 
mutate(BD_X1_FY_ort=mean(BD)) 

# vseti kullanarak vseti'ne işlemleri yap 
# vseti'ni X1 ve FY kullanarak grupla 
# vseti'nin sonuna BD'nin X1-FY gruplarına göre ortalamasını ekle

Şimdi grafikleri çıkartabiliriz. Önce toplamdan gidelim (ggplot2 ve ggplot2 cheat sheet)

ggplot(vseti, aes(x=X1, y=BD_X1_ort, fill=X1))+
geom_point(size=3,aes(colour=X1))+
labs(title="Grafigin Ana Basligi", x ="X1 Gruplari", y = "Ortalama Bagimli Degisken")+
theme(plot.title = element_text(hjust = 0.5))+
scale_y_continuous(breaks = seq(0,0.1,0.01))

#vseti'ni kullanarak BD~X1 grafiğini oluştur.
#noktaları göster, nokta büyüklüğü 3 olsun, X1'e göre renklendir
#eksenleri ve grafiğin adını yaz
#grafiğin adını ortala
#y eksenindeki ara değerleri ayarla

Şimdi yıllara göre grafikleri çıkartalım

ggplot(vseti, aes(x=FY, y=BD_X1_FY_ort, fill=X1))+
geom_line(aes(colour=X1, group=X1))+
geom_point(size=1,aes(colour=X1))+
labs(title="Grafigin Ana Basligi", x ="X1 Gruplari", y = "Yillara Gore Ortalama Bagimli Degisken")+
theme(plot.title = element_text(hjust = 0.5))+
scale_y_continuous(breaks = seq(0,0.1,0.01))+
scale_x_continuous(breaks = seq(2015,2017,1))

#vseti'ni kullanarak BD~FY grafiğini oluştur.
#çizgi grafik oluştur
#noktaları göster, nokta büyüklüğü 3 olsun, X1'e göre renklendir
#eksenleri ve grafiğin adını yaz
#grafiğin adını ortala
#y eksenindeki ara değerleri ayarla

Regresyon Modeli

Çalışmamızdaki temel varsayımımız (hipotez) aşağıdaki gibidir:

H0: X1 BD'yi etkilememektedir
H1: X1 BD'yi etkilemektedir (Beklenen işaret: )

Bu varsayımı test etmek modelimizi en küçük kareler regresyonu olarak kaydedelim. Modele birim ve yıl sayesinde sabit etkileri de test etmiş oluyoruz. Açıkçası birim ve yıl olmadan daha güzel sonuçlar elde edebilirsiniz ancak iyi bir dergide doğrudan kabul görmesi mümkün olmayacak. Modelin sonuçlarını görmek için ikinci komutu çalıştırın.

model<-lm(BD~X1+KD1+KD2+KD3+KD4+KD5+factor(ID)+factor(FY),data=vseti)
summary(model)

#factor() fonksiyonunun amacı her bir birim koduna ve yıla kukla değişken vermek yerine hepsini çalıştırmak. R her bir gruptan bir tanesini analize dahil etmeyecek

Standart Hatalar Düzeltilmiş Model

Doğrusal ile test ettiğimiz modelde bir sorun var. Modelin standart hatalarını birim ve yıllara göre kümeleyerek sonuçların istatistiki olarak sağlamlığını yapacağız.

model_SE<-coeftest(model,vcov=vcovHC(model,type="HC1",cluster=c(PDS$ID,PDS$FY)))

Modelin kantil regresyon ile test edilmesi

Varsayalım ki doğrusal regresyon için gerekli şartları sağlayamadık (Baltagi, 2011) ve modelimizi parametrik olmayan regresyon ile test etmemiz gerekiyor . Koenker ve Basset (1978) tarafından geliştirilen kantil regresyon bu amaçla kullanabileceğiniz yöntemlerden biri. Analiz için gereken paket bizzat regresyonu geliştiren Koenker’e ait. Modeli bir kere çalıştırmak ve kaydetmek yetmeyecek. Ne yapacağız peki? Toplamda 7 kere çalıştırmış olacağız. rq fonksiyonu içindeki tau hangi kantil değerinin test edileceğini göstermektedir.

model_0.05<-rq(BD~X1+KD1+KD2+KD3+KD4+KD5+factor(ID)+factor(FY),data=vseti, tau=0.05)
model_0.10<-rq(BD~X1+KD1+KD2+KD3+KD4+KD5+factor(ID)+factor(FY),data=vseti, tau=0.10)
model_0.25<-rq(BD~X1+KD1+KD2+KD3+KD4+KD5+factor(ID)+factor(FY),data=vseti, tau=0.25)
model_0.50<-rq(BD~X1+KD1+KD2+KD3+KD4+KD5+factor(ID)+factor(FY),data=vseti, tau=0.50)
model_0.75<-rq(BD~X1+KD1+KD2+KD3+KD4+KD5+factor(ID)+factor(FY),data=vseti, tau=0.75)
model_0.90<-rq(BD~X1+KD1+KD2+KD3+KD4+KD5+factor(ID)+factor(FY),data=vseti, tau=0.90)
model_0.95<-rq(BD~X1+KD1+KD2+KD3+KD4+KD5+factor(ID)+factor(FY),data=vseti, tau=0.95)

Kantil regresyon ve bootstrapping

En küçük kareler yönteminde standart hataları birim/yıl olarak kümeleyerek modelin sağlamlığını yapmıştık. Kantil regresyonda neden böyle bir şey yapmıyoruz? Kantil regresyonda standart hataların kümelenmesi için sonuçların yeniden önyüklenmesi gerekiyor (Hao ve Naiman, 2007). Bu işlem esas modelde yaptığımıza göre daha kolay.

summary.rq(model_0.05,se="boot")
summary.rq(model_0.10,se="boot")
summary.rq(model_0.25,se="boot")
summary.rq(model_0.50,se="boot")
summary.rq(model_0.75,se="boot")
summary.rq(model_0.90,se="boot")
summary.rq(model_0.95,se="boot")

Kantil regresyon için R2 hesaplamak

Koenker diyor ki

I don't much like R1, or R2 for that matter, so it isn't likely to be automatically provided in quantreg any time soon.
(R1 veya R2'den çok hoşlanmıyorum, bu yüzden de yakın bir zamanda quantreg'de otomatik olarak yer alması olası değil (Çeviri: GC)

Peki siz bu modellerin altında R² eklemek istiyorsanız ne yapmanız lazım? (Çözüm için Koenker ve Machado, 1999. R-squared in quantile regression, Dimitriy V. Masterov (https://stats.stackexchange.com/users/7071/dimitriy-v-masterov), R-squared in quantile regression, URL (version: 2014-12-17): https://stats.stackexchange.com/q/129251)

fit0 <- rq(BD~1,tau=0.95,data=vseti)
fit1<-rq(BD~X1+KD1+KD2+KD3+KD4+KD5+factor(ID)+factor(FY),data=vseti, tau=0.05)
rho <- function(u,tau=.5)u*(tau - (u < 0))
R1 <- 1 - fit1$rho/fit0$rho
R1

Bu formülü her bir kantil noktası için tek tek (0.05, 0.10, 0.25, 0.50, 0.75, 0.90 ve 0.95) hesaplamanız gerekiyor.

R Notebook

Sizinle paylaştığım kodları çalışmanızın (ödevler de dahil olmak üzere) analizine uygun olarak değiştirebilir, sonuçları tek bir dosyada raporlayabilirsiniz. Kodların ne işe yaradığını gösteremeyen bir not uygulaması yerine gerçek bir çalışma defteri isterseniz R Notebook işinizi görecektir. R Markdown kullanarak metni istediğiniz gibi biçimlendirebilirsiniz. Kodları R Notebook’ta nasıl kullanacaksınız?

```{r}
desk<-summary(vseti[,4:10])
desk
```

Yukarıdaki kod yığını (chunk) hem sizin aldığınız notları saklayacak hem de bu kodun ne anlama geldiğini merak ettiğinizde yığının en sağındaki oynat tuşuna basarak (“Run Current Chunk”) sonuçları da görmenizi sağlayacak.

Datacamp Career Track

Statistician with R (Tamamandı)
Quantitative Analyst with R (Devam ediyor)