7 Regresja logistyczna

7.1 Model

Regresja logistyczna (ang. logistic regression) jest techniką z rodziny klasyfikatorów liniowych z reprezentacją logistyczną, a formalnie należy do rodziny uogólnionych modeli liniowych (GLM). Stosowana jest wówczas, gdy zmienna wynikowa posiada dwa stany (sukces i porażka), kodowane najczęściej za pomocą 1 i 0. W tej metodzie modelowane jest warunkowe prawdopodobieństwo sukcesu za pomocą kombinacji liniowej predyktorów \(X\).

Ogólna postać modelu \[\begin{align} Y\sim &B(1, p)\\ p(X)=&\E(Y|X)=\frac{\exp(\beta X)}{1+\exp(\beta X)}, \end{align}\] gdzie \(B(1,p)\) jest rozkładem dwumianowym o prawdopodobieństwie sukcesu \(p\), a \(\beta X\) oznacza kombinację liniową parametrów modelu i wartości zmiennych niezależnych, przyjmując, że \(x_0=1\). Jako funkcji łączącej (czyli opisującej związek między kombinacją liniową predyktorów i prawdopodobieństwem sukcesu) użyto logitu. Pozwala on na wygodną interpretację wyników w terminach szans.

Szansą (ang. odds) nazywamy stosunek prawdopodobieństwa sukcesu do prawdopodobieństwa porażki \[\begin{equation} o = \frac{p}{1-p}. \end{equation}\]

Ponieważ będziemy przyjmowali, że \(p\in (0,1)\), to \(o\in (0,\infty)\), a jej logarytm należy do przedziału \((-\infty, \infty)\).

Zatem logarytm szansy jest kombinacją liniową predyktorów \[\begin{equation} \log\left[\frac{p(X)}{1-p(X)}\right]=\beta_0+\beta_1x_1+\ldots+\beta_dx_d. \end{equation}\]

7.2 Estymacja parametrów modelu

Estymacji parametrów modelu logistycznego dokonujemy za pomocą metody największej wiarogodności. Funkcja wiarogodności w tym przypadku przyjmuje postać \[\begin{equation} L(X_1,\ldots,X_n,\beta)=\prod_{i=1}^{n}p(X_i)^Y_i[1-p(X_i)]^{1-Y_i}, \end{equation}\] gdzie wektor \(\beta\) jest uwikłany w funkcji \(p(X_i)\). Maksymalizacji dokonujemy raczej po nałożeniu na funkcję wiarogodności logarytmu, bo to ułatwia szukanie ekstremum. \[\begin{equation} \log L(X_1,\ldots,X_n,\beta) = \sum_{i=1}^n(Y_i\log p(X_i)+(1-Y_i)\log(1-p(X_i))). \end{equation}\]

7.3 Interpretacja

Interpretacja (lat. ceteris paribus - “inne takie samo”) poszczególnych parametrów modelu jest następująca:

  • jeśli \(b_i>0\) - to zmienna \(x_i\) ma wpływ stymulujący pojawienie się sukcesu,
  • jeśli \(b_i<0\) - to zmienna \(x_i\) ma wpływ ograniczający pojawienie się sukcesu,
  • jeśli \(b_i=0\) - to zmienna \(x_i\) nie ma wpływu na pojawienie się sukcesu.

Iloraz szans (ang. odds ratio) stosuje się w przypadku porównywania dwóch klas obserwacji. Jest on jak sama nazwa wskazuje ilorazem szans zajścia sukcesu w obu klasach \[\begin{equation} OR = \frac{p_1}{1-p_1}\frac{1-p_2}{p_2}, \end{equation}\] gdzie \(p_i\) oznacza zajście sukcesu w \(i\)-tej klasie.

Interpretujemy go następująco:

  • jeśli \(OR>1\) - to w pierwszej grupie zajście sukcesu jest bardziej prawdopodobne,
  • jeśli \(OR<1\) - to w drugiej grupie zajście sukcesu jest bardziej prawdopodobne,
  • jeśli \(OR=1\) - to w obu grupach zajście sukcesu jest jednakowo prawdopodobne.

Przykład 7.1 Jako ilustrację działania regresji logistycznej użyjemy modelu dla danych ze zbioru Default pakietu ISLR.

library(ISLR)
head(Default)
##   default student   balance    income
## 1      No      No  729.5265 44361.625
## 2      No     Yes  817.1804 12106.135
## 3      No      No 1073.5492 31767.139
## 4      No      No  529.2506 35704.494
## 5      No      No  785.6559 38463.496
## 6      No     Yes  919.5885  7491.559

Zmienną zależną jest default, a pozostałe są predyktorami. najpierw dokonamy podziału próby na ucząca i testową, a następnie zbudujemy model.

set.seed(2019)
ind <- sample(1:nrow(Default), size = 2/3*nrow(Default))
dt.ucz <- Default[ind,]
dt.test <- Default[-ind,]
mod.logit <- glm(default~., dt.ucz, family = binomial("logit"))
summary(mod.logit)
## 
## Call:
## glm(formula = default ~ ., family = binomial("logit"), data = dt.ucz)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.128e+01  6.169e-01 -18.287   <2e-16 ***
## studentYes  -4.627e-01  2.862e-01  -1.617    0.106    
## balance      5.830e-03  2.883e-04  20.221   <2e-16 ***
## income       9.460e-06  9.833e-06   0.962    0.336    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1967.2  on 6665  degrees of freedom
## Residual deviance: 1046.5  on 6662  degrees of freedom
## AIC: 1054.5
## 
## Number of Fisher Scoring iterations: 8

Tylko income nie ma żadnego wpływu na prawdopodobieństwo stanu Yes zmiennej default. Zmienna balance wpływa stymulująco na prawdopodobieństwo pojawienia się sukcesu. Natomiast jeśli badana osoba jest studentem (studentYes), to ma wpływ ograniczający na pojawienie się sukcesu. Chcąc porównać dwie grupy obserwacji, przykładowo studentów z nie-studentami, możemy wykorzystać iloraz szans.

exp(cbind(OR = coef(mod.logit), confint(mod.logit))) %>% 
    kable(digits = 4)
OR 2.5 % 97.5 %
(Intercept) 0.0000 0.0000 0.0000
studentYes 0.6296 0.3598 1.1060
balance 1.0058 1.0053 1.0064
income 1.0000 1.0000 1.0000

Z powyższej tabeli wynika, że bycie studentem zmniejsza szanse na Yes w zmiennej default o około 37% (w stosunku do nie-studentów). Natomiast wzrost zmiennej balance przy zachowaniu pozostałych zmiennych na tym samym poziomie skutkuje wzrostem szans na Yes o około 0.6%.

Chcąc przeprowadzić predykcję na podstawie modelu dla ustalonych wartości cech (np. student = Yes, balance = $1000 i income = $40000) postępujemy następująco

dt.new <- data.frame(student = "Yes", balance = 1000, income = 40000)
predict(mod.logit, newdata = dt.new, type = "response")
##           1 
## 0.003931398

Otrzymany wynik jest oszacowanym prawdopodobieństwem warunkowym wystąpienia sukcesu (default = Yes). Widać zatem, że poziomy badanych cech sprzyjają raczej porażce.

Jeśli chcemy sprawdzić jakość klasyfikacji na zbiorze testowym, to musimy ustalić na jakim poziomie prawdopodobieństwa będziemy uznawać obserwację za sukces. W zależności od tego, na predykcji jakiego stanu zależy nam bardziej, możemy różnie dobierać ten próg (bez żadnych dodatkowych przesłanek najczęściej jest to 0.5).

pred <- predict(mod.logit, newdata = dt.test, type = "response")
pred.class <- ifelse(pred > 0.5, "Yes", "No")
(tab <- table(pred.class, dt.test$default))
##           
## pred.class   No  Yes
##        No  3211   71
##        Yes   15   37
(acc <- sum(diag(prop.table(tab))))
## [1] 0.9742052

Klasyfikacja na poziomie 97% wskazuje na dobre dopasowanie modelu.