В проекте используются данные по кредитным картам клиентов. Данные можно скачать отсюда: default of credit card clients Data Set
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
%matplotlib inline
%pylab inline
pylab.rcParams['figure.figsize'] = (16, 7)
mpl.rcParams['figure.figsize'] = (16,7)
Представлены данные за полгода (апрель - сентябрь 2005) по клиентам банка в Тайване.
Список переменных:
LIMIT_BAL - Величина выдаваемого кредита: включает как индивидуальный потребительский кредит, так и его / ее семейный (дополнительный) кредит.
SEX - Пол (1 = мужчина; 2 = женщина)
EDUCATION - Образование (1 = аспиратура; 2 = высшее образование; 3 = среднее образование; 4 = другое).
MARRIAGE - Семейное положение (1 = женат/замужем; 2 = одинокий; 3 = другое).
AGE - Возраст (в годах).
PAY_0 ... PAY_6 - История платежей. Отслеживались прошлые ежемесячные платежные записи (с Апреля по Сентябрь, 2005) следующим образом: PAY_0 = статус погашения в Сентябре, 2005, ..., PAY_6 = статус погашения в Апреле, 2005. Шкала измерения для статуса погашения: -1 = оплачен вовремя; 1 = оплата задержана на месяц; ... 8 = оплата задержана на 8 месяцев; 9 = оплата задержана на 9 месяцев и больше.
BILL_AMT1 ... BILL_AMT6 - Сумма выписки по счету в каждый из месяцев. BILL_AMT1 = Сентябрь, 2005; ... BILL_AMT6 = Апрель, 2005.
PAY_AMT1 ... PAY_AMT6 - Сумма предыдущего платежа. PAY_AMT1 = выплаченная сумма в Сентябре, 2005; ...
PAY_AMT6 = выплаченная сумма в Апреле, 2005.
Целевая переменная: default payment next month - неуплата в следующем месяце, бинарная (1 = да, 0 = нет).
data = pd.read_excel('default of credit card clients.xls', header=1)
data.head().T
data.info()
Как видно, в данных нет пропущенных значений, тип всех переменных int64.
Выведем основные харакетристики признаков.
data.describe(include='all').T
# для удобства переименуем целевую переменную
data = data.rename(columns = {'default payment next month':'target'})
ax = data['target'].hist(orientation='horizontal', figsize=(10,5))
ax.set_title("default payment next month distribution")
print('Distribution of target:')
data['target'].value_counts() / data.shape[0]
Выборка не сбалансирована, одного класса больше чем другого.
Исследуем категориальные признаки $SEX$ (пол), $EDUCATION$ (образование) и $MARRIAGE$ (семейное положение).
print(data['SEX'].value_counts())
print('***')
print(data['EDUCATION'].value_counts())
print('***')
print(data['MARRIAGE'].value_counts())
В переменной $EDUCATION$ кроме заявленных категорий (1 = аспиратура; 2 = высшее образование; 3 = среднее образование; 4 = другое) появились неизвестные категории 0, 5 и 6. Проведем дополнительное исследование этих параметров.
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
%%time
data2 = data[(data['EDUCATION'] == 0) | (data['EDUCATION'] == 5)| (data['EDUCATION'] == 6) |
(data['EDUCATION'] == 4)]
scaler = StandardScaler()
scaler.fit(data2)
print('Размер нового датасета: {}'.format(data2.shape))
tsne_model = TSNE(2, random_state=17)
data2_2d = tsne_model.fit_transform(scaler.transform(data2))
plt.scatter(data2_2d[:, 0], data2_2d[:, 1], color = ['red' if i == 0 else (
'blue' if i==5 else ('black' if i == 4 else 'green')) for i in data2['EDUCATION']])
Категория 0, выделенная красным, четко отделена от категорий 4,5 и 6. Значит можем объединить их в одну категорию 4.
data['EDUCATION'].replace([5, 6], 4, inplace=True)
Аналогично в признаке $MARRIAGE$ появилась неизвестная категория 0. Исследуем этот случай.
%%time
data2 = data[(data['MARRIAGE'] == 0) | (data['MARRIAGE'] == 3)]
scaler = StandardScaler()
scaler.fit(data2)
print('Размер нового датасета: {}'.format(data2.shape))
tsne_model = TSNE(2, random_state=17)
data2_2d = tsne_model.fit_transform(scaler.transform(data2))
plt.scatter(data2_2d[:, 0], data2_2d[:, 1], color = ['red' if i == 0 else 'blue' for i in data2['MARRIAGE']])
В целом, категория 0 выделена в отдельную группу, кроме нескольких точек.
Исследуем зависимость категориальных признаков $SEX$ (пол), $EDUCATION$ (образование), $PAY\_i, i = 1,...,6$ (статус платежа в $i$-том месяце) и $MARRIAGE$ (семейное положение) от целевой переменной $target$.
sns.barplot(x='SEX', y='target', data=data)
Доля неуплативших долг у мужчин ($SEX$=1), чуть больше, чем у женщин ($SEX$=2), но в целом существенной разницы нет.
sns.barplot(x='MARRIAGE', y='target', data=data)
Доля неуплативших долг в незаявленной категории ($MARRIAGE$ = 0) меньше, чем в других.
sns.barplot(x='EDUCATION', y='target', data=data)
Доля неуплативших долг людей со средним образованием ($EDUCATION$ = 3) больше, чем в остальных. Причем в незаявленной категории
($EDUCATION$ = 0) неуплативших нет.
data = data.rename(columns = {'PAY_0':'PAY_1'})
pays = ['PAY_%s' % i for i in range(1,7)]
for i, pay in enumerate(pays):
plt.subplot(3, 2, i + 1)
sns.barplot(x=pay,y='target', data=data)
Можно заметить, что чем больше задержка оплаты долга, тем вероятней его неуплата.
corr_matrix = data.drop('ID', axis=1).corr()
sns.heatmap(corr_matrix, cmap="Blues")
Видно, что признаки $BILL\_AMTi, i = 1,...6$ (сумма выписки по счету в $i$-том месяце) коррелируют между собой, как и признаки $PAY\_i, i = 1,...,6$ (статус платежа в $i$-том месяце)
bills = ['BILL_AMT%s' % i for i in range(1,7)]
plots = data[bills].hist(figsize=(12,10))
pay_amts = ['PAY_AMT%s' % i for i in range(1,7)]
plots = data[pay_amts].hist(figsize=(12,10), bins = 5)
Так как в выборке есть дисбаланс данных, то в качестве целевой метрики будем использовать roc_auc_score, так же, для контроля, будем использовать f1_score и accuracy.
from sklearn.metrics import accuracy_score, confusion_matrix, f1_score, roc_auc_score
Для построения модели будем использовать три разных классификатора: логистическую регрессию, случайный лес и градиентный бустинг.
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
def learning_model(model, X_train, y_train, X_test, y_test):
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
y_roc = model.predict_proba(X_test)[:,1]
score_roc = roc_auc_score(y_test, y_roc)
score_f1 = f1_score(y_test, y_pred, average='weighted')
score_acc = accuracy_score(y_test, y_pred)
print('roc_auc = {}'.format(score_roc))
print('f1_score = {}'.format(score_f1))
print('acuracy_score={}'.format(score_acc))
return y_roc, score_roc, score_f1, score_acc
Так как выборка не сбалансирована, то будем использовать дополнительную настройку для случайного леса и логистической регрессии - class_weight='balanced'
forest = RandomForestClassifier(n_estimators=50, random_state=17,
class_weight='balanced')
logit = LogisticRegression(random_state=17, class_weight= 'balanced')
boosting = GradientBoostingClassifier(random_state=17)
Удалим столбец $ID$, для дальнейшего исседования он не нужен. Признак $SEX$ перекодируем в значения 0 ($SEX$=2) и 1 ($SEX$=1) и переименуем этот признак в $male$ (0 – женщина, 1 – мужчина).
data = data.drop('ID', axis=1)
data['SEX'] = data['SEX'].apply(lambda x: 0 if x==2 else 1 )
data = data.rename(columns = {'SEX':'male'})
Можно заметить, что переменные $BILL\_AMTi, i = 1,...6$ сильно коррелируют между собой, создадим переменную $balance$, в которую запишем разницу между $BILL\_AMTi$ сумма выписки по счету в $i$-том месяце) и $PAY\_AMT\_i$ (выплаченная сумма в $i$-том месяце).
balances = ['balance_%s' % i for i in range(1,7)]
for i in range(6):
data[balances[i]] = data[bills[i]] - data[pay_amts[i]]
Скорее всего хорошие клиенты никогда (или редко) задерживают платеж, тогда у клиентов с низким средним значением $PAY\_i, i = 1,...,6$ вероятнее не будет неуплаты платежа в следующем месяце. Создадим бинарный признак $good\_client$, где 1 - среднее значение < 0, 0 - среднее значение > 0.
data['pay_mean'] = (data['PAY_1'] + data['PAY_2'] + data['PAY_3'] + data['PAY_4'] + data['PAY_5']\
+ data['PAY_6'] ) / 6
data['good_client'] = data['pay_mean'].apply(lambda x: 1 if x <= 0 else 0).astype('int')
data = data.drop('pay_mean', axis=1)
sns.barplot(x='good_client', y='target', data=data)
Судя по графику разница между "хорошими" клиентами и "плохими" есть, ее влияние на предсказание оценим ниже.
Из категориальных переменных $EDUCATION$, $MARRIAGE$ и $PAY\_i$ (статус погашения в $i$-том месяце) сделаем дамми-переменные
col_for_dum = ['EDUCATION', 'MARRIAGE'] + pays
data = pd.get_dummies(data, columns=col_for_dum )
data.head()
Разобьем выборку на тестовую и обучающую.
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(data.drop('target', axis=1), data['target'],
test_size=0.3,
random_state=17,
stratify = data['target'])
from sklearn.preprocessing import StandardScaler
Отдельно маштабируем признаки $PAY\_AMT\_i$ (выплаченная сумма в $i$-том месяце), $AGE$, $LIMIT\_BAL$ и новые колонки $balance\_i$ в тестовой и отложенной выборке.
scaler = StandardScaler()
X_train_scaled = pd.DataFrame(scaler.fit_transform(X_train), columns=X_train.columns)
X_test_scaled = pd.DataFrame(scaler.transform(X_test), columns=X_train.columns)
Построим предсказания на заданных классификаторах без настройки гиперпараметров.
print('RandomForest:')
learning_model(forest, X_train_scaled, y_train, X_test_scaled, y_test);
print('LogisticRegression:')
learning_model(logit, X_train_scaled, y_train, X_test_scaled, y_test);
print('GradientBoostingClassifier:')
learning_model(boosting, X_train_scaled, y_train, X_test_scaled, y_test);
Оценим важность признаков с помощью дерева решений.
from sklearn.tree import DecisionTreeClassifier
tree_sel = DecisionTreeClassifier(random_state=17)
tree_sel.fit(X_train_scaled, y_train)
features = pd.DataFrame(tree_sel.feature_importances_,
index=X_train.columns,
columns=['Importance'])
features.sort_values(by = 'Importance', ascending=False).head(20)
Наибольшую важность имеют признаки $good\_client$, $AGE$, $LIMIT\_BAL$, созданные признаки $balance\_i$ так же входят в 20 самых важных признаков.
f_sort = features.sort_values(by = 'Importance', ascending=False)
plt.plot(range(len(f_sort.Importance.tolist())),
f_sort.Importance.tolist())
Отберем только 45 самых важных признаков.
col = f_sort.iloc[:45].T.columns
X_train_new = X_train_scaled[col]
X_test_new = X_test_scaled[col]
from sklearn.model_selection import StratifiedShuffleSplit, GridSearchCV
sss = StratifiedShuffleSplit(n_splits=5, test_size=0.3, random_state=17)
logit_parameters = {'C': np.logspace(0, 1, 10)}
boosting_parameters = {'n_estimators': [100, 150],
'learning_rate': [0.01, 0.05, 0.5],
'max_depth': [5, 6, 7],
'max_features': [0.3, 0.5, 0.75, 1.0],
'min_samples_leaf': [4, 5, 6]}
forest_parameters = {'n_estimators': [50, 75, 100],
'max_depth': [3, 6, 9, 12],
'max_features': [0.5, 0.75, 1.0]}
grid_logit = GridSearchCV(logit, logit_parameters, n_jobs=-1, scoring ='roc_auc', cv=sss)
grid_boosting = GridSearchCV(boosting, boosting_parameters, n_jobs=-1, scoring ='roc_auc', cv=sss)
grid_forest = GridSearchCV(forest, forest_parameters, n_jobs=-1, scoring ='roc_auc', cv=sss)
y_score = []
roc_scores, f1_scores, acc_scores = [], [], []
for grid in (grid_logit, grid_boosting, grid_forest):
y, roc, f1, acc = learning_model(grid, X_train_new, y_train,
X_test_new , y_test);
y_score.append(y)
roc_scores.append(roc)
f1_scores.append(f1)
acc_scores.append(acc)
Лучший результат показал градиентный бустинг с $roc\_auc = 0.786$
max_feat = grid_boosting.best_params_['max_features']
min_leaf = grid_boosting.best_params_['min_samples_leaf']
n_est = grid_boosting.best_params_['n_estimators']
from sklearn.model_selection import learning_curve, validation_curve
def plot_learning_curve(estimator, title, X, y,label_y, ylim=None, cv=None, scoring='roc_auc',
n_jobs=1, train_sizes=np.linspace(.1, 1.0, 10)):
plt.figure()
plt.title(title)
if ylim is not None:
plt.ylim(*ylim)
plt.xlabel("Training examples")
plt.ylabel(label_y)
train_sizes, train_scores, test_scores = learning_curve(
estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes, scoring=scoring)
train_scores_mean = np.mean(train_scores, axis=1)
train_scores_std = np.std(train_scores, axis=1)
test_scores_mean = np.mean(test_scores, axis=1)
test_scores_std = np.std(test_scores, axis=1)
plt.grid()
plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
train_scores_mean + train_scores_std, alpha=0.1,
color="r")
plt.fill_between(train_sizes, test_scores_mean - test_scores_std,
test_scores_mean + test_scores_std, alpha=0.1, color="g")
plt.plot(train_sizes, train_scores_mean, 'o-', color="r",
label="Training score")
plt.plot(train_sizes, test_scores_mean, 'o-', color="g",
label="Cross-validation score")
plt.legend(loc="best")
return plt
plot_learning_curve(GradientBoostingClassifier(random_state=17, learning_rate=0.05, max_depth=6,
n_estimators=n_est,
min_samples_leaf=min_leaf, max_features=max_feat),
'Gradient Boosting',
X_train_new, y_train, scoring='roc_auc', cv=sss, n_jobs=-1, label_y='roc_auc')
plot_learning_curve(GradientBoostingClassifier(random_state=17, learning_rate=0.05, max_depth=6,
n_estimators=n_est,
min_samples_leaf=min_leaf, max_features=max_feat),
'Gradient Boosting',
X_train_new, y_train, scoring='f1', cv=sss, n_jobs=-1, label_y='f1')
plot_learning_curve(GradientBoostingClassifier(random_state=17, learning_rate=0.05, max_depth=6,
n_estimators=n_est,
min_samples_leaf=min_leaf, max_features=max_feat),
'Gradient Boosting',
X_train_new, y_train, scoring='accuracy', cv=sss, n_jobs=-1, label_y='accuracy')
def plot_validation_curve(estimator, X, y, cv_param_name, cv_param_values,ylim=None, cv=None,
scoring='roc_auc',
n_jobs=-1, train_sizes=np.linspace(.1, 1.0, 10)):
plt.figure()
if ylim is not None:
plt.ylim(*ylim)
plt.xlabel(cv_param_name)
plt.ylabel("Score")
train_sizes, train_scores, test_scores = learning_curve(
estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes, scoring=scoring)
val_train, val_test = validation_curve(estimator, X, y, cv_param_name,
cv_param_values, cv=cv,
scoring=scoring)
val_train_mean = np.mean(val_train, axis=1)
val_train_std = np.std(val_train, axis=1)
val_test_mean = np.mean(val_test, axis=1)
val_test_std = np.std(val_test, axis=1)
plt.grid()
plt.fill_between(cv_param_values, val_train_mean - val_train_std,
val_train_mean + val_train_std, alpha=0.1,
color="r")
plt.fill_between(cv_param_values, val_test_mean - val_test_std,
val_test_mean + val_test_std, alpha=0.1, color="g")
plt.plot(cv_param_values, val_train_mean, 'o-', color="r",
label="Training score")
plt.plot(cv_param_values, val_test_mean, 'o-', color="g",
label="Cross-validation score")
plt.legend(loc="best")
return plt
estimators = np.arange(25, 425, 25)
plot_validation_curve(GradientBoostingClassifier(random_state=17, learning_rate=0.05, max_depth=6,
min_samples_leaf =min_leaf, max_features=max_feat),
X_train_new, y_train,
cv_param_name='n_estimators',
cv_param_values= estimators,
scoring='roc_auc')
depth = np.arange(3, 25)
plot_validation_curve(GradientBoostingClassifier(random_state=17, learning_rate=0.05,
n_estimators=n_est,
min_samples_leaf =min_leaf, max_features=max_feat),
X_train_new, y_train,
cv_param_name='max_depth',
cv_param_values= depth,
scoring='roc_auc')
learn_rate = np.linspace(0.01, 1.0, 10)
plot_validation_curve(GradientBoostingClassifier(random_state=17, max_depth=6,
n_estimators=n_est,
min_samples_leaf =min_leaf, max_features=max_feat),
X_train_new, y_train,
cv_param_name='learning_rate',
cv_param_values= learn_rate,
scoring='roc_auc')
min_l = np.arange(3, 10)
plot_validation_curve(GradientBoostingClassifier(random_state=17, max_depth=6,
n_estimators=n_est,max_features=max_feat),
X_train_new, y_train,
cv_param_name='min_samples_leaf',
cv_param_values= min_l,
scoring='roc_auc')
labels = ['LogisticRegression', 'GradientBoosting', 'RandomForest']
from sklearn.metrics import roc_curve, auc
for i in range(3):
# Compute ROC curve and ROC area
fpr, tpr, _ = roc_curve(y_test, y_score[i])
roc_auc = auc(fpr, tpr)
# Plot of a ROC curve for a specific class
plt.figure(1)
plt.plot(fpr, tpr, label=labels[i]+', ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], 'k--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.legend(loc="lower right")
import itertools
def plot_confusion_matrix(cm, classes,
normalize=False,
title='Confusion matrix',
cmap=plt.cm.Blues):
"""
This function prints and plots the confusion matrix.
Normalization can be applied by setting `normalize=True`.
"""
plt.imshow(cm, interpolation='nearest', cmap=cmap)
plt.title(title)
plt.colorbar()
tick_marks = np.arange(len(classes))
plt.xticks(tick_marks, classes, rotation=45)
plt.yticks(tick_marks, classes)
if normalize:
cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
print("Normalized confusion matrix")
else:
print('Confusion matrix, without normalization')
#print(cm)
thresh = cm.max() / 2.
for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
plt.text(j, i, cm[i, j],
horizontalalignment="center",
color="white" if cm[i, j] > thresh else "black")
plt.tight_layout()
plt.ylabel('True label')
plt.xlabel('Predicted label')
font = {'size' : 15}
plt.rc('font', **font)
cnf_matrix = confusion_matrix(y_test, grid_logit.predict(X_test_new))
plt.figure(figsize=(8, 6))
plot_confusion_matrix(cnf_matrix, classes=['No default', 'default'], title=labels[0]+' Confusion matrix')
plt.show()
font = {'size' : 15}
plt.rc('font', **font)
cnf_matrix = confusion_matrix(y_test, grid_boosting.predict(X_test_new))
plt.figure(figsize=(8, 6))
plot_confusion_matrix(cnf_matrix, classes=['No default', 'default'], title=labels[1]+' Confusion matrix')
plt.show()
font = {'size' : 15}
plt.rc('font', **font)
cnf_matrix = confusion_matrix(y_test, grid_forest.predict(X_test_new))
plt.figure(figsize=(8, 6))
plot_confusion_matrix(cnf_matrix, classes=['No default', 'default'], title=labels[2]+' Confusion matrix')
plt.show()
scores = pd.DataFrame({'roc_auc_score': roc_scores, 'f1_score': f1_scores, 'accuracy': acc_scores},
index=labels)
scores
Построены модели предсказания оплаты долга клиентом в следующем месяце по информации о предыдущих шести месяцев. Модель градиентного бустинга дает примерно 79% точности на отложенных 30% выборки, но по матрицам ошибок видно, что градиентный бустинг сильно ошибается, предсказывая класс "нет уплаты" ("No default"), когда на самом деле неуплата ("default") происходит.
Были построены кривые обучения и валидационные кривые. Кривые обучения показывают, что при увеличении размера выборки, кривые будут сходиться. В данном случае, объем выборки недостаточен.
По валидационным кривым видно, что гиперпараметры, подобранные с помощью GridSearchCV являются оптимальными.