Regresja Logistyczna

Model regresji logistycznej jest szczególnym przypadkiem uogólnionego modelu liniowego. Znajduje zastosowanie, gdy zmienna zależna jest dychotomiczna, to znaczy przyjmuje tylko dwie wartości takie jak na przykład sukces lub porażka, wystąpienie lub brak pewnej jednostki chorobowej, kobieta lub mężczyzna. W zapisie matematycznym wartości te reprezentowane są jako 1 i 0.

Model

Niech \(Y\) będzie binarną zmienną objaśnianą, która z prawdopodobieństwem \(0<\pi(x)<1\) zależnym od wartości zmiennej objaśniającej typu ciągłego \(X\) przyjmuje wartość \(1\). Innymi słowy - prawdziwa jest równość

$\pi(x)=P(Y=1|X=x)=1-P(Y=0|X=x)$

Dla zjawisk najczęściej modelowanych przez statystyków funkcja \(\pi(x)\) jest monotoniczna. Co więcej, zazwyczaj obserwuje się, że taka sama zmiana wartości \(x\) powoduje gwałtowniejszą zmianę prawdopodopodobieństwa \(\pi(x)\), jesli przyjmuje ono wartości bliskie 0.5, niż gdy osiąga wartości bliskie skrajnym (0 lub 1). Wykres tej zależności może więc przpominać literę S. Do jej matematycznego opisu stosuje się funkcję logistyczną określoną dla \(-\infty<t<+\infty\) wzorem

$\textrm{logistic}(t)=\frac{\exp(t)}{1+\exp(t)}$

Przyjmuje ona wartości z przedziału \((0,1)\), a funkcja do niej odwrotna (nazywana logitową) dana jest wzorem

$\textrm{logit}(u)=\textrm{logistic}^{-1}(u)=\log\left(\frac{u}{1-u}\right)$

Dla przykładu rozważmy, jak zmienia się prawdopodobieństwo zakupu auta nowego zamiast używanego w zalezności od zarobków. W przypadku osoby zarabiającej rocznie 1.000.000$ wzrost przychodu o 50.000$ nie zmieni prawdopodobieństwa zakupu nowego auta (bliskiego 100%) tak znacznie, jak zrobiłby to w przypadku osoby zarabiającej 50.000$.

W modelu regresji logistycznej przyjmuje się, że

$\pi(x)=\frac{\exp(\alpha+\beta x)}{1+\exp(\alpha+\beta x)}$

Bezpośrednią konsekwencją tego założenia jest równość

$\textrm{logit}\left(\pi(x)\right)=\log\left(\frac{\pi(x)}{1-\pi(x)}\right)=\alpha+\beta x$

Występujące powyżej wyrażenie \(\frac{\pi(x)}{1-\pi(x)}\) nazywamy szansą (ang. odds). Określa ono stosunek prawdopodobieństwa sukcesu do prawdopodobieństwa porażki i jest monotonicznym przekształceniem prawdopodobieństwa na przedział \((0,+\infty)\). Oznacza to, że gdy rośnie (spada) prawdopodobieństwo, to jednocześnie rośnie (spada) szansa. Możemy stwierdzić, że logarytm szansy (lub logit prawdopodobieństwa) w sposób liniowy zalezy od wartości zmiennej objaśniającej.

Model regresji logistycznej może być stosowany także wtedy, gdy zmienna niezależna jest dychotomiczna.

Interpretacja współczynników

Poprzednia równość implikuje, że

$\frac{\pi(x)}{1-\pi(x)}=e^\alpha e^{\beta x}$

Wynika z tego, że wraz ze wzrostem wartości \(x\) o 1 jednostkę szansa wzrasta \(e^\beta\) razy. Zauważmy, że dla \(\beta=0\) prawdopodobieństwo przyjęcia przez zmienną wynikową \(Y\) wartości 1 nie zależy od zmiennej \(X\) (jest stałe). Jeśli \(\beta>0\), to \(e^\beta >1\), więc przyrost wartości \(x\) wiąże się ze wzrostem prawdopodobieństwa sukcesu. Dla \(\beta<0\) wiąże się ze spadkiem.

Powróćmy jeszcze do opisu wykresu zależności \(\pi(x)\). Obliczając pochodną otrzymujemy, że tangens kąta nachylenia stycznej do wykresu w punkcie \((x, \pi(x))\) wynosi

$\pi'(x)=\beta \pi(x)(1-\pi(x))$

Tak jak oczekiwaliśmy, wykres jest najbardziej stromy, gdy osiągana jest wartość \(\pi(x)=0.5\), co zachodzi dla \(x=- \frac{\alpha}{\beta}\) zwanego poziomem lub dawką medianową (ang. median effective level). Jeśli \(\beta<0\), to wykres przypomina lustrzane odbicie litery S.

Interpretację współczynników dla dychotomicznej zmiennej objaśniającej ułatwi nam pojęcie ilorazu szans (ang. odds ratio). Przypuśćmy, że badamy występowanie pewnego zjawiska w dwóch grupach. Niech dla przykładu prawdopodobieństwo wygrania w pewnej grze wynosi 0,8 dla kobiet i 0,5 dla mężczyzn. Zgodnie z przedstawioną wcześniej formułą szanse wynoszą odpowiednio \(\frac{0.8}{0.2}=4\) oraz \(\frac{0.5}{0.5}=1\). Iloraz szans w naszym przypadku wynosi więc \(\frac{4}{1}=4\). Możemy powiedzieć, że szansa wygranej jest 4 razy większa u kobiet niż u mężczyzn. W ogólności wartości powyżej jedynki oznaczają, że prawdopodobieństwo zajścia danego zjawiska jest większe w pierwszej z porównywanych grup.

W modelu z dychotomiczną zmienną niezależną szansa sukcesu dla poziomu referencyjnego, który został zakodowany jako 0, wynosi \(e^\alpha\). Wartość \(e^\beta\) jest natomiast ilorazem szans, czyli mówi, ile razy wzrasta szansa sukcesu dla poziomu zakodowanego jako 1 w stosunku do poziomu referencyjnego.

Założenia

Do tej pory uwzględnialiśmy tylko jedną zmienną niezależną. Regresję logistyczną można jednak stosować również dla większej liczby predyktorów, o ile są one nieskorelowane. W szczególności umożliwia to rozpatrywanie jakościowych zmiennych objaśniających o więcej niż dwóch poziomach. Dla wielu predyktorów model przyjmuje postać

$P(Y=1|X_1=x_1, \ldots, X_k=x_k)=\pi(x_1, \ldots, x_k)=1-P(Y=0|X_1=x_1, \ldots, X_k=x_k),$

$\pi(x_1, \ldots, x_k)= \frac{\exp(\beta_0+\beta_1 x_1 + \cdots + \beta_k x_k)}{1+\exp(\beta_0+\beta_1 x_1 + \cdots + \beta_k x_k)},$

gdzie \((X_1,\ldots, X_k)\) jest wektorem zmiennych objaśniającyh, a \(Y\) jest binarną zmienną objaśnianą.

Model regresji logistycznej zakłada spełnienie innych warunków niż model regresji liniowej. Nie obowiązują w nim na przykład założenia dotyczące normalności czy homoskedastyczności. Wymagane jest jednak, aby:

• obserwacje były od siebie niezależne;
• logit prawdopododobieństwa zależał w sposób liniowy od zmiennych objaśniających.

Na poprawność i dokładność wnioskowania wpływa również odpowiedni dobór zmiennych niezależnych. Uwzględnione powiny być tylko te, które wykazują istotny wpływ na zmienną \(Y\). Pominięcie którejś z nich skutkuje spadkiem jakości analizy. Zbyt mała liczba obserwacji powoduje ten sam efekt.

Przykład

Powyższe rozważania zilustrujemy teraz przykładem. Przy okazji zaprezenetujemy, jak przeprowadzić regresję logistyczną w pakiecie statystycznym R. Posłużymy się dostępną w sieci i niemal gotową do analizy bazą danych, którą zapiszemy w zmiennej baza.

baza <- read.csv("http://www.karlin.mff.cuni.cz/~pesta/prednasky/NMFM404/Data/binary.csv")

Wyświetlmy 10 pierwszych wierszy.

head(baza, 10)

##    admit gre  gpa rank
## 1      0 380 3.61    3
## 2      1 660 3.67    3
## 3      1 800 4.00    1
## 4      1 640 3.19    4
## 5      0 520 2.93    4
## 6      1 760 3.00    2
## 7      1 560 2.98    1
## 8      0 400 3.08    2
## 9      1 540 3.39    3
## 10     0 700 3.92    2

W pierwszej kolumnie (admit) zawarta jest informacja, czy badany student dostał się na studia drugiego stopnia (1 jeśli tak, 0 w przeciwnym przypadku). Zmienna gre (ang. Graduate Record Exam) odpowiada za liczbę punktów zdobytych w teście końcowym, a w kolumnie gpa (ang. grade point average) zapisana jest średnia ocen. Zmienna rank przyjmuje natomiast wartości 1, 2, 3 oraz 4 - tym niższą, im wyższy jest prestiż ukończonej wcześniej szkoły.

Przed przystąpieniem do analizy przyjrzyjmy się bliżej typom zmiennych.

str(baza)

## 'data.frame':    400 obs. of  4 variables:
##  $ admit: int  0 1 1 1 0 1 1 0 1 0 ...
##  $ gre  : int  380 660 800 640 520 760 560 400 540 700 ...
##  $ gpa  : num  3.61 3.67 4 3.19 2.93 3 2.98 3.08 3.39 3.92 ...
##  $ rank : int  3 3 1 4 4 2 1 2 3 2 ...

Widać, że wartości zmiennej rank zostały zakodowane jako liczbowe. Ponieważ nie mają one sensu liczbowego, lecz traktujemy je jako kategorie, dokonamy konwersji typu na czynnikowy.

baza$rank <- factor(baza$rank)

Przejdźmy do właściwej analizy, której wynik przechowamy w zmiennej reg_log.

reg_log <- glm(admit ~ gre + gpa + rank, data = baza, family = "binomial")

Funkcja summary pozwala na jego prezentację.

summary(reg_log)

## 
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = "binomial", 
##     data = baza)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6268  -0.8662  -0.6388   1.1490   2.0790  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.989979   1.139951  -3.500 0.000465 ***
## gre          0.002264   0.001094   2.070 0.038465 *  
## gpa          0.804038   0.331819   2.423 0.015388 *  
## rank2       -0.675443   0.316490  -2.134 0.032829 *  
## rank3       -1.340204   0.345306  -3.881 0.000104 ***
## rank4       -1.551464   0.417832  -3.713 0.000205 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.98  on 399  degrees of freedom
## Residual deviance: 458.52  on 394  degrees of freedom
## AIC: 470.52
## 
## Number of Fisher Scoring iterations: 4

Model regresji liniowej możemy zatem zapisać jako:

$\log\left( \frac{\pi(x)}{1-\pi(x)}\right)=-3,99+0,002x_{gre}+0,804x_{gpa}-0,675x_{rank2}-1,34x_{rank3}-1,551x_{rank4}$

dla \(x=(x_{gre}, x_{gpa}, x_{rank2}, x_{rank3}, x_{rank4})\), gdzie

• \(x_{gre}\) jest wynikiem uzyskanym w teście,
• \(x_{gpa}\) jest średnią ocen,
• \(x_{rank2}\), \(x_{rank3}\), \(x_{rank4}\) są binarnymi zmiennymi przyjmującymi wartośc 1 dla odpowiedniej kategorii obserwacji.

Kolejna kolumna zawiera błędy standardowe współczynników. Jeśli podzielimy przez nie wartości odpowiadających im współczynników, to otrzymamy wartości tzw. statystyki Walda. Znajdziemy je w następnej kolumnie. Przy założeniu, że współczynnik jest równy 0, czyli odpowiadająca mu zmienna nie ma wpływu na wynik, statystyka ta ma w przybliżeniu standardowy rozkład normalny. Możliwe jest więc wskazanie p-wartości (mamy je w kolejnej kolumnie) i ewentualne odrzucenie hipotezy o nieistotności czynnika. W naszym przykładzie wszystkie p-wartości są małe (<5%), zatem na wynik wpływa każda ze zmiennych.

Estymacja współczynników polega na maksymalizowaniu funkcji wiarygodności, czyli na takim ich dobraniu, by uzyskać największą wartość wyrażenia

$\prod_{i=1}^n\pi(x_i)^{y_i}(1-\pi(x_i))^{1-y_i},$

gdzie \(x_i\) jest wektorem wartości zmiennych objaśniających dla \(i\)-tej obserwacji, a \(y_i\) jest równe 1, jeśli \(i\)-ta obserwacja jest sukcesem. Dla porażki \(y_i\) przyjmuje wartość 0. W praktyce rozwiązuje się równoważny problem polegający na maksymalizacji logarytmu funkcji wiarygodności. W naszym przykładzie program wykorzystał do tego algorytm Fishera. Możemy łatwo odczytać zmaksymalizowaną wartość.

logLik(reg_log)

## 'log Lik.' -229.2587 (df=6)

Mnożąc ją razy -2, otrzymamy dewiancję, która także pozwala ocenić istotność zmiennych uwzględnonych w modelu.

-2*logLik(reg_log)

## 'log Lik.' 458.5175 (df=6)

Mogliśmy ją odczytać już wcześniej (Residual deviance). Była ona zestawiona z dewiancją modelu pozbawionego wszelkich predyktorów (Null deviance). Przy założeniu braku istotności zmiennych dla dużej liczby obserwacji różnica tych dwóch dewiancji ma rozkład \(\chi^2\) o liczbie stopni swobody równej liczbie zmiennych niezależnych w modelu, czyli w naszym przypadku 5. Poniżej obliczymy p-wartość.

with(reg_log, pchisq(null.deviance - deviance, df.null - df.residual, lower.tail = FALSE))

## [1] 7.578194e-08

Z dwóch różnych modeli lepiej dopasowany jest ten, który ma mniejszą dewiancję.

Korzystając z wyestymowanych współczynników możemy obliczyć prawdopodobieństwa dostania się na studia drugiego stopnia dla badanych studentów z bazy danych. W pierwszym kroku obliczmy logity prawdopodobieństw.

baza$log_odds <- predict(reg_log)
head(baza, 10)

##    admit gre  gpa rank    log_odds
## 1      0 380 3.61    3 -1.56712564
## 2      1 660 3.67    3 -0.88484417
## 3      1 800 4.00    1  1.03771175
## 4      1 640 3.19    4 -1.52733046
## 5      0 520 2.93    4 -2.00811132
## 6      1 760 3.00    2 -0.53234576
## 7      1 560 2.98    1 -0.32586874
## 8      0 400 3.08    2 -1.28321604
## 9      1 540 3.39    3 -1.38170577
## 10     0 700 3.92    2  0.07150324

Następnie przekształćmy je używając transformacji logistycznej.

baza$p_stwo <- exp(baza$log_odds)/(1+exp(baza$log_odds))
head(baza, 10)

##    admit gre  gpa rank    log_odds    p_stwo
## 1      0 380 3.61    3 -1.56712564 0.1726265
## 2      1 660 3.67    3 -0.88484417 0.2921750
## 3      1 800 4.00    1  1.03771175 0.7384082
## 4      1 640 3.19    4 -1.52733046 0.1783846
## 5      0 520 2.93    4 -2.00811132 0.1183539
## 6      1 760 3.00    2 -0.53234576 0.3699699
## 7      1 560 2.98    1 -0.32586874 0.4192462
## 8      0 400 3.08    2 -1.28321604 0.2170033
## 9      1 540 3.39    3 -1.38170577 0.2007352
## 10     0 700 3.92    2  0.07150324 0.5178682

Ten sam efekt można uzyskać, specyfikując parametr type w funkcji predict jako “response”.

baza$p_stwo2 <- predict(reg_log, type = "response")
head(baza, 10)

##    admit gre  gpa rank    log_odds    p_stwo   p_stwo2
## 1      0 380 3.61    3 -1.56712564 0.1726265 0.1726265
## 2      1 660 3.67    3 -0.88484417 0.2921750 0.2921750
## 3      1 800 4.00    1  1.03771175 0.7384082 0.7384082
## 4      1 640 3.19    4 -1.52733046 0.1783846 0.1783846
## 5      0 520 2.93    4 -2.00811132 0.1183539 0.1183539
## 6      1 760 3.00    2 -0.53234576 0.3699699 0.3699699
## 7      1 560 2.98    1 -0.32586874 0.4192462 0.4192462
## 8      0 400 3.08    2 -1.28321604 0.2170033 0.2170033
## 9      1 540 3.39    3 -1.38170577 0.2007352 0.2007352
## 10     0 700 3.92    2  0.07150324 0.5178682 0.5178682

Funkcji predict można także użyć do obliczania prawdopodobieństwa w modelu dla nowych danych. Rozpatrzmy dla przykładu zmiany prawdopodobieństwa w zależności od prestiżu ukończonej szkoły dla kogoś, kto uzyskał średni wynik w teście, a jego średnia ocen równa jest medianie. W tym celu stwórzmy najpierw odpowiednią bazę (ramkę) danych i wywołajmy funkcję predict. Zwróćmy przy tym uwagę, aby nazwy kolumn w nowej bazie danych były takie same jak nazwy predyktorów w modelu.

nowa_baza <- with(baza, data.frame(gre = mean(gre), gpa = median(gpa), rank = factor(1:4)))
nowa_baza$p_stwo <- predict(reg_log, newdata = nowa_baza, type = "response")
nowa_baza

##     gre   gpa rank    p_stwo
## 1 587.7 3.395    1 0.5176256
## 2 587.7 3.395    2 0.3532208
## 3 587.7 3.395    3 0.2193133
## 4 587.7 3.395    4 0.1852867