SPSS Örneği Lojistik Regresyon Analizi

# Veri setini oluşturma
veri <- data.frame(
  X = c(1,2,3,4,4,5,6,7,8,8,9),
  pass = c(0,0,0,0,0,0,1,1,0,1,1)
)

# Model
model <- glm(pass ~ X, data = veri, family = binomial)

# Sonuç
summary(model)
## 
## Call:
## glm(formula = pass ~ X, family = binomial, data = veri)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -6.3475     3.6336  -1.747   0.0807 .
## X             0.9843     0.5547   1.774   0.0760 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 14.4206  on 10  degrees of freedom
## Residual deviance:  7.6638  on  9  degrees of freedom
## AIC: 11.664
## 
## Number of Fisher Scoring iterations: 6
#odds
exp(coef(model))
## (Intercept)           X 
## 0.001751178 2.675995938

logit(p)=−6.3475+0.9843X denklemi elde edilir. X’teki 1 birimlik artış, başarılı olma durumunun log-odds’unu 0.9843 birim artırır. Bu artış pozitif, ilişki yönü olumlu. Bu artış, odds oranı cinsinden değerlendirildiğinde yaklaşık 2.68 katlık bir artışa karşılık gelmektedir.

# Her Gözlem için Tahmin Olasılıkları
veri$tahmin_olasilik <- predict(model, type = "response")
veri
##    X pass tahmin_olasilik
## 1  1    0     0.004664287
## 2  2    0     0.012384798
## 3  3    0     0.032467741
## 4  4    0     0.082399697
## 5  4    0     0.082399697
## 6  5    0     0.193744796
## 7  6    1     0.391374821
## 8  7    1     0.632459980
## 9  8    0     0.821582546
## 10 8    1     0.821582546
## 11 9    1     0.924939091

Modelden elde edilen tahmin olasılıkları incelendiğinde, X değişkeni arttıkça bireylerin başarılı olma olasılığının artış gösterdiği görülmektedir. Ancak bu artış doğrusal olmayıp lojistik eğri formundadır. Aynı X değerine sahip bireyler arasında (x=8 örneğinde)farklı sonuçların gözlenmesi, ilgili değişkenin tek başına başarıyı belirlemede yeterli olmadığını göstermektedir.Bu durum, başarı olgusunun olasılıksal bir yapıya sahip olduğunu ve lojistik regresyon modelinin bu belirsizliği modellemek amacıyla kullanıldığını ortaya koymaktadır.

#Sınıflama
veri$tahmin_sinif <- ifelse(veri$tahmin_olasilik >= 0.5, 1, 0)

table(Gercek = veri$pass, Tahmin = veri$tahmin_sinif)
##       Tahmin
## Gercek 0 1
##      0 6 1
##      1 1 3

Doğru tahminler 6 → 0 olup doğru şekilde 0 tahmin edilmiş (True Negative) 3 → 1 olup doğru şekilde 1 tahmin edilmiş (True Positive) Toplam doğru = 9

Hatalı tahminler 1 → 0 iken 1 tahmin edilmiş (False Positive) 1 → 1 iken 0 tahmin edilmiş (False Negative) Toplam hata = 2

Model Başarısı 9/11=0,818 %81,8 doğruluk

#Grafik
plot(veri$X, veri$pass, pch = 16)

x_seq <- seq(min(veri$X), max(veri$X), length.out = 100)
y_prob <- predict(model, newdata = data.frame(X = x_seq), type = "response")

lines(x_seq, y_prob, lwd = 2)

Lojistik regresyon eğrisi incelendiğinde, X değişkeni arttıkça bireylerin başarılı olma olasılığının S-şeklinde bir artış gösterdiği görülmektedir. Düşük X değerlerinde başarı olasılığı oldukça düşükken, orta değerlerde hızlı bir artış gözlenmekte ve yüksek X değerlerinde olasılık 1’e yaklaşmaktadır. Bununla birlikte, aynı X değerine sahip bireyler arasında farklı sonuçların gözlenmesi, modelin olasılıksal yapısından kaynaklanmaktadır.

Kalp Nakli Örneği

library(ggplot2)
library(openintro)
## Zorunlu paket yükleniyor: airports
## Zorunlu paket yükleniyor: cherryblossom
## Zorunlu paket yükleniyor: usdata
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
data("heart_transplant")
View(heart_transplant)

str(heart_transplant)
## tibble [103 × 8] (S3: tbl_df/tbl/data.frame)
##  $ id        : int [1:103] 15 43 61 75 6 42 54 38 85 2 ...
##  $ acceptyear: int [1:103] 68 70 71 72 68 70 71 70 73 68 ...
##  $ age       : int [1:103] 53 43 52 52 54 36 47 41 47 51 ...
##  $ survived  : Factor w/ 2 levels "alive","dead": 2 2 2 2 2 2 2 2 2 2 ...
##  $ survtime  : int [1:103] 1 2 2 2 3 3 3 5 5 6 ...
##  $ prior     : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ transplant: Factor w/ 2 levels "control","treatment": 1 1 1 1 1 1 1 2 1 1 ...
##  $ wait      : int [1:103] NA NA NA NA NA NA NA 5 NA NA ...
head(heart_transplant)
## # A tibble: 6 × 8
##      id acceptyear   age survived survtime prior transplant  wait
##   <int>      <int> <int> <fct>       <int> <fct> <fct>      <int>
## 1    15         68    53 dead            1 no    control       NA
## 2    43         70    43 dead            2 no    control       NA
## 3    61         71    52 dead            2 no    control       NA
## 4    75         72    52 dead            2 no    control       NA
## 5     6         68    54 dead            3 no    control       NA
## 6    42         70    36 dead            3 no    control       NA
# Yeni değişken oluşturma
heart_transplant <- heart_transplant %>%
  mutate(is_alive = ifelse(survived == "alive", 1, 0))

# Lojistik Regresyon
model <- glm(is_alive ~ age,
             data = heart_transplant,
             family = binomial)

summary(model)
## 
## Call:
## glm(formula = is_alive ~ age, family = binomial, data = heart_transplant)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  1.56438    1.01521   1.541   0.1233  
## age         -0.05847    0.02306  -2.535   0.0112 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 120.53  on 102  degrees of freedom
## Residual deviance: 113.67  on 101  degrees of freedom
## AIC: 117.67
## 
## Number of Fisher Scoring iterations: 4
# Odds Ratio
exp(coef(model))
## (Intercept)         age 
##   4.7797050   0.9432099
#Grafik
ggplot(heart_transplant, aes(x = age, y = is_alive)) +
  geom_jitter(width = 0, height = 0.05, alpha = 0.5) +
  stat_smooth(
    method = "glm",
    method.args = list(family = "binomial"),
    se = FALSE,
    color = "red"
  )
## `geom_smooth()` using formula = 'y ~ x'

Analizde bağımlı değişken “survived” değişkeni ikili (binary) formata dönüştürülmüş ve “alive” durumu 1, “dead” durumu 0 olarak kodlanarak “is_alive” değişkeni oluşturulmuştur.

Lojistik regresyon modeli şu şekilde elde edilmiştir:

logit(p)=1.564−0.058⋅age

Yaş değişkeninin katsayısı negatiftir (B = -0.058) ve istatistiksel olarak anlamlıdır (p = 0.011 < 0.05). Yaşta 1 birim artış, bireyin hayatta kalma durumuna ilişkin log-odds değerini 0.058 birim azaltır.

Yaştaki bir birimlik artış, bireyin hayatta kalma odds’unu yaklaşık 0.943 katına düşürmektedir.

(1−0.9432099)×100=5.67901

Yaştaki 1 birim artış, hayatta kalma odds’unu yaklaşık %5.7 azaltır.

Olasılık Ölçeğinde Grafik

#Tahmin Olasılıkları
library(broom)

heart_transplant_plus <- model %>%
  augment(type.predict = "response") %>%
  mutate(y_hat = .fitted)

ggplot(heart_transplant_plus, aes(x = age, y = y_hat)) +
  geom_point(alpha = 0.5) +
  geom_line(color = "red") +
  scale_y_continuous("Probability of being alive", limits = c(0, 1)) +
  labs(x = "Age")

Odds Ölçeğinde Grafik

heart_transplant_plus <- heart_transplant_plus %>%
  mutate(odds_hat = y_hat / (1 - y_hat))

ggplot(heart_transplant_plus, aes(x = age, y = odds_hat)) +
  geom_point(alpha = 0.5) +
  geom_line(color = "red") +
  labs(x = "Age", y = "Odds of being alive")

Log-Odds (Logit) Ölçeğinde Grafik

# log-odds hesaplama
heart_transplant_plus <- heart_transplant_plus %>%
  mutate(log_odds = log(odds_hat))

# log-odds grafiği
ggplot(heart_transplant_plus, aes(x = age, y = log_odds)) +
  geom_point(alpha = 0.5) +
  geom_line() +
  labs(x = "age", y = "Log(odds) of being alive")

Lojistik regresyon modelinden elde edilen tahminler olasılık, odds ve log-odds düzeylerinde incelenmiştir. İlk grafikte dikey eksende “hayatta kalma olasılığı”, ikinci grafikte “hayatta kalma odds’u”, üçüncü grafikte ise “log-odds” değerleri yer almaktadır. Olasılık grafiği, yaş arttıkça bireylerin hayatta kalma olasılığının azaldığını ve bu değerlerin 0 ile 1 arasında sınırlandığını göstermektedir. Örneğin genç yaşlarda hayatta kalma olasılığı yaklaşık 0.75 civarında iken, yaş ilerledikçe bu değer kademeli olarak düşmekte ve ileri yaşlarda 0.10 civarına kadar inmektedir. Odds, olasılığın doğrudan kendisi değil, olasılığın gerçekleşmeme durumuna oranıdır ve grafikte, yaşın düşük olduğu durumlarda hayatta kalma odds’unun daha yüksek, yaş arttıkça ise belirgin şekilde azaldığı görülmektedir. Örneğin olasılık 0.75 olduğunda buna karşılık gelen odds değeri 3’tür; bu da olayın gerçekleşme olasılığının gerçekleşmeme olasılığının 3 katı olduğunu ifade eder. Bu nedenle ikinci grafikte genç yaşlarda değerlerin 1’in üzerinde olması mümkündür. Yaş arttıkça odds değerinin hızlı biçimde düştüğü, ileri yaşlarda ise 1’in altına indiği ve bu durumun hayatta kalmanın hayatta kalmamaya göre daha düşük ihtimalli hale geldiğini gösterdiği görülmektedir. Log-odds grafiği incelendiğinde ise yaş ile hayatta kalma arasındaki ilişkinin doğrusal bir yapı sergilediği dikkat çekmektedir.

Kalp Nakli Yapılan Hastaların Sağkalım Durumlarının Kontrol Grubu ile Karşılaştırılması

Bu bölümde, kalp nakli yapılan hastaların sağkalım durumları yalnızca yaş değişkeni ile değil, aynı zamanda nakil durumu (transplant) dikkate alınarak incelenmiştir. Önceki analizlerde yaşın hayatta kalma üzerindeki etkisi ortaya konmuş; ancak çalışmanın temel amacı, kalp naklinin sağkalım üzerindeki etkisini değerlendirmektir. Bu doğrultuda lojistik regresyon modeline transplant değişkeni dahil edilerek, yaşın etkisi kontrol altında tutulmuş ve kalp naklinin hayatta kalma üzerindeki katkısı incelenmiştir. Böylece, nakil yapılan hastaların kontrol grubuna kıyasla daha uzun yaşayıp yaşamadığı istatistiksel olarak değerlendirilmiştir.

library(dplyr)


# model
model2 <- glm(is_alive ~ age + transplant,
              data = heart_transplant,
              family = binomial)

summary(model2)
## 
## Call:
## glm(formula = is_alive ~ age + transplant, family = binomial, 
##     data = heart_transplant)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)   
## (Intercept)          0.97311    1.07683   0.904  0.36616   
## age                 -0.07632    0.02551  -2.992  0.00277 **
## transplanttreatment  1.82316    0.66799   2.729  0.00635 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 120.53  on 102  degrees of freedom
## Residual deviance: 103.90  on 100  degrees of freedom
## AIC: 109.9
## 
## Number of Fisher Scoring iterations: 5
exp(coef(model2))
##         (Intercept)                 age transplanttreatment 
##           2.6461676           0.9265153           6.1914009

Model denklemi:

logit(p)=0.973−0.076⋅age+1.82⋅transplant

Yaş değişkeninin katsayısı negatiftir. Bu da yaş arttıkça hayatta kalma ihtimalinin azaldığını gösterir.

Odds ratio’ya baktığımızda yaş için değer yaklaşık 0.93 çıkmıştır. Bu, yaş 1 yıl arttığında hayatta kalma odds’unun yaklaşık %7 civarında azaldığı anlamına gelir.

Transplant değişkeninin katsayısı pozitiftir. Bu da kalp nakli yapılan hastaların hayatta kalma ihtimalinin arttığını gösterir.

Transplant için odds ratio yaklaşık 6.2’dir. Yani kalp nakli yapılan hastaların hayatta kalma odds’u, nakil yapılmayanlara göre yaklaşık 6 kat daha fazladır.

Her iki değişken de istatistiksel olarak anlamlıdır. Yani hem yaşın hem de nakil durumunun hayatta kalma üzerinde etkisi vardır

augment(model2, type.predict = "response")
## # A tibble: 103 × 9
##    is_alive   age transplant .fitted .resid   .hat .sigma  .cooksd .std.resid
##       <dbl> <int> <fct>        <dbl>  <dbl>  <dbl>  <dbl>    <dbl>      <dbl>
##  1        0    53 control     0.0443 -0.301 0.0219   1.02 0.000354     -0.304
##  2        0    43 control     0.0904 -0.435 0.0295   1.02 0.00104      -0.442
##  3        0    52 control     0.0476 -0.312 0.0225   1.02 0.000392     -0.316
##  4        0    52 control     0.0476 -0.312 0.0225   1.02 0.000392     -0.316
##  5        0    54 control     0.0412 -0.290 0.0213   1.02 0.000319     -0.293
##  6        0    36 control     0.145  -0.560 0.0403   1.02 0.00248      -0.571
##  7        0    47 control     0.0682 -0.376 0.0259   1.02 0.000666     -0.381
##  8        0    41 treatment   0.418  -1.04  0.0192   1.02 0.00477      -1.05 
##  9        0    47 control     0.0682 -0.376 0.0259   1.02 0.000666     -0.381
## 10        0    51 control     0.0512 -0.324 0.0231   1.02 0.000436     -0.328
## # ℹ 93 more rows

augment() fonksiyonu model nesnesinden yeni bir veri çerçevesi üretir ve tahmin edilen değerleri ekler.

type.predict = “response” yazıldığında .fitted değeri olasılık olarak gelir.

Örneğin, 53 yaşında ve nakil yapılmamış bir bireyin hayatta kalma olasılığı yaklaşık % 4.4 çıkmaktadır.

Dick Cheney için modele dayalı tahmin

ABD eski başkan yardımcısı Dick Cheney, 2012 yılında 71 yaşında kalp nakli olmuştur.

Peki modelimiz onun 5 yıllık sağkalım olasılığı için ne tahmin eder?

cheney <- data.frame(age = 71, transplant = "treatment")

augment(model2, newdata = cheney, type.predict = "response")
## # A tibble: 1 × 3
##     age transplant .fitted
##   <dbl> <chr>        <dbl>
## 1    71 treatment   0.0677

Model tahminine göre, 71 yaşında ve kalp nakli yapılmış bir bireyin 5 yıllık sağkalım olasılığı yaklaşık %6,7’dir.

mod_plus <- augment(model2, type.predict = "response") %>%
  mutate(alive_hat = round(.fitted))

mod_plus %>%
  select(is_alive, alive_hat) %>%
  table()
##         alive_hat
## is_alive  0  1
##        0 71  4
##        1 20  8

71 → ölü olup doğru şekilde 0 tahmin edilmiş (True Negative) 8 → hayatta olup doğru şekilde 1 tahmin edilmiş (True Positive)

Toplam doğru = 79

4 → ölü ama model “yaşar” demiş (False Positive) 20 → hayatta ama model “ölür” demiş (False Negative)

Toplam hata = 24

modelin 103 gözlemden 79’unu doğru sınıflandırdığı ve genel doğruluk oranının yaklaşık %77 olduğu görülmektedir.

Model, ölen bireyleri doğru tahmin etmede başarılı olmakla birlikte, hayatta kalan bireyleri tahmin etmede daha zayıf performans göstermektedir.

Öğrenme Günlüğü

Lojistik regresyon keyifli bir konuydu. ilk fırsatta çalışmaya başladım. Denklem ve katsayıların yorumlanması ile ilgili notlarımı aldım. Buradaki açıklamaları da pekişmesi için detaylı olarak yazmaya çalıştım. Verimli bir çalışma olduğunu düşünüyorum.