[Из песочницы] Соревнование Kaggle Home Credit Default Risk — анализ данных и простые предсказательные модели20.06.2018 06:02

На датафесте 2 в Минске Владимир Игловиков, инженер по машинному зрению в Lyft, совершенно замечательно объяснил, что лучший способ научиться Data Science — это участвовать в соревнованиях, запускать чужие решения, комбинировать их, добиваться результата и показывать свою работу. Собственно в рамках этой парадигмы я и решил посмотреть внимательнее на соревнование по оценке кредитного риска от Home Credit и объяснить (начинающим дата саентистам и прежде всего самому себе), как правильно анализировать подобные датасеты и строить под них модели.

(картинка отсюда)

Home Credit Group — группа банков и небанковских кредитных организаций, ведет операции в 11 странах (в том числе в России как ООО «Хоум Кредит энд Финанс Банк»). Цель соревнования — создать методику оценки кредитоспособности заемщиков, не имеющих кредитной истории. Что выглядит довольно благородно — заемщики этой категории часто не могут получить никакой кредит в банке и вынуждены обращаться к мошенникам и микрозаймам. Интересно, что заказчик не выставляет требований по прозрачности и интерпретируемости модели (как это обычно бывает в банках), можно использовать что угодно, хоть нейросети.

Обучающая выборка состоит из 300+ тыс. записей, признаков достаточно много — 122, среди них много категориальных (не числовых). Признаки довольно подробно описывают заемщика, вплоть до материала, из которого сделаны стены его жилища. Часть данных содержатся в 6 дополнительных таблицах (данные по кредитному бюро, балансу кредитной карты и предыдущим кредитам), эти данные нужно тоже как-то обработать и подгрузить к основным.

Соревнование выглядит как стандартная задача классификации (1 в поле TARGET означает любые сложности с платежами, 0 — отсутствие сложностей). Однако следует предсказывать не 0/1, а вероятность возникновения проблем (что, впрочем, довольно легко решают методы предсказания вероятностей predict_proba, которые есть у всех сложных моделей).

На первый взгляд датасет довольно стандартной для задач машинного обучения, организаторы предложили крупный приз в $70к, в итоге в соревновании на сегодня участвует уже больше 2600 команд, а битва идет за тысячные доли процентов. Однако, с другой стороны, такая популярность означает, что датасет исследован вдоль и поперек и создано много кернелов с хорошими EDA (Exploratory Data Analisys — исследование и анализ данных в сете, в том числе графический), Feature engineering’ом (работа с признаками) и с интересными моделями. (Кернел — это пример работы с датасетом, который может выложить любой желающий, чтобы показать свою работу другим кагглерам.)

Заслуживают внимания кернелы:

Для работы с данными обычно рекомендуется следующий план, которому мы и постараемся следовать.

Понимание проблемы и ознакомление с данными
Чистка данных и форматирование
EDA
Базовая модель
Улучшение модели
Интерпретация модели

В данном случае нужно взять поправку на то, что данные довольно обширны и сразу их можно и не осилить, есть смысл действовать поэтапно.

Начнем с импорта библиотек, которые нам понадобятся в анализе для работы с данными в виде таблиц, построения графиков и для работы с матрицами.

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
%matplotlib inline

Загрузим данные. Посмотрим, что у нас вообще есть. Такое расположение в каталоге »…/input/», кстати, связано с требованием по размещению своих кернелов на Kaggle.

import os
PATH="../input/"
print(os.listdir(PATH))

['application_test.csv', 'application_train.csv', 'bureau.csv', 'bureau_balance.csv', 'credit_card_balance.csv', 'HomeCredit_columns_description.csv', 'installments_payments.csv', 'POS_CASH_balance.csv', 'previous_application.csv']

Есть 8 таблиц с данными (не считая таблицы HomeCredit_columns_description.csv, в которой содержится описание полей), которые связаны между собой следующим образом:

application_train/application_test: Основные данные, заемщик идентифицируется по полю SK_ID_CURR
bureau: Данные по предыдущим займам в других кредитных организациях из кредитного бюро
bureau_balance: Ежемесячные данные по предыдущим кредитам по бюро. Каждая строка — месяц испльзования кредита
previous_application: Предыдущие заявки по кредитам в Home Credit, каждая имеет уникальное поле SK_ID_PREV
POS_CASH_BALANCE: Ежемесячные данные по кредитам в Home Creditс выдачей наличными и кредитам на покупки товаров
credit_card_balance: Ежемесячные данные по балансу кредитных карт в Home Credit
installments_payment: Платежная история предыдущих займов в Home Credit.

Сосредоточимся для начала на основном источнике данных и посмотрим, какую информацию из него можно извлечь и какие модели построить. Загрузим основные данные.

app_train = pd.read_csv (PATH + 'application_train.csv',)
app_test = pd.read_csv (PATH + 'application_test.csv',)
print («формат обучающей выборки:», app_train.shape)
print («формат тестовой выборки:», app_test.shape)
формат обучающей выборки: (307511, 122)
формат тестовой выборки: (48744, 121)

Итого у нас есть 307 тысяч записей и 122 признака в обучающей выборке и 49 тысяч записей и 121 признак в тестовой. Расхождение, очевидно, вызвано тем, что целевого признака TARGET в тестовой выборке нет, его-то мы и будем предсказывать.

Посмотрим на данные внимательнее

pd.set_option('display.max_columns', None) # иначе pandas не покажет все столбцы
app_train.head()

(показано первые 8 столбцов)

Довольно сложно смотреть данные в таком формате. Посмотрим на список столбцов:

app_train.info(max_cols=122) RangeIndex: 307511 entries, 0 to 307510 Data columns (total 122 columns): SK_ID_CURR 307511 non-null int64 TARGET 307511 non-null int64 NAME_CONTRACT_TYPE 307511 non-null object CODE_GENDER 307511 non-null object FLAG_OWN_CAR 307511 non-null object FLAG_OWN_REALTY 307511 non-null object CNT_CHILDREN 307511 non-null int64 AMT_INCOME_TOTAL 307511 non-null float64 AMT_CREDIT 307511 non-null float64 AMT_ANNUITY 307499 non-null float64 AMT_GOODS_PRICE 307233 non-null float64 NAME_TYPE_SUITE 306219 non-null object NAME_INCOME_TYPE 307511 non-null object NAME_EDUCATION_TYPE 307511 non-null object NAME_FAMILY_STATUS 307511 non-null object NAME_HOUSING_TYPE 307511 non-null object REGION_POPULATION_RELATIVE 307511 non-null float64 DAYS_BIRTH 307511 non-null int64 DAYS_EMPLOYED 307511 non-null int64 DAYS_REGISTRATION 307511 non-null float64 DAYS_ID_PUBLISH 307511 non-null int64 OWN_CAR_AGE 104582 non-null float64 FLAG_MOBIL 307511 non-null int64 FLAG_EMP_PHONE 307511 non-null int64 FLAG_WORK_PHONE 307511 non-null int64 FLAG_CONT_MOBILE 307511 non-null int64 FLAG_PHONE 307511 non-null int64 FLAG_EMAIL 307511 non-null int64 OCCUPATION_TYPE 211120 non-null object CNT_FAM_MEMBERS 307509 non-null float64 REGION_RATING_CLIENT 307511 non-null int64 REGION_RATING_CLIENT_W_CITY 307511 non-null int64 WEEKDAY_APPR_PROCESS_START 307511 non-null object HOUR_APPR_PROCESS_START 307511 non-null int64 REG_REGION_NOT_LIVE_REGION 307511 non-null int64 REG_REGION_NOT_WORK_REGION 307511 non-null int64 LIVE_REGION_NOT_WORK_REGION 307511 non-null int64 REG_CITY_NOT_LIVE_CITY 307511 non-null int64 REG_CITY_NOT_WORK_CITY 307511 non-null int64 LIVE_CITY_NOT_WORK_CITY 307511 non-null int64 ORGANIZATION_TYPE 307511 non-null object EXT_SOURCE_1 134133 non-null float64 EXT_SOURCE_2 306851 non-null float64 EXT_SOURCE_3 246546 non-null float64 APARTMENTS_AVG 151450 non-null float64 BASEMENTAREA_AVG 127568 non-null float64 YEARS_BEGINEXPLUATATION_AVG 157504 non-null float64 YEARS_BUILD_AVG 103023 non-null float64 COMMONAREA_AVG 92646 non-null float64 ELEVATORS_AVG 143620 non-null float64 ENTRANCES_AVG 152683 non-null float64 FLOORSMAX_AVG 154491 non-null float64 FLOORSMIN_AVG 98869 non-null float64 LANDAREA_AVG 124921 non-null float64 LIVINGAPARTMENTS_AVG 97312 non-null float64 LIVINGAREA_AVG 153161 non-null float64 NONLIVINGAPARTMENTS_AVG 93997 non-null float64 NONLIVINGAREA_AVG 137829 non-null float64 APARTMENTS_MODE 151450 non-null float64 BASEMENTAREA_MODE 127568 non-null float64 YEARS_BEGINEXPLUATATION_MODE 157504 non-null float64 YEARS_BUILD_MODE 103023 non-null float64 COMMONAREA_MODE 92646 non-null float64 ELEVATORS_MODE 143620 non-null float64 ENTRANCES_MODE 152683 non-null float64 FLOORSMAX_MODE 154491 non-null float64 FLOORSMIN_MODE 98869 non-null float64 LANDAREA_MODE 124921 non-null float64 LIVINGAPARTMENTS_MODE 97312 non-null float64 LIVINGAREA_MODE 153161 non-null float64 NONLIVINGAPARTMENTS_MODE 93997 non-null float64 NONLIVINGAREA_MODE 137829 non-null float64 APARTMENTS_MEDI 151450 non-null float64 BASEMENTAREA_MEDI 127568 non-null float64 YEARS_BEGINEXPLUATATION_MEDI 157504 non-null float64 YEARS_BUILD_MEDI 103023 non-null float64 COMMONAREA_MEDI 92646 non-null float64 ELEVATORS_MEDI 143620 non-null float64 ENTRANCES_MEDI 152683 non-null float64 FLOORSMAX_MEDI 154491 non-null float64 FLOORSMIN_MEDI 98869 non-null float64 LANDAREA_MEDI 124921 non-null float64 LIVINGAPARTMENTS_MEDI 97312 non-null float64 LIVINGAREA_MEDI 153161 non-null float64 NONLIVINGAPARTMENTS_MEDI 93997 non-null float64 NONLIVINGAREA_MEDI 137829 non-null float64 FONDKAPREMONT_MODE 97216 non-null object HOUSETYPE_MODE 153214 non-null object TOTALAREA_MODE 159080 non-null float64 WALLSMATERIAL_MODE 151170 non-null object EMERGENCYSTATE_MODE 161756 non-null object OBS_30_CNT_SOCIAL_CIRCLE 306490 non-null float64 DEF_30_CNT_SOCIAL_CIRCLE 306490 non-null float64 OBS_60_CNT_SOCIAL_CIRCLE 306490 non-null float64 DEF_60_CNT_SOCIAL_CIRCLE 306490 non-null float64 DAYS_LAST_PHONE_CHANGE 307510 non-null float64 FLAG_DOCUMENT_2 307511 non-null int64 FLAG_DOCUMENT_3 307511 non-null int64 FLAG_DOCUMENT_4 307511 non-null int64 FLAG_DOCUMENT_5 307511 non-null int64 FLAG_DOCUMENT_6 307511 non-null int64 FLAG_DOCUMENT_7 307511 non-null int64 FLAG_DOCUMENT_8 307511 non-null int64 FLAG_DOCUMENT_9 307511 non-null int64 FLAG_DOCUMENT_10 307511 non-null int64 FLAG_DOCUMENT_11 307511 non-null int64 FLAG_DOCUMENT_12 307511 non-null int64 FLAG_DOCUMENT_13 307511 non-null int64 FLAG_DOCUMENT_14 307511 non-null int64 FLAG_DOCUMENT_15 307511 non-null int64 FLAG_DOCUMENT_16 307511 non-null int64 FLAG_DOCUMENT_17 307511 non-null int64 FLAG_DOCUMENT_18 307511 non-null int64 FLAG_DOCUMENT_19 307511 non-null int64 FLAG_DOCUMENT_20 307511 non-null int64 FLAG_DOCUMENT_21 307511 non-null int64 AMT_REQ_CREDIT_BUREAU_HOUR 265992 non-null float64 AMT_REQ_CREDIT_BUREAU_DAY 265992 non-null float64 AMT_REQ_CREDIT_BUREAU_WEEK 265992 non-null float64 AMT_REQ_CREDIT_BUREAU_MON 265992 non-null float64 AMT_REQ_CREDIT_BUREAU_QRT 265992 non-null float64 AMT_REQ_CREDIT_BUREAU_YEAR 265992 non-null float64 dtypes: float64(65), int64(41), object(16) memory usage: 286.2+ MB

Напомню, подробные аннотации по полям — в файле HomeCredit_columns_description. Как видно из info, часть данных неполная и часть — категориальная, они отображаются как object. Большинство моделей с такими данными не работают, нам придется что-то с этим делать. На этом начальный анализ можно считать законченным, перейдем непосредственно к EDA

Exploratory Data Analysis или первичное исследование данных

В процессе EDA мы считаем основные статистики и рисуем графики, чтобы найти тренды, аномалии, паттерны и связи внутри данных. Цель EDA — узнать, что могут рассказать данные. Обычно анализ идет сверху вниз — от общего обзора к исследованию отдельных зон, которые привлекают внимание и могут представлять интерес. Впоследствии эти находки можно использовать в построении модели, выборе признаков для нее и в её интерпретации.

Распределение целевой переменной

app_train.TARGET.value_counts()

0 282686 1 24825 Name: TARGET, dtype: int64

plt.style.use('fivethirtyeight')
plt.rcParams["figure.figsize"] = [8,5]

plt.hist(app_train.TARGET)
plt.show()

Напомню, 1 означает проблемы любого рода с возвратом, 0 — отсутствие проблем. Как видно, в основном заемщики не имеют проблем с возвратом, доля проблемных около 8%. Это значит, что классы не сбалансированы и это, возможно, нужно будет учитывать при построении модели.

Исследование недостающих данных

Мы видели, что нехватка данных довольно существенна. Посмотрим более подробно, где и чего недостает.

# Функция для подсчета недостающих столбцов
def missing_values_table(df):
    
        # Всего недостает
        mis_val = df.isnull().sum()
        
        # Процент недостающих данных
        mis_val_percent = 100 * df.isnull().sum() / len(df)
        
        # Таблица с результатами
        mis_val_table = pd.concat([mis_val, mis_val_percent], axis=1)
        
        # Переименование столбцов
        mis_val_table_ren_columns = mis_val_table.rename(
        columns = {0 : 'Missing Values', 1 : '% of Total Values'})
        
        # Сортировка про процентажу
        mis_val_table_ren_columns = mis_val_table_ren_columns[
            mis_val_table_ren_columns.iloc[:,1] != 0].sort_values(
        '% of Total Values', ascending=False).round(1)
        
        # Инфо
        print ("В выбранном датафрейме " + str(df.shape[1]) + " столбцов.\n"      
            "Всего " + str(mis_val_table_ren_columns.shape[0]) +
              " столбцов с неполными данными.")
        
        # Возврат таблицы с данными
        return mis_val_table_ren_columns
    
missing_values = missing_values_table(app_train)
missing_values.head(10)

В выбранном датафрейме 122 столбцов. Всего 67 столбцов с неполными данными.

В графическом формате:

plt.style.use('seaborn-talk')

fig = plt.figure(figsize=(18,6))
miss_train = pd.DataFrame((app_train.isnull().sum())*100/app_train.shape[0]).reset_index()
miss_test = pd.DataFrame((app_test.isnull().sum())*100/app_test.shape[0]).reset_index()
miss_train["type"] = "тренировочная"
miss_test["type"]  =  "тестовая"
missing = pd.concat([miss_train,miss_test],axis=0)
ax = sns.pointplot("index",0,data=missing,hue="type")
plt.xticks(rotation =90,fontsize =7)
plt.title("Доля отсуствующих значений в данных")
plt.ylabel("Доля в %")
plt.xlabel("Столбцы")

На вопрос «что со всем этим делать» есть много ответов. Можно заполнить нулями, можно медианными значениями, можно просто удалить строки без нужной информации. Все зависит от модели, которую мы планируем использовать, так как некоторые совершенно отлично справляются с отсутствующими значениями. Пока запомним этот факт и оставим все как есть.

Типы столбцов и кодирование категориальных данных

Как мы помним. часть столбцов имеет тип object, то есть имеет не числовое значение, а отражает какую-то категорию. Посмотрим на эти столбцы более внимательно.

app_train.dtypes.value_counts()

float64 65 int64 41 object 16 dtype: int64

app_train.select_dtypes(include=[object]).apply(pd.Series.nunique, axis = 0)

NAME_CONTRACT_TYPE 2 CODE_GENDER 3 FLAG_OWN_CAR 2 FLAG_OWN_REALTY 2 NAME_TYPE_SUITE 7 NAME_INCOME_TYPE 8 NAME_EDUCATION_TYPE 5 NAME_FAMILY_STATUS 6 NAME_HOUSING_TYPE 6 OCCUPATION_TYPE 18 WEEKDAY_APPR_PROCESS_START 7 ORGANIZATION_TYPE 58 FONDKAPREMONT_MODE 4 HOUSETYPE_MODE 3 WALLSMATERIAL_MODE 7 EMERGENCYSTATE_MODE 2 dtype: int64

У нас есть 16 столбцов, в каждом из которых от 2 до 58 разных вариантов значений. В основном модели машинного обучения не могут ничего сделать с таким столбцами (кроме некоторых, например LightGBM или CatBoost). Так как мы планируем опробовать разные модели на датасете, то с этим надо что-то делать. Подхода здесь два:

Label Encoding — категориям присваиваются цифры 0, 1, 2 и так далее и записываются в тот же столбец
One-Hot-encoding — один столбец раскладывается на несколько по количеству вариантов и в этих столбцах отмечается, какой вариант у данной записи.

С Label Encoding есть небольшая проблема — он присваивает числовые значения, которые не имеют ничего общего с реальностью. Например, если мы имеем дело с числовым значением, то доход заемщика в 10000 однозначно больше и лучше, чем доход в 20000. Но можно ли сказать, что, например, один город лучше другого потому, что одному присвоено значение 100, а другому — 200?

One-Hot-encoding, с другой стороны, более безопасен, но может плодить «лишние» столбцы. Например, если мы закодируем тот же пол при помощи One-Hot, у нас получится два столбца, «пол мужской» и «пол женский», хотя хватило бы и одного, «мужчина ли это».

По хорошему для данного датасета надо бы кодировать признаки с низкой вариативностью при помощи Label Encoding, а все остальное — One-Hot, но для упрощения закодируем все по One-Hot. На скорость вычисления и результат это практически не повлияет. Сам процесс кодирования pandas делает очень просто

app_train = pd.get_dummies(app_train)
app_test = pd.get_dummies(app_test)

print('Training Features shape: ', app_train.shape)
print('Testing Features shape: ', app_test.shape)

Training Features shape: (307511, 246) Testing Features shape: (48744, 242)

Так как количество вариантов в столбцах выборок не равное, количество столбцов теперь не совпадает. Требуется выравнивание — нужно убрать из тренировочной выборки столбцы, которых нет в тестовой. Это делает метод align, нужно указать axis=1 (для столбцов).

#сохраним лейблы, их же нет в тестовой выборке и при выравнивании они потеряются. 
train_labels = app_train['TARGET']

# Выравнивание - сохранятся только столбцы. имеющиеся в обоих датафреймах
app_train, app_test = app_train.align(app_test, join = 'inner', axis = 1)

print('Формат тренировочной выборки: ', app_train.shape)
print('Формат тестовой выборки: ', app_test.shape)

# Add target back in to the data
app_train['TARGET'] = train_labels

Формат тренировочной выборки: (307511, 242) Формат тестовой выборки: (48744, 242)

Корреляция в данных

Хороший метод понять данные — рассчитать коэффициенты корреляции Пирсона для данных относительно целевого признака. Это не лучший метод показать релевантность признаков, но он прост и позволяет составить представление о данных. Интерпретировать коэффициенты можно следующим образом:

00-.19 «очень слабая»
20-.39 «слабая»
40-.59 «средняя»
60-.79 «сильная»
80–1.0 «очень сильная»

# Корреляция и сортировка
correlations = app_train.corr()['TARGET'].sort_values()

# Отображение
print('Наивысшая позитивная корреляция: \n', correlations.tail(15))
print('\nНаивысшая негативная корреляция: \n', correlations.head(15))

Наивысшая позитивная корреляция: DAYS_REGISTRATION 0.041975 OCCUPATION_TYPE_Laborers 0.043019 FLAG_DOCUMENT_3 0.044346 REG_CITY_NOT_LIVE_CITY 0.044395 FLAG_EMP_PHONE 0.045982 NAME_EDUCATION_TYPE_Secondary / secondary special 0.049824 REG_CITY_NOT_WORK_CITY 0.050994 DAYS_ID_PUBLISH 0.051457 CODE_GENDER_M 0.054713 DAYS_LAST_PHONE_CHANGE 0.055218 NAME_INCOME_TYPE_Working 0.057481 REGION_RATING_CLIENT 0.058899 REGION_RATING_CLIENT_W_CITY 0.060893 DAYS_BIRTH 0.078239 TARGET 1.000000 Name: TARGET, dtype: float64

Наивысшая негативная корреляция: EXT_SOURCE_3 -0.178919 EXT_SOURCE_2 -0.160472 EXT_SOURCE_1 -0.155317 NAME_EDUCATION_TYPE_Higher education -0.056593 CODE_GENDER_F -0.054704 NAME_INCOME_TYPE_Pensioner -0.046209 ORGANIZATION_TYPE_XNA -0.045987 DAYS_EMPLOYED -0.044932 FLOORSMAX_AVG -0.044003 FLOORSMAX_MEDI -0.043768 FLOORSMAX_MODE -0.043226 EMERGENCYSTATE_MODE_No -0.042201 HOUSETYPE_MODE_block of flats -0.040594 AMT_GOODS_PRICE -0.039645 REGION_POPULATION_RELATIVE -0.037227 Name: TARGET, dtype: float64

Таким образом, все данные слабо коррелируют с таргетом (кроме самого таргета, который, понятно, равен сам себе). Однако из данных выделяются возраст и некие «внешние источники данных». Вероятно, это некие дополнительные данные из других кредитных организаций. Забавно, что хотя цель декларируется как независимость от подобных данных в принятии кредитного решения, на деле мы будем основываться в первую очередь на них.

Возраст

Понятно, что чем старше клиент, тем выше вероятность возврата (до определенного предела, конечно). Но возраст почему-то указан в отрицательных днях до выдачи кредита, поэтому он положительно коррелирует с невозвратом (что выглядит несколько странно). Приведем его к положительному значению и посмотрим на корреляцию.

app_train['DAYS_BIRTH'] = abs(app_train['DAYS_BIRTH'])
app_train['DAYS_BIRTH'].corr(app_train['TARGET'])

-0.078239308309827088

Посмотрим на переменную внимательнее. Начнем с гистограммы.

# Гистограмма распределения возраста в годах, всего 25 столбцов
plt.hist(app_train['DAYS_BIRTH'] / 365, edgecolor = 'k', bins = 25)
plt.title('Age of Client'); plt.xlabel('Age (years)'); plt.ylabel('Count');

Сама по себе гистограмма распределения может сказать немного полезного, кроме того что мы не видим особых выбросов и все выглядит более-менее правдоподобно. Чтобы показать эффект влияния возраста на результат, можно построить график kernel density estimation (KDE) — распределение ядерной плотности, раскрашенный в цвета целевого признака. Он показывает распределение одной переменной и может быть истолкован как сглаженная гистограмма (рассчитывается как Гауссианское ядро по каждой точке, которое затем усредняется для сглаживания).

# KDE займов, выплаченных вовремя
sns.kdeplot(app_train.loc[app_train['TARGET'] == 0, 'DAYS_BIRTH'] / 365, label = 'target == 0')

# KDE проблемных займов
sns.kdeplot(app_train.loc[app_train['TARGET'] == 1, 'DAYS_BIRTH'] / 365, label = 'target == 1')

# Обозначения
plt.xlabel('Age (years)'); plt.ylabel('Density'); plt.title('Distribution of Ages');

Как видно, доля невозвратов выше для молодых людей и снижается с ростом возраста. Это не повод отказывать молодым людям в кредите всегда, такая «рекомендация» приведет лишь к потере доходов и рынка для банка. Это повод задуматься о более тщательном отслеживании таких кредитов, оценке и, возможно, даже каком-то финансовом образовании для молодых заемщиков.

Внешние источники

Посмотрим внимательнее на «внешние источники данных» EXT_SOURCE и их корреляцию.

ext_data = app_train[['TARGET', 'EXT_SOURCE_1', 'EXT_SOURCE_2', 'EXT_SOURCE_3', 'DAYS_BIRTH']]
ext_data_corrs = ext_data.corr()
ext_data_corrs

Также корреляцию удобно отображать при помощи heatmap

sns.heatmap(ext_data_corrs, cmap = plt.cm.RdYlBu_r, vmin = -0.25, annot = True, vmax = 0.6)
plt.title('Correlation Heatmap');

Как видно, все источники показывают негативную корреляцию с таргетом. Посмотрим на рапсределение KDE по каждому источнику.

plt.figure(figsize = (10, 12))

# итерация по источникам
for i, source in enumerate(['EXT_SOURCE_1', 'EXT_SOURCE_2', 'EXT_SOURCE_3']):
    
    # сабплот
    plt.subplot(3, 1, i + 1)
    # отрисовка качественных займов
    sns.kdeplot(app_train.loc[app_train['TARGET'] == 0, source], label = 'target == 0')
    # отрисовка дефолтных займов
    sns.kdeplot(app_train.loc[app_train['TARGET'] == 1, source], label = 'target == 1')
    
    # метки
    plt.title('Distribution of %s by Target Value' % source)
    plt.xlabel('%s' % source); plt.ylabel('Density');
    
plt.tight_layout(h_pad = 2.5)

Картина аналогична распределению по возрасту — с ростом показателя растет вероятность возврата кредита. Третий источник наиболее силен в этом плане. Хотя в абсолютном выражении корреляция с целевой переменной все еще в категории «очень низкая», источники внешних данных и возраст будут иметь наивысшее значение в построении модели.

Парный график

Для лучшего понимания взаимоотношений этих переменных можно построить парный график, в нем мы сможем увидеть взаимоотношения каждой пары и гистограмму распределения по диагонали. Выше диагонали можно показать диаграмму рассеяния, а ниже — 2d KDE.

#вынесем данные по возрасту в отдельный датафрейм
age_data = app_train[['TARGET', 'DAYS_BIRTH']]
age_data['YEARS_BIRTH'] = age_data['DAYS_BIRTH'] / 365

# копирование данных для графика
plot_data = ext_data.drop(labels = ['DAYS_BIRTH'], axis=1).copy()

# Добавим возраст
plot_data['YEARS_BIRTH'] = age_data['YEARS_BIRTH']

# Уберем все незаполненнные строки и ограничим таблицу в 100 тыс. строк
plot_data = plot_data.dropna().loc[:100000, :]

# Функиця для расчет корреляции
def corr_func(x, y, **kwargs):
    r = np.corrcoef(x, y)[0][1]
    ax = plt.gca()
    ax.annotate("r = {:.2f}".format(r),
                xy=(.2, .8), xycoords=ax.transAxes,
                size = 20)

# Создание объекта pairgrid object
grid = sns.PairGrid(data = plot_data, size = 3, diag_sharey=False,
                    hue = 'TARGET', 
                    vars = [x for x in list(plot_data.columns) if x != 'TARGET'])

# Сверху - скаттерплот
grid.map_upper(plt.scatter, alpha = 0.2)

# Диагональ - гистограмма
grid.map_diag(sns.kdeplot)

# Внизу - распределение плотности
grid.map_lower(sns.kdeplot, cmap = plt.cm.OrRd_r);

plt.suptitle('Ext Source and Age Features Pairs Plot', size = 32, y = 1.05);

Синим показаны возвратные кредиты, красным — невозвратные. Интерпретировать это все довольно сложно, но зато из этой картинки может выйти неплохой принт на майку или картина в музей современного искусства.

Исследование прочих признаков

Рассмотрим более подробно другие признаки и их зависимость от целевой переменной. Так как среди них много категориальных (а мы уже успели их закодировать), нам снова понадобятся исходные данные. Назовем их немного по-другому во избежание путаницы

application_train = pd.read_csv(PATH+"application_train.csv")
application_test = pd.read_csv(PATH+"application_test.csv")

Также нам понадобится пара функций для красивого отображения распределений и их влияния на целевую переменную. Большое за них спасибо автору вот этого кернела

def plot_stats(feature,label_rotation=False,horizontal_layout=True):
    temp = application_train[feature].value_counts()
    df1 = pd.DataFrame({feature: temp.index,'Количество займов': temp.values})

    # Расчет доли target=1 в категории
    cat_perc = application_train[[feature, 'TARGET']].groupby([feature],as_index=False).mean()
    cat_perc.sort_values(by='TARGET', ascending=False, inplace=True)
    
    if(horizontal_layout):
        fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12,6))
    else:
        fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(12,14))
    sns.set_color_codes("pastel")
    s = sns.barplot(ax=ax1, x = feature, y="Количество займов",data=df1)
    if(label_rotation):
        s.set_xticklabels(s.get_xticklabels(),rotation=90)
    
    s = sns.barplot(ax=ax2, x = feature, y='TARGET', order=cat_perc[feature], data=cat_perc)
    if(label_rotation):
        s.set_xticklabels(s.get_xticklabels(),rotation=90)
    plt.ylabel('Доля проблемных', fontsize=10)
    plt.tick_params(axis='both', which='major', labelsize=10)

    plt.show();

Итак, рассмотрим основные признаки колиентов

Тип займа

plot_stats('NAME_CONTRACT_TYPE')

Интересно, что револьверные кредиты (вероятно, овердрафты или что-то вроде того) составляют меньше 10% от общего количества займов. В то же время процент невозврата среди них гораздо выше. Хороший повод пересмотреть методику работы с этими займами, а может быть и отказаться от них вовсе.

Пол клиента

plot_stats('CODE_GENDER')

Женщин-клиентов почти вдвое больше мужчин, при этом мужчины показывают гораздо более высокий риск.

Владение машиной и недвижимостью

plot_stats('FLAG_OWN_CAR')
plot_stats('FLAG_OWN_REALTY')

Клиентов с машиной вдвое меньше «безлошадных». Риск по ним практически одинаковый, клиенты с машиной платят чуть лучше.

По недвижимости обратная картина — клиентов без нее вдвое меньше. Риск по владельцам недвижимости также чуть меньше.

Семейный статус

plot_stats('NAME_FAMILY_STATUS',True, True)

В то время как большинство клиентов состоит в браке, ниаболее рискованы клиенты в гражданском браке и одинокие. Вдовцы показывают минимальный риск.

Количество детей

plot_stats('CNT_CHILDREN')

Большинство клиентов бездетны. При этом клиенты с 9 и 11 детьми показывают полный невозврат

application_train.CNT_CHILDREN.value_counts()

0 215371 1 61119 2 26749 3 3717 4 429 5 84 6 21 7 7 14 3 19 2 12 2 10 2 9 2 8 2 11 1 Name: CNT_CHILDREN, dtype: int64

Как показывает подсчет значений, эти данные статистически незначимы — всего по 2 клиента обеих категорий. Однако, все 4 вышли в дефолт, равно как и половина клиентов с 6 детьми.

Количество членов семьи

plot_stats('CNT_FAM_MEMBERS',True)

Ситауция аналогична — чем меньше ртов, тем выше возвратность.

Тип дохода

plot_stats('NAME_INCOME_TYPE',False,False)

Матери-одиночки и безработные, скорее всего, уже отсекаются на этапе подачи заявки — их слишком мало в выборке. Но стабильно показывают проблемы.

Вид деятельности

plot_stats('OCCUPATION_TYPE',True, False)

application_train.OCCUPATION_TYPE.value_counts()

Laborers 55186 Sales staff 32102 Core staff 27570 Managers 21371 Drivers 18603 High skill tech staff 11380 Accountants 9813 Medicine staff 8537 Security staff 6721 Cooking staff 5946 Cleaning staff 4653 Private service staff 2652 Low-skill Laborers 2093 Waiters/barmen staff 1348 Secretaries 1305 Realty agents 751 HR staff 563 IT staff 526 Name: OCCUPATION_TYPE, dtype: int64

Здесь вызывают интерес водители и сотрудники безопасности, которые довольно многочисленны и выходят на проблемы чаще других категорий.

Образование

plot_stats('NAME_EDUCATION_TYPE',True)

Чем выше образование, тем лучше возвратность, очевидно.

Тип организации — работодателя

plot_stats('ORGANIZATION_TYPE',True, False)

Наивысший процент невозврата наблюдается в Transport: type 3 (16%), Industry: type 13 (13.5%), Industry: type 8 (12.5%) и в Restaurant (до 12%).

Распределение суммы кредитования

Рассмотрим распределение сумм кредитов и влияние их на возвратность

plt.figure(figsize=(12,5))
plt.title("Распределение AMT_CREDIT")
ax = sns.distplot(app_train["AMT_CREDIT"])

plt.figure(figsize=(12,5))

# KDE займов, выплаченных вовремя
sns.kdeplot(app_train.loc[app_train['TARGET'] == 0, 'AMT_CREDIT'], label = 'target == 0')

# KDE проблемных займов
sns.kdeplot(app_train.loc[app_train['TARGET'] == 1, 'AMT_CREDIT'], label = 'target == 1')

# Обозначения
plt.xlabel('Сумма кредитования'); plt.ylabel('Плотность'); plt.title('Суммы кредитования');

Как показывает график плотности, крепные суммы возвращаются чаще

Распределение по плотности проживания

plt.figure(figsize=(12,5))
plt.title("Распределение REGION_POPULATION_RELATIVE")
ax = sns.distplot(app_train["REGION_POPULATION_RELATIVE"])

plt.figure(figsize=(12,5))

# KDE займов, выплаченных вовремя
sns.kdeplot(app_train.loc[app_train['TARGET'] == 0, 'REGION_POPULATION_RELATIVE'], label = 'target == 0')

# KDE проблемных займов
sns.kdeplot(app_train.loc[app_train['TARGET'] == 1, 'REGION_POPULATION_RELATIVE'], label = 'target == 1')

# Обозначения
plt.xlabel('Плотность'); plt.ylabel('Плотность населения'); plt.title('Плотность населения');

Клиенты из более населенных регионов склонны лучше выплачивать кредит.

Таким образом, мы получили представление об основных признаках датасета и их влиянии на результат. Конкретно с перечисленными в этой статье мы ничего делать не будем, но они могут оказаться очень важными в дальнейшей работе.

Feature Engineering — преобразование признаков

Соревнования на Kaggle выигрываются преобразованием признаков — побеждает тот, кто смог создать самые полезные признаки из данных. По меньшей мере для структурированных данных выигрышные модели — это сейчас в основном разные варианты градиентного бустинга. Чаще всего эффективнее потратить время на преобразование признаков, чем на настройку гиперпараметров или подбор моделей. Модель все-таки может обучиться только по тем данным, которые ей переданы. Убедиться, что данные релевантны задаче — главная ответственность дата саентиста.

Процесс преобразования признаков может включать создание новых из имеющихся данных, выбор наиболее важных из имеющихся и т.д. Опробуем на этот раз полиномиальные признаки.

Полиномиальные признаки

Полиномиальный метод конструирования признаков заключается в то, что мы просто создаем признаки, которые являются степенью имеющихся признаков и их произведениями. В некоторых случаях такие сконструированные признаки могут иметь более сильную корреляцию с целевой переменной, чем их «родители». Хотя такие методы часто используются в статистических моделях, в машинном обучении они встречаются значительно реже. Впрочем. ничего не мешает нам их попробовать, тем более что Scikit-Learn имеет класс специально для этих целей — PolynomialFeatures — который создает полиномиальные признаки и их произведения, нужно указать лишь сами исходные признаки и максимальную степень, в которую их нужно возводить. Используем самые мощные по силе воздействия на результат 4 признака и степень 3, чтобы не слишком сильно усложнять модель и избежать оверфиттинга (перетренированности модели — её излишней подстройки под обучающую выборку).

# создадим новый датафрейм для полиномиальных признаков
poly_features = app_train[['EXT_SOURCE_1', 'EXT_SOURCE_2', 'EXT_SOURCE_3', 'DAYS_BIRTH', 'TARGET']]
poly_features_test = app_test[['EXT_SOURCE_1', 'EXT_SOURCE_2', 'EXT_SOURCE_3', 'DAYS_BIRTH']]

# обработаем отуствующие данные
from sklearn.preprocessing import Imputer
imputer = Imputer(strategy = 'median')

poly_target = poly_features['TARGET']

poly_features = poly_features.drop('TARGET', axis=1)

poly_features = imputer.fit_transform(poly_features)
poly_features_test = imputer.transform(poly_features_test)

from sklearn.preprocessing import PolynomialFeatures
                                  
# Создадим полиномиальный объект степени 3
poly_transformer = PolynomialFeatures(degree = 3)

# Тренировка полиномиальных признаков
poly_transformer.fit(poly_features)

# Трансформация признаков
poly_features = poly_transformer.transform(poly_features)
poly_features_test = poly_transformer.transform(poly_features_test)
print('Формат полиномиальных признаков: ', poly_features.shape)

Формат полиномиальных признаков: (307511, 35) Присвоить признакам имена можно при помощи метода get_feature_names

poly_transformer.get_feature_names(input_features = ['EXT_SOURCE_1', 'EXT_SOURCE_2', 'EXT_SOURCE_3', 'DAYS_BIRTH'])[:15]

['1', 'EXT_SOURCE_1', 'EXT_SOURCE_2', 'EXT_SOURCE_3', 'DAYS_BIRTH', 'EXT_SOURCE_1^2', 'EXT_SOURCE_1 EXT_SOURCE_2', 'EXT_SOURCE_1 EXT_SOURCE_3', 'EXT_SOURCE_1 DAYS_BIRTH', 'EXT_SOURCE_2^2', 'EXT_SOURCE_2 EXT_SOURCE_3', 'EXT_SOURCE_2 DAYS_BIRTH', 'EXT_SOURCE_3^2', 'EXT_SOURCE_3 DAYS_BIRTH', 'DAYS_BIRTH^2']

Итого 35 полиномиальных и производных признаков. Проверим их корреляцию с таргетом.

# Датафрейм для новых фич 
poly_features = pd.DataFrame(poly_features, 
                             columns = poly_transformer.get_feature_names(['EXT_SOURCE_1', 'EXT_SOURCE_2', 
                                                                           'EXT_SOURCE_3', 'DAYS_BIRTH']))

# Добавим таргет
poly_features['TARGET'] = poly_target

# рассчитаем корреляцию
poly_corrs = poly_features.corr()['TARGET'].sort_values()

# Отобразим признаки с наивысшей корреляцией
print(poly_corrs.head(10))
print(poly_corrs.tail(5))

EXT_SOURCE_2 EXT_SOURCE_3 -0.193939 EXT_SOURCE_1 EXT_SOURCE_2 EXT_SOURCE_3 -0.189605 EXT_SOURCE_2 EXT_SOURCE_3 DAYS_BIRTH -0.181283 EXT_SOURCE_2^2 EXT_SOURCE_3 -0.176428 EXT_SOURCE_2 EXT_SOURCE_3^2 -0.172282 EXT_SOURCE_1 EXT_SOURCE_2 -0.166625 EXT_SOURCE_1 EXT_SOURCE_3 -0.164065 EXT_SOURCE_2 -0.160295 EXT_SOURCE_2 DAYS_BIRTH -0.156873 EXT_SOURCE_1 EXT_SOURCE_2^2 -0.156867 Name: TARGET, dtype: float64 DAYS_BIRTH -0.078239 DAYS_BIRTH^2 -0.076672 DAYS_BIRTH^3 -0.074273 TARGET 1.000000 1 NaN Name: TARGET, dtype: float64

Итак, некоторые признаки показывают более высокую корреляцию, чем исходные. Есть смысл попробовать обучение с ними и без них (как и многое другое в машинном обучении, это можно выяснить экспериментально). Для этого создадим копию датафреймов и добавим туда новые фичи.

# загрузим тестовые признаки в датафрейм
poly_features_test = pd.DataFrame(poly_features_test, 
                                  columns = poly_transformer.get_feature_names(['EXT_SOURCE_1', 'EXT_SOURCE_2', 
                                                                                'EXT_SOURCE_3', 'DAYS_BIRTH']))

# объединим тренировочные датафреймы
poly_features['SK_ID_CURR'] = app_train['SK_ID_CURR']
app_train_poly = app_train.merge(poly_features, on = 'SK_ID_CURR', how = 'left')

# объединим тестовые датафреймы
poly_features_test['SK_ID_CURR'] = app_test['SK_ID_CURR']
app_test_poly = app_test.merge(poly_features_test, on = 'SK_ID_CURR', how = 'left')

# Выровняем датафреймы
app_train_poly, app_test_poly = app_train_poly.align(app_test_poly, join = 'inner', axis = 1)

# Посмотрим формат
print('Тренировочная выборка с полиномиальными признаками: ', app_train_poly.shape)
print('Тестовая выборка с полиномиальными признаками: ', app_test_poly.shape)

Тренировочная выборка с полиномиальными признаками: (307511, 277) Тестовая выборка с полиномиальными признаками: (48744, 277)

Тренировка модели

Базовый уровень

В расчетах нужно отталкиваться от какого-то базового уровня модели, ниже которого упасть уже нельзя. В нашем случае это могло бы быть 0,5 для всех тестовых клиентов — это показывает, что мы совершенно не представляем, вернет кредит клиент или нет. В нашем случае предварительная работа уже проведена и можно использовать более сложные модели.

Логистическая регрессия

Для расчета логистической регрессии нам нужно взять таблицы с закодированными категориальными признаками, заполнить недостающие данные и нормализовать их (привести к значениям от 0 до 1). Все это выполняет следующий код:

from sklearn.preprocessing import MinMaxScaler, Imputer

# Уберем таргет из тренировочных данных
if 'TARGET' in app_train:
    train = app_train.drop(labels = ['TARGET'], axis=1)
else:
    train = app_train.copy()
features = list(train.columns)

# копируем тестовые данные
test = app_test.copy()

# заполним недостающее по медиане
imputer = Imputer(strategy = 'median')

# Нормализация
scaler = MinMaxScaler(feature_range = (0, 1))

# заполнение тренировочной выборки
imputer.fit(train)

# Трансофрмация тренировочной и тестовой выборок
train = imputer.transform(train)
test = imputer.transform(app_test)

# то же самое с нормализацией
scaler.fit(train)
train = scaler.transform(train)
test = scaler.transform(test)

print('Формат тренировочной выборки: ', train.shape)
print('Формат тестовой выборки: ', test.shape)

Формат тренировочной выборки: (307511, 242) Формат тестовой выборки: (48744, 242)

Используем логистическую регрессию из Scikit-Learn как первую модель. Возьмем дефольную модель с одной поправкой — понизим параметр регуляризации C во избежание оверфиттинга. Синтаксис обычный — создаем модель, тренируем ее и пресказываем вероятность при помощи predict_proba (нам же нужна вероятность, а не 0/1)

from sklearn.linear_model import LogisticRegression

# Создаем модель
log_reg = LogisticRegression(C = 0.0001)

# Тр
    
            © Habrahabr.ru