R Dersleri - Ali Onur Gitmez
R Dersi 13
Ali Onur Gitmez
2024-07-15
Logistic Regression
Linear regresyonu dependent variablımızın continious olduğu veya ona yakınsadığı durumlarda kullanmıştık. Mesela bir kişinin maaşı veya yaşını tahmin etmek gibi analizler için linear regresyon en iyi seçenek oluyor. Ancak bunun dışında kalan veriler de var. Dependent variableın 0-1 değerlerini aldığı veya 1-5 arası gibi çok daha sınırlı değerler aldığı durumlarda linear regresyon kullanmak çok uygun olamayabilir. Bu gibi durumlarda logistic regresyon daha güçlü bir seçenek olarak ortaya çıkıyor. Mesela bir kişinin referendumda Evet veya Hayır yönünde nasıl oy kullanacağı, 3-4 sandiv arasından hangisini seçeceğini anlamak için bu modeli kullanırız.
Formülü bu şekilde gösteririz. z, bütün bağımsız değişkenlerimizin ağırlıklı toplamını temsil eder. Z’nin negatif olmasının sebebi, sigmoid fonksiyonunun şeklidir; bu fonksiyon büyüdükçe 1’e yaklaşırken, z küçüldükçe 0’a yaklaşır. Fonksiyonun paydasında “1 artı” terimi vardır çünkü sigmoid fonksiyonunun çıktısı her zaman 0 ile 1 arasında kalmalıdır. Üstte “1” olmasının sebebi ise bu normalizasyonu sağlamaktır; böylece çıktı, asla 1’i geçmez ve her zaman olasılık aralığında kalır. “1”, aynı zamanda sigmoid fonksiyonunun 0 ile 1 arasında düzgün ve sürekli bir geçiş yapmasını sağlar.
\(P(Y=1) = \frac{1}{1 + e^{-z}}\)
Logistic modellerde kullanacağımız bir kavram da odds. Bu bir olayın gerçekleşme ihtimalinin gerçekleşeme ihtimaline oranı olarak bulunur. Mesela bir olayın gerçekleşme ihtimali 0.8 ise gerçekleşmeme ihtimali 0.2dir. Bu bağlamda 0.8/0.2 ile log odds 4 olur. Lojistik regresyon ise bu değerin logunun alındığı log-odds dediğimiz bir kavram üstünden hesaplama yapar.
\(odds=(\frac{P(Y=1)}{1 - P(Y=1)})\)
Son olarak lojistik regresonla alakalı bir kavram da Maximum Likelihood Estimation (MLE). Bu lojistik regresyonun beta değerlerini yani coefficentları tahmin etmek için kullandığı modele verilen isim. Bu yöntem observe ettiğimiz datayı elde etmenin en mümkün olduğu coefficentları seçiyor. Bu yönüyle de linear regresyondan ayrılıyor ancak hala generalized linear models diye bir grubun altında inceliyoruz. İsminde linear var diye şaşırtmasın, model linear değil.
##Assumptions:
Lojistik regresyonda benzer assumptionlar vardır ancak modelin tipine de uyum sağlar. Mesela düz lojistik regresyonda outcome değişkeninin 0 veya 1 değerlerini alması gerekir. Aynı zamanda bağımlı değişkenimizin log oddları ile bağımsız değişkenlerimiz arasında linear bir ilişki olmalıdır. Bu modelin linear olduğu anlamına gelmez ancak log-odd dediğimiz logistic regresyon sonuçları ve bağımsız değişkenimiz arasındaki ilişki linear olmalı.
## Loading required package: carData## MLE of lambda Score Statistic (t) Pr(>|t|)
## 1.427 0.3199 0.7514
##
## iterations = 3Bu sonuçlara bakarak linearity assumptionını reddetmeye yeterli kanıt olmadığımı söyleyebiliriz.
Eğer istersek de modeli oluşturduktan sonra görsel olarak da buna bakabiliriz
model <- glm(am ~ mpg, data = mtcars, family = "binomial")
mtcars$logit <-
log(model$fitted.values / (1 - model$fitted.values))
plot(mtcars$mpg, mtcars$logit)
Potansiyel outlier değerleri görmek için Cook’s distance kullanabiliriz
cooksD <- cooks.distance(model)
plot(cooksD, type = "h", main =
"Cook's Distance for Influential Observations",
ylab = "Cook's Distance")
## Warning: NAs introduced by coercion## mpg cyl disp hp drat wt qsec vs am gear carb logit
## NA NA NA NA NA NA NA NA NA NA NA NA NA
## NA.1 NA NA NA NA NA NA NA NA NA NA NA NA
## NA.2 NA NA NA NA NA NA NA NA NA NA NA NA#Yorumlama
##
## Call:
## glm(formula = am ~ mpg, family = "binomial", data = mtcars)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.6035 2.3514 -2.808 0.00498 **
## mpg 0.3070 0.1148 2.673 0.00751 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 29.675 on 30 degrees of freedom
## AIC: 33.675
##
## Number of Fisher Scoring iterations: 5Bu modelin yorumlamasını linear regresyon gibi yapamıyoruz. Yorumlama basitçe mpg’de olan 1 birimlik değişim aracın manuel vites olma ihtimalini 0.307 log odds artıryor diyebiliriz. Ancak bu kimsenin anlayacağı bir yroumlama değil. Bizim için bir mantığı yok. Yine de bu ilişkinin pozitif olduğunu ve istatistiksel olarak anlamlı olduğunu bilebiliyoruz. Bu da null hipotezi reddetmemize yardımcı oluyor.
Peki düzgün yorumlama için ne yapabiliriz? Sonuçta ilişkinin anlamlı olması kadar gücü de bizim için önemli. Burada belli alterntiflerler karşımızıa çıkıyor.
Burada bir alternatif odds-ratio dediğimiz sistemi kullanmak. Genellikle tıp, biyoloji ve psikoloji alanında bu kullanılır. Buradan çıkan sonuç bize kısaca şunu söylüyor. Mpg’de olan her artışla beraber aracın manüel vitese sahip olma şansı %36 yükseliyor. Burada dikkat çekilecek bir nokta ihimal dememem. Odds dediğimiz olayın artmasından bahsediyor.
## (Intercept) mpg
## 0.001355579 1.359379288Her ne kadar bu güzel bir çözüm gibi gözükse bile önemli eksikleri var. Mesela hastanede hastalari 2 gruba ayirdigimizi ve 1 gruba ilac digerine ise placebo verdigimizi dusunelim. placebo grubunda hastlarin %50si olurken, ilac grubunda hastalarin %25i oluyor olsun. buradan odds ratio hesapladik. ikinci durumda ise yine ayni ornekten bakacak olursak placebo gurubundaki hastlarin %5i, ilac grubundaki hastlarin ise %2.5i oluyor olsun. burada da odds ratio hesaplayalim. simdi ikisinin de odds ratiosu ayni cikiyor. ancak bir yerde arada %25lik cok ciddi bir fark varken diger yerde %2.5lik bir fark var. Mesela bu örnekte eğer aracın manüel olma oddsu %50 olsaydı bunu 1:1 olarak belirtirdik. Budadaki 1.36’lık artış bunu 1.31:1 yapıyor. Yani bize bir bilgi vermesine rağmen tam istediğimiz gibi değil.
Marginal effects dediğimiz bir kavram ise buna alternatif. Diğer sayıların tabanına bağlı olarak bir yüzdelik artış yerine doğrudan bir artıştan bahsediyor ve böylece daha anlaşılır bir sonuç elde etmiş oluyoruz. Bu yöntem özellikle ekonomi alanında kullanılır ve benim de tercih ettiğim yöntem.
library(margins)
library(tidyverse)
marginal_effects <- margins(model)
marginal_effects_df <-
as.data.frame(marginal_effects)
summary(marginal_effects)## factor AME SE z p lower upper
## mpg 0.0465 0.0089 5.2419 0.0000 0.0291 0.0639Bu sonuçlar bize average marginal effects’in 0.047 olduğunu gösteriyor. Bu mpg’nin her 1 birimlik artışının aracın manuel olma ihimalini doğrudan %4.6 artırdığını söylüyor. Bu da odds ratio dedğimiz measure’a göre çok daha intuitive ve kafaya oturan cinsten. Bunu istersek plot ile de gösterebiliriz. Bu plot ise bize mpg’nin farklı noktalarının aracın manuel olma hitimalini ne kadar artıdığını gösteriyor.
ggplot(marginal_effects_df,
aes(x = mpg, y = dydx_mpg)) +
geom_line(color = "blue") +
geom_point(color = "red") +
labs(
title = "Average Marginal Effects of mpg
on the Probability of Manual Transmission",
x = "Miles per Gallon (mpg)",
y = "Marginal Effect on Probability"
) +
theme_minimal()
Bu örnek yine tek değişkenli bir regresyon örneğiydi. Hızlıca bir çok değişkenli örnek inceleyelim.
##
## Call:
## glm(formula = am ~ mpg + hp, family = "binomial", data = mtcars)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -33.60517 15.07672 -2.229 0.0258 *
## mpg 1.25961 0.56747 2.220 0.0264 *
## hp 0.05504 0.02692 2.045 0.0409 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 19.233 on 29 degrees of freedom
## AIC: 25.233
##
## Number of Fisher Scoring iterations: 7Burada yorumlarımız yine öncekinde yaptığımız gibi ancak hem yapabileceğimiz assumption testi hem de Average Marginal Effects değişiyor.
Birden çok bağımsız değişken olan veride multicolliniearity için de bakabiliriz. Burada baktığımız zaman multicolliniearity problemi varmış gibi gözüküyor. Ancak daha önceden de dediğim gibi o kadar da önemli bir problem olmayabilir.
## mpg hp
## 7.925516 7.925516Marginal effects noktasına geldiğimiz zaman iş biraz daha değişiyor.
marginal_effects_2 <- margins(model_2)
marginal_effects_2_df <-
as.data.frame(marginal_effects_2)
marginal_effects_2## Average marginal effects## glm(formula = am ~ mpg + hp, family = "binomial", data = mtcars)## mpg hp
## 0.1238 0.00541ggplot(marginal_effects_2_df,
aes(x = mpg, y = dydx_mpg)) +
geom_line(color = "blue") +
geom_point(color = "red") +
labs(
title = "Average Marginal Effects of mpg
on the Probability of Manual Transmission",
x = "Miles per Gallon (mpg)",
y = "Marginal Effect on Probability"
) +
theme_minimal()
ggplot(marginal_effects_2_df, aes(x = hp, y = dydx_hp)) +
geom_line(color = "green") +
geom_point(color = "orange") +
labs(
title = "Average Marginal Effects of hp
on the Probability of Manual Transmission",
x = "Horsepower (hp)",
y = "Marginal Effect on Probability"
) +
theme_minimal()
Burada Average Marginal Effectleri X değişkeni için Y değişkenini kontrol ederek yapıyoruz. Gördüğümüz üzere iki değişkeni kullandığımız zaman çok daha oynak bir yapıya sahip oldu.
Model karşılaştırma için ise AIC değerini kullanabiriz. Burada daha düşük değer modelin daha iyi olduğunu gösterir. Burada ikinci modelimizin daha iyi olduğunu görüyoruz. Aynı zamanda Null deviance ve Residual deviance arasındaki düşüş farkı da modelin hangisinin daha iyi olduğunu gösteriyor. Bu örnekte ikinci modelde düşüş daha yüksek olduğu için ikinci model daha iyi diyebiliriz.
##
## Call:
## glm(formula = am ~ mpg, family = "binomial", data = mtcars)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.6035 2.3514 -2.808 0.00498 **
## mpg 0.3070 0.1148 2.673 0.00751 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 29.675 on 30 degrees of freedom
## AIC: 33.675
##
## Number of Fisher Scoring iterations: 5##
## Call:
## glm(formula = am ~ mpg + hp, family = "binomial", data = mtcars)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -33.60517 15.07672 -2.229 0.0258 *
## mpg 1.25961 0.56747 2.220 0.0264 *
## hp 0.05504 0.02692 2.045 0.0409 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 43.230 on 31 degrees of freedom
## Residual deviance: 19.233 on 29 degrees of freedom
## AIC: 25.233
##
## Number of Fisher Scoring iterations: 7Modellerin kalitesini ROC Curve dediğimiz bir görsel ve AUC (Area Under the Curve) ile de görebiliriz. Bunlar bize modelimizin farklı sınıflar arasında ayrımları ne kadar başarılı yapabildiğini gösterir.
library(pROC)
pred_probs <- predict(model, type = "response")
roc_curve <- roc(mtcars$am, pred_probs)
## Area under the curve: 0.83AUC değeri 0 ve 1 arasındadır. 1’e yakın olan bir değer modelin tahmin gücünün iyi olduğunu gösterir. 0.5 ise modelin sallamadan farksız tahmin yaptğını gösterir. 0.83 değeri bizim modelimizin iyiye yakın olduğunu gösteriyor.
Multinomial Models
Ancak lojistik regresyon modellerinde sadece binomial değil ikiden fazla seçeneğin olduğu multinomial ve seçeneklerin sıralı olduğu ordinal analizleri de yapabiliriz. Örneklerle inceleyelim. Mesela veriye arabaların geldiği kıtaları koyalım. Asya, Avrupa ve Amerika olmak üzere. Böylece sıralı olmayan ama çoklu değeri olan bir değişken oluşturmuş olduk.
mtcars <- mtcars %>%
mutate(
brand = sapply(strsplit(rownames(mtcars),
" "), `[`, 1),
continent = case_when(
brand %in% c('Mazda', 'Datsun',
'Toyota', 'Honda') ~ 1,
brand %in% c('Hornet', 'Valiant', 'Duster',
'Cadillac', 'Lincoln', 'Chrysler',
'Dodge', 'AMC', 'Camaro', 'Pontiac',
'Ford') ~ 0,
brand %in% c('Merc', 'Fiat', 'Porsche', 'Lotus',
'Ferrari', 'Maserati', 'Volvo') ~ 2
)
)library(nnet)
mtcars$continent <- as.factor(mtcars$continent)
model_3 <- multinom(continent ~ mpg, data = mtcars)## # weights: 9 (4 variable)
## initial value 35.155593
## iter 10 value 24.795460
## final value 24.794146
## converged## Call:
## multinom(formula = continent ~ mpg, data = mtcars)
##
## Coefficients:
## (Intercept) mpg
## 1 -10.64156 0.5159800
## 2 -7.25556 0.4074629
##
## Std. Errors:
## (Intercept) mpg
## 1 3.527891 0.1818517
## 2 2.923730 0.1629697
##
## Residual Deviance: 49.58829
## AIC: 57.58829z_values <- summary(model_3)$coefficients /
summary(model_3)$standard.errors
p_values <- (1 - pnorm(abs(z_values), 0, 1)) * 2
p_values## (Intercept) mpg
## 1 0.002557889 0.004548736
## 2 0.013078985 0.012410995Bu sonuçlara bakarak mpg değeri hem Asya(1) hem de Avrupa(2) için istatistisel olarak anlamlı çıkmış. Bu da bize araçların mpg değerlerinin kıtaya göre değiştiğini gösteriyor. Estimation değerlerine baktığımız zaman da hem Avrupa hem de Asya Amerika’ya göre dah yüksek değerler almış. Bu da bu kıtalardan gelen araçların Amerika’dan gelenlere kıyasla daha iyi yakıt tüketimine sahip olduğunu gösteriyor. Burada dikkat edilmesi gereken bir nokta karşılaştırmalarımızı referans kategoriye göre yapıyoruz ve Avrupa ile Asya arasında yapmıyoruz. Farklı karşılaştırma yapmak istersek referans kategoriyi değiştirmemiz gerekir.
Marginal effectlerin hesaplaması R’da Stata’da olduğu kadar iyi ve güçlü değil. glm dışına çıktığımız zaman yapamıyoruz. O yüzden ne yazık ki bu noktada AME gösteremiyorum. Benzer bir durum assumption checkler için de geçerli. Klasik yöntemlerimizi burada kullanamıyoruz.
Ordinal Model
Son olaraksa ordinal logistic regression modeline bakacağız. Bu her ne kadar multinomial modele benzese bile yorumlanması farklı olduğu ve başka noktalarda da kullanılabileceği için işimize yarayan bir model olacak. Sadece binomial logistic gibi yaygın şekilde kullanılmadığı için hem assumption hem de AME olarak yine zayıf kalıyor. Yine de yorumlamasını öğrenip kullanmak önemli olacak. Bu yöntemi özellikle anket cevapları arasındakı farklara bakmak için kullanabiliriz. Mesela cevaplar katılıyorum, kısmen katılıyorum vb gibi gidiyorsa bu modeli uygulayabiliriz.
Önce değişkeni bi düzenleyelim.
library(MASS)
library(ordinal)
data(airquality)
airquality <- na.omit(airquality)
airquality$Ozone_cat <- cut(airquality$Ozone,
breaks =
quantile(airquality$Ozone,
probs = c(0, 0.33,
0.67, 1)),
labels = c("Low", "Medium",
"High"),
include.lowest = TRUE,
ordered_result = TRUE)##
## Re-fitting to get Hessian## Call:
## polr(formula = Ozone_cat ~ Wind + Temp, data = airquality, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## Wind -0.3529 0.08454 -4.175
## Temp 0.2182 0.03693 5.907
##
## Intercepts:
## Value Std. Error t value
## Low|Medium 12.1275 2.8612 4.2386
## Medium|High 14.9195 3.0371 4.9124
##
## Residual Deviance: 142.1852
## AIC: 150.1852Sonuçlara bakacak olursak rüzgarın ve sıcaklığın da etkisi istatistiksel olarak anlamlı gözüküyor. Yani her iki değer için de null hipotezi reddetmeyi başarıyoruz. Ancak rüzgarın etkisi negatif iken sıcaklığın etkisi pozitif. Şimdi bunların etkisine daha yakından bakalım:
## Wind Temp
## 0.7026448 1.2437995Rüzgar için odds ratio 1’in altında. Bu da değerin negatif olduğunun göstergesi. Sonuç bize rüzgardaki her birimlik artışın ozone kategorisinde bir üst kategoride bulunma oddunu %29.7 azaltıyor. Aynı şekilde sıcaklık 1 üstü ve bu da sıcaklktaki her birimlik artışın ozone kategorisinde bir üst kategoride bulunma oddunu %24.4 artıdığını gösteriyor.
Bu modelin önemli bir assumptionı da var. Her değer arasındaki ilişkinin eşit olduğunu varsayıyor. Bunu test etmek için:
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 9.72 2 0.01
## Wind 6.12 1 0.01
## Temp 8.61 1 0
## --------------------------------------------
##
## H0: Parallel Regression Assumption holdsSonuçlara bakacak olursak da assumption’ın hold etmediğini görüyoruz. H0 0.01 seviyesinde reddedilmiş görünüyor. Bu da bize modelimizin data için uygun olmadığının bir göstergesi. Bu durumlarda alternatif modeller kullanmak daha doğru olacaktır.
Linear Proability Model binomial logistic regressiona bir alternatif olarak öne çıkıyor ve bazı ekonometristler arasında çok popüler oldu. Bu model temel olarak logistic regresyonda kullandığımız 0 ve 1 bağımlı değişken sistemini linear regresyona uyarlama oluyor. Popüler olmasının temel sebebi doğrudan yorumlanabilmesi. Yani independent variablellarda olan 1 birimlik değişimin dependent variable’ın 0’dan 1’e çıkma ihtimalini % kaç artırdığına bakar. Özellikle ekonomistler arasında yayıldı ancak eksikleri de var. Mesela 0 altı 1 üstü gibi gerçek olması mümkün olmayan olasılık sonuçları da verebiliyor.
\(P(Y = 1 \mid X) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \varepsilon\)
##
## Call:
## lm(formula = am ~ mpg + hp + wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5189 -0.2470 -0.1013 0.2425 0.5876
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.195757 0.916432 0.214 0.8324
## mpg 0.036309 0.023984 1.514 0.1413
## hp 0.003892 0.001393 2.794 0.0093 **
## wt -0.338757 0.123810 -2.736 0.0107 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.335 on 28 degrees of freedom
## Multiple R-squared: 0.593, Adjusted R-squared: 0.5494
## F-statistic: 13.6 on 3 and 28 DF, p-value: 1.168e-05Bu sonuçlar her zaman Average Marginal Effect ile uyuşmaz çünkü LPM modeli marginal effecting predictorun her değeri için aynı olduğunu assume eder. Bu durum da 0-1 dışına çıkan değerlerin sebebidir. Bazı işler için kolay gibi gelse de bizim 0-1 değerlere logistic regresyon kullanma sebebimiz zaten bu tarz sorunların üstesinden gelmek. Bu durumda da LPM kullanmak garip oluyor elbet. Ancak basit ve anlaışır modeli kullanıp bu riski almayı tercih edenler de oluyor elbet. Tercih kişiye kalmış ama logistic regresyon bence daha iyi bir tercih.