Przewidywanie defaultu wśród posiadaczy kart kredytowych

Ostatni projekt, który publikowałem na blogu, dotyczył klasyfikacji wniosków o wydanie karty kredytowej. Z pomocą uczenia maszynowego starałem się w nim odzwierciedlić decyzje, jakie podejmują eksperci pracujący w banku. Dziś idę o krok dalej i pokażę, jak można przewidzieć zdarzenie niewypłacalności posiadaczy kart kredytowych.

---

1. Wstęp
2. Cel projektu
3. Założenia dotyczące projektu
4. Opis zbioru danych
5. Wczytanie i przygotowanie danych
6. Statystyka opisowa
7. Analiza pojedynczych zmiennych
8. Analiza zależności pomiędzy zmiennymi
9. Przygotowanie danych do modelowania
10. Model 1
11. Model 2
12. Finalna weryfikacja jakości modelu
13. Podsumowanie

---

Wstęp

Projekt ten można nazwać Credit Scoringiem, choć ten przypadek nieco różni się od klasycznej definicji. Skupię się bowiem na osobach posiadających kartę kredytową. Postaram się przewidzieć ich zachowania i wyznaczyć prawdopodobieństwo niepłacenia zobowiązań wobec banku.

Ten wpis skupia się na wybranym przykładzie. Warto mieć jednak  na uwadze,  że tego typu rozwiązania mają wiele zastosowań. Systemy scoringowe służą głównie do oceny ryzyka i wglądu w wypłacalnośc klientów. Te z kolei przekładają się na dalsze korzyści, takie jak:

• automatyzacja podejmowanych decyzji,
• ograniczanie ryzyka,
• optymaizacja kosztów.

Oprócz powyższych przykładów dosyć popularnym tematem jest tzw. risk based pricing. W skrócie polega on na dynamicznym dostosowaniu oferty (np. kredytu hipotecznego) do ryzyka związanego z podjęciem współpracy z danym klientem. To, co łączy te wszystkie zastosowania, to ich główny cel: zwiększenia dochodowości danej firmy.

Cel projektu

Głównym celem projektu jest zbudowanie modelu predykcyjnego, który pozwoli wyznaczyć PD (ang. Probability of Default), a więc prawdopodobieństwo wystąpienia niepożądanego dla banku zdarzenia – niewypłacalności klienta.

Założenia dotyczące projektu

Scoring kredytowy jest bardzo szerokim zagadnieniem, dlatego celowo narzucę w tym miejscu pewne ograniczenia. Użycie wszystkich technik w jednym projekcie, który publikuję na blogu, byłoby niezwykle czasochłonne i mało przystępne dla Ciebie drogi czytelniku 🙂

Chcę uniknąć sytuacji, w której publikuję jeden niezwykle obszerny artykuł, który porusza wiele aspektów danego zagadnienia, ale jednocześnie jest na tyle długi, że mało kto dociera do jego końca. Stawiam zatem na systematyczne poszerzeanie tematyki i publikowanie relatywnie „prostych” i przystępnych wpisów. Wybacz zatem, jeśli nie ujrzysz w tym wpisie wzmianki o kilku technikach charakterystycznych dla credit scoringu 🙂

Pozostałe założenia:

• Scoring kredytowy, to zagadnienie używane w sektorze finansowym, a więc model, który zbuduję musi być w pełni interpretowalny.
• W ostatnim projekcie użyłem drzewa decyzyjnego, dlatego dziś wybieram regresję logistyczną.
• Unikam wszystkich zbędnych transformacji danych, które mogą negatywnie wpłynąć na interpretowalność modelu.
• Na interpretowalność modelu negatywnie wpływa duża liczba zmiennych – będę się starać maksymalnie zmniejszyć liczbę zmiennych, przy jednoczesnym jak najmniejszym uszczerbku na jakości.
• W przypadku scoringu prawdopodobieństwo wystąpienia zdarzenia jest ważniejsze niż przypisanie obserwacji do danej klasy. Jako główną miarę jakości modelu użyję miary Gini.

Opis zbioru danych

Główne informacje o zbiorze:

• Default of Credit Card Clients Data Set – UCI.
• Autor: I-Cheng Yeh.
• Dodano w 2016.01.26.
• 24 zmienne, 30 000 obserwacji.

Opis zmiennych zawartych w zbiorze:

• X1: Suma_kredytow / Limit_na_karcie (danej osoby i najbliższej rodziny).
• X2: Plec (1 = mężczyzna; 2 = kobieta).
• X3: Edukacja (1 = wyzsze_pelne, link_2; 2 = wyzsze; 3 = srednie; 4 = inne).
• X4: Stan_cywilny (1 = w_zwiazku; 2 = kawaler_panna; 3 = inny).
• X5: Wiek (podany w latach).
• X6 – X11: Historia poprzednich płatności (od kwietnia 2005 to września 2005).
  • Historia została podana od najnowszych płatności do najstarszych, a więc: X6 – wrzesień, X11 – kwiecień 2005.
  • Zmienne przyjmują wartości (-1, 9), gdzie:
    • -1 – płatność w terminie,
    • 1,2,3…9 – liczba miesięcy opóźnienia.
• X12-X17: Suma_kwot_na_wyciagu (ang. Amount of bill statement).
  • X12 = wyciąg za wrzesień 2005, X13 = wyciąg za sierpień 2015, itd.
• X18-X23: Suma_poprzednich_platnosci.
  • X18 = suma płatności we wrześniu 2005, X19 = suma płatności w sierpniu 2015, itd.

Dysponuję danymi do września 2015, dlatego zakładam, że punkt obserwacji (punkt w czasie, od którego przewidujemy zachowania klientów), to październik 2005.

Autor niestety nie podał kilku znaczących informacji dotyczących:

• Okresu obserwacji (ang. outcome period), czyli okresu w czasie, w którym badamy zajście zdarzenia. Okres ten trwa od punktu obserwacji do punktu granicznego.
• Aplikacyjnej definicji defaultu, a więc zdarzenia, do które powinno dojść w okresie obserwacji, by można było mówić o default-cie. Wiedza ta mogłaby nieco ułatwić modelowanie.

Kolejnym ograniczeniem jest brak podziału zbioru. Oczywiście, mogę to zrobić samemu, lecz będzie to nieco przekłamane. W idealnym świecie powinienem dysponować jeszcze co najmniej jedną próbką testową: out-of-time. Zawierać ona powinna obserwacje i opisujące je zmienne z innego okresu niż obserwacje widoczne w zbiorze uczącym i walidującym. Pozwoliłoby to na symulację produkcyjnego działania modelu, na „świeżych” danych. O to jednak trudno w darmowych, ogólnodostępnych zbiorach. Muszę się zatem ograniczyć do tego, co jest dostępne 🙂

Wczytanie i przygotowanie danych

Wczytuję wszystkie niezbędne biblioteki. W tym wpisie będzie mniej sklearna, a więcej statsmodels. Czemu? Powodów jest kilka:

• w moim odczuciu regresja logistyczna z biblioteki statsmodels jest bardziej przyjazna niż ta z sklearn,
• regresja logistyczna z biblioteki statsmodels daje bardziej rozbudowane podsumowanie modelu i przypomina standard znany z implementacji logrega z R,
• z mojego doświadczenia wynika, że w przypadku instytucji finansowych, to R i statsmodels są standardem.

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import scipy.stats
import statsmodels.formula.api as sm
from sklearn.model_selection import train_test_split
from sklearn.metrics import auc, roc_curve
import statsmodels.api as sm_api
from statsmodels.stats.outliers_influence import variance_inflation_factor as vif

Dane są zapisane w pliku Excel, dlatego w pierwszej kolejności sprawdzam dostępne w pliku zakładki, a dopiero później go wczytuję.

ef = pd.ExcelFile('data/credit_card_default.xls')
print(ef.sheet_names)

out: ['Data']

df = ef.parse('Data')
df.head()

	X1	X2	X3	X4	X5	X6	X7	X8	X9	X10	X11	X12	X13	X14	X15	X16	X17	X18	X19	X20	X21	X22	X23	Y
ID	LIMIT_BAL	SEX	EDUCATION	MARRIAGE	AGE	PAY_0	PAY_2	PAY_3	PAY_4	PAY_5	PAY_6	BILL_AMT1	BILL_AMT2	BILL_AMT3	BILL_AMT4	BILL_AMT5	BILL_AMT6	PAY_AMT1	PAY_AMT2	PAY_AMT3	PAY_AMT4	PAY_AMT5	PAY_AMT6	default payment next month
1	20000	2	2	1	24	2	2	-1	-1	-2	-2	3913	3102	689	0	0	0	0	689	0	0	0	0	1
2	120000	2	2	2	26	-1	2	0	0	0	2	2682	1725	2682	3272	3455	3261	0	1000	1000	1000	0	2000	1
3	90000	2	2	2	34	0	0	0	0	0	0	29239	14027	13559	14331	14948	15549	1518	1500	1000	1000	1000	5000	0
4	50000	2	2	1	37	0	0	0	0	0	0	46990	48233	49291	28314	28959	29547	2000	2019	1200	1100	1069	1000	0

Jak się okazuje, zbiór ma nietypową strukturę i zawiera nagłówek w dwóch pierwszych wierszach. Nanoszę zatem małą korektę, zmieniam nazwy zmienych na polskie i usuwam zmienną „id”.

df = ef.parse('Data', skiprows=1, names = ['id', 'lim_kredytu', 'plec', 'wyksztalcenie', 'stan_cywilny', 'wiek', 'opozn_plat_wrz', 'opozn_plat_sie', 'opozn_plat_lip', 'opozn_plat_cze', 'opozn_plat_maj', 'opozn_plat_kwi', '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', 'y'])
df.drop('id', axis = 1, inplace = True)
df.head()

	lim_kredytu	plec	wyksztalcenie	stan_cywilny	wiek	opozn_plat_wrz	opozn_plat_sie	opozn_plat_lip	opozn_plat_cze	opozn_plat_maj	opozn_plat_kwi	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	y
0	20000	2	2	1	24	2	2	-1	-1	-2	-2	3913	3102	689	0	0	0	0	689	0	0	0	0	1
1	120000	2	2	2	26	-1	2	0	0	0	2	2682	1725	2682	3272	3455	3261	0	1000	1000	1000	0	2000	1
2	90000	2	2	2	34	0	0	0	0	0	0	29239	14027	13559	14331	14948	15549	1518	1500	1000	1000	1000	5000	0
3	50000	2	2	1	37	0	0	0	0	0	0	46990	48233	49291	28314	28959	29547	2000	2019	1200	1100	1069	1000	0
4	50000	1	2	1	57	-1	0	-1	0	0	0	8617	5670	35835	20940	19146	19131	2000	36681	10000	9000	689	679	0

Ok. Teraz całość wygląda prawidłowo. Sprawdzam rozmiar zbioru.

print('Zbiór zawiera {} obserwacji, opisanych z pomocą {} zmiennych.'.format(df.shape[0], df.shape[1]))

out: Zbiór zawiera 30000 obserwacji, opisanych z pomocą 24 zmiennych.

Przed przejściem do modelowania należy również sprawdzić zbiór pod kątem ewentualnych braków danych.

df.isnull().any().any()

out: False

Brak braków 🙂

Ostatnim elementem, który powinienem sprawdzić, są typy poszczególnych zmiennych.

print('Typy poszczególnych zmiennych:')
df.dtypes

lim_kredytu	int64
plec	int64
wyksztalcenie	int64
stan_cywilny	int64
wiek	int64
opozn_plat_wrz	int64
opozn_plat_sie	int64
opozn_plat_lip	int64
opozn_plat_cze	int64
opozn_plat_maj	int64
opozn_plat_kwi	int64
kwota_wyciagu_wrz	int64
kwota_wyciagu_sie	int64
kwota_wyciagu_lip	int64
kwota_wyciagu_cze	int64
kwota_wyciagu_maj	int64
kwota_wyciagu_kwi	int64
platnosc_wrz	int64
platnosc_sie	int64
platnosc_lip	int64
platnosc_cze	int64
platnosc_maj	int64
platnosc_kwi	int64
y	int64

Teoretycznie mamy same liczby całkowite. W praktyce natomiast zmienne: „plec”, „wyksztalcenie” i „stan_cywilny”, to zmienne kategoryczne. Zmieniam zatem ich wartości.

df.plec.replace([1,2], [0, 1], inplace = True)
df.wyksztalcenie.replace([0, 1, 2, 3, 4, 5, 6], ['nieznane', 'wyzsze_pelne', 'wyzsze', 'srednie', 'inne', 'nieznane', 'nieznane'], inplace = True)
df.stan_cywilny.replace([0, 1, 2, 3], ['nieznany', 'w_zwiazku', 'kawaler_panna', 'inny'], inplace = True)

Statystyka opisowa

Zmienne numeryczne

df.describe()

	lim_kredytu	plec	wiek	opozn_plat_wrz	opozn_plat_sie	opozn_plat_lip	opozn_plat_cze	opozn_plat_maj	opozn_plat_kwi	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	y
count	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0	30000.0
mean	167484.32266666667	0.6037333333333333	35.4855	-0.0167	-0.13376666666666667	-0.1662	-0.22066666666666668	-0.2662	-0.2911	51223.3309	49179.07516666667	47013.1548	43262.94896666666	40311.40096666667	38871.7604	5663.5805	5921.1635	5225.6815	4826.076866666666	4799.387633333334	5215.502566666667	0.2212
std	129747.66156719506	0.4891291960904071	9.217904068090183	1.1238015279973212	1.19718597303439	1.1968675684467378	1.1691386224022984	1.1331874060026166	1.1499876256077741	73635.86057552874	71173.76878252918	69349.38742703729	64332.85613391704	60797.155770264195	59554.10753674454	16563.280354026534	23040.870402054872	17606.96146980426	15666.159744031342	15278.30567914539	17777.465775434066	0.4150618056909671
min	10000.0	0.0	21.0	-2.0	-2.0	-2.0	-2.0	-2.0	-2.0	-165580.0	-69777.0	-157264.0	-170000.0	-81334.0	-339603.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
25%	50000.0	0.0	28.0	-1.0	-1.0	-1.0	-1.0	-1.0	-1.0	3558.75	2984.75	2666.25	2326.75	1763.0	1256.0	1000.0	833.0	390.0	296.0	252.5	117.75	0.0
50%	140000.0	1.0	34.0	0.0	0.0	0.0	0.0	0.0	0.0	22381.5	21200.0	20088.5	19052.0	18104.5	17071.0	2100.0	2009.0	1800.0	1500.0	1500.0	1500.0	0.0
75%	240000.0	1.0	41.0	0.0	0.0	0.0	0.0	0.0	0.0	67091.0	64006.25	60164.75	54506.0	50190.5	49198.25	5006.0	5000.0	4505.0	4013.25	4031.5	4000.0	0.0
max	1000000.0	1.0	79.0	8.0	8.0	8.0	8.0	8.0	8.0	964511.0	983931.0	1664089.0	891586.0	927171.0	961664.0	873552.0	1684259.0	896040.0	621000.0	426529.0	528666.0	1.0

Pierwsze wnioski:

• W kilku przypadkach mamy znaczne różnice pomiędzy wartościami średnimi, a medianą – oznaka skośności rozkładu.
• W kolejnych krokach będę badać korelację. Korelacja Pearsona zakłada normalność rozkładów zmiennych, będzie konieczne zatem użycie korelacji rang Spearmana.

Zmienne kategoryczne

df.select_dtypes(include = ['object']).describe()

	wyksztalcenie	stan_cywilny
count	30000	30000
unique	5	4
top	wyzsze	kawaler_panna
freq	14030	15964

Skromniejsze podsumowanie niż w przypadku zmiennych numerycznych. Nie wiele z niego wynika, poza tym, że najwięcej jest osób o stanie cywilnym kawaler/panna. Większość posiadaczy kart kredytowych to osoby o wykształceniu wyższym.

Analiza pojedynczych zmiennych

Analiza zmiennej celu

df['y'].value_counts(normalize = True)

plt.figure(figsize=(10,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.countplot(df['y'], palette = ['#eb6c6a', '#f0918f']).set(title = 'Wykres częstości - "y"', xlabel = 'Default', ylabel = 'liczba obserwacji')
plt.show()

out: 0 0.779
     1 0.221

Defaulty stanowią ok. 22% wszystkich obserwacj. Nie jest źle, choć mogłoby być nieco lepiej.

Analiza zmiennych numerycznych

W tym punkcie zbadam zmienne numeryczne. Skupię się na analizie wykresów pod kątem ewentualnych anomalii i wykrycia wartości odstających. Dodatkowo zbadam rozkład zmiennych z użyciem biblioteki scipy.

Zmienna „lim_kredytu”

# Histogram
plt.figure(figsize=(13,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.distplot(df['lim_kredytu'], kde = False, bins = 30, color = '#eb6c6a').set(title = 'Histogram - "lim_lredytu"', xlabel = 'limit na karcie kredytowej', ylabel = 'liczba obserwacji')
plt.show()

# Wykres gęstości
plt.figure(figsize=(13,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.kdeplot(df['lim_kredytu'], shade = True, color = '#eb6c6a').set(title = 'Wykres gęstości - "lim_lredytu"', xlabel = 'limit na karcie kredytowej', ylabel = '')
plt.show()

# Wykres pudełkowy
plt.figure(figsize=(13,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.boxplot(df['lim_kredytu'], color = '#eb6c6a').set(title = 'Wykres pudełkowy - "lim_lredytu"', xlabel = 'limit na karcie kredytowej', ylabel = 'liczba obserwacji')
plt.show()

# Test na normalność rozkładu
# Zakładany poziom istotności alfa = 0.05.
if(scipy.stats.normaltest(df['lim_kredytu'])[1] < 0.05):
    print('Odrzucam hipotezę zerową i przyjmuję hipotezę alternatywną: zmienna nie pochodzi z rozkładu normalnego.')
else:
    print('Przyjmuję hipotezę zerową. Zmienna pochodzi z rozkładu normalnego.')

out: Odrzucam hipotezę zerową i przyjmuję hipotezę alternatywną: zmienna nie pochodzi z rozkładu normalnego.

Zmienna „wiek”

# Histogram
plt.figure(figsize=(13,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.distplot(df['wiek'], kde = False, bins = 30, color = '#eb6c6a').set(title = 'Histogram - "wiek"', xlabel = 'wiek', ylabel = 'liczba obserwacji')
plt.show()

# Wykres gęstości
plt.figure(figsize=(13,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.kdeplot(df['wiek'], shade = True, color = '#eb6c6a').set(title = 'Wykres gęstości - "wiek"', xlabel = 'wiek', ylabel = '')
plt.show()

# Wykres pudełkowy
plt.figure(figsize=(13,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.boxplot(df['wiek'], color = '#eb6c6a').set(title = 'Wykres pudełkowy - "wiek"', xlabel = 'wiek', ylabel = '')
plt.show()

# Test na normalnośc rozkładu
if(scipy.stats.normaltest(df['wiek'])[1] < 0.05):
    print('Odrzucam hipotezę zerową i przyjmuję hipotezę alternatywną: zmienna nie pochodzi z rozkładu normalnego.')
else:
    print('Przyjmuję hipotezę zerową. Zmienna pochodzi z rozkładu normalnego.')

out: Odrzucam hipotezę zerową i przyjmuję hipotezę alternatywną: zmienna nie pochodzi z rozkładu normalnego.

Opóźnienia płatności

Pomijam. Nie ma tu mowy o żadnych wartościach odstających. Przekładam tu racjonalne myślenie ponad statystykę 🙂

Kwota widoczna na wyciągu z rachunku

W tym przypadku pochylę się nad całą grupą zmiennych. By wykonać wizualizację danych na jednym wykresie, zmienię nieco ich strukturę.

mdf_kwota = pd.melt(df, value_vars=['kwota_wyciagu_kwi', 'kwota_wyciagu_maj', 'kwota_wyciagu_cze', 'kwota_wyciagu_lip', 'kwota_wyciagu_sie', 'kwota_wyciagu_wrz'], var_name=['miesiac'], value_name = 'kwota_wyciagu') 
mdf_kwota.head()

	miesiac	kwota_wyciagu
0	kwota_wyciagu_kwi	0
1	kwota_wyciagu_kwi	3261
2	kwota_wyciagu_kwi	15549
3	kwota_wyciagu_kwi	29547
4	kwota_wyciagu_kwi	19131

plt.figure(figsize=(13,13))
sns.set(font_scale=1.4, style="whitegrid")
sns.boxplot(data = mdf_kwota, y = 'miesiac', x = 'kwota_wyciagu', color = '#eb6c6a').set(title = 'Kwoty widoczna na wyciągu dla poszczególnych miesięcy')
plt.show()

Widać całą masę outlierów. Nie znaczy to oczywiście, że muszę je wszystkie usuwać. Kropki symbolizują wartości odstające, zostały wyznaczone metodą 1.5 odstępu międzykwartylowego. Raczej nie będę z niej korzystać. Powód jest prosty: oznaczałaby ona usunięcie znacznej części obserwacji ze zbioru, a tego bym wolał uniknąć.

Kwota płatności w danym miesiącu

Postępuję podobnie jak w przypadku kwoty widocznej na wyciągu.

mdf_platnosc = pd.melt(df, value_vars=['platnosc_kwi', 'platnosc_maj', 'platnosc_cze', 'platnosc_lip', 'platnosc_sie', 'platnosc_wrz'], var_name=['miesiac'], value_name = 'platnosc')
mdf_platnosc.head()

Jeszcze więcej outlierów niż w poprzednim przypadku. Ten wykres dobitnie ukazuje, dlaczego metoda 1.5 IQR nie powinna być stosowana w tym przypadku 🙂

Analiza zmiennych kategorycznych

Zmienna „plec”

print('Rozkład zmiennej "plec":')
print(df['plec'].value_counts(normalize = True))

plt.figure(figsize=(10,7))
sns.set(font_scale=1.4)
sns.countplot(df['plec'], palette = 'Blues_d', order = df['plec'].value_counts().index).set(title = 'Wykres częstości - "plec"', xlabel = 'plec', ylabel = 'liczba obserwacji')
plt.show()

out: 
Rozkład zmiennej "plec":
     1 0.604
     0 0.396

Mamy nieco więcej kobiet, choć różnica nie jest ogromna.

Zmienna „wyksztalcenie”

print('Rozkład zmiennej "wyksztalcenie"')
print(df['wyksztalcenie'].value_counts(normalize = True))

plt.figure(figsize=(10,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.countplot(df['wyksztalcenie'], palette = ['#eb6c6a', '#f0918f', '#f2a3a2', '#f5b5b4', '#f7c8c7']).set(title = 'Wykres częstości - "wyksztalcenie"', xlabel = 'wyksztalcenie', ylabel = 'liczba obserwacji')
plt.show()

out: Rozkład zmiennej "wyksztalcenie": wyzsze 0.468 wyzsze_pelne 0.353 srednie 0.164 nieznane 0.011 inne 0.004

Zmienna „stan_cywilny”

print('Rozkład zmiennej "stan_cywilny"')
print(df['stan_cywilny'].value_counts(normalize = True))

plt.figure(figsize=(10,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.countplot(df['stan_cywilny'], palette = ['#eb6c6a', '#f0918f', '#f2a3a2', '#f5b5b4', '#f7c8c7'], order = df['stan_cywilny'].value_counts().index).set(title = 'Wykres częstości - "stan_cywilny"', xlabel = 'stan_cywilny', ylabel = 'liczba obserwacji')
plt.show()

out: 
Rozkład zmiennej "stan_cywilny":
    kawaler_panna   0.532
    w_zwiazku       0.455
    inny            0.011
    nieznany        0.002

Analiza zależności pomiędzy zmiennymi

Analiza zależności zmiennych objaśniających ze zmienną celu

Zmienne numeryczne a zmienna celu

W przypadku zmiennych numerycznych użyję analizy korelacji Spearmana, ponieważ:

• nie ma ona założenia co do normalności rozkładu, a zgodnie z przeprowadzonymi testami odrzuciłem tę hipotezę w przypadku wszystkich badanych zmiennych,
• jako że korelacja Spearmana jest tak naprawdę korelacją rang, to mogę jej użyć do badania zależności pomiędzy zmiennymi numerycznymi i zmienną celu, która jest zmienną dychotomiczną.

corr_num_y = pd.DataFrame(scipy.stats.spearmanr(df.drop(['plec', 'wyksztalcenie', 'stan_cywilny'], axis = 1))[0],
columns = df.drop(['plec', 'wyksztalcenie', 'stan_cywilny'], axis = 1).columns,
index = df.drop(['plec', 'wyksztalcenie', 'stan_cywilny'], axis = 1).columns).iloc[-1,:-1]
corr_num_y = pd.DataFrame(corr_num_y)
corr_num_y.reset_index(inplace=True)

# wykres słupkowy współczynnika korelacji
plt.figure(figsize=(10,7))
sns.set(font_scale=1.4)
sns.barplot(data = corr_num_y.sort_values('y', ascending=False), x = 'y', y = 'index', palette = 'Blues_d').set(title = 'Wykres słupkowy - współczynnik korelacji \n zmiennych numerycznych ze zmienną celu', xlabel = 'współczynnik korelacji', ylabel = 'nazwa zmiennej')
plt.show()

By mieć lepszy obraz tego, jak duża jest zależność pomiędzy zmiennymi, wizualizuję wartości bezwzględne.

corr_num_y['abs']= corr_num_y['y'].abs()

plt.figure(figsize=(10,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.barplot(data = corr_num_y.sort_values('abs', ascending=False), x = 'abs', y = 'index', color = '#eb6c6a').set(title = 'Wykres słupkowy - współczynnik korelacji \n zmiennych numerycznych ze zmienną celu \n (wartości bezwzględne).', xlabel = 'współczynnik korelacji', ylabel = 'nazwa zmiennej')
plt.show()

Wnioski:

• im bliżej października, tym większa korelacja poszczególnym grup zmiennych ze zmienną celu:
• opóźnienia w płatnościach (korelacja pozytywna),
• kwoty na wyciągach z konta (korelacja negatywna),
• wysokość płatności (korelacja negatywna),
• lim_kredytu jest negatywnie skorelowany ze zmienną celu, a więc osoby o większym limicie rzadziej mają problemy ze spłatą zobowiązań w październiku,
• zmienna wiek jest nieznacznie skorelowana ze zmienną celu.

Przyjrzyjmy się nieco bliżej wspomnianej zmiennej „lim_kredytu”, by potwierdzić powyższe wyniki.

ct = pd.crosstab(df.lim_kredytu > df.lim_kredytu.mean(), df.y)
ct.index = ['limit_ponizej_sredniej', 'limit_powyzej_sredniej']
ct['% of default'] = round(ct[1]/(ct[0]+ct[1])*100,2)
print(ct)

	0	1	% of default
limit_ponizej_sredniej	12448	4646	27.18
limit_powyzej_sredniej	10916	1990	15.42

W przypadku kredytów o limicie powyżej średniej, odsetek defaultów jest niemal dwukrotnie mniejszy niż w przypadku kredytów o limicie poniżej średniej. Potwierdza to zatem negatywną korelację ze zmienną celu.

Zmienne kategoryczne a zmienna celu

Do badania zależności pomiędzy poszczególnymi zmiennymi kategorycznymi użyję współczynnika V Craméra. Jest on jedną z metod służących do badania zależności pomiędzy zmiennymi kategorycznymi.

Niestety nie udało mi się znaleźć żadnej biblioteki z gotową metodą wyznaczania współczynnika V Craméra, dlatego sam definiuję funkcje, które wyznaczą macierz zależności.

def CramersV(tab):
    a = scipy.stats.chi2_contingency(tab)[0]/sum(tab.sum())
    b = min(tab.shape[0]-1, tab.shape[1]-1,)
    return(np.sqrt(a/b))

def CalculateCrammersV(tab):
    ret = []
    for m in tab:
        row = []
        for n in tab:
            cross_tab = pd.crosstab(tab[m].values,tab[n].values)
            row.append(CramersV(cross_tab))
        ret.append(row)
    return pd.DataFrame(ret, columns=tab.columns, index=tab.columns)

y	1.000
wyksztalcenie	0.073
plec	0.040
stan_cywilny	0.034

cat_y = pd.DataFrame(CalculateCrammersV(df[['y','wyksztalcenie', 'plec', 'stan_cywilny']]).iloc[0,1:])
cat_y.reset_index(inplace = True)

plt.figure(figsize=(10,7))
sns.set(font_scale=1.4, style="whitegrid")
sns.barplot(data = cat_y, y = 'y', x = 'index', palette = ['#eb6c6a', '#f0918f', '#f2a3a2', '#f5b5b4', '#f7c8c7']).set(title = 'Wykres słupkowy - współczynnik V Crammera', xlabel = 'zmienne kategoryczne', ylabel = 'Crammer\'s V')
plt.show()

Wnioski:

• największa zależność występuje pomiędzy zmienną celu a zmienną „wyksztalcenie”,
• najmniejszy poziom zależności występuje pomiędzy zmienną celu a zmienną”stan_cywilny”,
• we wszystkich trzech przypadkach poziom zależności jest relatywnie niski.

Analiza korelacji pomiędzy zmiennymi numerycznymi

corr_num = pd.DataFrame(scipy.stats.spearmanr(df.drop(['plec', 'wyksztalcenie', 'stan_cywilny'], axis = 1))[0],
columns = df.drop(['plec', 'wyksztalcenie', 'stan_cywilny'], axis = 1).columns,
index = df.drop(['plec', 'wyksztalcenie', 'stan_cywilny'], axis = 1).columns)

plt.figure(figsize=(15,6))
sns.set(font_scale=1)
sns.heatmap(corr_num.abs(), cmap="Reds", linewidths=.5).set(title='Heatmap-a współczynnika korelacji rang Spearmana')
plt.show()

Nie sposób umieścić całej macierzy na jednej stronie, więc pozostanę przy heatmapie 🙂

Widać, że kwota wyciągu w danym miesiącu jest skorelowana z płatnością w miesiącu kolejnym (np. kwota wyciągu za sierpień jest skorelowana z płatnością za wrzesień), co oczywiście ma sens 🙂 Im więcej się zadłużamy, tym więcej spłacamy. Bez tego na kredytach zarabialiby tylko specjaliści od windykacji 😉

Analiza zależności pomiędzy zmiennymi kategorycznymi

crammer = CalculateCrammersV(df[['y','wyksztalcenie', 'plec', 'stan_cywilny']])

plt.figure(figsize=(15,6))
sns.set(font_scale=1.4)
sns.heatmap(crammer, cmap="Reds", linewidths=.5).set(title='Heatmap-a współczynnika zależności Crammera')
plt.show()

Analiza zależności pomiędzy zmiennymi numerycznymi a kategorycznymi

W przypadku zmiennych numerycznych i kategorycznych nie mówimy już o korelacji, lecz o poziomie zależności. Miarą, która jest najbliższa współczynnikowi korelacji, w tym przypadku jest współczynnik korelacji wielorakiej (współczynnik R w modelach regresji) – użyję jej do badania zależności pomiędzy zmiennymi kategorycznymi (nominalnymi i porządkowymi) i numerycznymi (liczby całkowite i zmiennoprzecinkowe). Analogiczną miarą może być Eta-squared.

# dzielę zmienne na dwie grupy
cat_cols = ['plec', 'wyksztalcenie', 'stan_cywilny', 'y']
num_cols = df.drop(cat_cols, axis = 1).columns

# w pętli buduję kolejne modele regresji liniowej
cols = []
for cat in cat_cols:
    rows = []
    for num in num_cols:
        formula = num + '~' +cat
        model = sm.ols(formula=formula,data=df)
        rows.append(np.sqrt(model.fit().rsquared))
    cols.append(rows)
corr_num_cat = pd.DataFrame(cols, index = cat_cols, columns = num_cols)

# wykres zależności
plt.figure(figsize=(15,6))
sns.set(font_scale=1.4)
sns.heatmap(corr_num_cat, cmap="Reds", linewidths=.5).set(title='Heatmap-a współczynnika korelacji wielorakiej')
plt.show()

Przygotowanie danych do modelowania

Usunięcie wartości odstających

Zanim przystąpię do modelowania muszę usunąć wartości odstające. Reguła 1.5 IQR nie wchodzi w grę – sprawdziłem to i usuwała ona ok. 35% obserwacji. Nie mogę również zastosować reguły 3sigma, gdyż zakłada ona normalny rozkład zmiennych. Usuwam zatem 5% skrajnych wartości.

lower_bound = df.quantile(0.025)
upper_bound = df.quantile(0.975)
num_of_outl = (df[upper_bound.index] > upper_bound).sum()
outliers = pd.DataFrame({'upper_bound':upper_bound, 'num_of_outliers_U':num_of_outl})
outliers.drop(['y', 'wiek', 'plec', 'opozn_plat_wrz', 'opozn_plat_sie', 'opozn_plat_lip', 'opozn_plat_cze', 'opozn_plat_maj', 'opozn_plat_kwi'], inplace = True, axis = 0)

Z grupy zmiennych analizowanych pod kątem outlierów usunąłem m.in.”y” i „wiek”. Ciężko w przypadku nich mówić o odstających wartościach 🙂

for row in outliers.iterrows():
    df = df[df[row[0]] <= row[1]['upper_bound']]
    print(df.shape)

out: (25737, 24)

Ostatecznie w zbiorze pozostało 25737 obserwacij.

Zmiana kodowania zmiennych

Zmienne „y” i „plec” już przyjmują wartości binarne. Mamy zatem dwie zmienne, których kodowanie trzeba zmienić: „wyksztalcenie” i „stan_cywilny”.

df = pd.get_dummies(df)

„Pułapka dummy coding” – w modelach liniowych trzeba uważać na problem współliniowości. Umieszczając wszystkie kategorie danej zmiennej w modelu mamy idealną współliniowość. Muszę zatem usunąć po jednej kategorii z tych zmiennych, które zostały zapisane w więcej niż jednej kolumnie.

df.drop(axis = 1, labels = ['stan_cywilny_inny', 'wyksztalcenie_inne'], inplace = True)

Podział zbioru

Dzielę zbiór na 2 części:

• Zbiór treningowy – na nim wykonam walidację krzyżową. Celem sprawdzianu krzyżowego będzie dobór zmiennych użytych w modelu.
• Zbiór testowy – finalny sprawdzian dla zbudowanego modelu. Sprawdzę jak zbudowany algorytm performuje na danych, których nigdy wcześniej nie widział.

Podział: 80% – zbiór treningowy i walidacyjny, 20% – zbiór testowy, z uwzględnieniem stratyfikacji.

y = df.y.copy()
df.drop('y', axis = 1, inplace = True)
x_tr, x_te, y_tr, y_te = train_test_split(sm_api.add_constant(df).rename(columns = {'const':'intercept'}), y, test_size = 0.2, random_state = 2062018, stratify = y)
x_tr, x_va, y_tr, y_va = train_test_split(x_tr, y_tr, test_size = 0.25, random_state = 2062018, stratify = y_tr)

Model 1

Opis modelu:

• dobór zmiennych: brak.

Będzie to mój punk odniesienia. Do wyniku uzyskanego przez ten model będę nawiązywać w kolejnych punktach.

Buduję funkcję, która pozwoli mi w prosty sposób ocenić model. Podczas budowy kolejnych modeli będę jedynie wybierać zmienne do modelowania zbioru.

def build_model(selected_features, plots):
    """Logistic regression model evaluation.

    Parameters:
    -----------
    selected_features : list, list of selected featured
    
    plots : bool, decission whether to print plots
   
    Returns:
    --------
    gini_score: float, gini score
    """
    model = sm.Logit(y_tr, x_tr[selected_features]).fit(disp = 0)
    pre = model.predict(x_va[selected_features])
    fpr, tpr, thresholds = roc_curve(y_va, pre)
    gini_score = 2 * auc(fpr, tpr) - 1
    print('Gini score: {}.'.format(round(gini_score,3)))
    
    if plots == True:
        # ROC CURVE
        plt.figure()
        lw = 2
        plt.plot(fpr, tpr, color='blue',
                 lw=lw, label='ROC curve (area under curve = %0.3f)' % auc(fpr, tpr))
        plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
        plt.xlim([0.0, 1.0])
        plt.ylim([0.0, 1.05])
        plt.xlabel('False Positive Rate')
        plt.ylabel('True Positive Rate')
        plt.title('ROC curve')
        plt.legend(loc="lower right")
        plt.show()

        # DENSITY PLOT
        pre = pd.DataFrame(pre, columns = ['prob'])
        pre['pred_class'] = 1
        pre.loc[pre.prob < 0.5,'pred_class'] = 0
        pre['real_class'] = y_va

        sns.distplot(pre[pre.real_class == 0]['prob'], label = 'negatives').set(title = 'Density plot', xlabel = 'probability')
        sns.distplot(pre[pre.real_class == 1]['prob'], label = 'positives').set(title = 'Density plot', xlabel = 'probability')
        plt.legend()
        plt.show()
    
    return gini_score, model

out: Gini score: 0.496.

Wynik w okolicy 0.496 można uznać za zadowalający 🙂 To, co może martwić, to liczba obecnych zmiennych. Postaram się ją nieco zmniejszyć.

Model 2

Opis modelu:

• dobór zmiennych: forward selection.

W modelu numer 2 decyduję się na dobór zmiennych oparty o metodę krokową: forward selection. Nie ma jej dostępnej w żadnej znanej mi bibliotece Python, dlatego użyję funkcji, którą kiedyś znalazłem na jednym z blogów.

def forward_selection(data, response):
    """Linear model designed by forward selection.

    Parameters:
    -----------
    data : pandas DataFrame with all possible predictors and response

    response: string, name of response column in data

    feature_limit: int, stop criterion, maximum number of features

    Returns:
    --------
    best_features: an "optimal" fitted statsmodels linear model
    with an intercept
    selected by forward selection
    evaluated by Gini score
    """
    if 'intercept' in data.columns:
        data.drop('intercept', axis = 1, inplace = True)
    remaining = set(data.columns)
    remaining.remove(response)
    best_features = []
    current_score, best_new_score = 0.0, 0.0
    while remaining and current_score == best_new_score:
        scores_with_candidates = []
        for candidate in remaining:
            formula = "{} ~ {} + 1".format(response, ' + '.join(best_features + [candidate]))
            _pre = sm.logit(formula, data).fit(disp = 0).predict(x_va)
            fpr, tpr, thresholds = roc_curve(y_va, _pre)
            score = 2 * auc(fpr, tpr) - 1
            scores_with_candidates.append((score, candidate))
        scores_with_candidates.sort()
        best_new_score, best_candidate = scores_with_candidates.pop()
        if current_score < best_new_score:
            remaining.remove(best_candidate)
            best_features.append(best_candidate)
            current_score = best_new_score
    formula = "{} ~ {} + 1".format(response, ' + '.join(best_features))
    best_features.append('intercept')
    return best_features

# forward selection
selected_features = forward_selection(pd.concat([x_tr, y_tr], axis = 1), 'y')

# Liczba wybranych zmiennych
print('Liczba zmiennych: {}.'.format(len(selected_features)))

# weryfikacja jakości modelu
score, model = build_model(selected_features, True)
model.summary()

out:
Liczba zmiennych: 16.
Gini score: 0.496.`

Udało się utrzymać wynik przy znacznie mniejszej liczbie zmiennych 🙂 Rzućmy jednak jeszcze okiem na współczynniki modelu, które są widoczne w kolumnie „coef”. Jak je interpretować?

Otóż są to logarytmy szans. W prosty sposób możemy je przekształcić na szanse. Jeśli dodatkowo odejmiemy od nich wartość 1, to uzyskamy dane proste do interpretacji. Zmienne o współczynnikach ujemnych zmniejszają ryzyko wystąpienia defaultu. Zmienne o współczynnikach dodatnich zwiększają ryzyko.

(np.exp(model.params).sort_values(ascending = False)-1).round(5)

Nazwa zmiennej	Szanse
opozn_plat_wrz	0.68347
stan_cywilny_w_zwiazku	0.17987
opozn_plat_maj	0.06940
opozn_plat_sie	0.06235
opozn_plat_cze	0.05945
opozn_plat_lip	0.05824
wiek	0.00173
lim_kredytu	-0.00000
platnosc_maj	-0.00003
platnosc_lip	-0.00003
platnosc_wrz	-0.00003
platnosc_sie	-0.00005
plec	-0.07712
intercept	-0.61647
wyksztalcenie_nieznane	-0.67651
stan_cywilny_nieznany	-0.83840

Dla przykładu każdy miesiąc opóźnienia płatności (stan na wrzesień) zwiększa szanse na wystąpienie problemów z płatnością zobowiązania o ponad 68%. Fakt, że osoba jest w związku, zwiększa szanse defaultu o blisko 18%.

Często spotykam się z opinią, że współczynniki modelu świadczą o ważności zmiennej. Nic bardziej mylnego 🙂 Dopiero współczynniki wymnożone przez wartości poszczególnych zmiennych mówią nam o tym, jak ważna jest dana zmienna. Zmienne binarne o dużych wartościach szans, mogą mieć niekiedy mniejszy wpływ niż zmienne ciągłe (które mogą przyjmować duże wartości, np. saldo konta) o niskich wartościach szans. Liczy się więc iloczyn obu wartości.

Finalna weryfikacja jakości modelu

Używam uzyskanych w poprzednim kroku zmiennych do zbudowania finalnego modelu.

model = sm.Logit(y_tr, x_tr[selected_features]).fit(disp = 0)
pred = model.predict(x_te[selected_features])
fpr, tpr, thresholds = roc_curve(y_te, pred)
gini_score = 2 * auc(fpr, tpr) - 1
print('Gini score: {}.'.format(round(gini_score,3)))

out: Gini score: 0.467.

Przy zaledwie 16 zmiennych i braku znacznych przekształceń zbioru uznaje ten wynik za zadowalający. W kolejnych wpisach będę nawiązywać do tego projektu i postaram się przedstawić metody, które pomagają w doborze zmiennych i są powszechnie wykorzystywane przy problemach scoringu, np. analiza współczynników (VIF, IV, WoE) i analiza wykresów separacji gęstości prawdopodobieństw obu klas.

Podsumowanie

Dziękuję Ci za dobrnięcie do końca 🙂 Jeśli masz jakiekolwiek uwagi, pytania, lub spostrzeżenia, to proszę, podziel się nimi w komentarzu na social media, lub skontaktuj się ze mną.

Poniżej znajduje się lista plików użytych, wspomnianych, bądź też utworzonych w projekcie:

• Kod do pobrania z GitHub.
• Zbiór danych z którego korzystałem.
• Moje pozostałe projekty.

photo: Pixaby.com (StockSnap)