1. Titanic utasai és a Logisztikus regresszió alapfeladata

Töltsük be a Moodle-n található titanic.csv adattáblát, ami , ami a Titanic 891 véletlenszerűen kiválasztott utasáról tartalmaz tartalmaz 12 változót:

  • PassengerId: Az utas azonosító sorszáma
  • Survived: Túlélte-e az utas a jégheggyel történő találkozást? (1: igen 0: nem)
  • Pclass: A kabin osztálya (1. 2. és 3. osztály)
  • Name: Az utas neve
  • Sex: Az utas neme
  • Age: Az utas életkora években
  • SibSp: Az utassal együtt utazó testvérek vagy házastársak száma
  • Parch: Az utassal együtt utazó szülők vagy gyermekek száma
  • Ticket: Az utas jegyének sorszáma
  • Fare: Fizetett videldíj angol fontban
  • Cabin: A kabin sorszáma
  • Embarked: Melyik állomáson szállt fel a hajén? (S: Southampton Q: Queenstown C: Cherbourg)

A továbbiakban viszont csak 4 változóval fogunk foglalkozni, szóval válasszuk ki ezeket: Survived, Sex, Age, Fare. A többi változót eldobhatjuk a titanic data frame-ből:

setwd("E:/Egyetem Tanársegéd/1. félév/Öko. I/R jegyzet")

library(readr)
## Warning: package 'readr' was built under R version 4.0.3
titanic <- read_csv("titanic.csv")
## 
## -- Column specification --------------------------------------------------------
## cols(
##   PassengerId = col_double(),
##   Survived = col_double(),
##   Pclass = col_double(),
##   Name = col_character(),
##   Sex = col_character(),
##   Age = col_double(),
##   SibSp = col_double(),
##   Parch = col_double(),
##   Ticket = col_character(),
##   Fare = col_double(),
##   Cabin = col_character(),
##   Embarked = col_character()
## )
titanic <- titanic[,c("Survived", "Sex", "Age", "Fare")]
str(titanic)
## tibble [891 x 4] (S3: tbl_df/tbl/data.frame)
##  $ Survived: num [1:891] 0 1 1 1 0 0 0 0 1 1 ...
##  $ Sex     : chr [1:891] "male" "female" "female" "female" ...
##  $ Age     : num [1:891] 22 38 26 35 35 NA 54 2 27 14 ...
##  $ Fare    : num [1:891] 7.25 71.28 7.92 53.1 8.05 ...

Az str függvény eredménye alapján érdemes lenne még az elemzéshez factor típusúvá konvertálni a Sex változót character típusból és a Survived változót numeric típusból (hiszen a Survived is nominális mérési skála):

titanic$Sex <- as.factor(titanic$Sex)
titanic$Survived <- as.factor(titanic$Survived)

summary(titanic)
##  Survived     Sex           Age             Fare       
##  0:549    female:314   Min.   : 0.42   Min.   :  0.00  
##  1:342    male  :577   1st Qu.:20.12   1st Qu.:  7.91  
##                        Median :28.00   Median : 14.45  
##                        Mean   :29.70   Mean   : 32.20  
##                        3rd Qu.:38.00   3rd Qu.: 31.00  
##                        Max.   :80.00   Max.   :512.33  
##                        NA's   :177

Ezzel akkor így szépen rendben vagyunk, a summary függvény minden változóra a mérési skálájának megfelelő leíró statisztikákat adja vissza (nominálisokra eloszlás, numerikus változókra átlag, kvantilisek, min-max).

Láthatjuk, hogy az Age változóban vannak hiányzó értékek. Most csúnyák leszünk, és ezeket az átlagos korrak pótoljuk:

titanic$Age[is.na(titanic$Age)] <- mean(titanic$Age, na.rm = TRUE)

summary(titanic)
##  Survived     Sex           Age             Fare       
##  0:549    female:314   Min.   : 0.42   Min.   :  0.00  
##  1:342    male  :577   1st Qu.:22.00   1st Qu.:  7.91  
##                        Median :29.70   Median : 14.45  
##                        Mean   :29.70   Mean   : 32.20  
##                        3rd Qu.:35.00   3rd Qu.: 31.00  
##                        Max.   :80.00   Max.   :512.33

Oké, nincs semmi gáz már az Age változóban sem! :)

A logisztikus regresszió alapfeladata az lenne, hogy a binárisan adott eredményváltozót, tehát azt, hogy az adott utas túlélte-e a hajókirándulást (Survived) próbáljuk megbecsülni az utas neme, kora és fizetett viteldíja alapján, mint magyarázóváltozók!

2. A Logit transzformáció

Először egyszerűsítünk, és csak a két numerikus változónk (Age és Fare) lesz a két magyarázóváltozónk: \(x_1\): Age és \(x_2\): Fare

Ha egy regresszióban az eredményváltozónk egy nominális változó, 2 lehetséges értékkel, mint a Survived, akkor a regressziós egyenlet sima OLS becslése NEM fog működni, mivel az \(SSE\) minimalizálásával az \(\hat{y}\) (becsült \(y\)) értékeink bárhova eshetnek a (-∞, ∞) intervallumon, de nekünk csak a 0 és 1 értékek lesznek jók \(\hat{y}\)-nak!

A nagy ötlet: Alakítsuk át a bináris \(y\) értékeket, hogy a (-∞, ∞) tartományon értelmezhetőek legyenek. Majd ezekere az átalakított értékekre írjunk fel egy sima lineáris regressziós egyenletet \(\beta_j\)-kel meg mindennel!

Az átalakítás lépései:

  1. A \(P(y=1)=p\) (túlélési) valószínűségek folytonosak már a \((0, 1)\) intervallumon. Emiatt próbáljuk meg ezeket becsülni a tényleges \(y\)-ok helyett!
  2. Az odds hányadosok: \(odds = \frac{p}{1-p}\) (mennyivel valószínűbb, hogy egy utas túlél, mint hogy meghal) már a [0, ∞) tartományon mozognak…
  3. … és a \(ln(odds)=ln(\frac{p}{1-p})\) értékek már a kívánt (-∞, ∞) tartományba esnek! Az \(ln(odds)\) értékeket hívjuk logitoknak!

A lineáris egynletünk tehát a \(ln(odds)\)-okat fogja becsülni: \(ln(\frac{p}{1-p})=\beta_1x_1 + \beta_2x_2 + \beta_0\).

A logit becslésekből majd vissza tudjuk számolni a \(p = P(y=1)\) túlélési valószínűségeket: \(p=\frac{exp(logit)}{1+exp(logit)}\)

Ha \(p>0.5\), akkor pedig \(\hat{y}=1\) (túléltünk), egyébként \(\hat{y}=0\) (meghaltunk).

Először, mint az OLS regresszió elején adunk valami kezdeti random értékeket minden \(\beta_j\)-ra:

Beta1 <- -0.05
Beta2 <- 0.05
Beta0 <- -0.5

Majd számoljuk ki a logitokat mind a 891 utasra a regressziós egyenletből:

titanic$Logit <- Beta1 * titanic$Age + Beta2 * titanic$Fare + Beta0
head(titanic)
## # A tibble: 6 x 5
##   Survived Sex      Age  Fare  Logit
##   <fct>    <fct>  <dbl> <dbl>  <dbl>
## 1 0        male    22    7.25 -1.24 
## 2 1        female  38   71.3   1.16 
## 3 1        female  26    7.92 -1.40 
## 4 1        female  35   53.1   0.405
## 5 0        male    35    8.05 -1.85 
## 6 0        male    29.7  8.46 -1.56

A logitokból számoljuk ki a \(P(y=1) = p\) valószínűségeket:

titanic$Probability <- exp(titanic$Logit) / (1+ exp(titanic$Logit))
head(titanic)
## # A tibble: 6 x 6
##   Survived Sex      Age  Fare  Logit Probability
##   <fct>    <fct>  <dbl> <dbl>  <dbl>       <dbl>
## 1 0        male    22    7.25 -1.24        0.225
## 2 1        female  38   71.3   1.16        0.762
## 3 1        female  26    7.92 -1.40        0.197
## 4 1        female  35   53.1   0.405       0.600
## 5 0        male    35    8.05 -1.85        0.136
## 6 0        male    29.7  8.46 -1.56        0.173

Végül adjuk meg a \(\hat{y}\)-kat:

titanic$y_kalap <- 0 # alapból azt mondom mindenki meghal
titanic$y_kalap[titanic$Probability > 0.5] <- 1 # akinek a modell szerint 50%-nál nagyobb a túlélési valószínűsége, nála felülírom a 0-t egy 1-gyel
head(titanic)
## # A tibble: 6 x 7
##   Survived Sex      Age  Fare  Logit Probability y_kalap
##   <fct>    <fct>  <dbl> <dbl>  <dbl>       <dbl>   <dbl>
## 1 0        male    22    7.25 -1.24        0.225       0
## 2 1        female  38   71.3   1.16        0.762       1
## 3 1        female  26    7.92 -1.40        0.197       0
## 4 1        female  35   53.1   0.405       0.600       1
## 5 0        male    35    8.05 -1.85        0.136       0
## 6 0        male    29.7  8.46 -1.56        0.173       0

Láthatjuk, hogy nem olyan rossz a modellünk az első 6 sor alapján: csak a 3. sorban mondta tévesen az \(\hat{y}\) becslés azt, hogy utasunk meghal (\(\hat{y}=0\)) ahelyett, ami a valóság, hogy túlélte a találkát a jégheggyel (\(y=1\)). A maradék 5 sorban jók vagyunk, mert \(\hat{y}=y\).

3. A Maximum Likelihood becslés

Eddig minden remekül halad, de a mi \(\beta_j\)-ink teljesen randomak. Olyanok kellenek, amik minimalizálják avalmilyen értelemben az eltéréseket \(\hat{y}\) és a valós \(y\)-ok között! De hogyan találjuk meg a \(\beta_j\)-kat, amik minimalizálják a becslési hibát?

Egy úgynevezett Maximum Likelihood módszert (ML röviden) használunk! Ez azt jelenti, hogy olyan \(\beta_j\)-kat keresünk amik maximalizálják annak a valószínűségét (likelihood-ját), hogy egy véletlen mintavételből épp az előttunk lévő, a titanic data frame-ben lakó, mintát kapjuk!!!

Oké, de hogyan mérjük ezt a valószínűséget?

Először is, meghatározzuk, hogy milyen valószínűséggel találunk meg minden utast olyan \(\beta_j\)-k mellett, amiket épp megadtunk.

Ezt kiszámolni elég könnyű, hiszen a valóságban \(y=1\)-es utasok \(P(y=1) = p\) valószínűséggel fordulnak elő a világban. Míg a valóságban \(y=0\)-s utasok \(P(y=0) = 1-p\) valószínűséggel fordulnak elő a világban. Ezt az információt meg tudjuk adni egy új oszlopában a data frame-ünknek:

titanic$Likelihood <- 0 # kezdetben mindenki 0 valszínűséggel fordulhat elő
titanic$Likelihood[titanic$Survived == 1] <- titanic$Probability[titanic$Survived == 1] # p megy az y=1 esetekhez
titanic$Likelihood[titanic$Survived == 0] <- 1 - titanic$Probability[titanic$Survived == 0] # 1-p megy az y=0 esetekhez

head(titanic)
## # A tibble: 6 x 8
##   Survived Sex      Age  Fare  Logit Probability y_kalap Likelihood
##   <fct>    <fct>  <dbl> <dbl>  <dbl>       <dbl>   <dbl>      <dbl>
## 1 0        male    22    7.25 -1.24        0.225       0      0.775
## 2 1        female  38   71.3   1.16        0.762       1      0.762
## 3 1        female  26    7.92 -1.40        0.197       0      0.197
## 4 1        female  35   53.1   0.405       0.600       1      0.600
## 5 0        male    35    8.05 -1.85        0.136       0      0.864
## 6 0        male    29.7  8.46 -1.56        0.173       0      0.827

Ebben a Likelihood nevű oszlopban található akkor, hogy a jelenlegi random választott \(\beta_j\) értékeink mellett egy konkrét egyed bekövetkezési valószínűsége mekkora egy véletlen mintavétel esetén.

Például az első sorban található elhunyt (nem túlélő) 22 éves férfi, aki 7.25 fontot fizetett az útért 77.51%-os valószínűsége van megjelennie a mintánkban a jelenlegi random választott \(\beta_j\) értékeink mellett!

Vagy a harmadik sorban található túlélő 26 éves nő, aki 7.925 fontot fizetett az útért 19.72%-os valószínűsége van megjelennie a mintánkban a jelenlegi random választott \(\beta_j\) értékeink mellett!

Második lépésben feltesszük, hogy a 891 elemű mintánk FÜGGETLEN megfigyelésekből áll, a teljes minta bekövetkezési valószínűsége az új, Likelihood oszlopban lévő értékek szorzata

Azaz \(\prod_{y=1}{p} \times \prod_{y=0}{1-p} = L\) … ez itt annak Likelihood-ja a jelenlegi \(\beta_j\)-k mellett, hogy a mi 891 utasunkat húzzuk ki, ha a Titanic utasaiból 891 elemű mintát veszünk

prod(titanic$Likelihood)
## [1] 5.023418e-273

Ínnye! Ez gyakorlatilag egy nagy kövér 0-ra jön ki (ugye \(10^{-273}\) a nagységrend). :( Ez nem meglepő, hiszen 891 db 0 és 1 közötti számot szoroztunk össze…nem véletlen, hogy numerikusan annyira 0 közel vagyunk, hogy a gépállat az értéket konkrétan 0-nak is veszi! Emiatt inkább a \(-ln(L) = \sum_{y=1}{ln(p)}+\sum_{y=0}{ln(1-p)}\) értéket MINIMALIZÁLJUK majd:

-sum(log(titanic$Likelihood))
## [1] 626.9916

Akkor hát keressük is meg mik az a \(\beta_j\)-k, amikkel a lehető legkisebb \(-ln(L)\) értéket kapjuk! A 3. gyakorlathoz hasonlóan megcsinálunk egy külön R függvényt a negatív log-likelihoodra, azaz a \(-ln(L)=\)-re, ami a \(l(\beta_0, \beta_1,\beta_2)\) Alakot ölti. Tehát, megadja a negatív log-likelihoodot \(\beta\)-k függvényében:

# a függvény megadása
neg_loglikelihood <- function(x) {
  titanic$Logit <- x[1] * titanic$Age + x[2] * titanic$Fare + x[3]
  titanic$Probability <- exp(titanic$Logit) / (1+ exp(titanic$Logit))
  
  titanic$Likelihood <- 0 # kezdetben mindenki 0 valszínűséggel fordulhat elő
  titanic$Likelihood[titanic$Survived == 1] <- titanic$Probability[titanic$Survived == 1] # p megy az y=1 esetekhez
  titanic$Likelihood[titanic$Survived == 0] <- 1 - titanic$Probability[titanic$Survived == 0] # 1-p megy az y=0 esetekhez
  
  -sum(log(titanic$Likelihood))
}

# a függvény használata úgy, hogy mindhárom Béta1=-0.05 || Béta2=0.05 || Béta0=-0.5
neg_loglikelihood(c(-0.05,0.05,-0.5))
## [1] 626.9916

Láss csodát, ugyan azt kapjuk negatív log-likelihoodra, mint az előbb! :)

Na, ezt az neg_loglikelihood függvényt felhasználva meg tudjuk kerestetni gépállattal hol vannak azok a \(\beta\)-k, amik mellett a legkisebb negatív log-likelihoodot, azaz a legnagyobb minta bekövetkezési valószínűséget kapjuk:

eredmeny <- optim(c(-0.05,0.05,-0.5),neg_loglikelihood) # az optimalizálást (hiba minimalizálást) az eddig is használt random béta értékekből indítjuk

Az újonnan létrehozott, eredmeny objektumból ki tudjuk olvasni, hogy mik lettek a \(\beta\)-k, és meg tudjuk nézni azt is, hogy mi az a legkisebb negatív log-likelihood érték, amit el tudtunk érni:

eredmeny$par # a Béták
## [1] -0.01738088  0.01619721 -0.45787001
eredmeny$value # a minimalizált -ln(L) érték
## [1] 553.9448

Ezzel a végső regressziós egyenletünk: \(logit(Survived)=-0.017 \times Age + 0.016 \times Fare - 0.46\). Azaz a \(\beta\)-k maximum likelihood becslése.

4. Logisztikus regresszió beépített R függvénnyel

Szerencsére, a logisztikus (vagy röviden csak logit) regressziónak is van olyan beépített R függvénye, mint az OLS regressziónak az lm. A Logit regresszió R függvénye a glm névre hallgat, és az eredményét ugyan úgy a summary függvénnyel lehet meglesni, mint a sima lm függvény esetében csináltuk:

logisztikus <- glm(Survived ~ Age + Fare, family = binomial(link='logit'), data = titanic) # itt a family paraméterben tudjuk konkretizálni, hogy az eredményváltozónk bináris, és, hogy ezt logit transzformációval kezelje a gépállat
summary(logisztikus)
## 
## Call:
## glm(formula = Survived ~ Age + Fare, family = binomial(link = "logit"), 
##     data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6464  -0.9008  -0.8294   1.2662   1.7993  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.457668   0.181886  -2.516  0.01186 *  
## Age         -0.017387   0.005666  -3.069  0.00215 ** 
## Fare         0.016195   0.002316   6.994 2.68e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.7  on 890  degrees of freedom
## Residual deviance: 1107.9  on 888  degrees of freedom
## AIC: 1113.9
## 
## Number of Fisher Scoring iterations: 5

Szuper, a \(\beta\)-k ugyan azok, mint amik nekünk is kijöttek manuálisan:

eredmeny$par
## [1] -0.01738088  0.01619721 -0.45787001

Láthatjuk, hogy a regressziós eredménytábla olyan, mint a lineáris esetben: minden \(\beta\)-hoz tartozik egy standard hiba, amiből lehet próbafüggvényt és p-értéket számolni a szokásos \(H_0: \beta_j=0\) és \(H_1: \beta_j \neq 0\) hipotézisvizsgálathoz. Ez alapján azt mondhatjuk, hogy a Fare és Age magyarázóváltozók hatása a mintán kívüli világban (sokaságban) is szignifikáns még \(\alpha=1\%\) esetén is, de a tengelymetszet már nem szignifikáns 1%-on, csak 5%-os \(\alpha\) mellett! Plusz, mivel a Fare változó p-értéke kisebb, mint az Age p-értéke, így itt is azt mondhatjuk, hogy a sokságban a viteldíj jobban befolyásolja a túlélést, mint az utas kora!

4.1. Marginális hatások Logisztikus regresszióban

Node, hogyan is kell érteni a magyarázóváltozóink margináláis hatását az eredményváltozóra? Mit jelentenek azok a \(\beta\)-k, amikről itt az előbb eldöntöttük nagyon ügyesen, hogy hatásuk szignifikáns a túlélésre mintán kívüli megfigyelsekre is?

Nos, a \(\beta\)-k önmagukban nem igazán jól értelmezhetők logisztikus regresszióban. Hiszen ezek a modell egyenletéből adódóan ceteris paribus +1 egység \(x\) növekedésének hatását az eredményváltozó \(logit\) transzformáltjára mutatják ki.
De ha \(e\)-re emeljük az együtthatókat, akkor már az \(e^{logit(y)}\)-ra, azaz az \(odds\) hányadosokra gyakorolt marginális hatásait kapjuk meg a magyarázóváltozóinknak!

Szóval tegyük is ezt meg:

exp(logisztikus$coefficients)
## (Intercept)         Age        Fare 
##   0.6327578   0.9827629   1.0163264

Értelmezzük akkor a két együtthatót (a tengelymetszettel ugyan úgy nem törődünk igazából, mint ahogy az OLS modellek \(\beta_0\)-val sem törődtünk értelmezési szempontból):

  • \(\beta_{Age}=0.98\): Ha minden más magyarázóváltozó változatlansága mellett egy utas életkora 1 évvel nő, akkor az ő túlélési ODDSai 0.98-szorosra változnak, azaz 2%-kal csökkennek.
  • \(\beta_{Fare}=1.016\): Ha minden más magyarázóváltozó változatlansága mellett egy utas viteldíja 1 fonttal nő, akkor az ő túlélési ODDSai 1.016-szorosra változnak, azaz 1.6%-kal nőnek.

Teljesen logikusak az előjelek:

  • Ugye először a gyerekeket mentették ki, szóval nyilván minél idősebb vagy, változatlan fizetett viteldíj mellett, annál kisebb a túlélési odds-od (a túlélési és halálozási valószínűség hányadosa).
  • Utána pedig a nők mellett a minél pénzesebb, gazdag embereket (akik sokat fizettek a jó kabinokért) rakták a csónakokba. Így logikus, hogy minél többet fizettél az útért, változatlan kor mellett, annál inkább nő a túlélési odds-od. Ugye DiCaprio karaktere a filmben nemcsak férfi volt, hanem szegény is! :)

5. Két logisztikus modell összehasonlítása

Csináljunk most egy olyan logisztikus regressziót, amiben ott van magyarázóváltozóként az utas neme (Sex) is:

logisztikus_bovebb <- glm(Survived ~ Age + Fare + Sex, family = binomial(link='logit'), data = titanic)
summary(logisztikus_bovebb)
## 
## Call:
## glm(formula = Survived ~ Age + Fare + Sex, family = binomial(link = "logit"), 
##     data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3396  -0.6159  -0.5786   0.8000   2.0599  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.921416   0.230665   3.995 6.48e-05 ***
## Age         -0.010303   0.006538  -1.576    0.115    
## Fare         0.011673   0.002339   4.990 6.03e-07 ***
## Sexmale     -2.399166   0.171021 -14.029  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.7  on 890  degrees of freedom
## Residual deviance:  881.8  on 887  degrees of freedom
## AIC: 889.8
## 
## Number of Fisher Scoring iterations: 5

Láthatjuk, hogy a Sex egy Férfi-e? (1: igen 0: nem) dummy változóval lett bevonva a modellbe a Sexmale változónévből. Ha \(e\)-re emelem ezt a \(\beta\)-t, akkor értelmezhetővé is válik a nem marginális hatása:

exp(logisztikus_bovebb$coefficients[4]) # a Sex változó bétája a 4. az együtthatók sorában
##    Sexmale 
## 0.09079361

Változatlan viteldíj és életkor mellett a férfiaknak kb. 91%-kal kisebb a túlélési odds-uk, mint a nőké! Kellemetlen! :( Sőt, a Sex bevonásával az Age változó nem-szignifikánsá vált minden szokásos \(\alpha\)-n. Úgy tűnik, hogy az életkor csak a férfi/nő hatást közvetítette eddig. Azonos nemű és viteldíjú embereket véve az életkornak nincs szignifikáns hatása a túlésére.

De honnan tudjuk, hogy a Sex változóval bővített modell tényleg jobban illeszkedik az eredményváltozóra, mint az anélküli modell? Hát az információs kritériumokat szerencsére a logisztikus regresszióban is lehet használni! Ugyanis, az \(SSE\) funkcióját a logitban a negatív log-likelihood veszi át! Hiszen itt a negatív log-likelihood (\(-ln(L)\)) minimalizálásával kaptuk meg a regresszió \(\beta_j\) együtthatóit! Innentől kezdve pl. egy OLS regressziós \(AIC=ln(SSE)+2 \times k\) képlet logisztikus regresszióban \(AIC=-ln(L)+2 \times k\)-vá válik. Egy \(BIC=ln(SSE) + k \times ln(n)\) képletből pedig a logitban \(BIC=-ln(L)+k \times ln(n)\) lesz. Ebben az esetben is \(k\) természetesen a modell paramétereinek, azaz \(\beta\)-inak száma.

Itt is igaz az, mint OLS-nél, hogy mivel \(-ln(L)\)-ből (tehát lényegében a modell teljes hibájából) indulnak ki, és arra pakolnak rá valami \(k\)-tól függő büntetőtagot, így arra jók, hogy ha 2 modellre kiszámítjuk őket, akkor a kisebb \(IC\)-jű modellt jobbnak tudjuk tekinteni!

Szóval akkor hasonlítsuk is össze \(IC\)-k alapján a két logit modellünket:

AIC(logisztikus, logisztikus_bovebb)
##                    df       AIC
## logisztikus         3 1113.8897
## logisztikus_bovebb  4  889.7983
BIC(logisztikus, logisztikus_bovebb)
##                    df       BIC
## logisztikus         3 1128.2667
## logisztikus_bovebb  4  908.9676

Még a magyarázóváltozókszámát legszigorúbban büntető \(BIC\) is a bővebb modell mellett teszi le voksát, így valóban megérte bevenni a modellbe a Sex változót! :)

Ezen a ponton persze még kísérletezhetünk az Age kihagyásával a modellből, de most inkább egy másik irány felé vesszük az utunk: egy logisztikus regressziós modell teljesítményének önálló értékelése felé.

6. A logisztikus regresszió teljesítményének kiértékelése

Az előző fejezetben láttuk tehát, hogy lehet összehasonlítani az információs kritériumok segítségével két logit modell becslési pontosságát. De hogyan tudjuk azt megmondani, hogy egy logit modell önmagában mennyire jó? Hogy hány %-ban magyaráztuk meg a modellünkkel eredményváltozó alakulását? Tehát, valami \(R^2\)-szerű dolgot szeretnénk logisztikus regresszióra is!

Első körben nézzük meg, hogy a glm függvénnyel felépített logisztikus_bovebb modellből hogyan tudunk becsléseket, azaz \(\hat{y}\)-kat berakni a titanic nevű data frame-be.

Kipucoljuk a manuális maximum likelihood becsléshez készített oszlopainkat:

titanic[,5:8] <- NULL
head(titanic)
## # A tibble: 6 x 4
##   Survived Sex      Age  Fare
##   <fct>    <fct>  <dbl> <dbl>
## 1 0        male    22    7.25
## 2 1        female  38   71.3 
## 3 1        female  26    7.92
## 4 1        female  35   53.1 
## 5 0        male    35    8.05
## 6 0        male    29.7  8.46

Majd alkalmazzuk a predict függvényt a logisztikus_bovebb modell R objektumán:

titanic$Probability_Bovebb <- predict(logisztikus_bovebb, newdata = titanic, type = "response")
head(titanic)
## # A tibble: 6 x 5
##   Survived Sex      Age  Fare Probability_Bovebb
##   <fct>    <fct>  <dbl> <dbl>              <dbl>
## 1 0        male    22    7.25              0.165
## 2 1        female  38   71.3               0.796
## 3 1        female  26    7.92              0.678
## 4 1        female  35   53.1               0.765
## 5 0        male    35    8.05              0.149
## 6 0        male    29.7  8.46              0.156

A type = "response" beállítással a predict függvény minden egyedre a \(P(y=1)\) valószínűségeket rakja vissza a táblába. Ha nem állítottuk volna be ezt a paramétert, akkor a logitokat rakta volna vissza.

Láthatjuk, hogy a Sex-et is tartalmazó logit modellünk alapján a 3. sorban a becsült túlélési valószínűség kerekítve \(67.8\%\).

Ebből a \(P(y=1)\) oszlopból kiindulva megadhatjuk a \(\hat{y}\)-kat is! Ahogy korábban is csináltuk, most is mehet az a szabály, hogy \(P(y=1)>0.5 \rightarrow \hat{y}=1\), és \(\hat{y}=0\) egyébként. Ekkor az \(50\%=0.5\).öt egy úgynevezett cut-off értéknek hívjuk, és igazából nem muszáj 50%-nak választani! Az 50% a logikus alapbeállítás, de a következő alkalommal látni fogunk olyan eseteket, ahol megéri már százalékot választani cut-offnak ahhoz, hogy \(\hat{y}=1\) becsléseket adjunk a \(P(y=1)\) értékek alapján.
De akkor most maradunk a \(c=0.5\)-ös cut-offnál:

titanic$y_kalap <- 0 # alapból azt mondom mindenki meghal
titanic$y_kalap[titanic$Probability_Bovebb > 0.5] <- 1 # akinek a modell szerint 50%-nál nagyobb a túlélési valószínűsége, nála felülírom a 0-t egy 1-gyel
head(titanic)
## # A tibble: 6 x 6
##   Survived Sex      Age  Fare Probability_Bovebb y_kalap
##   <fct>    <fct>  <dbl> <dbl>              <dbl>   <dbl>
## 1 0        male    22    7.25              0.165       0
## 2 1        female  38   71.3               0.796       1
## 3 1        female  26    7.92              0.678       1
## 4 1        female  35   53.1               0.765       1
## 5 0        male    35    8.05              0.149       0
## 6 0        male    29.7  8.46              0.156       0

Elég menő a modell első ránézésre, mert az első 6 sorban a becslésünk egyezik az valós Survived (azaz \(y\)) értékekkel!

De lássuk az összképet! Nézzünk megy egy gyakorisági kereszttáblát a Survived és y_kalap változók szerint! Ezt hívja a szaknyelv konfúziós vagy klasszifikációs mátrixnak!

Ezt megcsinálhatom simán a table függvénnyel. Most jól fog jönni nekünk, ha a kereszttáblába beírja az R az oszlopneveket is, így a table függvénybe egy data frame-et rakunk, aminek két oszlopa a y_kalap és a Survived :

table(titanic[,c("y_kalap", "Survived")])
##        Survived
## y_kalap   0   1
##       0 462 106
##       1  87 236

Lássuk hát a különböző gyakoriságok értelmezését a kereszttáblában:

  • \(462\): 462 utast becsült a logit modellünk úgy halottnak, hogy a valóságban is meghaltak.
  • \(236\): 236 utast becsült a logit modellünk úgy túlélőnek, hogy a valóságban is túlélték az utazást.
  • \(106\): 106 utast becsült a logit modellünk tévesen halottnak, hiszen a valóságban túlélték az utazást.
  • \(87\): 87 utast becsült a logit modellünk tévesen túlélőnek, hiszen a valóságban meghaltak.

Ezekből a gyakoriságokból pedig ötféle pontossági mutatót is tudunk számolni a logit modellünk kiértékeléséhez:

  • accuracy: \(\frac{462+236}{462+236+106+87}=0.7834=78.3\%\) Az összes utas 78.34%-át osztályozta helyesen a modell.
  • precision(1): A modell által élőnek prediktált utasok közül \(\frac{236}{236+87}=73.07\%\) élte tényleg túl az utat.
  • precision(0): A modell által halottnak prediktált utasok közül \(\frac{462}{462+106}=81.34\%\) lett tényleg halott.
  • recall(1): A modell által élőnek prediktált utasok között van az összes valójában túlélő utas \(\frac{236}{236+106}=69.01\%\)-a.
  • recall(0): A modell által halottnak prediktált utasok között van az összes valójában halott utas \(\frac{462}{462+87}=84.15\%\)-a.

Azt, hogy az 5 lehetséges kiértékelési mutató közül melyiket alkalmazzuk alapvetően az üzleti környezet dönti el.
Például, ha egy telekom cég kér meg minket, hogy keressünk neki potenciális ügyfeleket, akik mondjuk telefonos megkeresés hatására szerződést kötnének a vállalattal, akkor csinálhatjuk azt, hogy korábbi hasonló marketing kampány adatok alapján egy olyan bináris eredményváltozóra építünk logit modellt, ahol 1 = szerződéskötés és 0 = telefon lecsapása.
Egy ilyen logit modell kiértékelése során, ha azt mondja a marketing igazgató, hogy most nincs sok pénze a kampányra, akkor a modellt a precision(1) mutató alapján érdemes értékelni: ekkor ugye épp az a cél, hogy akiket 1-nek, azaz megkeresendőnek prediktálunk, azoknak minél nagyobb %-a kössön a céggel tényleg szerződést. Nincs sok pénz, ne legyen sok elpazarolt hívás.
De ha a főni épp azt mondja, hogy most van sok zseton a hívásokra, és találjunk neki minél több potenciális ügyfelet, akkor a modellt épp a recall(1) mutató szerint kell értékelni. Hiszen ekkor a lehető legnagyobb %-ban azonosítani akarunk a modellel minden \(y=1\) egyedet. Még ha ennek az is az ára, hogy valamivel több téves megkeresés is lesz a hívások között.

7. Ritka eseményre épített logisztikus regresszió - Csődmodellezés

Töltsük be a Moodle-n található csod.xlsx adattáblát, ami 1969 magyar vállalatról tartalmaz 8 számviteli jellegű változót:

  • csod: A vállalat csődbe ment-e (1: igen; 0: nem)
  • eszkfseb100: eszközök forgási sebessége = nettó árbevétel / eszközök teljes értéke
  • stokear100: saját tőke aránya = saját tőke értéke / források teljes értéke
  • bonitas100: bonitás = kötelezettségek értéke / saját tőke értéke
  • jovedelm100: eszközjövedelmezőség = (adózott eredmény + értékcsökkenési leírás) / eszközök teljes értéke
  • feszkar100: forgó eszközök aránya = forgóeszközök értéke / eszközök teljes értéke
  • likvid100: likviditás = (forgóeszközök értéke -rövid lejáratú kötelezettségek értéke) / eszközök teljes értéke

Készítsünk a csod változóra egy logisztikus regressziós modellt, ahol a data frame-ben lévő összes többi változót használjuk fel magyarázóváltozóként. Biztos, ami biztos alapon a modell építése előtt az eredményváltozónkat konvertáljuk át factor típusúvá:

library(readxl)
csod <- read_excel("csod.xlsx")

csod$csod <- as.factor(csod$csod)
str(csod)
## tibble [1,969 x 7] (S3: tbl_df/tbl/data.frame)
##  $ csod       : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ eszkfseb100: num [1:1969] 29.64 125.77 219.11 6.09 4.59 ...
##  $ stokear100 : num [1:1969] 99.9 98.9 99.8 57.9 16.3 ...
##  $ bonitas100 : num [1:1969] 0.124 1.114 0.158 72.687 512.34 ...
##  $ jovedelm100: num [1:1969] 15.23 -22.69 32.22 -5.48 -2.62 ...
##  $ feszkar100 : num [1:1969] 9.98 100 88.71 1.24 92.81 ...
##  $ likvid100  : num [1:1969] 9.85 98.9 88.55 1.07 92.76 ...
csod_logit <- glm(csod ~ ., family = binomial(link='logit'), data = csod)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(csod_logit)
## 
## Call:
## glm(formula = csod ~ ., family = binomial(link = "logit"), data = csod)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0316  -0.2365  -0.1288  -0.0762   3.6420  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.1360175  0.8974953  -3.494 0.000476 ***
## eszkfseb100 -0.0034575  0.0012475  -2.772 0.005579 ** 
## stokear100  -0.0462510  0.0103454  -4.471 7.80e-06 ***
## bonitas100  -0.0012246  0.0003241  -3.778 0.000158 ***
## jovedelm100 -0.0106813  0.0024883  -4.293 1.77e-05 ***
## feszkar100   0.0425958  0.0106822   3.988 6.68e-05 ***
## likvid100   -0.0213708  0.0101194  -2.112 0.034697 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 831.34  on 1968  degrees of freedom
## Residual deviance: 421.26  on 1962  degrees of freedom
## AIC: 435.26
## 
## Number of Fisher Scoring iterations: 15

Láthatjuk, hogy minden magyarázóváltozónk szignifikáns minden szokásos \(\alpha\)-n, kivéve a likvid100, de még az is szignifikánsnak vehető \(\alpha=5\%\) mellett. Most nem végzünk modellszelekciót az \(IC\)-k, alapján, hanem élből ezt a modellt próbáljuk meg értékelni a klasszifikációs mátrix alapján. A cut-off most is legyen az első körben logikus 50%:

csod$CsodVsz <- predict(csod_logit, newdata = csod, type = "response") # csődvalószínűség becslése

csod$BecsultCsod <- 0 # kezdetben mindenki fizetőképes marad
csod$BecsultCsod[csod$CsodVsz > 0.5] <- 1 # akinek 50%-nál magasabb a csődvalószínűsége, azt csődösnek jelöljük

table(csod[,c("csod", "BecsultCsod")]) # klasszifikációs mátrix
##     BecsultCsod
## csod    0    1
##    0 1862    0
##    1   53   54

Nézzünk meg a klasszifikációs mátrixból az accuracy mutatót: \(\frac{1862+54}{1862+54+53}=97.3\%\). Tyű, ez a 97%-os pontosság nem semmi! NODE: pesszimista közgazdászok vagyunk! Nézzünk csak rá a csődesemény (\(y=1\)) precision és recall mutatójára!

\(precision(1)=\frac{54}{54+0}=100\%\). Tehát, a modell által csődosnek becsült cégekből mind csődös a valóságban is! Ez elég jónak néz ki még mindig!

\(recall(1)=\frac{54}{54+53}=50.5\%\). Szóval, az összes csődös cégnek csak kicsivel több, mint a felét becsüli tényleg csődösnek a logit modellünk! Ez azért már GázGéza!
Ha jobban belegondolunk, egy banknak az a fontos, hogy a lehető legtöbb tényleg csődös céget azonosítsuk, hiszen az nagyobb probléma, ha a bank egy olyan cégnek ad hitelt, ami utána csődbemegy, mint ha egy olyannak nem ad hitelt, aki amúgy nem megy csődbe!

Szóval, az a nagy problémánk jelenleg, hogy az \(y=1\) esetekből nagyon kevésnek (kicsivel több, mint 50%-nak) lesz a modellünk alapján \(\hat{y}=1\) becslése!
Ha meglessük a klasszifikációs mátrixunkat, akkkor látjuk is szívünk nagy fájdalmának az okát: összesen \(53+54=107\) db csődös vállalat van abban az \(1969\) elemű mintában. Tehát a csődesemény (\(y=1\)) egy kimondottan ritka esemény, relatív gyakorisága mindössze \(\frac{107}{1969}=5.4\%\). Ezt úgy is szoktuk mondani kultúrnyelven, hogy a bináris eredméynváltozónk kiegyensúlyozatlan.

Mindez pedig azt eredményezi, hogy a becsült csődvalószínűségeink, azaz a \(P(y=1)\) valószínűségek, eloszlása elég szépen jobbra elnyúló lesz:

hist(csod$CsodVsz)

Láthatjuk, hogy alig akad olyan becsült \(P(y=1)\) valószínűség, ami 50%-nál nagyobb. Szóval, elég ritkásan fogja a modell osztogatni 50%-os cut-off mellett az \(\hat{y}=1\) becsléseket, ahogy látjuk is a klasszifikációs mátrixban. Valami hasonló amúgy ez a probléma, mint a heteroszkedaszticitás volt az OLS modellben: az eredményváltozó egy adott értéktartományából kevés megfigyelésünk van, és ezen az értéktartományon jobban hibáznak a regressziónk \(\hat{y}\) becslései.

A probléma egy kézzelfekvő megoldása az lehet, hogy hát ha tudjuk, hogy a mintánkban az egyedek 5.4%-a csődös cég, akkor vegyünk minden céget csődösnek, aminél a modellünk szerint a csődvalószínűség 5.4%-nél nagyobb. Magyarul, minden céget csődösnek mondunk, ahol a becsült csődvalószínűség az “átlagosnál” nagyobb.
Szakszóval meg azt mondjuk, hogy a cut-off értéket levisszük 50%-ról 5.4%-ra! Nézzük is meg, hogy így mi történik a klasszifikációs mátrixunkkal így:

csod$BecsultCsod <- 0 # kezdetben mindenki fizetőképes marad
csod$BecsultCsod[csod$CsodVsz > 0.054] <- 1 # akinek 5.4%-nál magasabb a csődvalószínűsége, azt csődösnek jelöljük

table(csod[,c("csod", "BecsultCsod")]) # klasszifikációs mátrix
##     BecsultCsod
## csod    0    1
##    0 1578  284
##    1   23   84

Láthatjuk, hogy kiegyensúlyozottabban viselkedik az eredményváltozó mindkét kategóriájában a modellünk!

A teljes osztályozási pontosság, az accuracy itt \(\frac{1578+84}{1578+84+23+284}=94.6\%\)-ra romlott a 97%-ról. Némileg rosszabb a dolog, de egyáltalán nem vészes!

Ellenben, \(recall(1)=\frac{84}{84+23}=78.5\%\)! Tehát az összes valójában csődös vállalat 78.5%-kát helyesen azonosította a modell!

Azért, ne örüljünk ennyire, mert ezt a 78.5%-os recall(1)-et úgy csináltunk, hogy a másik irányon nagyot rontottunk: \(percision(1)=\frac{84}{84+284}=22.8\%\). Ínnye! A modell által csődösnek becsült vállalatoknak csak durván 23%-a csődös a valóságban is.
Tehát oké, hogy a valójában csődös vállalatoknak nagy %-át megtalálja a modell, de ezt úgy teszi, hogy javarészt vaktában lövöldözik! Ez azért nem oké, mert így már számottevő csődmentes cégnek nem adna hitelt a bank, amiért már fájhat egy pénzéhes bankigazgató szíve!

Az egyensúlyt biztosító cut-off megtalálásában tud nekünk segítséget nyújtani a ROC-görbe!

8. A ROC görbe

A ROC görbe használatához először is ki kell jelölnünk egy úgynevezett pozitív osztályt. Ez a bináris eredményváltozónk azon értéke a kettőből, amit igazán pontosan meg akarunk becsülni az elemzésünk során.
Értelemszerűen, mi most a csődös vállalatokat szeretnénk minél pontosabban beazonosítani, azaz nagy %-ukat megtalálni, de anélkül, hogy túl sok téves riasztásunk legyen. Tehát, jelen példában a pozitív osztályunk az \(y=1\).

Ezek után két fogalmat kell bevezetni.

  1. A pozitív osztály True Positive Rate-je, más szóval a modell szenzitivitása. Ez gyakorlatilag a pozitív osztály recall mutatójának álneve. Tehát, azt adja meg, hogy a valóságban az eredményváltozó pozitív osztályába tartozó egyedek hány %-át azonosítja helyesen a modellünk.
  2. A negatív osztály True Positive Rate-je, más szóval a modell specificitása. Ez gyakorlatilag a negatív osztály recall mutatójának álneve. Tehát, azt adja meg, hogy a valóságban az eredményváltozó negatív osztályába tartozó egyedek hány %-át azonosítja helyesen a modellünk.

A ROC görbe amúgy igazából az 1 - specificitást használja. Tehát, azt az arányt, ami azt mutatja meg, hogy a valóságban az eredményváltozó negatív osztályába tartozó egyedek hány %-a lett tévesen a pozitív osztályba sorolva. Ez így lényegében a pozitív osztály False Positive Rate-je!

Számoljuk ki a jelenlegi modellünk klasszifikációs mátrixából a fenti 2+1 mutatót:

table(csod[,c("csod", "BecsultCsod")]) # klasszifikációs mátrix
##     BecsultCsod
## csod    0    1
##    0 1578  284
##    1   23   84

Tehát, szenzitivitás = True Positive Rate = \(TPR=\frac{84}{84+23}=78.5\%\)

Aztán, specificitás = \(\frac{1578}{1578+284}=84.7\%\)

Végül pedig 1 - specificitás = False Positive Rate = \(FPR = 1-0.847=15.3\%\) vagy \(\frac{284}{284+1578}\)

Szóval, végső soron egy olyan cut-off érték után kutatunk, aminél a TPR a lehető legnagyobb, míg a FPR a lehető legkisebb! A ROC görbe pedig pont azt csinálja, hogy sok-sok cut-off érték mellett végigszámolja a klasszifikációs mátrixot és vele a pozitív osztály TPR és FPR értékét.
Az eredményeket pedig egy kooridnáta rendszerben ábrázolja, ahol x koordináták = FPR és y koordináták = TPR a sok-sok cut-off érték mellett!

8.1. ROC görbe működése R-ben

Lássunk hát akkor egy ilyen ROC görbét az R-ben is! Ehhez egy pROC névre hallgató csomagot kell telepíteni.

install.packages("pROC")
library(pROC) # Szokásos módon ne törődjünk az esetleges Warningokkal! :)
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

Egy ilyen ROC görbét a csomag pROC függvényével hazhatunk létre. Értelemszerűen kelleni fognak hozzá a valós \(y\) értékek, és a logit modellból becsült \(P(y=1)\) valószínűségek. Az eredményeket elmentjük egy ROC_adatok című objektumba:

ROC_adatok <- roc(csod$csod, csod$CsodVsz)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases

Ebben a ROC_adatok objektumban (ami, ha megfigyeljük egy lista típusú cucc) megvannak a megvizsgált cut-off értékek a thresholds elemben. Ezeket a modellünk \(P(y=1)\) valószínűségek eloszlásának percentilisei alapján alapján határozza meg gépállat:

summary(ROC_adatok$thresholds)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##     -Inf 0.003662 0.010336          0.036219      Inf

Láthatjuk, hogy a vizsgált cut-offok mediánja kb. 1%-nál van.

A sensitivities és specificities elemekben értelemszerűen A vizsgált cut-offokhoz tartozó szenzitivitásokat és specificitásokat találhatjuk.

Ha megnézzük, hogy nagyjából hányadik cut-off lehet a nálunk most beállított \(5.4\%\), akkor megleshetjük, hogy kb. ugyan azt kapjuk-e szenzitivitásra és specificitásra, amiket mi is számolgattunk a klasszifikációs mátrixból:

ROC_adatok$thresholds[round(ROC_adatok$thresholds,3)==0.054] # 3 tizedesre kerekítjük a vizsgált cut-off-okat, hogy legyen egyezés
## [1] 0.05368371 0.05406985 0.05408672 0.05420627 0.05434217 0.05449725

Talán a legközelebb a mi 0.054-ünkhöz a 0.05406985 van. Nézzük meg ez hányadik elem, és nézzük meg a megfelelő értékeket a sensitivities és specificities elemekből:

which(round(ROC_adatok$thresholds,3)==0.054) # kereken 1600. elem a 0.05406985-nek megfelelő
## [1] 1599 1600 1601 1602 1603 1604
ROC_adatok$sensitivities[1600]
## [1] 0.7850467
ROC_adatok$specificities[1600]
## [1] 0.8480129

Láthatjuk, hogy szenzitivitás = 78.5% és specificitás = 84.8%, ahogy mi is számoltuk itt a 8. fejezet legelején.

Na, ezek után rajzoljuk ki ezt a cuki ROC görbét:

plot(ROC_adatok)

Kicsit fura az ábra. Ahogy írtam, az y tengelyen a TPR, azaz a szenzitivitás található. Ellenben az x tengelyre a gép a specificitás-t írja, noha az 1 - specificitást kéne ide írni.
Igen, csak az x tengely skálázása fordított!! Tehát, nem 0-tól, hanem 1-től indul…szóval valójában mégis az 1 - specificitás van az x tengelyen csak gépállat ilyen idiótán írja!. Such is life, mondá a bölcs kínai, ezzel együtt kell élni.

Érdemes még megfigyelni, hogy az ábrára be van húzva egy 45 fokos egyenes is. Ez egy olyan logit modellenek a ROC görbéje lenne, aminél minden egyedre \(P(y=1)=0.5\) lenne a becslés. Tehát ez az a modell, ami random találgatja, hogy a megfigyelés most éppen \(y=1\) vagy \(y=0\)-é.
Láthatjuk, hogy a mi modellünk ennél az állapotnál lényegesen jobb, mert minden cut-off (kivéve persze a kereken 0 és 1 értékeket) mellett magasabb TPR-t (szenzitivitás) és alacsonyabb FPR-t (1 - specificitás) produkál.
Mivel ez a ROC görbe igazából egy egységnégyzetben mozog (mindkét tengelyen egy \([0,1]\) intervallumon mozgó arányszámot ábrázolunk), így tudjuk, hogy az egyenes alatti terület \(0.5\).

A logit modellünkhöz tartozó ROC görbe alatti terület ennél a \(0.5\)-nél tehát tuti nagyobb. De konkrétan mennyi? Ezt integrálással ki is lehet számítani, és az R ki is számítja nekünk, a ROC_adatok objektum auc elemében épp ez a terület van benne. Az AUC rövidítés az Area Under the ROC Curve kódja:

ROC_adatok$auc
## Area under the curve: 0.9021

A ROC görbe alatti terület tehát \(0.9021\), ami nagyon király! Miért is? Hát a legjobb, amit csinálni tudunk, hogy 1 = 100% TPR-t érünk el, 0% FPR mellett! Magyarul felköltözünk a koordináta rendszerben található egységnégyzet bal felső sarkába, azaz a \((0,1)\) pontba. Egy ilyen optimális ROC görbe pedig akkor egy egységnyi oldalú négyzet, aminek területe kereken \(1\). Ehhez az ideális állapothoz pedig elég közel vagyunk, ha a ROC görbe alatti terület \(0.9\)!

Sőt, konkrétan azt is megtudhatjuk, hogy melyik ROC-on vizsgált cut-off esetén vagyunk a legközelebb ehhez az ideális \((0,1)\) ponthoz! Egyszerűen kiszámoljuk a ROC görbe minden x-y kooridináta párjának és a \((0,1)\) pontnak az euklideszi távolságát, és megkeressük hogy a legkisebb távolságot melyik cut-offhoz tartozik.
Szerencsére, mindezt tudja beépítve a pROC csomag:

legjobb_cutoff <- coords(ROC_adatok, "best", ret = "threshold", best.method = "closest.topleft")
## Warning in coords.roc(ROC_adatok, "best", ret = "threshold", best.method =
## "closest.topleft"): The 'transpose' argument to FALSE by default since pROC
## 1.16. Set transpose = TRUE explicitly to revert to the previous behavior,
## or transpose = TRUE to silence this warning. Type help(coords_transpose) for
## additional information.
# árulkodó az utolsó paraméter értéke, hiszen a diagram topfleft pontjához legközelebbi pontjához eső cut-off-ot kutatjuk :)

legjobb_cutoff
##    threshold
## 1 0.04340297

Tehát, a legjobb cut-offunk a ROC görbe alapján a kb. \(4.3\%\). Ez kb. 1%-ponttal alacsonyabb, mint az az \(5.4\%\), amit az \(y=1\) értékek mintabeli relatív gyakorisága alapján mondtunk.
De hát lessük meg, hogy alakul ezzel a cut-offal a klasszifikációs mátrixunk:

csod$BecsultCsod <- 0 # kezdetben mindenki fizetőképes marad
csod$BecsultCsod[csod$CsodVsz > legjobb_cutoff[1,1]] <- 1 # akinek kb. 4.3%-nál magasabb a csődvalószínűsége, azt csődösnek jelöljük
# figyeljünk rá, hogy a legjobb_cutoff egy 1 soros és 1 oszlopos dataframe objektum
# ezért kell rá ilyen kacifántosan, [1,1] módon hivatkozni

table(csod[,c("csod", "BecsultCsod")]) # klasszifikációs mátrix
##     BecsultCsod
## csod    0    1
##    0 1524  338
##    1   16   91

Olybá tűnik, hogy a ROC görbe alapján még megengedhetünk némi veszteséget a csődös (\(y=1\)) osztály precision mutatójából: \(\frac{91}{91+338}=21.2\%\)-ra csökken az értéke.
Ez azért van a ROC görbe szerint, mert a csődös osztály recall-ját sikerül még szebbé javítani: \(\frac{91}{91+16}=85.0\%\)-ra. És a False Positive Rate is elég tolrálható: \(\frac{338}{1524+338}=18.15\%\). Hiszen, jó sok nem csődös vállalat van. Ennyiből megengedheti magának a bank, hogy 18%-nak nem ad hitelt, ha így a valódi csődösök 85%-át helyesen azonosítja!

Az optimális cut-off választás egy alternatív módja az úngynevezett Youden-szabály. Itt azt a cut-offot választjuk a ROC görbén megvizsgált cut-off értékek közül, minél a szenzitivitás + specificitás összeg maximális.

A coords függvénnyel meg tudjuk nézni azt is, hogy mi a Youden-szabály szerint optimális cut-off:

coords(ROC_adatok, "best", ret = "threshold", best.method = "youden")
## Warning in coords.roc(ROC_adatok, "best", ret = "threshold", best.method =
## "youden"): The 'transpose' argument to FALSE by default since pROC 1.16. Set
## transpose = TRUE explicitly to revert to the previous behavior, or transpose
## = TRUE to silence this warning. Type help(coords_transpose) for additional
## information.
##    threshold
## 1 0.04340297

Ugyan az a \(4.3\%\), ami a bal felső sarokhoz legközelebb eső pont keresése során is kijött. Tök jó, két módszer is ugyan azt a cut-offot ajánlja, szóval olyan rossz csak nem lehet! :)

8.2. Cut-Off és az Asszimetrikus Döntés

De nem feltétlenül mindig olyan egyszerű a feladatunk, hogy egy sima ROC görbével el tudjuk dönteni azt, hogy mi az optimális cut-off értékünk!

Tegyük fel, hogy a bankigazgatónk azt mondja, hogy a tapasztalatai szerint egy csődös ügyfélen 10-szer magasabb költsége van átlagosan, mint amennyi bevételt termel neki egy rendesen törlesztő, nem csődös ügyfél. Azaz, ha a Logit modell helyesen azonosít egy csődmentes ügyfelet azon a bank 1 egységnyi pénzt nyer, míg ha a modell elhibázik egy valóságban csődös ügyfelet, azon 10 egység pénzt veszít a bank.
Szóval, ha a helyesen nem csődösnek azonosított ügyfelek számát True Negative-nak, azaz \(TN\)-nek jelöljük, míg a tévesen csődmentesnek ítélt ügyfelek számát False Negative-nak, azaz \(FN\)-nek jelöljük, akkor a bankunk hasznossági függvénye: \[U(TN,FN)=\alpha \times TN + \beta \times FN\] Ahol \(\alpha=1\) és \(\beta=-10\).

Nézzük meg, hogy a jelenlegi ROC görbe alapú \(c=4.3\%\)-os cur-offunk mekkora hasznosságot eredményez a klasszifikációs mátrix alapján:

table(csod[,c("csod", "BecsultCsod")]) # klasszifikációs mátrix
##     BecsultCsod
## csod    0    1
##    0 1524  338
##    1   16   91

A jelenlegi hasznosságunk tehát: \(U(1524, 16)= 1 \times 1524 - 10 \times 16 = 1364\). Lássuk, találunk-e ennél jobb hasznosságot eredméynező cut-offot!

Ehhez először írjunk egy saját függvényt, ami egy előre megadott cut-off mellett kiszámolja nekünk a hasznosságunkat \(\alpha=1\) és \(\beta=-10\) mellett:

# a függvény megadása
hasznossag <- function(x) {
  csod$BecsultCsod <- 0 # kezdetben mindenki fizetőképes marad
  csod$BecsultCsod[csod$CsodVsz > x] <- 1 # akinek kb. x-nél magasabb a csődvalószínűsége, azt csődösnek jelöljük
  klassz_matrix <- table(csod[,c("csod", "BecsultCsod")]) # klasszifikációs mátrix elmentése
  # a True negative = TN értékek a kalsszifikációs mátrix 1. sorában és 1. oszlopában,
  # a False Negative = FN értékek a 2. sorában és 1. oszlopában találhatók
  hasznossag <- 1*klassz_matrix[1,1] - 10*klassz_matrix[2,1]
  return(hasznossag)
}

# a függvény használata úgy, hogy az aktuális cut-off = a ROC görbe által javasolt legjobb cut-off, a 4.3%
hasznossag(legjobb_cutoff[1,1])
## [1] 1364

Tök jó, ugyan az a hasznosságfüggvény értékünk jött ki, mint amit az előbb “kézzel” is kiszámoltunk \(4.3\%\)-os cut-off mellett: \(1364\).

Na, akkor most eresszük rá a hasznossági függvényünkre az R optimize függvényét, ami szépen megadja ennek a hasznosság értéknek a maximumát eredményező, optimális cut-offot. A beállításoknál a második paraméterben a c(0,1)-gyel megadjuk, hogy csak a \([0,1]\) intervallumban keressen cut-off értékeket a gépállat, hiszen más értéktartománynak a csődvalószínűségek körében nincs értelme. A maximum = TRUE beállítással pedig jelezzük, hogy a hasznossag függvényt maximalizálni és NEM minimalizálni szeretnénk:

eredmeny <- optimize(hasznossag, c(0,1), maximum = TRUE)

Az újonnan létrehozott, eredmeny objektumból ki tudjuk olvasni, hogy mi lett az optimális cut-off, és meg tudjuk nézni azt is, hogy mi az a legnagyobb hasznosság érték, amit el tudtunk érni:

eredmeny$maximum # az optimális cut-off
## [1] 0.09458189
eredmeny$objective # a maximális hasznosság éték
## [1] 1437

Ezzel láthatjuk, hogy a jelenlegi preferenciáink mellett a maximális hasznosságot a kb. \(c=9.46\%\)-os cut-off értékel érjük el, ami egy jó \(5\) százalkponttal magasabb, mint a szimmetrikus preferenciák esetén a ROC görbe által javasolt \(4.3\%\).
A \(c=9.46\%\)-os cut-offal pedig \(U=1437\)-es hasznosság értéket lehet elérni, a \(4.3\%\)-as cut-off \(U=1364\) hasznosságához képest.

9. Pszeudo R-négyzet mutató

Láthattuk, hogy a logisztikus regresszió becslési pontosságának megállapítása messze nem egy egyszerű feladat, és rengeteg eszköz van rá (klasszifikációs mátrix, ROC görbe, miegymás) rá.
Viszont, ha nagyon szeretnénk, akkor tudunk egy olyan mutatót is adni a logit modellek becslési pontosságának értékelésére, mint amilyen az OLS modell R-négyzet értéke is volt!

Az ötlet ugyan az, amiről már az előző gyakorlaton is volt szó. Mivel a modell negatív log-likelihood értékének (\(-ln(L)\)) minimalizálásával nyerjük a \(\beta_j\)-et, így lényegében ez az érték hasonlóan jellemzi a logisztikus regressziós modell összes hibáját, mint az \(SSE\) az OLS modell esetében.
Láttuk is, hogy az \(IC\) képletekben le is lehet cserélni gyakorlatilag az \(SSE\)-t \(-ln(L)\)-re.

Egy logit modell \(-2ln(L)\) értékét R-ből lekérhetjük a deviance függvénnyel (amivel OLS modellben az \(SSE\)-t lekértük).

Neg_lnL_targymodell <- deviance(csod_logit)
Neg_lnL_targymodell
## [1] 421.2593

A \(-2\)-nek lesz majd jelentősége, szóval ne osszuk az eredményt 2-vel! A lényegen ez nem változtat! Ezt a \(-2ln(L)\) értéket fel tudjuk használni egy \(R^2\)-szerű pontossági mutató megalkotására a logisztikus regresszióhoz! Ez a mutató lesz a McFadden-féle pszeudo R-négyzet érték!

Induljunk ki a sima OLS esetén használt \(R^2\) képletből:

Azt akkor már többször is megideologizáltuk, hogy logisztikus regresszióban \(ESS=SSE \approx -ln(L)\). De mi lesz a \(TSS=SST\) megfelelője?

Ugye az \(SST\)-re úgy gondolunk folytonos számértékű eredményváltozó esetében, mint az eredményváltozóban lévő összes megmagyarázható információra. De egy olyan értelmezést is adtunk neki, hogy \(SST\) egy olyan modell \(SSE\)-je, amiben nincs magyarázóváltozó. Azaz, \(SST\) nem más, mint a nullmodell \(SSE\)-je.

Na, ezt a második megközelítést lehet örökíteni a logisztikus regresszióra! Mivel a magyarázóváltozók nélküli logit modell \(-2ln(L)\)-jét minden további nélkül ki tudjuk számolni.

Az előző gyakorlaton eljutottunk odáig képletek szintjén, hogy: \(-ln(L) = \sum_{y=1}{ln(p)}+\sum_{y=0}{ln(1-p)}\). Szóval csak ennek a képletnek az eredményét kell majd beszorozni \(2\)-vel.

Namost, egy nullmodellben, ahol nincs magyarázóváltozónk, a \(p=P(y=1)\) valószínűség nem más, mint a megfigyelt mintánkban a \(y=1\) értékek relatív gyakorisága! Stat. II. ismétlés: valószínűség legjobb becslése a (kedvező esetek / összes eset) elven számolt relatív gyakoriság!

Ezt a relatív gyakoriságot a leíró statisztikák alapján simán ki tudjuk számolni:

summary(csod$csod)
##    0    1 
## 1862  107
p_nullmodell <- 107/(107+1862)
p_nullmodell
## [1] 0.05434231

Ez ugye a már ismert \(5.4\%\). A leíró statisztikákból pedig azt is tudjuk, hogy \(107\) db \(y=1\) elemet kell venn a \(-ln(L)\) szummájában \(ln(p)\)-vel, és a maradék \(1862\) db \(y=0\) elemet pedig \(ln(1-p)\)-vel. Ezzel akkor meg is van \(-2ln(L)\) a nullmodellre:

Neg_lnL_nullmodell <- -2*(107*log(p_nullmodell) + 1862*log(1-p_nullmodell))
Neg_lnL_nullmodell
## [1] 831.3419

Miután megvan \(SST\) kistesója is a logit modellben, simán megalkothatjuk a pszeudo R-négyzet mutatónk képletét, ami McFadden elképzelése szerint a fentebb látott \(SST\)-t és \(SSE\)-t használó \(R^2\) formulát tükrözi le a kiszámolt \(-2ln(L)\) értékekkel:

Ezt a képletet a részeredményeinkből könnyen ki tudjuk számolni R-ben:

pszeudo_Rnegyzet <- (Neg_lnL_nullmodell - Neg_lnL_targymodell) / Neg_lnL_nullmodell
pszeudo_Rnegyzet
## [1] 0.4932779

Tehát, a McFadden-féle pszeudo \(R^2\) szerint a modellünk magyarázóereje \(49.3\%\). Ez valami olyasmit jelent, hogy a magyarázóváltozókat használó logit modellünk 49.3%-kal ad jobb becsléseket \(P(y=1)\)-re, mint a magyarázóváltozók nélküli nullmodell.

Természetesen az értéket beépített R függvénnyel is ki tudjuk számolni. A DescTools névre hallgató csomag kell hozzá:

install.packages("DescTools")
library(DescTools) # Szokásos módon ne törődjünk az esetleges Warningokkal! :)
## Warning: package 'DescTools' was built under R version 4.0.5

A függvény pedig PseudoR2 névre hallgat, és egy paraméterben kell megadni, hogy a McFadden-féle verziót akarjuk látni:

PseudoR2(csod_logit, "McFadden")
##  McFadden 
## 0.4932779

Kijött ugyan az a \(49.3\%\), ait mi is számoltunk. Hurrá! :)

Érdemes még megfigyelni, hogy a \(-2ln(L_{targy})\) és \(-2ln(L_{null})\) értékek benne voltak az alap logisztikus regresszió summary() függvénnyel kapott összefoglaló táblájában is:

summary(csod_logit)
## 
## Call:
## glm(formula = csod ~ ., family = binomial(link = "logit"), data = csod)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0316  -0.2365  -0.1288  -0.0762   3.6420  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.1360175  0.8974953  -3.494 0.000476 ***
## eszkfseb100 -0.0034575  0.0012475  -2.772 0.005579 ** 
## stokear100  -0.0462510  0.0103454  -4.471 7.80e-06 ***
## bonitas100  -0.0012246  0.0003241  -3.778 0.000158 ***
## jovedelm100 -0.0106813  0.0024883  -4.293 1.77e-05 ***
## feszkar100   0.0425958  0.0106822   3.988 6.68e-05 ***
## likvid100   -0.0213708  0.0101194  -2.112 0.034697 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 831.34  on 1968  degrees of freedom
## Residual deviance: 421.26  on 1962  degrees of freedom
## AIC: 435.26
## 
## Number of Fisher Scoring iterations: 15

A \(-2ln(L_{targy})\) itt Residual deviance néven, míg a \(-2ln(L_{null})\) pedig Null deviance néven fut az összefoglaló tábla alján. :)

10. A logisztikus regresszió illeszkedésének globális tesztelése

Zárszóként, nézzük meg mire is kellett ezekben a \(-2ln(L)\) értékekben az a \(2\)-es szorzó.

Gondolom többekben felmerült, hogy van-e olyan hipotézisvizsgálat logisztikus regresszióban, mint az \(ANOVA\) a sima OLS modellben, ahol \(H_0: \forall \beta_j=0\). Nos, a válasz az, hogy van, és \(\chi^2\)-féle illeszkedésvizsgálat a neve.

Szóval, ennél a próbánál \(H_0: \forall \beta_j=0\) és \(H_1: \exists j:\beta_j \neq0\) a mintán kívüli világban. Pont, mint a sima OLS modell esetében a globális F-próbában.

A hipotézisvizsgálat prübafüggvénye, és eloszlása igaz \(H_0\) esetén sok-sok mintában:

A képletben \(m\) a modellben lévő magyarázóváltozók száma.

Na, és ezért kell az a \(2\) szorzó a \(-ln(L)\)-ekhez! Mert ennek a próbafüggvénynek, csak a \(2\) szorzóval lesz \(\chi_m^2\) eloszlása igaz \(H_0\) esetén!

Számoljuk is ki a próbafüggvényt:

Probafv <- Neg_lnL_nullmodell - Neg_lnL_targymodell
Probafv
## [1] 410.0826

A szabdságfokunk \(m\) a magyarázóváltozók száma. Ezt simán megkapjuk úgy, hogy \(\beta\) együtthatók száma - 1 (a tengelymetszet miatt):

m <- length(csod_logit$coefficients) - 1
m
## [1] 6

Aztán jöhet is a \(\chi_m^2\) eloszlásból a p-érrték. A próba jobboldali, mint az OLS globális F-próbája, szóval kell az eloszlásfüggvényhez az “1-”:

1 - pchisq(Probafv,m)
## [1] 0

Nagyon jó, a p-érték gyakorlatilag 0, tehát minden szokásos \(\alpha\)-nál kisebb. Szóval egyértelműen elfogadhatjuk a \(H_1\)-et, miszerint a logit modellünkben van legalább 1 \(\beta_j\), ami a sokaságban (mintán kívüli világban) sem nulla!

Ezt a \(\chi^2\)-próbát úgy lehet kiszámoltatni az R-el, ha a nullmodellt legyártjuk egy külön R objektumba, és azt eresztjük össze az aktuális tárgymodellünkkel el anova függvényen belül:

logit_nullmodell <- glm(csod ~ 1, family = binomial (link = "logit"), data = csod)
# figyeljük meg, hogy 1-gyel jelezzük, hogy nem akarunk magyarázóváltozókat a modellben!

anova(logit_nullmodell, csod_logit)
## Analysis of Deviance Table
## 
## Model 1: csod ~ 1
## Model 2: csod ~ eszkfseb100 + stokear100 + bonitas100 + jovedelm100 + 
##     feszkar100 + likvid100
##   Resid. Df Resid. Dev Df Deviance
## 1      1968     831.34            
## 2      1962     421.26  6   410.08

Láthatjuk, hogy a próbafüggvény stimmel azzal a \(410.08\)-cal, amit mi számoltunk!
“ÖrömÉsBódogságVan”! :)