5. HAFTANIN TEKRARI

Bağımlı değişkenin iki kategorili olduğu durumlar için lojistik regresyon konusunu işledik. Doğrusal regresyonun neden bu tür verilerde yetersiz kaldığını, olasılık-odds-logit dönüşümlerini, R’da glm() ile model kurmayı, odds oranı yorumlamayı ve model uyum değerlendirmesini öğrendik.

Doğrusal Regresyonun Yetersizliği

İkili bağımlı değişkene (0-1) doğrusal regresyon uygulandığında iki temel sorun ortaya çıkıyor:

  1. tahmin edilen olasılıklar [0,1] aralığının dışına çıkabiliyor,

  2. hata terimleri normal dağılmıyor ve varyanslar heterojen oluyor.

Cox ve Wermuth (1992) makalesinde de belirtildiği gibi, olasılıklar %20-%80 aralığında kaldığı sürece doğrusal regresyon idare eder ama uç olasılıklarda ciddi sorunlar yaratıyor. Bu nedenle ikili bağımlı değişkenlerde lojistik regresyon tercih edilmeliymiş.

Olasılık, Odds ve Log-Odds (Logit)

Derste en çok üstünde durduğumuz kısımdı. Lojistik regresyonun mantığını kavramak için üç ölçek arasındaki dönüşüm söz konusu:

  • Olasılık (p): [0,1] arasında, sezgisel ama doğrusal değil.

  • Odds = p/(1-p): [0, +∞) arasında. Bir olayın gerçekleşme şansının gerçekleşmeme şansına oranı. Ama bu sadece pozitif ama bize -/+ lazım o yüzden log alacağız.

  • Log-odds (logit) = ln(odds): (-∞, +∞) arasında. Bu dönüşümle model doğrusal hale geliyor ve bağımsız değişkenlerle doğrusal denklem kurulabiliyor.

Derste Pass verisi üzerinde bu üç ölçeği grafik üzerinde ayrı ayrı inceledik. Olasılık ölçeğinde S şeklinde (sigmoid) eğri, odds ölçeğinde üstel artış, log-odds ölçeğinde ise doğrusal bir ilişki gördük.

glm() ile Lojistik Regresyon

R’da lojistik regresyon glm(formula, family=binomial,data) ile kuruluyor. family =binomial argümanı varsayılan olarak logit bağlantı fonksiyonunu kullanıyor. Katsayılar log-odds cinsinden verildiği için yorumlama kolaylığı açısından exp(coef()) ile odds oranına dönüştürülüyor. Bu bir aslında gerekçesi olmayan sadece yorumlanması için yapılan bir işlemmiş.

library(haven)
pass <- read_sav("Pass.sav")
mod_pass <- glm(pass ~ X, family = binomial, data = pass)
summary(mod_pass)
exp(coef(mod_pass))

pass verisinde OR = 2.676 bulundu, yani ön test puanındaki 1 birimlik artış geçme oddsunu yaklaşık 2.68 kat artırıyor.

Model Uyumu ve Testler

Model Uyumu ve diğer Testler

  • Null vs Residual Deviance: Null deviance sadece sabit terimli model, residual deviance yordayıcı eklenmiş model. Fark ki-kare dağılımına uyuyor; fark kritik değerden büyükse model anlamlı.

  • Nagelkerke R²: [0,1] arasında; doğrusal regresyondaki R²’nin yaklaşık bir karşılığı.

  • Hosmer-Lemeshow testi: H0 model veriye uygundur. p > .05 isteniyor.

  • Wald testi: Bu test ile her bir katsayının anlamlılığı ayrıca test edilebiliyor. Küçük örneklemlerde z testinden daha güvenilir.

  • Sınıflandırma tablosu: Tahmin edilen olasılıklar .50 eşiğiyle 0/1’e dönüştürülüp gerçek değerlerle karşılaştırılıyor. Doğru sınıflandırma yüzdesi hesaplanıyor.

FARKLI VERİ SETİ İLE UYGULAMA: birthwt

Chatten destek aldığım kısım veri seti seçmek içindir ama uygulamaları kendim yaptım: Derste öğrendiğimiz lojistik regresyon adımlarını MASS paketindeki birthwt veri setiyle uygulayacağım. Bu veri seti 189 anneye ait bilgileri içeriyor. Bağımlı değişken düşük doğum ağırlığı durumu (low: 0 = normal, 1 = düşük ağırlık), yordayıcılar ise annenin yaşı (age), sigara içme durumu (smoke: 0/1) ve hamilelik öncesi kilosu (lwt, pound cinsinden).

library(MASS)
library(dplyr)
library(ggplot2)
library(broom)
library(knitr)
library(ResourceSelection)
library(DescTools)
library(aod)
data(birthwt)
head(birthwt)
##    low age lwt race smoke ptl ht ui ftv  bwt
## 85   0  19 182    2     0   0  0  1   0 2523
## 86   0  33 155    3     0   0  0  0   3 2551
## 87   0  20 105    1     1   0  0  0   1 2557
## 88   0  21 108    1     1   0  0  1   2 2594
## 89   0  18 107    1     1   0  0  1   0 2600
## 91   0  21 124    3     0   0  0  0   0 2622
str(birthwt)
## 'data.frame':    189 obs. of  10 variables:
##  $ low  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
##  $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
##  $ race : int  2 3 1 1 1 3 1 3 1 1 ...
##  $ smoke: int  0 0 1 1 1 0 0 0 1 1 ...
##  $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ht   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ui   : int  1 0 0 1 1 0 0 0 0 0 ...
##  $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
##  $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...

Veri setinde 189 gözlem ve 10 değişken var.

birthwt$smoke <- factor(birthwt$smoke, levels = c(0, 1), labels = c("Hayır", "Evet"))
table(birthwt$low)
## 
##   0   1 
## 130  59

Düşük doğum ağırlığı olan bebek sayısı 59, normal olan 130.

birthwt %>%
  group_by(smoke) %>%
  summarise(n = n(),
            dusuk = sum(low),
            oran = round(mean(low), 2)) %>%
  kable()
smoke n dusuk oran
Hayır 115 29 0.25
Evet 74 30 0.41

Doğrusal vs Lojistik Regresyon Karşılaştırması

Derste pass verisinde gördüğümüz gibi, doğrusal regresyon ikili bağımlı değişkende [0,1] dışına çıkabiliyor. Bunu burada anne yaşı üzerinden gösteriyorum

ggplot(birthwt, aes(x = age, y = low)) +
  geom_jitter(width = 0, height = 0.05, alpha = 0.4) +
  geom_smooth(method = "lm", se = FALSE, color = "lightblue") +
  geom_smooth(method = "glm", se = FALSE, color = "lightpink",
              method.args = list(family = "binomial")) +
  theme_minimal() +
  labs(title = "Doğrusal (mavi) vs LR (pembe)",
       x = "Anne yaşı", y = "Düşük doğum ağırlığı")

Mavi doğrusal çizgi belirli yaş değerlerinde 0’ın altına çıkabilirken, pembe lojistik eğri her zaman [0,1] aralığında kalıyor.

Lojistik Regresyon Modeli

model_lr <- glm(low ~ age + smoke + lwt, family = binomial, data = birthwt)
summary(model_lr)
## 
## Call:
## glm(formula = low ~ age + smoke + lwt, family = binomial, data = birthwt)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  1.36823    1.01426    1.35    0.177  
## age         -0.03899    0.03273   -1.19    0.233  
## smokeEvet    0.67076    0.32588    2.06    0.040 *
## lwt         -0.01214    0.00613   -1.98    0.048 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 222.88  on 185  degrees of freedom
## AIC: 230.9
## 
## Number of Fisher Scoring iterations: 4

Odds

exp(coef(model_lr)) %>% round(3)
## (Intercept)         age   smokeEvet         lwt 
##       3.928       0.962       1.956       0.988
tidy(model_lr, conf.int = TRUE, exponentiate = TRUE) %>%
  mutate(across(where(is.numeric), ~round(., 3))) %>%
  kable(col.names = c("Değişken", "OR", "SH", "z", "p", "GA Alt", "GA Üst"))
Değişken OR SH z p GA Alt GA Üst
(Intercept) 3.928 1.014 1.35 0.177 0.565 30.690
age 0.962 0.033 -1.19 0.233 0.900 1.024
smokeEvet 1.956 0.326 2.06 0.040 1.033 3.720
lwt 0.988 0.006 -1.98 0.048 0.975 0.999

OR > 1 ise risk artırıcı, OR < 1 ise koruyucu faktör. Sigara içmenin düşük doğum ağırlığı oddsunu ne kadar artırdığına ve anne kilosunun koruyucu etkisine dikkat edilmeli.

Model Uyumu

Nagelkerke R²

PseudoR2(model_lr, which = "Nagelkerke")
## Nagelkerke 
##     0.0851

Hosmer-Lemeshow Testi

y_tahmin <- fitted(model_lr)
y_gozlem <- birthwt$low
hl <- hoslem.test(y_gozlem, y_tahmin, g = 10)
hl
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  y_gozlem, y_tahmin
## X-squared = 7, df = 8, p-value = 0.5

p > .05 ise model veriye iyi uyum sağlıyor

Wald Testleri

wald.test(b = coef(model_lr), Sigma = vcov(model_lr), Terms = 3)
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 4.2, df = 1, P(> X2) = 0.04
wald.test(b = coef(model_lr), Sigma = vcov(model_lr), Terms = 4)
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 3.9, df = 1, P(> X2) = 0.048

Sınıflandırma tablosu

birthwt$tahmin_prob <- predict(model_lr, type = "response")
birthwt$tahmin_sinif <- ifelse(birthwt$tahmin_prob >= 0.5, 1, 0)

sinif_tablo <- table(Gercek = birthwt$low,
                      Tahmin = birthwt$tahmin_sinif)
sinif_tablo
##       Tahmin
## Gercek   0   1
##      0 123   7
##      1  53   6
dogruluk <- sum(diag(sinif_tablo)) / sum(sinif_tablo)
cat("Doğru sınıflandırma yüzdesi:", round(dogruluk * 100, 1), "%\n")
## Doğru sınıflandırma yüzdesi: 68.3 %

Üç Ölçekte Görselleştirme

Derste pass verisinde olasılık, odds ve log-odds grafiklerini çizdik. Burada da anne kilosu üzerinden benzer görselleri deniyoprum

model_basit <- glm(low ~ lwt, family = binomial, data = birthwt)

pred_veri <- data.frame(lwt = seq(80, 250, by = 1))
pred_veri$prob <- predict(model_basit, newdata = pred_veri, type = "response")
pred_veri$odds <- pred_veri$prob / (1 - pred_veri$prob)
pred_veri$logodds <- log(pred_veri$odds)
ggplot(pred_veri, aes(x = lwt, y = prob)) +
  geom_line(color = "lightblue", linewidth = 1) +
  ylim(0, 1) +
  theme_minimal() +
  labs(title = "Olasılık ölçeği", x = "Annekilosu", y = "P(Düşük doğum ağırlığı)")

ggplot(pred_veri, aes(x = lwt, y = odds)) +
  geom_line(color = "lightpink", linewidth = 1) +
  theme_minimal() +
  labs(title = "Odds", x = "Anne kilosu", y = "Odds")

ggplot(pred_veri, aes(x = lwt, y = logodds)) +
  geom_line(color = "orange", linewidth = 1) +
  theme_minimal() +
  labs(title = "Logit",
       x = "Anne Kilosu", y = "Log-Odds")

Log-odds ölçeğinde ilişkinin doğrusallaştığı açıkça gördüm

YANSITMA

LG neden logit dönüşümüne ihtiyaç duyduğunu anladım. İlk bakışta olasılıktan oddsa sonda da logodds a karmaşık görünüyor ama aslında bir sebebi varmış. olasılık [0,1]’e sıkışmış, odds üst sınırı kaldırıyor [0,+∞), log-odds ise alt sınırı da kaldırarak (-∞,+∞) aralığına taşıyıp doğrusal modelleme imkanı sağlıyor. Odds yorumlama konusunda da bir şeyler öğrendim. OR= 2.68 demek olasılık 2.68 kat artıyor değil, odds 2.68 kat artıyor demek. Bu ayrım başta kafa karıştırıcı olsa da, derste yapılan örneklerle somutlaştıu. Model uyumu değerlendirme kısmında Hosmer-Lemeshow (Her seferinde bu testin nasıl yazıldığına bakmam gerekiyor…) testinin mantığını kavramakta zorlandım. yol analizindeki kikare testinde olduğu gibi, burada da H0’ın reddedilmemesi isteniyor, yani p >.05 olması iyi bir şey. Bu, daha önce öğrendiğimiz YEM uyum testleriyle paralel bir mantıktı aslında. Sınıflandırma tablosu ve doğruluk yüzdesi hesaplama kısmı pratik açıdan en anlaşılır kısmıydı kesinlikle. Modelin genel doğruluğunun yanıltıcı olabileceğini (örn.: dengesiz sınıflarda çoğunluk sınıfını sürekli tahmin ederek yüksek doğruluk elde edilebileceğini) unutmamak lazımmış. Şimdiye kadar sürekli bağımlı değişkenlerle çalışırken kategorik bağımlı değişkenlerle de analiz yapabilmeye başladık. Bu, ölçme değerlendirme alanında sıkça karşılaşılan bir durum (geçti/kaldı, doğru/yanlış gibi) olduğundan benim için doğrudan uygulanabilir bir konu.