4. HAFTANIN TEKRARI

Bu hafta çoklu regresyonun devamı olarak kategorik değişkenlerle regresyon (dummy), düzenleyicilik (moderation) analizi, Anscombe’un dörtlüsü ve nedensel dörtlü konularını işledik. İİlk dönemde işlediğimizSimpson Paradoksuna da değindik.

Dummy ile Kategorik Değişkenler

Regresyon modellerinde kategorik değişkenleri kullanabilmek için dummy -yapay- değişkenler oluşturmak gerekiyor. k düzeye sahip bir kategorik değişken için k-1 tane dummy değişken oluşturulur..

ÖNEMLİ NOT: Referans grubu olarak seçilen kategoriye ait bireyler, tüm dummy değişkenlerde 0 değerini alır.

Derste “Performansd1” veri setinde medeni durum üç düzeydi, iki dummy değişken oluyor: D1 (Evli) ve D2 (Bekar). Referans grup “Diğer” olduğunda, kesişim katsayısı (b0) diğer grubunun yordanan puanını, b1 evliler ile diğerleri arasındaki farkı, b2 bekarlar ile diğerleri arasındaki farkı veriyordu. R’da lm() fonksiyonu kategorik değişkeni otomatik olarak dummy’ye çevirir ya da fastDummies::dummy_cols() ile elle oluşturulabilirmiş.

#DENEME İÇİN YAPTIK AMA KNİT EDERKEN SORUN OLMAMASI İÇİN eval=FALSE yazdık chunkın setupına. Bunu chatten öğrendim hoşuma gitti.

library(fastDummies)

dataf <- dummy_cols(Performansd1, select_columns = 'Medeni')
model_dummy <- lm(Performans ~Medeni_1 + Medeni_2, data = dataf)

model_2 <- lm(Performans ~ Medeni, data = Performansd1)

Moderation Analizi

Düzenleyici etki, bir bağımsız değişken ile bağımlı değişken arasındaki ilişkinin üçüncü bir değişkenin düzeyine göre değişmesidir. Regresyon denklemi: \(Y = a + b_1X + b_2Z + b_3(X \times Z)\)

Derste iki örnek incelendi:

Örnek 1 (Hassles.sav): Sorunlar, sosyal destek ve psikolojik belirtiler. Etkileşim terimi anlamlı çıktı, düşük sosyal destekteki bireylerde sorun arttıkça belirtiler dramatik artarken, yüksek destektekilerde bu artış çok daha yavaştı.

Örnek 2 (duzenleyici.Rds): Uyku saati, kahve tüketimi ve derse dikkat. etkileşim anlamlıydı. Cohen f² ile etki büyüklüğü hesaplandı.

#DERS DENEMESİ

Xc <- c(scale(X, center = TRUE, scale = FALSE))
Zc <- c(scale(Z, center = TRUE, scale = FALSE))

fitMod <- lm(Y ~ Xc + Zc + Xc * Zc)

library(effectsize)
cohens_f_squared(fitMod)

library(rockchalk)
plotSlopes(fitMod, plotx = "Xc", modx = "Zc", modxVals = "std.dev")

Anscombe’un Dörtlüsü

Dört farklı veri seti aynı ortalama, standart sapma, korelasyon ve regresyon katsayılarına sahip olmalarına rağmen tamamen farklı örüntüler gösteriyordu: doğrusal ilişki, eğrisel ilişki, aykırı değerli doğrusal ilişki ve etkili değerli ilişkisizlik. istatistikler tek başına yetmez, verileri mutlaka görselleştirmek gerekir.

Nedensel Dörtlü

Collider (seçilim yanlılığı, düzeltme yapılmamalı),

Confounder (sahte ilişki, düzeltme yapılmalı),

Mediator (mekanizma, amaca göre düzeltilir veya düzeltilmez),

M-Bias (ölçülemeyen faktörler nedeniyle düzeltme zararlı olabilir).

Simpson Paradoksu

Berkeley kabul oranları örneğinde, genel olarak kadınların kabul oranı erkeklerden düşük görünürken, bölüm bazında bakıldığında çoğu bölümünde kadınlar daha yüksek oranda kabul almıştı. Sebebi: kadınlar kabul oranı düşük bölümlere daha fazla başvuruyorlardı. Bir ilişkiyi incelerken önemli bir değişkeni göz ardı etmek sonuçları tamamen yanıltabilir.

FARKLI VERİ SETİ İLE UYGULAMA: iris

Bu hafta derste öğrendiğimiz dummy ve düzenleyicilik analizini R’ın iris veri seti üzerinde uygulayacağım. Bu veri seti 3 tür bitkinin yaprak ölçümlerini içermektedir - 3 düzeyli kategorik değişken dummy kodlama için ideal. Bu veri setini daha önce kullanmıştık ve datacampta da egzersizlerde vardı.

library(dplyr)
library(ggplot2)
library(broom)
library(knitr)
library(fastDummies)
head(iris, 5)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
table(iris$Species)
## 
##     setosa versicolor  virginica 
##         50         50         50

Dummy Denemesi

iris_d <- dummy_cols(iris, select_columns = "Species")
head(iris_d[, c("Species", "Species_setosa", "Species_versicolor", "Species_virginica")])
##   Species Species_setosa Species_versicolor Species_virginica
## 1  setosa              1                  0                 0
## 2  setosa              1                  0                 0
## 3  setosa              1                  0                 0
## 4  setosa              1                  0                 0
## 5  setosa              1                  0                 0
## 6  setosa              1                  0                 0

Referans grubu olarak “setosa”yı seçerek Sepal.Length’i yordayan bir regresyon modeli

model_dummy <- lm(Sepal.Length ~ Species_versicolor + Species_virginica, data = iris_d)
tidy(model_dummy) %>% kable(digits = 3)
term estimate std.error statistic p.value
(Intercept) 5.01 0.073 68.76 0
Species_versicolor 0.93 0.103 9.03 0
Species_virginica 1.58 0.103 15.37 0 
model_oto <- lm(Sepal.Length ~ Species, data = iris)
tidy(model_oto) %>% kable(digits = 3)
term estimate std.error statistic p.value
(Intercept) 5.01 0.073 68.76 0
Speciesversicolor 0.93 0.103 9.03 0
Speciesvirginica 1.58 0.103 15.37 0

Kesişim katsayısı setosa grubunun ortalama Sepal.Length’ini veriyor. Versicolor katsayısı, versicolor türünün setosa’ya göre ne kadar farklı olduğunu, virginica katsayısı da virginica’nın setosa’ya göre farkını gösteriyor.

iris %>%
  group_by(Species) %>%
  summarise(ort = mean(Sepal.Length)) %>%
  kable(digits = 3)
Species ort
setosa 5.01
versicolor 5.94
virginica 6.59

Gerçek ortalamalara bakıldığında katsayıları doğrulandı. setosa ortalaması kesişim değerine ve farklar da eğim katsayılarına eşit.

Moderation Analizi

Petal.Length’in Sepal.Length üzerindeki etkisinin türe göre farklılaşıp farklılaşmadığını test edeceğim. Burada Species kategorik bir düzenleyici olarak modele giriyor.

Derste öğrendiğimiz gibi, önce sürekli değişkeni sabitlicez

iris$Petal.Length_c <- c(scale(iris$Petal.Length, center = TRUE, scale = FALSE))

Modeller farklı farklı

model_ana <- lm(Sepal.Length ~ Petal.Length_c + Species, data = iris)
glance(model_ana)$adj.r.squared
## [1] 0.833
model_mod <- lm(Sepal.Length ~ Petal.Length_c * Species, data = iris)
summary(model_mod)
## 
## Call:
## lm(formula = Sepal.Length ~ Petal.Length_c * Species, data = iris)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7348 -0.2279 -0.0313  0.2437  0.9361 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         6.251      0.637    9.81   <2e-16 ***
## Petal.Length_c                      0.542      0.277    1.96    0.052 .  
## Speciesversicolor                  -0.731      0.641   -1.14    0.256    
## Speciesvirginica                   -1.449      0.658   -2.20    0.029 *  
## Petal.Length_c:Speciesversicolor    0.286      0.295    0.97    0.334    
## Petal.Length_c:Speciesvirginica     0.453      0.290    1.56    0.120    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.336 on 144 degrees of freedom
## Multiple R-squared:  0.84,   Adjusted R-squared:  0.835 
## F-statistic:  152 on 5 and 144 DF,  p-value: <2e-16
tidy(model_mod) %>% kable(digits = 3)
term estimate std.error statistic p.value
(Intercept) 6.251 0.637 9.810 0.000
Petal.Length_c 0.542 0.277 1.959 0.052
Speciesversicolor -0.731 0.641 -1.140 0.256
Speciesvirginica -1.449 0.658 -2.203 0.029
Petal.Length_c:Speciesversicolor 0.286 0.295 0.969 0.334
Petal.Length_c:Speciesvirginica 0.453 0.290 1.563 0.120

R² karşılaştırması

cat("Ana etki R²:", round(glance(model_ana)$adj.r.squared, 3), "\n")
## Ana etki R²: 0.833
cat("Etkileşim R²:", round(glance(model_mod)$adj.r.squared, 3), "\n")
## Etkileşim R²: 0.835

Düzeltilmiş R² değerini karşılaştırıyoruz çünkü derste öğrendiğimiz gibi, değişken eklendikçe R² her zaman artar ama düzeltilmiş R² sadece yeni değişken gerçekten bilgi sağlıyorsa artar.

Etki büyüklüğü

library(effectsize)
cohens_f_squared(model_mod)
## # Effect Size for ANOVA (Type I)
## 
## Parameter              | Cohen's f2 (partial) |      95% CI
## -----------------------------------------------------------
## Petal.Length_c         |                 4.76 | [3.72, Inf]
## Species                |                 0.48 | [0.29, Inf]
## Petal.Length_c:Species |                 0.02 | [0.00, Inf]
## 
## - One-sided CIs: upper bound fixed at [Inf].

Görselleştirme

ggplot(iris, aes(x = Petal.Length, y = Sepal.Length, color = Species)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal() +
  labs(title = "Petal.Length Etkisinin Türe Göre Farklılaşması",
       x = "Petal Uzunluğu", y = "Sepal Uzunluğu")

Grafikte her tür için regresyon doğrusunun eğiminin farklı olup olmadığı görülebilir. Eğimler farklıysa düzenleyicilik etkisi var demek.

Kısmi Regresyonlar

Derste car::avPlots() fonksiyonuyla kısmi regresyon grafikleri incelenmişti:

library(car)
avPlots(model_mod)

Derste Anscombe’un dörtlüsü ile görselleştirmenin önemini gördük. Bunu iris verisi üzerinde şöyle gösterebiliriz: sadece korelasyona baktığımızda tüm türlerde Petal.Length ve Sepal.Length arasında pozitif ilişki var ama bu ilişkinin gücü ve doğası türe göre farklılaşıyor.

iris %>%
  group_by(Species) %>%
  summarise(r = cor(Petal.Length, Sepal.Length),
            n = n()) %>%
  kable(digits = 3)
Species r n
setosa 0.267 50
versicolor 0.754 50
virginica 0.864 50 
ggplot(iris, aes(x = Petal.Length, y = Sepal.Length)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~Species, scales = "free") +
  theme_minimal() +
  labs(title = "Türlere Göre Petal-Sepal İlişkisi")

Tüm veriyi tek bir regresyonla modellersek güçlü bir ilişki görürüz. Ancak alt gruplara bakıldığında ilişkinin yapısı farklılaşıyor. Bu, derste Simpson Paradoksu yani.

YANSITMA

Causal quartet konusunu çok zor kavradım. Aynı veri üzerinde collider, confounder, mediator ve M-bias yapılarının hepsi aynı korelasyon değerlerini üretmesine rağmen doğru analiz stratejisi tamamen farklılaşıyordu. Özellikle mediator ve confounder ayrımı benim için önemliydi: ikisinde de üçüncü bir değişken söz konusu ama birinde düzeltme gerekiyor , diğerinde ise düzeltme yaparsanız asıl etki sıfırlanıyor. Dummy konusu ise geçen dönem jamovide yaptığım analizlerden aşinaydım ama R’da lm() fonksiyonunun otomatik olarak referans grubu seçip dummy’leri oluşturması, fastDummies paketiyle elle kontrol edilmesi ve referans grubu değiştirmenin katsayıları nasıl etkilediğini görmek pekiştirici - ve aslında pratik oldu jamovide kolaydı ama spssde karıştırıcı oluyordu- oldu. Anscombe’un dörtlüsünde de aynı ortalama, aynı korelasyon, aynı regresyon denklemi farklı örüntülere sahip olabiliyor. Analiz yapmadan önce mutlaka veriyi görselleştirmemiz gerekişyor.Hazır paketlerden faydalandığımız için dersin teorik kısmı ağırlaşmış olsa bile kodlama kısmı biraaz daha kolaylaştı gibi. Yani teorik kısmı kavramak hala güç onu anlayıp kodlamak hala çok zor ama çok basitleştiilmiş veri seti ve regresyon bağlantıları bilindiğinde değişkenler DOĞRU etiketlendiğinde (bağımnlı, bağımsız aracı vb.) R üzerinde işlemi yapmak daha kolaydı ilk dönemki fonksiyonlara ve apply ailesine göre.