Это третья часть статьи, описывающая модель машинного обучения для прогнозирования дефолтов по кредитам. Первая часть описывает бизнес-проблему, обрисовывает в общих чертах модель и обсуждает ее влияние на бизнес. Вторая часть — это глубокое погружение в подготовку данных. Эта часть представляет собой подробное описание моделирования и оценки модели. Он основан на этой записной книжке Kaggle.
Моделирование
У нас есть довольно большой (560 тыс. наблюдений) набор табличных данных. Таким образом, я использую модель XGBoost. Во-первых, я подогнал базовую модель XGBoost с выбором гиперпараметров по умолчанию. Затем я подгоню другую модель XGBoost с гиперпараметрами, настроенную с помощью Optuna.
time1 = time.time()
xgb = XGBClassifier(tree_method = 'gpu_hist')
xgb.fit(X_train, y_train)
precision_t, recall_t, threshold = \
precision_recall_curve(y_train, xgb.predict_proba(X_train)[:, 1])
auc_precision_recall_train = auc(recall_t, precision_t)
temp = recall_t[(recall_t>0.095)&(recall_t<0.105)]
temp = temp[int(len(temp)/2)]
indexx = ((np.where(recall_t==temp)))[0][0]
r10prec_train = precision_t[indexx]
precision_t, recall_t, threshold = \
precision_recall_curve(y_test, xgb.predict_proba(X_test)[:, 1])
auc_precision_recall_test = auc(recall_t, precision_t)
temp = recall_t[(recall_t>0.095)&(recall_t<0.105)]
temp = temp[int(len(temp)/2)]
indexx = ((np.where(recall_t==temp)))[0][0]
r10prec_test = precision_t[indexx]
display('Train Accuracy: ', accuracy_score(y_train,xgb.predict(X_train)))
display('F1 score: ', f1_score(y_train,xgb.predict(X_train)))
display('ROCAUC: ', roc_auc_score(y_train,xgb.predict(X_train)))
display('PRAUC: ', auc_precision_recall_train)
display('R10P: ', r10prec_train)
# Performance evaluation:
display('Test Accuracy: ', accuracy_score(y_test,xgb.predict(X_test)))
display('F1 score: ', f1_score(y_test,xgb.predict(X_test)))
display('ROCAUC: ', roc_auc_score(y_test,xgb.predict(X_test)))
display('PRAUC: ', auc_precision_recall_test)
display('R10P: ', r10prec_test)
display(time.time()-time1)
fig, ax = plt.subplots()
ax.plot(recall_t, precision_t, color='purple')
ax.set_title('Precision-Recall Curve, test')
ax.set_ylabel('Precision')
ax.set_xlabel('Recall')
ax.set_ylim(bottom=0, top=1.02)
plt.show()
Эта модель XGBoost имеет 57% (55%) ROCAUC в обучающей (тестовой) выборке. Точность при 10% отзыве составляет 74% (58%). Это означает, что хотя модель плохо работает для большинства кредитов, есть некоторые кредиты, которые она может надежно идентифицировать как рискованные. В большинстве оценочных показателей разница между показателями обучающей выборки и тестовой выборки незначительна. Таким образом, доказательств переоснащения мало.
Код ниже использует Optuna для настройки гиперпараметров. Во-первых, я определяю пространство поиска гиперпараметров. Затем Optuna выполняет итеративный поиск в этом пространстве, чтобы максимизировать показатель ROCAUC.
### Fit XGBoost using Optuna hyperparameter optimization ###
time1 = time.time()
def objective(trial, cv_runs=1, n_splits=2, n_jobs=-1):
cv_regularizer=0.01
# Usually values between 0.1 and 0.2 work fine.
params = {
"tree_method": 'gpu_hist',
"verbosity": 0, # 0 (silent) - 3 (debug)
"n_estimators": trial.suggest_int("n_estimators", 500, 1500),
"max_depth": trial.suggest_int("max_depth", 2, 6),
"learning_rate": trial.suggest_uniform("learning_rate", 0.005, 0.2),
"colsample_bytree": trial.suggest_uniform("colsample_bytree", 0.1, 0.95),
"subsample": trial.suggest_uniform("subsample", 0.5, 0.95),
"alpha": trial.suggest_loguniform("alpha", 0.1, 10.0),
"lambda": trial.suggest_loguniform("lambda", 0.1, 150.0),
"gamma": trial.suggest_loguniform("gamma", 1e-10, 10.0),
"min_child_weight": trial.suggest_loguniform("min_child_weight", 0.1, 10),
"n_jobs": n_jobs,
"scale_pos_weight": trial.suggest_uniform("scale_pos_weight", 1, 3),
}
# usually it makes sense to resrtict hyperparameter space from some solutions which Optuna will find
# e.g., for tmx-joined data only (downsampled tmx), optuna keeps selecting depths of 2 and 3.
# for my purposes (smooth left side of prc, close to 1), those solutions are no good.
temp_out = []
for i in range(cv_runs):
X = X_train
y = y_train
model = XGBClassifier(**params)
rkf = KFold(n_splits=n_splits, shuffle=True)
X_values = X.values
y_values = y.values
y_pred = np.zeros_like(y_values)
y_pred_train = np.zeros_like(y_values)
for train_index, test_index in rkf.split(X_values):
X_A, X_B = X_values[train_index, :], X_values[test_index, :]
y_A, y_B = y_values[train_index], y_values[test_index]
model.fit(X_A, y_A, eval_set=[(X_B, y_B)], verbose = False)
y_pred[test_index] += model.predict_proba(X_B)[:, 1]
#y_pred_train[train_index] += model.predict_proba(X_A)[:, 1]
#precision_t, recall_t, threshold = precision_recall_curve(y_train, y_pred_train)
#score_train = auc(recall_t, precision_t)
precision_t, recall_t, threshold = precision_recall_curve(y_train, y_pred)
score_test = auc(recall_t, precision_t)
#score_train = roc_auc_score(y_train, y_pred_train)
#score_test = roc_auc_score(y_train, y_pred)
#overfit = score_train-score_test
#temp_out.append(score_test-cv_regularizer*overfit)
temp_out.append(score_test)
return (np.mean(temp_out))
study = optuna.create_study(direction="maximize")
study.optimize(objective, n_trials=50)
study = optuna.create_study(direction="maximize")
study.enqueue_trial(starting_hyperparameters)
study.optimize(objective, n_trials=50)
print('Total time for hypermarameter optimization ', time.time()-time1)
hp = study.best_params
for key, value in hp.items():
print(f"{key:>20s} : {value}")
print(f"{'best objective value':>20s} : {study.best_value}")
optuna_hyperpars = study.best_params
optuna_hyperpars['tree_method']='gpu_hist'
optuna_hyperpars['scale_pos_weight']=1
Я использую словарь «starting_hyperprameters» для сохранения гиперпараметров, на которых сошлась Optuna. Затем я подгоняю модель XGBoost с этими гиперпараметрами.
starting_hyperparameters = {'n_estimators': 1186,
'max_depth': 4,
'learning_rate': 0.04582452146963614,
'colsample_bytree': 0.600410879432218,
'subsample': 0.5561312005431341,
'alpha': 5.0981168922187665,
'lambda': 4.090612170344411,
'gamma': 0.0029390495081042726,
'min_child_weight': 0.21653462492301606,
'scale_pos_weight': 1,
'tree_method': 'gpu_hist'}
optuna_hyperpars = starting_hyperparameters
optuna_xgb = XGBClassifier(**optuna_hyperpars, seed=8)
optuna_xgb.fit(X_train, y_train)
precision_t, recall_t, threshold = \
precision_recall_curve(y_train, optuna_xgb.predict_proba(X_train)[:, 1])
auc_precision_recall_train = auc(recall_t, precision_t)
temp = recall_t[(recall_t>0.095)&(recall_t<0.105)]
temp = temp[int(len(temp)/2)]
indexx = ((np.where(recall_t==temp)))[0][0]
r10prec_train = precision_t[indexx]
fig, ax = plt.subplots()
ax.plot(recall_t, precision_t, color='purple')
ax.set_title('Precision-Recall Curve, train')
ax.set_ylabel('Precision')
ax.set_xlabel('Recall')
ax.set_ylim(bottom=0, top=1.02)
plt.show()
precision_t, recall_t, threshold = \
precision_recall_curve(y_test, optuna_xgb.predict_proba(X_test)[:, 1])
auc_precision_recall_test = auc(recall_t, precision_t)
temp = recall_t[(recall_t>0.095)&(recall_t<0.105)]
temp = temp[int(len(temp)/2)]
indexx = ((np.where(recall_t==temp)))[0][0]
r10prec_test = precision_t[indexx]
fig, ax = plt.subplots()
ax.plot(recall_t, precision_t, color='purple')
ax.set_title('Precision-Recall Curve, test')
ax.set_ylabel('Precision')
ax.set_xlabel('Recall')
ax.set_ylim(bottom=0, top=1.02)
plt.show()
print('Train Accuracy: ', accuracy_score(y_train,optuna_xgb.predict(X_train)))
print('F1 score: ', f1_score(y_train,optuna_xgb.predict(X_train)))
print('ROCAUC: ', roc_auc_score(y_train,optuna_xgb.predict(X_train)))
print('PRAUC: ', auc_precision_recall_train)
print('R10P: ', r10prec_train)
# Performance evaluation:
print('Test Accuracy: ', accuracy_score(y_test,optuna_xgb.predict(X_test)))
print('F1 score: ', f1_score(y_test,optuna_xgb.predict(X_test)))
print('ROCAUC: ', roc_auc_score(y_test,optuna_xgb.predict(X_test)))
print('PRAUC: ', auc_precision_recall_test)
print('R10P: ', r10prec_test)
print('Time to do hyperparameter optimization: ', time.time()-time1)
Производительность этой модели аналогична базовой модели XGBoost, описанной выше. Эта модель обеспечивает лучшую точность при 10% отзыве (60% против 58%), но имеет несколько худшие показатели F1 и ROCAUC. Это означает, что от оптимизации гиперпараметров не так много пользы. Этот результат не очень удивителен, поскольку известно, что XGBoost представляет собой семейство моделей, устойчивых к выбору гиперпараметров.
Интерпретация модели
Отправной точкой в интерпретации модели является важность признаков. В приведенном ниже коде используется пакет SHAP для отображения 30 наиболее важных функций.
# template here:
# https://www.kaggle.com/code/kaanboke/catboost-lightgbm-xgboost-explained-by-shap/notebook
explainerxgbc = shap.TreeExplainer(optuna_xgb)
shap_values_XGBoost_test = explainerxgbc.shap_values(X_test)
shap_values_XGBoost_train = explainerxgbc.shap_values(X_train)
vals = np.abs(shap_values_XGBoost_test).mean(0)
feature_names = X_test.columns
feature_importance = pd.DataFrame(list(zip(feature_names, vals)),
columns=['col_name','feature_importance_vals'])
feature_importance.sort_values(by=['feature_importance_vals'],
ascending=False, inplace=True)
#display(feature_importance)
shap.summary_plot(shap_values_XGBoost_test,
X_test, plot_type="bar", plot_size=(8,8), max_display=20)
shap.summary_plot(shap_values_XGBoost_train,
X_train,plot_type="dot", plot_size=(8,8), max_display=20)
На приведенном выше графике показаны наиболее важные функции и их связь с целью. Подземный переход – самая важная характеристика. Высокие значения грунтового основания (например, E, F, G) означают более низкое качество кредита. График показывает, что высокие значения грунтового основания имеют высокие значения SHAP. Другими словами, кредиты более низкого качества имеют более высокую вероятность дефолта. Точно так же высокие коэффициенты LTI и DTI связаны с более высоким риском дефолта. Все эти результаты интуитивно понятны.
Чтобы получить представление о внутренней работе модели, рассмотрим два случая. Во-первых, посмотрите на кредит, который модель считает безопасным:
indx = 1
#fig = plt.subplots(figsize=(2,2),dpi=200)
ax_2 = shap.decision_plot(explainerxgbc.expected_value,
shap_values_XGBoost_test[indx],
X_test.iloc[indx],link= "logit")
shap.initjs()
shap.force_plot(explainerxgbc.expected_value,
shap_values_XGBoost_test[indx],
X_test.iloc[[indx]],link= "logit")
Модель оценивает риск дефолта по этому кредиту в 7%. На двух рисунках ниже показано, какие функции использует модель, чтобы прийти к такому выводу.
Кредит имеет sub_grade A2. Целевая кодировка sub_grade — 0,048. Это значение имеет большое понижающее влияние на прогнозируемую вероятность дефолта. Процентная ставка составляет 6%, что ниже среднего 13%. График показывает, что такая процентная ставка соответствует уменьшенному риску дефолта. Точно так же высокий балл FICO заставляет модель воспринимать ссуду как менее рискованную.
Теперь давайте посмотрим на кредит, который модель считает рискованным:
indx = 10
#fig = plt.subplots(figsize=(2,2),dpi=200)
ax_2= shap.decision_plot(explainerxgbc.expected_value,
shap_values_XGBoost_test[indx],
X_test.iloc[indx],link= "logit")
shap.initjs()
shap.force_plot(explainerxgbc.expected_value,
shap_values_XGBoost_test[indx],
X_test.iloc[[indx]],link= "logit")
Прогнозируемая вероятность дефолта составляет 52%. Это намного выше, чем средний уровень дефолтов в 20%. Subgrade E5 имеет целевое кодирование 0,423. Это единственный наиболее важный фактор высокой прогнозируемой вероятности дефолта. Модель воспринимает этот кредит как более рискованный, поскольку его LTI составляет 0,42. Процентная ставка до 25% является еще одной причиной высокого прогнозируемого риска дефолта.
Влияние на бизнес
Этот проект прекрасно иллюстрирует, как машинное обучение может создать ценность, если мы тщательно обдумаем стоящую бизнес-проблему, а не сосредоточимся на нескольких узких показателях оценки модели.
Если вы рассматриваете этот проект исключительно как упражнение по интеллектуальному анализу данных/прогнозированию, то, возможно, это провал. Точность модели составляет 81%, что едва превышает 80% по сравнению с тривиальной моделью, всегда предсказывающей отсутствие дефолта. Оценка F1 составляет всего 18%, что не впечатляет. ROCAUC составляет 55%, что чуть лучше, чем случайное предположение.
Однако приведенная выше модель прогнозирования может создать ценность при правильном использовании. Позвольте мне показать вам, как это сделать.
Приведенная выше кривая точности-отзыва показывает, что модель на самом деле очень хорошо предсказывает дефолты примерно для 20% просроченных кредитов. Например, его точность при 10% отзыве составляет 60%. Это означает, что если мы используем модель только для прогнозирования дефолта по тем кредитам, в которых она наиболее уверена, мы можем избежать по крайней мере 10% кредитных убытков, отказываясь при этом от относительно небольшого количества выгодных возможностей кредитования.
Чтобы количественно оценить ценность, создаваемую моделью, я рассматриваю проблему, с которой сталкиваются все инвесторы Lending Club как группа, стремящиеся максимизировать свою инвестиционную прибыль. Они могут сделать это, избегая x% самых рискованных кредитов, определенных моей моделью. Чтобы сделать этот анализ простым и доступным для более широкой аудитории, я абстрагируюсь от временного дисконтирования, сложных процентов и расчетов текущей стоимости. Все приведенные ниже расчеты предполагают простые проценты.
Я покажу, что инвесторы выиграют, если будут избегать 10% самых рискованных кредитов и вместо этого вкладывать эти деньги в безрисковый актив (казначейские облигации). Значение, созданное моделью, представляет собой разницу в доходах между двумя стратегиями. Первый заключается в том, чтобы инвестировать в кредиты, кроме 10% наиболее рискованных, а оставшиеся 10% средств инвестировать в безрисковые активы. Второй — это инвестиции во все кредиты LendingClub. Выявляя самые рискованные кредиты, эта модель позволяет инвесторам экономить деньги.
Неясно заранее, создаст ли такая модель ценность. Более рискованные кредиты обычно имеют более высокую процентную ставку. Избегая таких кредитов, мы теряем процентный доход. Lending Club использовал собственную модель риска для установления процентных ставок. Таким образом, моя модель может приносить инвесторам прибыль тогда и только тогда, когда для некоторых кредитов она предсказывает дефолты лучше, чем модель Lending Club.
Приведенный ниже код рассчитывает прибыль от отказа от 10% самых рискованных кредитов. Он учитывает скорость возмещения и процентный доход от кредитов.
recoveries['total_recovery'] = recoveries.total_rec_prncp + recoveries.total_rec_int + recoveries.recoveries
recoveries['tot_recov_rp'] = recoveries.total_recovery/recoveries.loan_amnt
recoveries['tot_recov_rt'] = recoveries.total_recovery/(recoveries.loan_amnt*((recoveries.int_rate/100+1)**3))
# in a few cases the recoveries seem to exceed total proceeds.
# this may be the result of ignoring time discounting and compound interest in my calculations
# ior those may be rarer 60-months loans. so i do not adjust such cases.
display(recoveries.describe(), recoveries.head())
X_test = X_test_0.copy()
test = X_test[['loan_amnt', 'int_rate']]
test.int_rate = test.int_rate/100+1
test['y_pred'] = optuna_xgb.predict_proba(X_test)[:,1]
test['id'] = test00.id
test['y'] = y_test
test = pd.merge(test, recoveries[['id', 'loan_amnt', 'int_rate', 'total_recovery']], on='id', how = 'left')
display(recoveries.head(), test.head())
display(test.loc[test.y==0].count(), test.loc[test.y==1].count())
# select riskiest loans by calculating decision threshold, giving 10% recall:
desired_recall = 0.1
temp = recall_t[(recall_t>(desired_recall-0.001))&(recall_t<(desired_recall+0.001))]
temp = temp[int(len(temp)/2)]
indexx = ((np.where(recall_t==temp)))[0][0]
r10threshold = threshold[indexx]
p90risk = r10threshold
test['p90risk'] = (test.y_pred>=p90risk).astype(int)
display(test.shape, p90risk)
risky_loans = test[test.p90risk==1]
risky_loans.loc[risky_loans.total_recovery.isnull(),'total_recovery']=\
risky_loans.loan_amnt_x*(risky_loans.int_rate_x**3)
# when the loan is repaid, I calculate total return and save it into total_recovery column.
display(risky_loans.head(), risky_loans.shape)
proceeds = risky_loans.total_recovery.sum()
proceeds_tnotes = (risky_loans.loan_amnt_x.sum())*(1.02**3)
print("Investors' total returns from investing in risky loans", int(proceeds))
print("Investors' total returns from investing in T-notes", int(proceeds_tnotes))
print("Investors' savings: $", int(proceeds_tnotes-proceeds))
estimated_savings = (proceeds_tnotes-proceeds)*3.3*4
print("Estimated investors' savings from all LendingClub loans: $", estimated_savings)
# In this particular sample, the savings are maximized at 8.5% recall threshold.
# Then the saving for investors will be $4.9M in test set or 64.6M for all loans.
print('Total time to run the script: ', time.time() - time0)
Эта модель создает ценность для инвесторов, предоставляющих ссуды P2P. Используя его, чтобы избежать самых рискованных кредитов, инвесторы сэкономили бы 60 млн долларов США.
На этом построение модели заканчивается. В следующей части мы обсудим, как запустить его в производство.
