Przed kilkoma tygodniami zainicjowałem na blogu cykl: „zaawansowana analityka”. Czas więc wytoczyć duże działa. Użyję ich do spekulacji nad sensownością połączenia XGBoost + scoring kredytowy z uwzględnieniem pewnych założeń. Zapraszam! 🙂
Pod koniec marca rozpocząłem na kanale Data Science Plus omawianie różnych podejść do poprawy modelu regresji logistycznej. Założyłem pełną interpretowalność, dlatego w kolejnych odcinkach (DPS 002, DSP 003, DSP 004) omawiałem relatywnie proste techniki, jak manipulacja punktem odcięcia, czy kategoryzacja zmiennych.
Co istotne, cały czas skupiam się na aspekcie finansowym, jaki niesie potencjalne wdrożenie budowanego modelu. To zysk banku jest tu kluczową metryką. Jako metrykę pomocniczą używam Gini. Jeśli nie widziałeś wspomnianych odcinków, to zachęcam Cię do ich obejrzenia.
A co by było, gdyby...¶
Puśćmy wodze wyobraźni. Co by było, gdyby XGBoost był w 100% interpretowalny? Albo, gdyby powstał jego całkowicie interpretowalny odpowiednik? A może już istnieją metody pozwalające na jego wystarczającą interpretację? Na jaki zysk moglibyśmy liczyć? Jak poprawiłaby się rentowność instytucji używających dziś regresji logistycznej? Czy dobry model regresji logistycznej daje dużo gorsze wyniki od modelu zbudowanego z użyciem XGBoost? Na te i inne pytania będę starał się odpowiedzieć w tym i kilku kolejnych wpisach. 🙂
Kilka wyjaśnień na początek¶
Czemu XGBoost, a nie inna drzewiasta implementacja korzystająca z Gradient Boostingu?
Bo jeśli chodzi o dane zapisane w formie tabelarycznej, to na Kaggle niepodzielnie króluje ensembling (źródła w sekcji "Linki" na końcu wpisu), którego najpopularniejszym reprezentantem jest właśnie XGBoost (choć coraz większą popularność zdobywa LightGBM). Nie chcę wchodzić w rozważania, czy dla tego konkretnego problemu lepiej sprawdzi się XGBoost, GradientBoostingClassifier, czy LightGBM. Zakładam, że XGBoost jest godnym reprezentantem metod drzewiastych w starciu z "dziadkiem" (historia sięgająca XIX wieku) o nazwie regresja logistyczna.Czy będę "żyłować" model?
Tylko trochę. Mocno ogranicza mnie czas. Cały wpis (tekst + kod) powstaje w dwa wieczory jednego majowego deszczowego i zimnego weekendu. 😉Czy algorytm, który wygra, jest bezsprzecznie lepszy?
Oczywiście nie. Jest to nic innego jak luźne porównanie wyników dwóch modeli korzystających z różnych algorytmów. 😉 Gdyby regresja dawała lepsze wyniki niż ensembling, to królowałaby na Kaggle. A gdyby ensembling był w 100% interpretowalny i stabilny, to korzystałby z jego dobrodziejstwa każdy bank. Proste. 🙂Czemu pomijasz EDA?
Znam omawiany zbiór bardzo dobrze. Przyglądałem mu się bliżej np. podczas tego projektu.
Wystarczy tego pisania, czas przejść do modelowania. Zaczynamy! 🙂
Nie wiesz jak rozpocząć analizę danych? Sprawdź ten wpis, w których pokauję 5 prostych kroków rozpoczynających analizę danych.
1. Import bibliotek.¶
from boruta import BorutaPy
from creme_de_la_creme_of_ML import mean_encoding
import numpy as np
import pandas as pd
from scipy.stats import randint, uniform
from sklearn.metrics import roc_auc_score, classification_report, confusion_matrix, accuracy_score, recall_score, precision_score, roc_curve
from sklearn.model_selection import RandomizedSearchCV, cross_val_score
from sklearn.tree import DecisionTreeClassifier
from xgboost import XGBClassifier
2. Wczytanie zbioru danych.¶
x_tr = pd.read_csv('data/x_tr.csv', index_col = 0).drop(columns = 'intercept')
x_te = pd.read_csv('data/x_te.csv', index_col = 0).drop(columns = 'intercept')
x_va = pd.read_csv('data/x_va.csv', index_col = 0).drop(columns = 'intercept')
y_tr = pd.read_csv('data/y_tr.csv', header = None, index_col = 0)
y_te = pd.read_csv('data/y_te.csv', header = None, index_col = 0)
y_va = pd.read_csv('data/y_va.csv', header = None, index_col = 0)
Zbiorów: x_tr, x_va użyję do rozwijania modelu. Finalny test przeprowadzę na zbiorze x_te.
3. Modelowanie.¶
3.1. Bazowy model.¶
Charakterystyka:
- dobór zmiennych: brak,
- dobór parametrów: brak,
- operacje na zbiorze: brak.
model_1 = XGBClassifier()
model_1.fit(x_tr, y_tr.values.ravel())
pred_1 = model_1.predict_proba(x_va)
gini_1 = 2 * roc_auc_score(y_va, pred_1[:, 1]) - 1
print('Gini score:', gini_1.round(3))
Podejrzyjmy istotność zmiennych, by sprawdzić:
- Czy jest zgodna z intuicją?
- Czy nie ma żadnych zbędnych zmiennych?
pd.DataFrame(model_1.feature_importances_, index = x_tr.columns).sort_values(by = 0, ascending = False)
3.2. Model 2.¶
Charakterystyka:
- dobór zmiennych: brak,
- dobór parametrów: brak,
- operacje na zbiorze: część zmiennych z kodowaniem MeanEncoding.
Po raz perwszy na łamach bloga będę korzystać z MeanEncodingu. Całą ideę opiszę w jednym z kolejnych wpisów, dlatego pominę teraz jakiekolwiek wyjaśnienia tej metody. 🙂
Sprawdzę jeszcze, które zmienne warto poddać procesowi zmiany kodowania. Chcę się skupić jedynie na zmiennych kategorycznych. Zmienne liczbowe, takie jak "lim_kredytu" na ten moment pomijam.
x_tr.head()
zmienne_me = ['plec', 'opozn_plat_wrz', 'opozn_plat_sie', 'opozn_plat_lip', 'opozn_plat_cze', 'opozn_plat_maj', 'opozn_plat_kwi',
'wyksztalcenie_nieznane', 'wyksztalcenie_srednie', 'wyksztalcenie_wyzsze', 'wyksztalcenie_wyzsze_pelne', 'stan_cywilny_kawaler_panna', 'stan_cywilny_nieznany', 'stan_cywilny_w_zwiazku']
dfs = mean_encoding([x_tr, x_va, x_te], y_tr, zmienne_me, 1.0)
model_2 = XGBClassifier()
model_2.fit(dfs[0], y_tr.values.ravel())
pred_2 = model_2.predict_proba(dfs[1])
gini_2 = 2 * roc_auc_score(y_va, pred_2[:, 1]) - 1
print('Gini score:', gini_2.round(3))
A więc lekki spadek. Spróbuję więc zmienić kodowanie pozostałych zmiennych numerycznych.
3.3. Model 3.¶
Charakterystyka:
- dobór zmiennych: brak,
- dobór parametrów: brak,
- operacje na zbiorze: część zmiennych z kategoryzacją pojedynczym drzewem.
def categorize_data(x, y):
"""Tree based data categorization.
Parameters:
-----------
x : Pandas Series, feature to categorize
y : Pandas Series, target variable
Returns:
--------
dt_model : sklearn.tree.DecisionTreeClassifier, fitted decision tree model
categories : numpy.ndarray, list of categories
"""
dt_model = DecisionTreeClassifier(max_depth = 3)
dt_model.fit(x.values.reshape(-1, 1), y)
categories = dt_model.apply(x.values.reshape(-1, 1))
return dt_model, categories
zmienne_do_kategoryzacji = ['lim_kredytu', 'wiek', 'kwota_wyciagu_wrz', 'kwota_wyciagu_sie', 'kwota_wyciagu_lip', 'kwota_wyciagu_cze',
'kwota_wyciagu_maj', 'kwota_wyciagu_kwi', 'platnosc_wrz', 'platnosc_sie', 'platnosc_lip', 'platnosc_cze',
'platnosc_maj', 'platnosc_kwi']
for zmienna in zmienne_do_kategoryzacji:
model, kategorie = categorize_data(x_tr[zmienna], y_tr)
x_tr = x_tr.assign(_kat = kategorie)
x_va = x_va.assign(_kat = model.apply(x_va[zmienna].values.reshape(-1, 1)))
x_te = x_te.assign(_kat = model.apply(x_te[zmienna].values.reshape(-1, 1)))
nazwa_zmiennej = zmienna + '_kat'
x_tr.rename(columns = {'_kat' : nazwa_zmiennej}, inplace = True)
x_va.rename(columns = {'_kat' : nazwa_zmiennej}, inplace = True)
x_te.rename(columns = {'_kat' : nazwa_zmiennej}, inplace = True)
x_tr.drop(columns = zmienna, inplace = True)
x_va.drop(columns = zmienna, inplace = True)
x_te.drop(columns = zmienna, inplace = True)
model_3 = XGBClassifier()
model_3.fit(x_tr, y_tr.values.ravel())
pred_3 = model_3.predict_proba(x_va)
gini_3 = 2 * roc_auc_score(y_va, pred_3[:, 1]) - 1
print('Gini score:', gini_3.round(3))
Ponownie pogorszenie pierwotnego wyniku. Sprawdzę, jak zmienne kategoryzowane sprawdzą się w połączeniu z MeanEncoding.
Szukasz czegoś prostszego? Sprawdź ten projekt: Klasyfikacja wniosków o wydanie karty kredytowej
3.4. Model 4.¶
Charakterystyka:
- dobór zmiennych: brak,
- dobór parametrów: brak,
- operacje na zbiorze: część zmiennych z kategoryzacją pojedynczym drzewem, a nastepnie MeanEncoding.
zmienne_me = x_tr.columns.tolist()
dfs = mean_encoding([x_tr, x_va], y_tr, zmienne_me, smooth=1)
model_4 = XGBClassifier()
model_4.fit(dfs[0], y_tr.values.ravel())
pred_4 = model_4.predict_proba(dfs[1])
gini_4 = 2 * roc_auc_score(y_va, pred_4[:, 1]) - 1
print('Gini score:', gini_4.round(3))
Niestety jeszcze jedno pogorszenie wyniku. Za kulisami testowałem inne wartości parametru smooth i niewiele to zmieniało. Skoro pierwszy, bazowy model poradził sobie z zadaniem najlepiej, to istnieje szansa, że zbiór jest relatywnie prosty i nie ma sensu wyciagać na niego armaty.
Sprawdzę jeszcze jak zachowa się model z innymi parametrami XGBoost.
3.5. Model 5.¶
Charakterystyka:
- dobór zmiennych: brak,
- dobór parametrów: RandomizedSearchCV,
- operacje na zbiorze: część zmiennych z kategoryzacją pojedynczym drzewem.
parametry_xgb = {
'n_estimators' : randint(50, 500),
'max_depth' : randint(2, 10),
'learning_rate' : uniform(loc = 0, scale = 1),
'gamma' : uniform(loc = 0, scale = 1),
'min_child_weight' : randint(1, 5),
'subsample': uniform(loc = 0, scale = 1),
'colsample_bytree': uniform(loc = 0, scale = 1),
'reg_alpha' : uniform(loc = 0, scale = 1),
'reg_lambda' : uniform(loc = 0, scale = 1),
'booster' : ['gbtree', 'gblinear', 'dart'],
'n_jobs' : [2]
}
rs = RandomizedSearchCV(XGBClassifier(), parametry_xgb, n_iter = 20, n_jobs = 2, cv = 5)
rs.fit(x_tr, y_tr.values.ravel())
wybrane_parametry = rs.best_params_
wybrane_parametry = {'booster': 'gbtree',
'colsample_bytree': 0.5349522340521695,
'gamma': 0.3349906283498246,
'learning_rate': 0.13815586195893292,
'max_depth': 2,
'min_child_weight': 2,
'n_estimators': 54,
'n_jobs': 2,
'reg_alpha': 0.4989236903863318,
'reg_lambda': 0.34314930043475145,
'subsample': 0.6568640048717483}
model_5 = XGBClassifier(**wybrane_parametry)
model_5.fit(x_tr, y_tr.values.ravel())
pred_5 = model_5.predict_proba(x_va)
gini_5 = 2 * roc_auc_score(y_va, pred_5[:, 1]) - 1
print('Gini score:', gini_5.round(3))
Wreszcie jest znaczna poprawa. 🙂 Sprawdzę jeszcze istotność zmiennych - być może w zbiorze są teraz jakieś nadmiarowe zmienne.
pd.DataFrame(model_5.feature_importances_, index = x_tr.columns).sort_values(by = 0, ascending = False)
Jak się okazuje 3 zmienne są nieistotne dla modelu - jak widzisz, wykształcenie wyższe i wyższe pełne okazują się bezwartościowe. 😉 Usuwam je i jeszcze raz sprawdzam jakość.
do_usuniecia = ['kwota_wyciagu_lip_kat', 'wyksztalcenie_wyzsze_pelne', 'wyksztalcenie_wyzsze']
model_5.fit(x_tr.drop(columns = do_usuniecia), y_tr.values.ravel())
pred_5 = model_5.predict_proba(x_va.drop(columns = do_usuniecia))
gini_5 = 2 * roc_auc_score(y_va, pred_5[:, 1]) - 1
print('Gini score:', gini_5.round(3))
3.6. Model 6.¶
Charakterystyka:
- dobór zmiennych: Boruta,
- dobór parametrów: RandomizedSearchCV,
- operacje na zbiorze: część zmiennych z kategoryzacją pojedynczym drzewem.
Skoro usunięcie kilku zmiennych dało niewielką poprawę, to warto chyba przetestować nieco bardziej zaawansowaną metodę doboru zmiennych niż feature importance, czy forward selection. 🙂 Decyduję się na sprawdzony algorytm: Boruta. Podobnie jak w przypadku MeanEncoding nie będę go teraz omawiać. Zrobię to w oddzielnym artykule, bo jest o czym pisać. 😉
model_6 = XGBClassifier(**wybrane_parametry)
feat_selector = BorutaPy(model_6, n_estimators='auto', random_state=1)
feat_selector.fit(x_tr.values, y_tr.values.ravel())
Pobieram listę zmiennych.
finalne_zmienne = x_tr.columns[feat_selector.support_]
Sprawdzam ranking zmiennych.
pd.DataFrame(feat_selector.ranking_, index = x_tr.columns, columns = ['ranking_zmiennych']).sort_values('ranking_zmiennych', ascending = True)
Uczę ponownie model na finalnym zestawie zmiennych.
model_6.fit(x_tr[finalne_zmienne], y_tr.values.ravel())
pred_6 = model_6.predict_proba(x_va[finalne_zmienne])
gini_6 = 2 * roc_auc_score(y_va, pred_6[:, 1]) - 1
print('Gini score:', gini_6.round(3))
I jest kolejna poprawa wyniku! 🙂 Sprawdzę jeszcze zachowanie modelu na zbiorze testowym, którego wcześniej nie używałem.
pred_te = model_6.predict_proba(x_te[finalne_zmienne])
gini_te = 2 * roc_auc_score(y_te, pred_te[:, 1]) - 1
print('Gini score:', gini_te.round(3))
Stabilność jest więcej niż zadowalająca. 🙂 Ten mini-projekt uznaję za mały sukces, bo początki nie były łatwe. Kategoryzacja zmiennych i MeanEncoding nic nie wniosły, a wręcz pogarszały bazowy wynik. Dopiero dobór parametrów modelu w połączeniu z zaawansowanym algorytmem doboru zmiennych znacząco poprawiły wynik.
pd.DataFrame({'model_1':['brak','brak','brak', 0.55],
'model_2':['brak','brak','mean_encoding', 0.549],
'model_3':['brak','brak','kategoryzacja_drzewem', 0.543],
'model_4':['brak','brak','kategoryzacja_drzewem + mean_encoding', 0.536],
'model_5':['feature_importance','RandomizedSearchCV','kategoryzacja_drzewem', 0.579],
'model_6':['Boruta','RandomizedSearchCV','kategoryzacja_drzewem', 0.584]}, index = ['metoda_doboru_zmiennych', 'metoda_doboru_parametrow', 'operacje_na_zbiorze', 'gini_va']).T
4. Dobór punktu odcięcia.¶
Zgodnie ze schematem obranym w tym wpisie spróbuję dobrać punkt odcięcia tak, by maksymalizować zysk banku.
punkty_odciecia = np.arange(0, 1.01, 0.01)
wyniki = {}
for punkt in punkty_odciecia:
tn, fp, fn, tp = confusion_matrix(y_te, pred_te[:, 1]>punkt).ravel()
koszt = tn * -15 + fn * 75
wyniki[punkt] = koszt
# Wyniki zapisuję w serii danych.
wyniki = pd.Series(wyniki)
najnizszy_koszt = wyniki.iloc[wyniki.argmin()]
najlepszy_punkt_odciecia = wyniki.iloc[wyniki.argmin():wyniki.argmin()+1].index[0]
print('Najmniejszy koszt wynosi: {}. Odpowiadający mu punkt odcięcia wynosi: {}.'.format(najnizszy_koszt, najlepszy_punkt_odciecia))
5. Szacowanie rentowności modelu na zbiorze testowym.¶
Przyjmuję dokładnie takie same założenia, jak te omawiane w tym miejscu. W skrócie koszty kolejnych decyzji wyglądają następująco:
- True Negative - udzielamy kredytu, a klient go spłaca z odsetkami. Zakładam, że średni koszt = -10 (bank zyskuje 10%).
- True Positive - wykrycie nierzetelnego klienta sprawia, że nic nie stracimy, lecz też nic nie zyskujemy. Koszt = 0.
- False Negative - najgorszy z rozpatrywanych scenariuszy. Zakładam średni koszt = 75 (bank traci 75% z wartości oczekiwanej kredytu).
- False Positive - prócz niewykorzystanych szans koszt = 0.
def podsumowanie_modelu(punkt_odciecia):
tn, fp, fn, tp = confusion_matrix(y_te, pred_te[:,1]>punkt_odciecia).ravel()
specyficznosc = tn / (tn+fp)
czulosc = tp / (tp + fn)
return tn, fp, fn, tp
tn, fp, fn, tp = podsumowanie_modelu(najlepszy_punkt_odciecia)
liczba_kredytow = tn + fn
print('Osiągnięta rentowność: {}%'.format(np.abs((tn * -15 + fn * 75)/liczba_kredytow).round(2)))
Mała uwaga: porównując osiągniętą rentowność z wynikiem uzyskanym przez regresję logistyczną, zauważymy spadek zwrotu z zainwestowanego kapitału (6.31% vs 5.28%). To dlatego, że przy niższych punktach odcięcia podejmujemy trafniejsze decyzje, ale akceptujemy mniej wniosków. Celem jest tu maksymalizacja zysku, dlatego wybrałem taki, a nie inny punkt odcięcia.
6. XGBoost vs regresja logistyczna - porównanie wyników.¶
Poniżej znajduje się tabela porównująca wyniki uzyskane przez oba modele:
- XGBoost - model numer 6 z tego projektu,
- regresja logistyczna - model 3.4. - najlepszy logit jaki do tej pory udało mi się zbudować dla tego zbioru.
pd.DataFrame({'Regresja logistyczna [model 3.4]' : [0.566, '268200 zł', '6.31%', '16920 zł'],
'XGBoost [model 6]' : [0.58, '330500 zł', '5.28%', '17445 zł']},
index = ['gini_te', 'zainwestowany_kapitał', 'stopa_zwrotu', 'zysk']).T
Podsumowanie¶
Mam nadzieję, że powyższy projekt przypadł Ci do gustu. Ja miałem sporo frajdy, walcząc o poprawę wyniku. 🙂 Proszę, daj znać, co o nim myślisz. Które z technik przypadku Ci do gustu, a które zamieniłbyś na coś innego? Napisz w komentarzu i podziel się swoim zdaniem ze mną i innymi czytelnikami. 🙂
Bardzo ciekawy artykuł. Fajnie poczytać coś o XGB.
Cieszę się, że Ci się podobał. Z czasem będzie nieco więcej wpisów o XGB i LightGBM. 🙂
Co to za biblioteka: creme_de_la_creme_of_ML?
Chciałbym dowiedzieć się w jaki sposób zmienne zostały zakodowane poprzez użycie MeanEncoding.
Czy możesz napisać o tym coś więcej ?
Sebastian, dziękuję za komentarz. 🙂
Zapomniałem o tym wspomnieć we wpisie. Jest to moja prywatna biblioteka, którą kiedyś żartobliwie w ten sposób nazwałem. Mam w niej przygotowane klasy i metody, których niegdyś nigdzie nie mogłem znaleźć w Python, a których to regularnie używam.
To temat na zupełnie odrębny wpis. Komentarz byłby zbyt obszerny. Poinformuję Cię, gdy opublikuję post w którym to wyjaśniam.