W tym artykule przedstawię sposoby estymacji algorytmu regresji logistycznej do klasyfikacji oraz oceny tej klasyfikacji przy użyciu Pythona i dostępnych dla tego języka biblioteki scikit-learn. Dodatkowo przeprowadzimy estymację w języku R oraz dokonamy interpretacji uzyskanych wyników. Jeżeli nie czytałeś poprzednich artykułów to zapraszam: część 1 oraz część 2.
Regresję będziemy przeprowadzać na zbiorze danych Titanic: cleaned data który pobrałem ze strony Kaggle.com. Jest to zmodyfikowana wersja bardzo popularnego zbioru danych zawierającego informacje dotyczące pasażerów Titanica. Zmienną celu, którą będziemy starali się przewidzieć jest informacja czy pasażer przeżył (1) czy nie przeżył (0) katastrofy. Zbiór danych jest oczyszczony z dużej ilości wartości Null oraz utworzono dodatkową zmienną wielkość rodziny. Poza tym zbiór danych składa się z 14 zmiennych.
Zanim przystąpimy do samej estymacji przeprowadzimy bardzo krótką analizę eksploracyjną zbioru danych oraz przygotujemy go w celu wykorzystania do dalszej pracy.
Krótka analiza eksploracyjna
Importujemy potrzebne nam pakiety:
import pandas as pd import numpy as np import matplotlib.pyplot as plt import os
oraz wczytujemy pobrany zbiór danych:
titanic = pd.read_csv('train_clean.csv')
titanic.head()I naszym oczom ukazuje się wczytany zbiór.
Na co możemy zwrócić od razu uwagę to wartości NaN w kolumnie Cabin, z którymi będziemy musieli sobie jakoś poradzić. Co więcej, zbiór danych zawiera zmienne kategoryczne takie jak płeć oraz Embarked, które oznacza port w którym osoba była zaokrętowana. Takie zmienne będziemy musieli transformować do zmiennych numerycznych ze względu na fakt, że algorytm regresji logistycznej z pakietu Scikit-learn przyjmuje jako dane wejściowe jedynie dane o formacie numerycznym.
Możemy zauważyć również zmienne, które odrazu możemy wyrzucić z zbioru danych ponieważ, że tak napiszę na „zdrowy rozsądek” możemy stwierdzić, że nie mają wartości objaśniającej dla naszej zmiennej celu – przeżycie. Taką zmienną jest imię i nazwisko paseżera. Pozbędziemy się również zmiennej id danego pasażera oraz numer biletu, z tego samego powodu co imię. Można by rozważać czy numer biletu nie przenosi ważnych informacji o kabinie czy w której klasie podróżował pasażer itd. Jednakże takie informacje mamy zapisane w zmiennych Cabin oraz Pclass. Ostatnią zmienną, której się pozbędziemy na tym etapie to zmienna Title. Wskazuje ona na płeć osoby podróżującej jednakże ze względu na to że mamy już zmienną która płeć określa nie potrzebujemy tej w w zbiorze.
Przed pozbyciem się zmiennych nieinformatywych ze zbioru, sprawdzimy jeszcze jego wymiary oraz wartości NaN z kolumny Cabin:
print(titanic.shape)
print('\n')
print(titanic.isna().sum())Otrzymujemy poniższy widok:
Mamy 891 obserwacji oraz 14 zmiennych. Jedynie zmienna Cabin posiada 687 wartości NaN. W stosunku do ilości obserwacji, zmienna ta posiada za dużo wartości NaN i również przenosi zbyt mało informacji dla nas, więc również i jej się pozbędziemy. Należy pamiętać, że to co tutaj robimy to tylko przykładowa estymacja regresji logistycznej. Gdy mierzymy się z rzeczywistymi problemami usuwanie zmiennych ze zbioru danych jest ostatecznością i musimy naprawdę rozważyć czy ta zmienna nam się nie przyda.
Pozbywamy się zmiennych ze zbioru danych:
titanic.drop(columns=['Cabin', 'Name', 'PassengerId', 'Ticket', 'Title'],
inplace=True)Następnie dokonamy binaryzacji zmiennej Sex:
titanic['Sex'] = titanic['Sex'].apply((lambda x: 1 if x=='male' else 0))
Oraz stworzymy dummy variables na podstawie zmiennej Embarked ponieważ przyjmuje 3 wartości i nie możemy po prostu dokonać zwykłego nadania wartości 0, 1, 2, dla kategorii tej zmiennej:
titanic = pd.get_dummies(titanic)
Na tak przygotowanym zbiorze danych możemy zaczynać pracę.
Estymacja w pythonie przy użyciu scikit-learn
Zaimportujmy potrzebne klasy i funkcje do przeprowadzenia regresji:
from sklearn.linear_model import LogisticRegression from sklearn.model_selection import train_test_split
Oprócz zaimportowania samej klasy do szacowania regresji logistycznej, zaimportowałem również klasę train_test_split, dzięki której zbiór danych podzielę na zbiór uczony oraz testowy w celu oceny jego jakości. Zanim dokonamy estymacji musimy przekształcić nasze dane z formatu dataframe z pakietu pandas na typ array z pakietu numpy, ponieważ taki format przyjmuje scikit-learn.
#macierz X bez zmiennej celu X = titanic.loc[:, titanic.columns != 'Survived'].values #wektor y zmienna celu y = titanic['Survived'].values
oraz dokonamy podziału na zbiór uczony oraz testowy w proporcjach 70:30:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.30)
Mając tak podzielony zbiór możemy przystąpić do estymacji regresji logistycznej, inicjujemy instancję klasy LogisticRegression:
lr = LogisticRegression(penalty='l2',
C=1.0,
fit_intercept=True,
solver='lbfgs',
max_iter=1000,
n_jobs=-1)Pakiet scikit-learn dla regresji logistycznej dostarcza nam wiele parametrów, których wartości możemy dobierać w celu lepszego dopasowania naszego modelu do danych. Zazwyczaj dokonuje się tego przy użyciu metody tzw. strojenia hiperparametrycznego (hyperparameter tuning) – czyli wskazania pewnych wartości parametrów a następnie estymowanie modelu dla danego parametru i porównanie jego jakości z pozostałymi modelami (z innymi wartościami parametrów). W ten sposób można dobrać optymalne wartości paramaterów, należy uważać jednakże na możliwość przetrenowania modelu. Najlepiej wydzielić sobie bardzo małą próbkę testową np. 10% całego zbioru i wykorzystać ją do ostatecznego przetestowania jakości modelu – już po strojeniu hiperparametrycznym.
Wszystkie parametry dostępne są w dokumentacji pakietu, tutaj skupię się na najważniejszych.
- C – parametr ten określa odwrotną „siłę” regularyzacji, dla kogoś kto nie wie czym jest regularyzacja – w skrócie jest to najczęściej nałożenie kary na estymowane wartości parametrów w celu ograniczenia nadmiernego dopasowania modelu do danych, w tym wypadku parametr C przyjmuje wartości od 0 do inf, im niższe wartości tym siła regularyzacji jest większa im wyższe tym mniejsza my zostawiamy wartość docelową czyli 1.0 brak regularyzacji, ten parametr jest wskazywany jako hiperparametr przez twórców scikita,
- penalty – parametr określający jaki rodzaj kary ma być nałożony przy estymacji parametrów, na pewno na blogu pojawi się temat regularyzacji więc nie będę szczegółowo opisywał czym się różnią poszczególne rodzaje kar, do wyboru mamy kary: l1, l2 i elasticnet – wybieram domyślną czyli l2,
- fit_intercept – parametr określa czy chcemy dodać wyraz wolny do naszego modelu – wybieramy, że chcemy go dodać
- solver – jest to wybór algorytmu optymalizacyjnego czyli szukającego maksimum naszej funkcji wiarygodności w modelu, dokładny opis czym się różnią poszczególne solvery i jakie opcje możemy z nimi wybierać również znajduje się w dokumentacji – nie jest to aż tak istotne w przypadku estymacji tego modelu
- max_iter – za pomocą tego parametru określamy liczbę iteracji naszego algorytmu szukającego maksimum funkcji wiarygodności – najlepiej wybrać większą liczbę iteracji żeby być pewnym, że algorytm zbiegł do maksimum, jeżeli nie osiągnie zbieżności zostaniemy o tym poinformowani,
- n_jobs – za pomocą tego parametru ustalamy ile rdzeni jednocześnie ma przetwarzać algorytm na naszym komputerze, wartość -1 oznacza, że mają być użyte wszystkie dostępne.
Okej, zdefiniowaliśmy nasz model więc załadujmy nasze dane uczące – nic prostszego:
lr.fit(X_train, y_train)
Otrzymaliśmy trochę rozjechany obrazek ale pokazujący jak zdefiniowaliśmy nasz model, dodatkowo nie wyskoczyły żadne ostrzeżenia to znaczy że algorytm osiągnął zbieżność. Teraz zobaczmy jaką dokładność dopasowania osiągnęliśmy na zbiorze testowym. W tym celu, musimy dokonać predykcji na podstawie macierzy X zbioru testowego. A następnie otrzymane przewidziane wartości porównamy do prawdziwych wartości z wektora y_test.
from sklearn.metrics import accuracy_score #importujemy miarę dokładności y_pred = lr.predict(X_test) #dokonujemy predykcji na zbiorze testowym print(accuracy_score(y_test,y_pred)) #printujemy wyniki
Otrzymaliśmy wynik w przybliżeniu 0.77 dokładności predykcji. Nie tak źle jednakże zanim będziemy się cieszyć musimy sprawdzić rozkład wartości naszej zmiennej celu. Miara dokładności jest narażona na pewne błędy lub może bardziej niedopatrzenia. Mianowicie, jeżeli rozkład klas jest bardzo nierównomierny w przypadku modelowanej zmiennej to wyniki miary dokładności mogą być zafałszowane. Wyobraźmy sobie sytuację gdy modelujemy oszustwa na kartach kredytowych. W przypadku takiego zbioru danych zwykle rozkład wartości to 99% obserwacji należących do klasy 0 czyli brak oszustwa natomiast 1% do klasy 1 czyli oszustwo. Gdy estymujemy nasz model i osiągnie on dokładność predykcji na poziomie 99% to nie będzie to powód do cieszenia się bo równie dobrze możemy każdą obserwację klasyfikować jako brak oszustwa i osiągniemy taki sam poziom dokładności a użyteczność modelu jest zerowa. Dalej pokażę jak sobie z tym radzić za pomocą macierzy pomyłek oraz miar precyzji i swoistości oraz krzywej ROC. Ale na razie sprawdźmy procentowy rozkład zmiennej celu:
print(titanic.Survived.value_counts()) #wartości rozkładu zmiennej
print("\n")
#procent obserwacji kategorii 1
print("{:.3}".format(titanic.Survived.value_counts()[1]/titanic.shape[0]*100))
Jak widać rozkład jest nierównomierny. Około 38.4% obserwacji przynależy do klasy 1. W tym momencie musimy użyć miar precision oraz recall żeby lepiej ocenić jakość naszego modelu. Importujemy więc miary precision, recall, f1 oraz macierz pomyłek:
from sklearn.metrics import precision_score, recall_score, f1_score from sklearn.metrics import confusion_matrix
Następnie obliczymy macierz i ładnie ją wyprintujemy:
tn, fp, fn, tp = confusion_matrix(y_true=y_test, y_pred=y_pred).ravel()
print("True Positives: \t{}\t False Positives: \
\t{}\t \nFalse Negatives: \t{}\t \
True Negatives: \t{}\t".format(tp, fp, fn, tn))
Otrzymując macierz pomyłek mamy już lepszy pogląd na ocenę jakości naszego modelu. Charakterystyczna jest dosyć duża liczna False Negatives czyli obserwacji, które przynależą rzeczywiście do klasy Pozytywnej (1) a nasz model zaklasyfikował je jako negatywne. Widzimy, że nasz model dosyć słabo radzi sobie z przewidywaniem osób które przeżyły na podstawie danego wektora cech. Zobaczmy na miary precyzji, czułości i f1:
print('Precyzja: %.3f' % precision_score(y_true=y_test, y_pred=y_pred))
print('Czułość: %.3f' % recall_score(y_true=y_test, y_pred=y_pred))
print('F1: %.3f' % f1_score(y_true=y_test, y_pred=y_pred))
Miara czułości wskazuje nam w jakim procencie obserwacje rzeczywiście przynależące do klasy 1 (osoba, która przeżyła) zostały przewidziane jako takie obserwacje. W naszym przypadku wynik jest dosyć niski 65.7% obserwacji należących do klasy 1 rzeczywiście zostało tak zaklasyfikowanych.
Z kolei precyzja wskazuje na procent obserwacji zaklasyfikowanych jako pozytywne i które również zostały potwierdzone jako obserwacje pozytywne. Miara ta pozwala nam zobrazować na ile możemy ufać przewidywaniom klasy pozytywnej przez nasz model. Wynik 78.9% nie jest najgorszy.
Miara F1 stanowi połączenie miar precyzji i czułości. Jest to wygodna forma porównania jakości kilku modeli ponieważ otrzymujemy jedną zsyntetyzowaną wartość opartą o miary precyzji i czułości.
Ostatnią metodą sprawdzenia jakości modeli, którą chciałbym pokazać jest krzywa ROC i miara pola pod tą krzywą czyli AUC. Narzędzie to jest świetne do wizualnej oceny jakości modelu oraz porównywania jakości kilku modeli. W takim razie zobrazujemy sobie krzywą dla naszego modelu.
#potrzebne pakiety
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, auc
from scipy import interp
import seaborn as sns
#obliczenie prawdopodobieństw modelu
proba = lr.predict_proba(X_test)
#wygladzanie krzwej
from scipy.ndimage.filters import gaussian_filter1d
#styl
sns.set(font_scale=1.6)
sns.set_style("darkgrid")
plt.figure(figsize=(10, 10))
#ROC i AUC wykres
fpr, tpr, threshold = roc_curve(y_test, proba[:,1])
roc_auc = auc(fpr, tpr)
plt.plot(gaussian_filter1d(fpr, sigma=0.9),
gaussian_filter1d(tpr, sigma=0.9),
lw=3,color='tomato',
label='%s (AUC = %0.2f)'
% ('Regresja Logistyczna', roc_auc))
plt.xlim([-0.02, 1.0])
plt.ylim([0, 1.02])
plt.xlabel('Odsetek fałszywie pozytywnych')
plt.ylabel('Odsetek prawdziwie pozytywnych')
plt.legend(loc="lower right");
Ogólna zasada kciuka w przypadku interpretacji krzywej ROC jest taka, że im model jest lepszy tym krzywa będzie zbliżała się do lewego górnego rogu. Dodatkowo w legendzie umieściłem wartość miary AUC czyli pola pod krzywą. Dla idealnego modelu wartość ta będzie równa 1. W szczególnym przypadku, gdy mamy binarną zmienną celu oraz zbiór ma równy rozkład klas (tzn. 50% obserwacji należy do jednej klasy) to możemy sobie wyobrazić prosty model losowy, którego zadaniem byłoby przewidywanie tylko jednej klasy dla całego zbioru danych. Krzywa ROC takiego modelu byłaby prostą przechodzącą od punktu (0,0) do punktu (1,1). W przypadku analizowania takich zbiorów często umieszcza się taką linię na wykresach ROC, żeby dodatkowo umieścić punkt odniesienia.
Regresja logistyczna w języku R wraz z interpretacją parametrów
Dalej przeprowadzimy estymację modelu z wykorzystaniem języka R i dostępnych dla niego bibliotek. Użyjemy tego samego zbioru, dlatego zbiór w którym pozbyłem się zmiennych itd zapisałem i wykorzystałem od razu do estymacji w R. Dodatkowo zmienną family_size zamieniłem na zmienną binarną – czyli czy dana osoba posiada dzieci czy nie posiada.
#pakiet zapewniający dodatką funkcjonalność dla modeli statystycznych
install.packages("aod")
library(aod)
#wczytanie zbioru danych
titanic <- read.csv("titanic_blog.csv")
#funkcja określająca czy posiada dziecko czy nie
have_children <- function(x){
if(x==0){
return(0)
} else {
return(1)
}
}
#aplikowanie funkcji
titanic$Family_Size <- sapply(titanic$Family_Size, have_children)
#estymacja modelu
logit <- glm(Survived ~., data=titanic, family="binomial")
#podsumowanie modelu
summary(logit)
#przedziały ufności
confint(logit)
#szanse
exp(coef(logit))Po kolei opiszę co oznacza dana linijka kodu. Zainstalowałem oraz wczytałem dodatkowy pakiet „aod” zapewnia on dodatkowe narzędzia do modeli statystycznych takich jak test walda czy np. przedziały ufności które tutaj wykorzystałem („confint”). Następnie utworzyłem prostą funkcje, którą wykorzystałem do transformacji zmiennej porządkowej family_size na zmienną binarną. Dalej, wykorzystałem wbudowaną w R funkcję glm, którą wykorzystuje się do budowania modeli statystycznych. Za pomocą family=”binomial” wskazuje, że ma być to model logitowy. Dalej otrzymujemy podsumowanie naszego modelu.
Gwiazdki przy danych zmiennych wskazują nam statystyczną istotność oszacowania danego parametru przy danym poziomie istotności. Widzimy również wartość kryterium informacyjnego Akaike – AIC, miara ta pomaga w przypadku porównania jakości dwóch modeli.
Pakiet R standardowo jako oszacowane wartości parametrów w kolumnie Estimate zwraca logarytmy szans. Jeżeli chcielibyśmy otrzymać proste do interpretacji iloraze szans musimy po prostu użyć funkcji exp na wartościach oszacowanych parametrów naszych zmiennych. Przykład ze zmienną Age:
exp(coef(logit))['Age']
Wynik ten możemy zinterpretować następująco: Osoba o rok starsza ma o 0.96 razy mniejsze szanse na przeżycia niż osoba od niej o rok młodsza. Na podstawie takiego wyniku możemy np. wnioskować, że osoby młodsze mogły być ewakuowane pierwsze itd. Podobnie możemy interpretować pozostałe zmienne, oraz stosować inne miary jak ryzyko, które opisałem w pierwszym artykule. Mając oszacowane logarytmy szans należy jedynie dokonywać przekształceń matematycznych, żeby otrzymać miarę, której potrzebujemy.
Bardzo często wartości logarytmów szans podaje się z przedziałami ufności, ma to na celu rozpoznanie, które zmienne są istotne oraz w jakich przedziałach się znajdują. Aby otrzymać takie ilorazy z przedziałami zastosujemy następującą komendę:
exp(cbind(OR = coef(logit), confint(logit)))
Otrzymaliśmy Odds Ratio wraz z 95% przedziałem ufności. Przypominam zasadę kciuka, jeżeli przedział zawiera 1 to oszacowany parametr nie jest istotny statystycznie na wybranym poziomie istotności. W ten sposób otrzymujemy ładną i spójną tabelę parametrów gotowych do interpretacji.
Była to ostatnia część z serii artykułów poświęconych algorytmowi regresji logistycznej. Oczywiście, nie wszystkie zagadnienia z zakresu tego modelu zostały tutaj poruszone a jedynie główny opis teoretyczny jak i praktyczne zastosowanie. Bardziej dociekliwym tematu gorąco polecam pozycję Logistic Regression – Kleinbaum.