3. HAFTANIN TEKRARI

Üçüncü haftamızda iki büyük konu işledik: çoklu regresyon modelleri ve yol analizi.

Çoklu regresyonda basit modelden etkileşim modeline kadar ilerlerken, yol analizinde ise lavaan paketiyle birden fazla regresyon denklemini bir arada test etmeyi öğrendik.

Basit ve Çoklu Regresyon

İlk olarak moderndive paketindeki evals veri seti üzerinden basit doğrusal regresyon modeli kurduk. Güzellik puanı (bty_avg) bağımsız değişken, öğretim değerlendirme puanı (score) bağımlı değişken olarak modele almıştık. lm() fonksiyonu ile model kurmuştuk, get_regression_table() ile sonuçlara baktık. Güzellik puanındaki 1 birimlik artışın değerlendirme puanını ortalama 0.067 birim artırdığı tespit ettik.. Ancak R² değeri sadece %3.5 çıktı yani güzellik puanı, ders değerlendirmesindeki varyansın çok küçük bir kısmını açıklıyor. Ardından cinsiyet (gender) değişkeni eklenerek iki farklı çoklu regresyon modeli kuruldu:

  • Paralel eğimler modeli (score ~ bty_avg + gender): Her iki cinsiyet için güzelliğin etkisi aynı varsayılır, sadece kesişim noktaları farklıdır. Erkek hocaların kadınlardan ortalama 0.172 puan daha yüksek aldığı görüldü.

  • Etkileşim modeli (score ~ bty_avg * gender): Güzelliğin etkisinin cinsiyete göre değişip değişmediği test ettik derste. Etkileşim terimi anlamlı çıktı (p = .015), yani cinsiyet bir moderatör olarak güzelliğin etkisini farklılaştırıyordu. R² %3’ten %7’ye yükseldi.

library(moderndive)
library(ggplot2)

model_basit <- lm(score ~ bty_avg, data = evals)
summary(model_basit)

model_paralel <- lm(score ~ bty_avg + gender, data = evals)
get_regression_table(model_paralel)

model_etkilesim <- lm(score ~ bty_avg * gender, data = evals)
get_regression_table(model_etkilesim)

Çoklu Bağlantı Problemi

ISLR paketindeki Credit veri setinde kredi kartı borcu “Balance” incelendi. Hatta derste bu veri setini çağrımakta ve balance ne olduğunu bulamamıştık neyseki kayıtta burasını daha kolay takip edbildim. Limit ve Income değişkenleri ayrı ayrı modele alındığında ikisi de borcu pozitif yönde yorduyordu. Ancak birlikte modele konulduğunda Income’ın katsayısı negatife döndü. Bunun nedeni iki bağımsız değişken arasındaki yüksek korelasyondu (r = .79). Bu durum çoklu bağlantı problemidir, bağımsız değişkenler arası korelasyon .80’i aştığında katsayılar gerçek etkisinden sapar.

library(ISLR)
data(Credit)


lm(Balance ~ Limit, data = Credit)

lm(Balance ~ Income, data = Credit)


lm(Balance ~ Limit + Income, data = Credit)

cor(Credit$Income, Credit$Limit) # 0.79

Path

DERS NOT: Yol analizi, birden fazla çoklu regresyon modelinin lavaan paketiyle bir arada test edilmesidir.

Derste “hastalık faktörleri” örneği üzerinde çalışıldı: egzersiz ve dayanıklılık >form ve stres >hastalık.

Modeldeki temel kavramlar: dışsal değişkenler (egzersiz, dayanıklılık) modelde nedenleri gösterilmeyen değişkenler; içsel değişkenler (form, stres, hastalık) nedenleri modelde açıklanan değişkenler. İçsel değişkenlerin açıklanamayan kısmı bozukluk (disturbance) olarak adlandırılır ve gizil değişken olarak ele alınır.

Model-veri uyumunu değerlendirmek için Ki-kare testi, RMSEA, CFI ve SRMR indeksleri inceledik. İlk modelde uyum zayıf çıktı (RMSEA = .168). Modifikasyon indeksleri incelenerek form -> stres yolu eklendi, ardından anlamsız yollar çıkarılarak daha tutumlu bir model elde ettik.Son modelde uyum indeksleri kabul edilebilir düzeye ulaştık.

library(lavaan)
library(semPlot)

yol_model <- 'stres ~ egzersiz + dayaniklilik
              hastalik ~ egzersiz + dayaniklilik + form + stres
              form ~ egzersiz + dayaniklilik
              egzersiz ~~ dayaniklilik'
yol_fit <- sem(yol_model, veri)
fitMeasures(yol_fit, c("chisq","df","pvalue","cfi","rmsea","srmr"))

yol_model_v2 <- 'stres ~ dayaniklilik
                 hastalik ~ form + stres
                 form ~ egzersiz + dayaniklilik
                 stres ~ form
                 egzersiz ~~ dayaniklilik'

yol_fit_v2 <- sem(yol_model_v2, veri)
fitMeasures(yol_fit_v2, c("rmsea","cfi","srmr"))

FARKLI VERİ SETİ İLE UYGULAMA: mtcars

Derste öğrendiğimiz basit regresyon, çoklu regresyon, çoklu bağlantı kontrolü ve yol analizi adımlarını R’ın kendi mtcars veri seti üzerinde uygulayacağım. Derslerde ve ödevelrde sıklıkla kullandığım için tanıdık bir veri seti olduğu için tercih ettim.

library(dplyr)
library(ggplot2)
library(broom)
library(lavaan)
library(semPlot)
library(knitr)

head(mtcars)
##                    mpg cyl disp  hp drat   wt qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.62 16.5  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.88 17.0  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.32 18.6  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.21 19.4  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.44 17.0  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.46 20.2  1  0    3    1
str(mtcars)
## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...
mtcars$am <- factor(mtcars$am, levels = c(0, 1), labels = c("otomatik", "möanuel"))

Basit Doğrusal Regresyon

Araç ağırlığının (wt) yakıt verimliliği (mpg) üzerindeki etkisini inceliyorum önce:

model_basit <- lm(mpg ~ wt, data = mtcars)
summary(model_basit)
## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.543 -2.365 -0.125  1.410  6.873 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   37.285      1.878   19.86  < 2e-16 ***
## wt            -5.344      0.559   -9.56  1.3e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.05 on 30 degrees of freedom
## Multiple R-squared:  0.753,  Adjusted R-squared:  0.745 
## F-statistic: 91.4 on 1 and 30 DF,  p-value: 1.29e-10
tidy(model_basit) %>% kable(digits = 3)
term estimate std.error statistic p.value
(Intercept) 37.28 1.878 19.86 0
wt -5.34 0.559 -9.56 0 
ggplot(mtcars, aes(x = wt, y = mpg)) +
  geom_point(alpha = 1, size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "pink") +
  theme_minimal() +
  labs(title = "Araç ağırlığı ile yakıt verimliliği",
       x = "Ağırlık", y = "Yakıt Verimliliği")

Çoklu Regresyon: Paralel Eğimler ve Etkileşim

Vites tipinin moderatör etkisini test edeceğim.

SORUM: Ağırlığın mpg üzerindeki etkisi otomatik ve manuel vitesli araçlarda farklılaşıyor mu?

Paralel Eğimler Modeli

data(mtcars)

model_paralel <- lm(mpg ~ wt + am, data = mtcars)

summary(model_paralel)
## 
## Call:
## lm(formula = mpg ~ wt + am, data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.530 -2.362 -0.132  1.403  6.878 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.3216     3.0546   12.22  5.8e-13 ***
## wt           -5.3528     0.7882   -6.79  1.9e-07 ***
## am           -0.0236     1.5456   -0.02     0.99    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.1 on 29 degrees of freedom
## Multiple R-squared:  0.753,  Adjusted R-squared:  0.736 
## F-statistic: 44.2 on 2 and 29 DF,  p-value: 1.58e-09
tidy(model_paralel) %>% kable(digits = 3)
term estimate std.error statistic p.value
(Intercept) 37.322 3.055 12.218 0.000
wt -5.353 0.788 -6.791 0.000
am -0.024 1.546 -0.015 0.988 
ggplot(mtcars, aes(x = wt, y = mpg, color = am)) +
  geom_point(alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal() +
  labs(title = "Paralel eğimler",
       x = "Ağırlık", y = "mpg", color = "Vites")

Etkileşim Modeli

data(mtcars)
model_etkilesim <- lm(mpg ~ wt * am, data = mtcars)
summary(model_etkilesim)
## 
## Call:
## lm(formula = mpg ~ wt * am, data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.600 -1.545 -0.533  0.901  6.091 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   31.416      3.020   10.40  4.0e-11 ***
## wt            -3.786      0.786   -4.82  4.6e-05 ***
## am            14.878      4.264    3.49   0.0016 ** 
## wt:am         -5.298      1.445   -3.67   0.0010 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.59 on 28 degrees of freedom
## Multiple R-squared:  0.833,  Adjusted R-squared:  0.815 
## F-statistic: 46.6 on 3 and 28 DF,  p-value: 5.21e-11
tidy(model_etkilesim) %>% kable(digits = 3)
term estimate std.error statistic p.value
(Intercept) 31.42 3.020 10.40 0.000
wt -3.79 0.786 -4.82 0.000
am 14.88 4.264 3.49 0.002
wt:am -5.30 1.445 -3.67 0.001 
cat("Basit model R²:", round(summary(model_basit)$r.squared, 3), "\n")
## Basit model R²: 0.753
cat("Paralel model R²:", round(summary(model_paralel)$r.squared, 3), "\n")
## Paralel model R²: 0.753
cat("Etkileşim modeli R²:", round(summary(model_etkilesim)$r.squared, 3), "\n")
## Etkileşim modeli R²: 0.833
ggplot(mtcars, aes(x = wt, y = mpg, color = am)) +
  geom_point(alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal() +
  labs(title = "Etkileşim Modeli",
       x = "Ağırlık", y = "mpg", color = "Vites")

Grafik incelendiğinde, manuel vitesli araçlarda ağırlık arttıkça yakıt verimliliğinin daha hızlı düştüğünü anlayabiliyorum.

Çoklu Bağlantı Kontrolü

mtcars_num <- mtcars %>% select(mpg, wt, hp, disp)
cor(mtcars_num) %>% round(2) %>% kable()
mpg wt hp disp
mpg 1.00 -0.87 -0.78 -0.85
wt -0.87 1.00 0.66 0.89
hp -0.78 0.66 1.00 0.79
disp -0.85 0.89 0.79 1.00

wt ve disp arasında çok yüksek korelasyon bulunuyor. Derste öğrendiğimiz gibi, bu iki değişkeni aynı modele koyduğumuzda bakmamız gerekiyor:

cat("wt tek başına:\n")
## wt tek başına:
tidy(lm(mpg ~ wt, data = mtcars))[2, c("term","estimate","p.value")] %>% kable(digits = 3)
term estimate p.value
wt -5.34 0 
cat("\ndisp tek başına:\n")
## 
## disp tek başına:
tidy(lm(mpg ~ disp, data = mtcars))[2, c("term","estimate","p.value")] %>% kable(digits = 3)
term estimate p.value
disp -0.041 0 
cat("\nBirlikte:\n")
## 
## Birlikte:
tidy(lm(mpg ~ wt + disp, data = mtcars))[2:3, c("term","estimate","p.value")] %>% kable(digits = 3)
term estimate p.value
wt -3.351 0.007
disp -0.018 0.064

Tek başlarına anlamlı olan her iki değişken birlikte modele alındığında, disp değişkeninin anlamlılığı düşüyor. Bu, derste Credit veri setinde gördüğümüz çoklu bağlantı probleminin bir örneği aslında. wt ve disp arasındaki korelasyon olduğundan, bu iki değişkeni birlikte kullanmamak gerek.

Yol Analizi

Derste illness verisi üzerinde yaptığımız gibi, mtcars verileriyle bir yol analizi modeli kuracağım. (Aşağıdaki hipotezleri test etmem için chat bana yazdı ödev olarak.)

“CHAT:

  • Beygir gücü (hp) ve ağırlık (wt), çeyrek mil süresini (qsec) etkiler.
  • Beygir gücü (hp) ve ağırlık (wt), yakıt tüketimini (mpg) etkiler.
  • Çeyrek mil süresi (qsec), yakıt tüketimini (mpg) etkiler.”

modelde hp ve wt dışsal değişkenler, qsec arabulucu değişken, mpg ise bağımlı değişken konumundadır.

mtcars_yol <- mtcars %>% select(mpg, hp, wt, qsec) %>% mutate(across(everything(), as.numeric))

Modelin Tanımlanması

yol_model <- '
  qsec ~ hp + wt
  mpg  ~ hp + wt + qsec
  hp ~~ wt
'

yol_fit <- sem(yol_model, data = mtcars_yol)
summary(yol_fit, standardized = TRUE)
## lavaan 0.6-20 ended normally after 33 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        10
## 
##   Number of observations                            32
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   qsec ~                                                                
##     hp               -0.027    0.004   -7.560    0.000   -0.027   -1.048
##     wt                0.942    0.253    3.720    0.000    0.942    0.516
##   mpg ~                                                                 
##     hp               -0.018    0.014   -1.272    0.203   -0.018   -0.203
##     wt               -4.359    0.704   -6.191    0.000   -4.359   -0.708
##     qsec              0.511    0.411    1.243    0.214    0.511    0.151
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   hp ~~                                                                 
##     wt               42.812   13.757    3.112    0.002   42.812    0.659
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .qsec              1.076    0.269    4.000    0.000    1.076    0.348
##    .mpg               5.814    1.454    4.000    0.000    5.814    0.165
##     hp             4553.965 1138.491    4.000    0.000 4553.965    1.000
##     wt                0.927    0.232    4.000    0.000    0.927    1.000

Model-Veri Uyum İndeksleri

fitMeasures(yol_fit, c("chisq", "df", "pvalue",
                        "cfi", "tli",
                        "rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
                        "srmr"))
##          chisq             df         pvalue            cfi            tli 
##              0              0             NA              1              1 
##          rmsea rmsea.ci.lower rmsea.ci.upper           srmr 
##              0              0              0              0

Uyum indekslerini değerlendirelim: Ki-kare testinin anlamlı olmamasını (p > .05), RMSEA’nın .08’den küçük olmasını, CFI’nın .90’dan büyük olmasını ve SRMR’nin .08’den küçük olmasını bekliyoruz.

Artık Kovaryans Matrisi

resid(yol_fit, type = "normalized")
## $type
## [1] "normalized"
## 
## $cov
##      qsec mpg hp wt
## qsec    0          
## mpg     0   0      
## hp      0   0  0   
## wt      0   0  0  0

DERS NOTU: Standartlaştırılmış artık değerlerin mutlak değerce 1.96’dan küçük olması istenir.

Modifikasyon İndeksleri

modindices(yol_fit, sort = TRUE)
## [1] lhs      op       rhs      mi       epc      sepc.lv  sepc.all sepc.nox
## <0 rows> (or 0-length row.names)

Yol Diyagramı

p <- semPaths(yol_fit, whatLabels = "est",
         sizeMan = 10,
         edge.label.cex = 1.15,
         style = "ram",
         layout = "spring",
         nCharNodes = 0, nCharEdges = 0)

p2 <- semptools::mark_sig(p, yol_fit)
plot(p2)

Yol Katsayılarının Yorumu

parameterEstimates(yol_fit, standardized = TRUE) %>%
  filter(op == "~") %>%
  select(lhs, rhs, est, std.all, pvalue) %>%
  mutate(pvalue = ifelse(pvalue < .001, "< .001", round(pvalue, 3))) %>%
  kable(digits = 3, col.names = c("Bağımlı", "Bağımsız", "B", "β", "p"))
Bağımlı Bağımsız B β p
qsec hp -0.027 -1.048 < .001
qsec wt 0.942 0.516 < .001
mpg hp -0.018 -0.203 0.203
mpg wt -4.359 -0.708 < .001
mpg qsec 0.511 0.151 0.214

Standardize edilmiş katsayılar incelendiğinde, hangi değişkenin mpg’yi en güçlü şekilde yordadığı ve qsec’in arabulucu rolü hakkında yorum yapılabilir oldu.

YANSITMA

Ödevleri haftalar öncesinde yapmaya başlamam arağmen ne yapacağımı analamadığım için yapamamıştım daha doğrusu yayıonlama formatına getirmedim ve hoşuma gitmemişti. Hocamızın mesajı sonrası içeriği anladıktan sonra yaptıklarımı aşamalandırmak ve başlıklandırmnak daha kolay oldu.

Geçen haftalarda öğrendiğimiz veri temizliğinin üzerine regresyon analizi ve yol analizi gibi daha ileri düzey teknikleri kullandık. Benim için çoklu bağlantı probleminin nasıl katsayıları tamamen çarpıtabildiğini görmek işlevsel oldu. Çünkü bunu raporlasak ve yorumlasak bile bu kada rüstüne düşmemiştim. Credit veri setinde Income’ın tek başına pozitif olan etkisi, Limit ile birlikte modele konduğunda negatife dönüyordu, bu Simpson Paradoksu’na benzer bir durum. Ders notunda vurgulanan “bağımsız değişkenler arası korelasyon .80’i aştığında aynı modele almayın” kuralını artık unutmam mümkün değil.

Yol analizi konusunda, aslında her bir denklemin ayrı bir çoklu regresyon olduğunu kavramak işleri kolaylaştırdı. Fark, bu regresyonların birbirine bağlı bir sistem olarak ele alınması ve model veri uyumunun bütünsel olarak değerlendirilmesi. Modifikasyon indekslerini kullanarak modeli revize ettik. Hocanın “istatistik sizi oraya yönlendirse bile, kuramsal olarak mantıklı değilse o yolu eklemeyin” uyarısı çok değerli bir ilke ama sosyal bilimlerde sayısal değelerden bağımsız olduğunda da bu ilke manipülasyona biraz açık olduğunu alanyazındaki yayınlardan sıklıkla anlayabiliyorum…

Etkileşim modeli ile paralel eğimler arasındaki fark da ilginçti. Başta anlamamlandıramamıştım. Bir kategorik değişkenin sürekli değişkenle çarpımını modele eklemek, moderasyon etkisini test etmenin en doğrudan yoluymuş. mtcars verisinde vites tipinin ağırlık ile yakıt ilişkisini nasıl modere ettiğini görmek, bu kavramı somutlaştırdı. Son olarak da ödevlerimi zamanında yapacağım, kaçınarak ilerleyemeceğimi kabul ediyorum.