Сравнение алгоритмов градиентного бустинга или история знает только первых…

Всем привет ! Данная статья написана по итогам обучения на курсе Otus ML Basic и в ней я проведу сравнение алгоритмов градиентного бустинга. Почему бустинг, спросите вы ? Понятно, что нейронные сети интереснее, но не всегда их применение целесообразно и есть задачи для которых классические методы машинного обучения являются лучшим выбором. Бустинг является одним из наиболее эффективных классических алгоритмов и поскольку существуют различные его реализации, то мы проведем сравнение, чтобы понять, кто из них демонстрирует лучшие результаты. Познакомимся с участниками турнира, чьи реализации алгоритма градиентного бустинга будут участвовать в сравнении:

• Sklearn;

• XGBoost;

• LightGBM;

• Catboost;

Напомню, что бустинг реализует идею построения "сильной" модели на основе композиции базовых алгоритмов (как правило, деревьев решений), точность предсказания которых может быть лишь немногим выше случайного угадывания. Общий подход к реализации выглядит следующим образом:

• строим алгоритмы последовательно;

• каждый следующий строится на ошибках предыдущего;

• решение принимается методом взвешенного голосования;

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

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

from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score
from sklearn.preprocessing import LabelEncoder, StandardScaler

from sklearn.ensemble import GradientBoostingClassifier
from xgboost import XGBClassifier
from lightgbm import LGBMClassifier
from catboost import CatBoostClassifier

from sklearn.model_selection import GridSearchCV

import warnings
warnings.filterwarnings("ignore")

Загрузим датасет и посмотрим на основные параметры

df = pd.read_csv('data/WA_Fn-UseC_-Telco-Customer-Churn.csv', index_col=0)
df.info()

    Index: 7043 entries, 7590-VHVEG to 3186-AJIEK
    Data columns (total 20 columns):
     #   Column            Non-Null Count  Dtype  
    ---  ------            --------------  -----  
     0   gender            7043 non-null   object 
     1   SeniorCitizen     7043 non-null   int64  
     2   Partner           7043 non-null   object 
     3   Dependents        7043 non-null   object 
     4   tenure            7043 non-null   int64  
     5   PhoneService      7043 non-null   object 
     6   MultipleLines     7043 non-null   object 
     7   InternetService   7043 non-null   object 
     8   OnlineSecurity    7043 non-null   object 
     9   OnlineBackup      7043 non-null   object 
     10  DeviceProtection  7043 non-null   object 
     11  TechSupport       7043 non-null   object 
     12  StreamingTV       7043 non-null   object 
     13  StreamingMovies   7043 non-null   object 
     14  Contract          7043 non-null   object 
     15  PaperlessBilling  7043 non-null   object 
     16  PaymentMethod     7043 non-null   object 
     17  MonthlyCharges    7043 non-null   float64
     18  TotalCharges      7043 non-null   object 
     19  Churn             7043 non-null   object 
    dtypes: float64(1), int64(2), object(17)
    memory usage: 1.1+ MB

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

Предварительная обработка данных

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

df.isna().sum()

    gender              0
    SeniorCitizen       0
    Partner             0
    Dependents          0
    tenure              0
    PhoneService        0
    MultipleLines       0
    InternetService     0
    OnlineSecurity      0
    OnlineBackup        0
    DeviceProtection    0
    TechSupport         0
    StreamingTV         0
    StreamingMovies     0
    Contract            0
    PaperlessBilling    0
    PaymentMethod       0
    MonthlyCharges      0
    TotalCharges        0
    Churn               0
    dtype: int64

Пропусков в данных нет и это хорошо, но есть 17 категориальных признаков, которые необходимо привести к числовому виду:

    Index(['gender', 'Partner', 'Dependents', 'PhoneService', 'MultipleLines',
           'InternetService', 'OnlineSecurity', 'OnlineBackup', 'DeviceProtection',
           'TechSupport', 'StreamingTV', 'StreamingMovies', 'Contract',
           'PaperlessBilling', 'PaymentMethod', 'TotalCharges', 'Churn'],
          dtype='object')

Посмотрим на количество уникальных значений для каждого из атрибутов

df.nunique()

    gender                 2
    SeniorCitizen          2
    Partner                2
    Dependents             2
    tenure                73
    PhoneService           2
    MultipleLines          3
    InternetService        3
    OnlineSecurity         3
    OnlineBackup           3
    DeviceProtection       3
    TechSupport            3
    StreamingTV            3
    StreamingMovies        3
    Contract               3
    PaperlessBilling       2
    PaymentMethod          4
    MonthlyCharges      1585
    TotalCharges        6531
    Churn                  2
    dtype: int64

Начнем с замены бинарных категориальных признаков значениями 1/0

bin_cat_cols_list = []
for index, value in df.nunique().items():    
  if value == 2:        
    bin_cat_cols_list.append(index)        
    print(f"Index : {index}, Value : {value}") 

    Index : gender, Value : 2
    Index : SeniorCitizen, Value : 2
    Index : Partner, Value : 2
    Index : Dependents, Value : 2
    Index : PhoneService, Value : 2
    Index : PaperlessBilling, Value : 2
    Index : Churn, Value : 2

for col in bin_cat_cols_list:
    print(col, df[col].unique())

    gender ['Female' 'Male']
    SeniorCitizen [0 1]
    Partner ['Yes' 'No']
    Dependents ['No' 'Yes']
    PhoneService ['No' 'Yes']
    PaperlessBilling ['Yes' 'No']
    Churn ['No' 'Yes']

bin_cat_cols_list.remove('SeniorCitizen')

for col in bin_cat_cols_list:
    print(col, df[col].unique())

    gender ['Female' 'Male']
    Partner ['Yes' 'No']
    Dependents ['No' 'Yes']
    PhoneService ['No' 'Yes']
    PaperlessBilling ['Yes' 'No']
    Churn ['No' 'Yes']

g_dict = {'Female':0, 'Male':1}
df['gender'] = df['gender'].map(g_dict)

yn_dict = {'Yes':1, 'No':0}
for col in bin_cat_cols_list[1:]:
    df[col] = df[col].map(yn_dict)

Посмотрим, что у нас получилось по итогам преобразования бинарных атрибутов

df.info()

    Index: 7043 entries, 7590-VHVEG to 3186-AJIEK
    Data columns (total 20 columns):
     #   Column            Non-Null Count  Dtype  
    ---  ------            --------------  -----  
     0   gender            7043 non-null   int64  
     1   SeniorCitizen     7043 non-null   int64  
     2   Partner           7043 non-null   int64  
     3   Dependents        7043 non-null   int64  
     4   tenure            7043 non-null   int64  
     5   PhoneService      7043 non-null   int64  
     6   MultipleLines     7043 non-null   object 
     7   InternetService   7043 non-null   object 
     8   OnlineSecurity    7043 non-null   object 
     9   OnlineBackup      7043 non-null   object 
     10  DeviceProtection  7043 non-null   object 
     11  TechSupport       7043 non-null   object 
     12  StreamingTV       7043 non-null   object 
     13  StreamingMovies   7043 non-null   object 
     14  Contract          7043 non-null   object 
     15  PaperlessBilling  7043 non-null   int64  
     16  PaymentMethod     7043 non-null   object 
     17  MonthlyCharges    7043 non-null   float64
     18  TotalCharges      7043 non-null   object 
     19  Churn             7043 non-null   int64  
    dtypes: float64(1), int64(8), object(11)
    memory usage: 1.1+ MB

Поработаем с оставшимися 11 категориальными признаками и начнем с приведения TotalCharges к типу float

df['TotalCharges'] = pd.to_numeric(df['TotalCharges'], errors='coerce')
df.isna().sum()

    gender               0
    SeniorCitizen        0
    Partner              0
    Dependents           0
    tenure               0
    PhoneService         0
    MultipleLines        0
    InternetService      0
    OnlineSecurity       0
    OnlineBackup         0
    DeviceProtection     0
    TechSupport          0
    StreamingTV          0
    StreamingMovies      0
    Contract             0
    PaperlessBilling     0
    PaymentMethod        0
    MonthlyCharges       0
    TotalCharges        11
    Churn                0
    dtype: int64

Видим, что есть 11 пропущенных значений в Total Charges, записей немного, поэтому, просто удалим их из набора данных

df.dropna(inplace = True)

Оставшиеся категориальные признаки преобразуем с использованием LabelEncoder пакета sklearn

obj_cols = df.select_dtypes(include='object').columns

for col in obj_cols:
    print(col, df[col].unique())

    MultipleLines ['No phone service' 'No' 'Yes']
    InternetService ['DSL' 'Fiber optic' 'No']
    OnlineSecurity ['No' 'Yes' 'No internet service']
    OnlineBackup ['Yes' 'No' 'No internet service']
    DeviceProtection ['No' 'Yes' 'No internet service']
    TechSupport ['No' 'Yes' 'No internet service']
    StreamingTV ['No' 'Yes' 'No internet service']
    StreamingMovies ['No' 'Yes' 'No internet service']
    Contract ['Month-to-month' 'One year' 'Two year']
    PaymentMethod ['Electronic check' 'Mailed check' 'Bank transfer (automatic)'
     'Credit card (automatic)']

label_encoder = LabelEncoder()

for col in obj_cols:
    df[col] = label_encoder.fit_transform(df[col])

for col in obj_cols:
    print(col, df[col].unique())

    MultipleLines [1 0 2]
    InternetService [0 1 2]
    OnlineSecurity [0 2 1]
    OnlineBackup [2 0 1]
    DeviceProtection [0 2 1]
    TechSupport [0 2 1]
    StreamingTV [0 2 1]
    StreamingMovies [0 2 1]
    Contract [0 1 2]
    PaymentMethod [2 3 0 1]

Проверим, что у нас получилось после всех преобразований:

df.info()

    Index: 7032 entries, 7590-VHVEG to 3186-AJIEK
    Data columns (total 20 columns):
     #   Column            Non-Null Count  Dtype  
    ---  ------            --------------  -----  
     0   gender            7032 non-null   int64  
     1   SeniorCitizen     7032 non-null   int64  
     2   Partner           7032 non-null   int64  
     3   Dependents        7032 non-null   int64  
     4   tenure            7032 non-null   int64  
     5   PhoneService      7032 non-null   int64  
     6   MultipleLines     7032 non-null   int64  
     7   InternetService   7032 non-null   int64  
     8   OnlineSecurity    7032 non-null   int64  
     9   OnlineBackup      7032 non-null   int64  
     10  DeviceProtection  7032 non-null   int64  
     11  TechSupport       7032 non-null   int64  
     12  StreamingTV       7032 non-null   int64  
     13  StreamingMovies   7032 non-null   int64  
     14  Contract          7032 non-null   int64  
     15  PaperlessBilling  7032 non-null   int64  
     16  PaymentMethod     7032 non-null   int64  
     17  MonthlyCharges    7032 non-null   float64
     18  TotalCharges      7032 non-null   float64
     19  Churn             7032 non-null   int64  
    dtypes: float64(2), int64(18)
    memory usage: 1.4+ MB

Видим, что все признаки теперь числовые и наш датасет готов к дальнейшей работе. Но прежде чем погрузиться в создание моделей градиентного бустинга, проведем разведочный анализ данных, aka Exploratory Data Analysis, он же EDA

EDA

Посмотрим на корреляцию в данных

fig, ax = plt.subplots(figsize=(12, 6))
sns.heatmap(df.corr(), annot=True, linewidths=.5, fmt= '.1f',ax=ax);

Визуализируем парные зависимости выбранных признаков

sns.pairplot(data=df[['tenure','Contract','MonthlyCharges','TotalCharges','Churn']], hue='Churn')
plt.show;

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

Корреляция целевой переменной с другими признаками

Посмотрим на корреляцию оттока (Churn) с другими признаками

plt.figure(figsize=(8,6))
df.corr()['Churn'].sort_values(ascending = False).plot(kind='bar')

Числовые признаки и целевая переменная

Посмотрим на распределение некоторых числовых признаков в разрезе целевой переменной

fig = plt.subplots(nrows = 1,ncols = 3,figsize = (20,7))

plt.subplot(1,3,1)
ax = sns.kdeplot(df.MonthlyCharges[(df["Churn"] == 0)], color='#008080', fill= True, alpha=.7, linewidth=0)
ax = sns.kdeplot(df.MonthlyCharges[(df["Churn"] == 1)], color='#FF6347', fill= True, alpha=.7, linewidth=0)
ax.legend(["Not Churn","Churn"],loc='upper right')
ax.set_ylabel('Density')
ax.set_xlabel('Monthly Charges')
ax.set_title('Distribution of Monthly Charges by Churn')

plt.subplot(1,3,2)
ax = sns.kdeplot(df.TotalCharges[(df["Churn"] == 0)], color='#008080', fill= True, alpha=.7, linewidth=0)
ax = sns.kdeplot(df.TotalCharges[(df["Churn"] == 1)], color='#FF6347', fill= True, alpha=.7, linewidth=0)
ax.legend(["Not Churn","Churn"],loc='upper right')
ax.set_ylabel('Density')
ax.set_xlabel('Total Charges')
ax.set_title('Distribution of Total Charges by Churn')

plt.subplot(1,3,3)
ax = sns.kdeplot(df.tenure[(df["Churn"] == 0)], color='#008080', fill= True, alpha=.7, linewidth=0)
ax = sns.kdeplot(df.tenure[(df["Churn"] == 1)], color='#FF6347', fill= True, alpha=.7, linewidth=0)
ax.legend(["Not Churn","Churn"],loc='upper right')
ax.set_ylabel('Density')
ax.set_xlabel('Tenure')
ax.set_title('Distribution of Tenure by Churn')

plt.show();

Полученные диаграммы позволяют сделать несколько выводов:

• Диаграмма распределения ежемесячных платежей (Monthly Charges) показывает, что к оттоку склонны клиенты с большими суммами платежей, возможно, неожиданные счета за роуминг влияют на лояльность клиентов;

• Среди клиентов с большой общей суммой счетов (Total Charges) выше доля лояльных клиентов;

• Распределение по времени контракта (Tenure) демонстрирует лучшее разделение по целевой переменной - лояльные клиенты имеют давние контракты, в то время как новые клиенты наиболее склонны к оттоку;

И также посмотрим на распределение целевой переменной:

palette = ['#008080','#FF6347', '#E50000', '#D2691E']

l1 = list(df['Churn'].value_counts())
pie_values = [l1[0] / sum(l1) * 100, l1[1] / sum(l1) * 100]

fig = plt.subplots(nrows = 1,ncols = 2,figsize = (20,7))

plt.subplot(1,2,1)
plt.pie(pie_values,labels = ['Not-Churn Customers','Churn Customers'], 
        autopct = '%1.2f%%',
        explode = (0.1,0),
        colors = palette,
        wedgeprops = {'edgecolor': 'black','linewidth': 1, 'antialiased' : True})
plt.title('Churn and Not-Churn Customers %');

plt.subplot(1,2,2)
ax = sns.countplot(data = df, 
                   x='Churn',
                   palette = palette, 
                   edgecolor = 'black')
for i in ax.containers:
    ax.bar_label(i,)
ax.set_xticklabels(['Not-Churn Customers','Churn Customers'])
    
plt.title('Churn and Not-Churn Customers')
plt.show()

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

Разделение и масштабирование

Разделение набора данных на обучающую и тестовую выборки

Как обычно, перед обучением модели, нам необходимо разделить датасет на обучающую (train) и тестовую (test) выборки. Используем для этого функцию train_test_split пакета sklearn и не забудем про параметр stratify, учитывая несбалансированность набора данных:

X = df.drop('Churn', axis=1)
y = df['Churn']

X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, test_size = 0.3, random_state = 13)

Масштабирование признаков

Масштабирование числовых признаков набора данных является общим требованием для многих моделей машинного обучения и мы воспользуемся функцией StandardScaler модуля preprocessing пакета sklearn. StandardScaler центрирует значения признаков относительно нуля, вычитая среднее значение каждого признака, а затем масштабирует их, деля на стандартное отклонение:

где это среднее значение для признака из обучающего набора, - стандартное отклонение

numeric_columns = ['tenure', 'MonthlyCharges', 'TotalCharges']
std_scaler = StandardScaler()

X_train[numeric_columns] = std_scaler.fit_transform(X_train[numeric_columns])
X_test[numeric_columns]= std_scaler.transform(X_test[numeric_columns])

Теперь мы готовы строить наши модели.

Sklearn

Открывать турнир будет реализация градиентного бустинга от sklearn и для удобства определим функцию, возвращающую значения метрик, по которым мы будем оценивать и сравнивать модели. Поскольку у нас задача бинарной классификации, то будем использовать соответствующие метрики, в частности - accuracy, precision, recall, f1-score и ROC-AUC. Сравнивать, так уж сравнивать )

def quality(true_y, prediction_y):
    """
    Evaluates and returns the following metrics: Accuracy, Precision, Recall, F1-score, AUC
    """
    accuracy = round(accuracy_score(true_y, prediction_y), 3)
    precision = round(precision_score(true_y, prediction_y), 3)
    recall = round(recall_score(true_y, prediction_y), 3)
    f1 = round(f1_score(true_y, prediction_y), 3)
    auc = round(roc_auc_score(true_y, prediction_y), 3)
    print(f" Accuracy: {accuracy}")
    print(f"Precision: {precision}")
    print(f"   Recall: {recall}")
    print(f" F1-score: {f1}")
    print(f"      AUC: {auc}")
    return [accuracy, precision, recall, f1, auc]

Измеренные значения метрик будем складывать в словарь

results = {}

В качестве точки отсчета запустим классификатор без настройки, со значениями гиперпараметров по умолчанию

# first run with default parameters
sgb_clf = GradientBoostingClassifier(random_state=13)

sgb_clf.fit(X_train, y_train)
y_pred = sgb_clf.predict(X_test)

results['Sklearn'] = quality(y_test, y_pred)

     Accuracy: 0.807
    Precision: 0.678
       Recall: 0.522
     F1-score: 0.59
          AUC: 0.716

Посмотрим на первые полученные результаты

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

Теперь у нас есть baseline и пора улучшить результаты.

Sklearn: кривая валидации для количества деревьев

Первый гиперпараметр который мы попробуем настроить это n_estimators или количество выполняемых этапов бустинга - в нашем случае, числа деревьев решений (decision trees), используемых в качестве базового алгоритма. И для начала построим кривую валидации, отображающую зависимость результатов (по метрике ROC-AUC) от количества деревьев решений

n_trees = [1, 3, 5, 10, 50, 100, 200, 300, 400, 500]
quals_train = []
quals_test = []
for n in n_trees:
    clf = GradientBoostingClassifier(n_estimators=n, random_state=13)
    clf.fit(X_train, y_train)
    q_train = roc_auc_score(y_train, clf.predict(X_train))
    q_test = roc_auc_score(y_test, clf.predict(X_test))

    quals_train.append(q_train)
    quals_test.append(q_test)

plt.figure(figsize=(8, 5))
plt.plot(n_trees, quals_train, marker='.', label='train')
plt.plot(n_trees, quals_test, marker='.', label='test')
plt.xlabel('Number of trees')
plt.ylabel('AUC-ROC')
plt.title('Sklearn GB Validation Curve')
plt.legend()

plt.show();

Если отсортировать результаты в порядке убывания значения выбранной метрики

sorted(list(zip(quals_test, n_trees)), reverse=True)

    [(0.7162731634117349, 100),
     (0.7099715876725712, 200),
     (0.7046383786215936, 400),
     (0.7040699019205077, 300),
     (0.7035014252194217, 500),
     (0.7016757404293955, 50),
     (0.6225124828967916, 10),
     (0.5, 5),
     (0.5, 3),
     (0.5, 1)]

то увидим, что лучшие результаты достигаются на 100 деревьях

Sklearn: кривые валидации для скорости обучения (learning rate)

Посмотрим, как влияет гиперпараметр learning rate на качество алгоритма и склонность к переобучению. Для построения кривых валидации воспользуемся методом staged_predict, позволяющим получать результаты на каждом этапе бустинга, по мере добавления очередного дерева решений:

for learning_rate in [1, 0.5, 0.3, 0.2, 0.1]:

    gbm = GradientBoostingClassifier(n_estimators=150, learning_rate=learning_rate, random_state=13).fit(X_train, y_train)
    
    test_deviance = np.zeros((gbm.n_estimators,), dtype=np.float64)
    for i, y_pred in enumerate(gbm.staged_predict(X_test)):
        test_deviance[i] = roc_auc_score(y_test, y_pred)
    
    train_deviance = np.zeros((gbm.n_estimators,), dtype=np.float64)
    for i, y_pred in enumerate(gbm.staged_predict(X_train)):
        train_deviance[i] = roc_auc_score(y_train, y_pred)

    plt.figure()
    plt.plot(test_deviance, 'r', linewidth=2)
    plt.plot(train_deviance, 'g', linewidth=2)
    plt.legend(['test', 'train'])
    
    plt.title('GBM lr=%.1f, test roc-auc=%.3f, best_est=%d' % (learning_rate, test_deviance.max(), test_deviance.argmax()+1))
    plt.xlabel('Number of trees')
    plt.ylabel('Metric')

Видим, что максимальное значение метрики ROC-AUC достигается при learning rate (lr) равном 0.1 и количестве этапов бустинга (n_estimators) равном 79.

Запустим классификатор sklearn с максимизирующими значение ROC-AUC параметрами lr=0.1, n_estimators=79

sgb_clf = GradientBoostingClassifier(n_estimators=79, learning_rate=0.1, random_state=13)

sgb_clf.fit(X_train, y_train)
y_pred = sgb_clf.predict(X_test)

results['Sklearn-VC'] = quality(y_test, y_pred)

     Accuracy: 0.809
    Precision: 0.684
       Recall: 0.528
     F1-score: 0.596
          AUC: 0.72

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

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

Sklearn: настройка гиперпараметров по сетке с кроссвалидацией

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

# Define Gradient Boosting classifier with default parameters
clf = GradientBoostingClassifier(random_state=13)

# Estimate grid of the classifier hyperparameters
parameters = {'n_estimators':[10,50,80,150],
              'max_depth':[1,2,3,5],
              'learning_rate':[1,0.5,0.3,0.2,0.1]
             }

# Define GridSearch parameters
gs = GridSearchCV(clf,                 # Classifier object to optimize
                  parameters,          # Grid of the hyperparameters
                  scoring='roc_auc',   # Classification quality metric to optimize
                  cv=5                 # Number of folds in KFolds cross-validation
                 )

# Run Grid Search optimization
gs.fit(X_train, y_train)
gs.best_params_

    {'learning_rate': 0.2, 'max_depth': 1, 'n_estimators': 150}

pred_gs = gs.predict(X_test)

results['Sklearn-GS'] = quality(y_test, pred_gs)

     Accuracy: 0.808
    Precision: 0.681
       Recall: 0.522
     F1-score: 0.591
          AUC: 0.717

Итоговые результаты алгоритма sklearn:

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

После оптимизации параметров с использованием GridSearch метрика ROC-AUC чуть лучше чем при использовании параметров по умолчанию, но несколько хуже результатов, полученных с подобранными на кривых валидации параметрами.

XGBoost

Переходим к тестированию реализации алгоритма градиентного бустинга пакета xgboost и начнем с параметров по умолчанию.

xgb_clf = XGBClassifier(random_state=13)

xgb_clf.fit(X_train, y_train)
y_pred = xgb_clf.predict(X_test)

results['XGBoost'] = quality(y_test, y_pred)

     Accuracy: 0.786
    Precision: 0.621
       Recall: 0.503
     F1-score: 0.556
          AUC: 0.696

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

И у нас есть baseline для xgboost...

XGBoost: кривая валидации для количества деревьев

Построим кривую валидации для настройки параметра n_estimators

n_trees = [1, 3, 5, 10, 50, 100, 200, 300, 400, 500]
quals_train = []
quals_test = []
for n in n_trees:
    clf = XGBClassifier(n_estimators=n, random_state=13)
    clf.fit(X_train, y_train)
    q_train = roc_auc_score(y_train, clf.predict(X_train))
    q_test = roc_auc_score(y_test, clf.predict(X_test))

    quals_train.append(q_train)
    quals_test.append(q_test)

plt.figure(figsize=(8, 5))
plt.plot(n_trees, quals_train, marker='.', label='train')
plt.plot(n_trees, quals_test, marker='.', label='test')
plt.xlabel('Number of trees')
plt.ylabel('AUC-ROC')
plt.title('XGBoost Validation Curve')
plt.legend()

plt.show()

Отсортируем в порядке убывания значения выбранной метрики

sorted(list(zip(quals_test, n_trees)), reverse=True)

    [(0.7062379385699934, 10),
     (0.7055296442187416, 50),
     (0.6958172082730621, 100),
     (0.6847365156520968, 200),
     (0.6816628288735529, 300),
     (0.6809401499903911, 500),
     (0.6796489944061432, 400),
     (0.6699756843872593, 5),
     (0.6428717739810286, 3),
     (0.5, 1)]

Лучший результат получаем для n_estimators = 10

XGBoost: кривые валидации для скорости обучения

Посмотрим, как влияет параметр learning rate на качество алгоритма и склонность к переобучению

for learning_rate in [1, 0.5, 0.3, 0.2, 0.1]:

    xgb = XGBClassifier(n_estimators=150, learning_rate=learning_rate, random_state=13, verbose=-1).fit(X_train, y_train)

    test_deviance = np.zeros((xgb.n_estimators,), dtype=np.float64)
    for i in range(xgb.n_estimators):
        y_pred_test = xgb.predict(X_test, iteration_range=(0,i))
        test_deviance[i] = roc_auc_score(y_test, y_pred_test)

    train_deviance = np.zeros((xgb.n_estimators,), dtype=np.float64)
    for i in range(xgb.n_estimators):
        y_pred_train = xgb.predict(X_train, iteration_range=(0,i))
        train_deviance[i] = roc_auc_score(y_train, y_pred_train)

    plt.figure()
    plt.plot(test_deviance[1:], 'r', linewidth=2)
    plt.plot(train_deviance[1:], 'g', linewidth=2)
    plt.legend(['test', 'train'])
    
    plt.title('XGBoost lr=%.1f, test roc-auc=%.3f, best_est=%d' % (learning_rate, test_deviance.max(), test_deviance.argmax()))
    plt.xlabel('Number of trees')
    plt.ylabel('Metric')

Максимальное значение метрики ROC-AUC достигается при learning rate (lr) равном 0.5 и количестве этапов бустинга (n_estimators) равном 10.

Запустим модель с найденными оптимальными параметрами lr=0.5, n_estimators=10

xgb_clf = XGBClassifier(n_estimators=10, learning_rate=0.5, random_state=13)
xgb_clf.fit(X_train, y_train)
y_pred = xgb_clf.predict(X_test)

results['XGBoost-VC'] = quality(y_test, y_pred)

     Accuracy: 0.803
    Precision: 0.661
       Recall: 0.531
     F1-score: 0.589
          AUC: 0.716

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

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

XGBoost: настройка гиперпараметров по сетке с кроссвалидацией

Посмотрим, какие результаты нам удастся получить после поиска по сетке с использованием GridSearchCV

# Define Gradient Boosting classifier with default parameters
clf = XGBClassifier(random_state=13)

# Estimate grid of the classifier hyperparameters
parameters = {'n_estimators':[10,50,100],
              'max_depth':[1,2,3,5],
              'learning_rate':[1,0.5,0.3]
             }

# Define GridSearch parameters
gs = GridSearchCV(clf,                 # Classifier object to optimize
                  parameters,          # Grid of the hyperparameters
                  scoring='roc_auc',   # Classification quality metric to optimize
                  cv=5                 # Number of folds in KFolds cross-validation
                 )

# Run Grid Search optimization
gs.fit(X_train, y_train)
gs.best_params_

    {'learning_rate': 0.5, 'max_depth': 1, 'n_estimators': 50}

pred_gs = gs.predict(X_test)

results['XGBoost-GS'] = quality(y_test, pred_gs)

     Accuracy: 0.806
    Precision: 0.668
       Recall: 0.535
     F1-score: 0.594
          AUC: 0.719

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

Видим, что с использованием GridSearchCV результаты еще улучшились и по метрике ROC-AUC xgboost вышел на второе промежуточное место.

LightGBM

Третий участник - реализация алгоритма градиентного бустинга пакета LightGBM и, как обычно, первый запуск "из коробки", со значениями гиперпараметров по умолчанию

lgbm_clf = LGBMClassifier(verbose=-1, random_state=13)

lgbm_clf.fit(X_train, y_train)
y_pred = lgbm_clf.predict(X_test)

results['LightGBM'] = quality(y_test, y_pred)

     Accuracy: 0.794
    Precision: 0.643
       Recall: 0.504
     F1-score: 0.565
          AUC: 0.702

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

Ну что ж, неплохо для начала...

LightGBM: кривая валидации для количества деревьев

Построим кривую валидации для гиперпараметра n_estimators

n_trees = [1, 3, 5, 10, 50, 100, 200, 300, 400, 500]
quals_train = []
quals_test = []
for n in n_trees:
    clf = LGBMClassifier(n_estimators=n, verbose=-1, random_state=13)
    clf.fit(X_train, y_train)
    q_train = roc_auc_score(y_train, clf.predict(X_train))
    q_test = roc_auc_score(y_test, clf.predict(X_test))

    quals_train.append(q_train)
    quals_test.append(q_test)

plt.figure(figsize=(8, 5))
plt.plot(n_trees, quals_train, marker='.', label='train')
plt.plot(n_trees, quals_test, marker='.', label='test')
plt.xlabel('Number of trees')
plt.ylabel('AUC-ROC')
plt.title('LightGBM Validation Curve')
plt.legend()

plt.show();

Отсортируем по убыванию значения выбранной метрики

sorted(list(zip(quals_test, n_trees)), reverse=True)

    [(0.7075290941542413, 50),
     (0.7015503073111397, 100),
     (0.6927435214945183, 200),
     (0.6884557802227645, 300),
     (0.6865961479374308, 500),
     (0.682139819951691, 400),
     (0.6565635468343097, 10),
     (0.5301488281209544, 5),
     (0.5, 3),
     (0.5, 1)]

Лучший результат достигается для числа деревьев равного 50, но мы еще не настраивали learning rate...

LightGBM: кривые валидации для скорости обучения

Посмотрим, как влияет параметр learning_rate на качество алгоритма и склонность к переобучению

for learning_rate in [1, 0.5, 0.3, 0.2, 0.1]:

    lgb = LGBMClassifier(n_estimators=150, learning_rate=learning_rate, random_state=13, verbose=-1).fit(X_train, y_train)

    test_deviance = np.zeros((lgb.n_estimators,), dtype=np.float64)
    for i in range(lgb.n_estimators):
        y_pred_test = lgb.predict(X_test, num_iteration=i)
        test_deviance[i] = roc_auc_score(y_test, y_pred_test)

    train_deviance = np.zeros((lgb.n_estimators,), dtype=np.float64)
    for i in range(lgb.n_estimators):
        y_pred_train = lgb.predict(X_train, num_iteration=i)
        train_deviance[i] = roc_auc_score(y_train, y_pred_train)

    plt.figure()
    plt.plot(test_deviance[1:], 'r', linewidth=2)
    plt.plot(train_deviance[1:], 'g', linewidth=2)
    plt.legend(['test', 'train'])
    
    plt.title('LightGBM lr=%.1f, test roc-auc=%.3f, best_est=%d' % (learning_rate, test_deviance.max(), test_deviance.argmax()))
    plt.xlabel('Number of trees')
    plt.ylabel('Metric')

Максимальное значение метрики ROC-AUC достигается с параметром learning rate равным 0.3 и n_estimators равным 12.

Запустим модель с найденными оптимальными значениями гиперпараметров lr=0.3, n_estimators=12

lgbm = LGBMClassifier(n_estimators=12, learning_rate=0.3, verbose=-1, random_state=13)
lgbm.fit(X_train, y_train)
y_pred = lgbm.predict(X_test)

results['LightGBM-VC'] = quality(y_test, y_pred)

     Accuracy: 0.803
    Precision: 0.664
       Recall: 0.524
     F1-score: 0.586
          AUC: 0.714

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

Видим, что результаты улучшились.

LightGBM: настройка гиперпараметров с использованием GridSearchCV

Проведем настройку гиперпараметров поиском по сетке с использованием GridSearchCV

# Define Gradient Boosting classifier with default parameters
clf = LGBMClassifier(verbose=-1, random_state=13)

# Estimate grid of the classifier hyperparameters
parameters = {'n_estimators':[10,50,100,150],
              'max_depth':[1,2,3,5],
              'learning_rate':[1,0.5,0.3,0.2,0.1]
             }

# Define GridSearch parameters
gs = GridSearchCV(clf,                 # Classifier object to optimize
                  parameters,          # Grid of the hyperparameters
                  scoring='roc_auc',   # Classification quality metric to optimize
                  cv=5                 # Number of folds in KFolds cross-validation
                 )

# Run Grid Search optimization
gs.fit(X_train, y_train)
gs.best_params_

    {'learning_rate': 0.2, 'max_depth': 1, 'n_estimators': 150}

pred_gs = gs.predict(X_test)
results['LightGBM-GS'] = quality(y_test, pred_gs)

     Accuracy: 0.806
    Precision: 0.669
       Recall: 0.533
     F1-score: 0.593
          AUC: 0.719

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

После настройки на GridSearchCV результаты LightGBM по метрике ROC-AUC сравнялись с xgboost - плотная борьба...

CatBoost

На десерт протестируем реализацию алгоритма бустинга пакета catboost от Yandex и для начала оценим метрики "из коробки", то есть, со значениями гиперпараметров по умолчанию

catboost = CatBoostClassifier(logging_level='Silent', random_state=13)

catboost.fit(X_train, y_train)
pred = catboost.predict(X_test)

results['Catboost'] = quality(y_test, pred)

     Accuracy: 0.799
    Precision: 0.656
       Recall: 0.51
     F1-score: 0.574
          AUC: 0.706

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

И у нас есть первый результат catboost, который мы попробуем улучшить

Catboost: кривая валидации для количества деревьев

Как обычно, начнем настройку с кривой валидации для количества деревьев (n_estimators)

n_trees = [1, 3, 5, 10, 50, 100, 200, 300, 400, 500]
quals_train = []
quals_test = []
for n in n_trees:
    clf = CatBoostClassifier(iterations=n, logging_level='Silent', random_state=13)
    clf.fit(X_train, y_train)
    q_train = roc_auc_score(y_train, clf.predict(X_train))
    q_test = roc_auc_score(y_test, clf.predict(X_test))

    quals_train.append(q_train)
    quals_test.append(q_test)

plt.figure(figsize=(8, 5))
plt.plot(n_trees, quals_train, marker='.', label='train')
plt.plot(n_trees, quals_test, marker='.', label='test')
plt.xlabel('Number of trees')
plt.ylabel('AUC-ROC')
plt.title('Catboost Validation Curve')
plt.legend()

plt.show();

Отсортируем результаты в порядке убывания выбранной метрики

sorted(list(zip(quals_test, n_trees)), reverse=True)

    [(0.7219435458906844, 100),
     (0.7210522802935365, 10),
     (0.7113398443478572, 300),
     (0.7110170554517953, 50),
     (0.7071292041671413, 400),
     (0.7028414628953876, 500),
     (0.6999363628308299, 200),
     (0.69259449774393, 5),
     (0.680159932979589, 1),
     (0.6744320123729989, 3)]

Лучший результат достигается для 100 деревьев

Catboost: кривые валидации для скорости обучения

Продолжим настройку и посмотрим, как гиперпараметр learning rate влияет на качество алгоритма и склонность к переобучению

n_iterations = 150
for learning_rate in [1, 0.5, 0.3, 0.2, 0.1]:

    cbt = CatBoostClassifier(iterations=n_iterations, learning_rate=learning_rate, logging_level='Silent', random_state=13).fit(X_train, y_train)
    
    test_deviance = np.zeros((n_iterations,), dtype=np.float64)
    for i, y_pred in enumerate(cbt.staged_predict(X_test, prediction_type='Class', ntree_start=0, ntree_end=i)):
        test_deviance[i] = roc_auc_score(y_test, y_pred)
    
    train_deviance = np.zeros((n_iterations,), dtype=np.float64)
    for i, y_pred in enumerate(cbt.staged_predict(X_train, prediction_type='Class', ntree_start=0, ntree_end=i)):
        train_deviance[i] = roc_auc_score(y_train, y_pred)

    plt.figure()
    plt.plot(test_deviance, 'r', linewidth=2)
    plt.plot(train_deviance, 'g', linewidth=2)
    plt.legend(['test', 'train'])
    
    plt.title('Catboost lr=%.1f, test roc-auc=%.3f, best_est=%d' % (learning_rate, test_deviance.max(), test_deviance.argmax()+1))
    plt.xlabel('Number of trees')
    plt.ylabel('Metric')

Максимальное значение метрики ROC-AUC достигается при learning rate (lr) равном 0.1 и количестве этапов бустинга (n_estimators) равном 98.

Запустим модель с найденным оптимальным набором гиперпараметров lr=0.1, n_estimators=98

catboost = CatBoostClassifier(iterations=98, learning_rate=0.1, logging_level='Silent', random_state=13)
catboost.fit(X_train, y_train)
pred = catboost.predict(X_test)

results['Catboost-VC'] = quality(y_test, pred)

     Accuracy: 0.81
    Precision: 0.677
       Recall: 0.545
     F1-score: 0.604
          AUC: 0.726

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

И у нас смена лидера - catboost вырывается вперед !

Catboost: настройка гиперпараметров по сетке с кроссвалидацией

Проведем завершающую настройку поиском по сетке с использованием GridSearchCV

# Define Gradient Boosting classifier with default parameters
clf = CatBoostClassifier(logging_level='Silent', random_state=13)

# Estimate grid of the classifier hyperparameters
parameters = {'n_estimators':[10,50,100,150],
              'max_depth':[1,2,3,5],
              'learning_rate':[1,0.5,0.3,0.2,0.1]
             }

# Define GridSearch parameters
gs = GridSearchCV(clf,                 # Classifier object to optimize
                  parameters,          # Grid of the hyperparameters
                  scoring='roc_auc',   # Classification quality metric to optimize
                  cv=5                 # Number of folds in KFolds cross-validation
                 )

# Run Grid Search optimization
gs.fit(X_train, y_train)
gs.best_params_

    {'learning_rate': 0.1, 'max_depth': 3, 'n_estimators': 100}

pred_gs = gs.predict(X_test)
results['CatBoost-GS'] = quality(y_test, pred_gs)

     Accuracy: 0.808
    Precision: 0.676
       Recall: 0.533
     F1-score: 0.596
          AUC: 0.72

pd.DataFrame(results, index = ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']).T

И на GridSearchCV catboost показывает результаты чуть хуже...

Результаты и выводы

Итоговая таблица лидеров

Отсортируем итоговую турнирную таблицу по убыванию метрики ROC-AUC

И чемпионом становится catboost !

Визуализация результатов

Как известно, одна картинка стоит тысячи слов, поэтому визуализируем полученные результаты

plt.figure(figsize=(15, 6))
x = np.arange(5)

for key, value in results.items():
    plt.plot(x, results[key], marker='x', label=key);
    
plt.xticks(x, ['Accuracy', 'Precision', 'Recall', 'F1-score', 'AUC']);
plt.ylim(0.49, 0.82)
plt.legend(prop ={'size': 10});

Выводы

• Из коробки на первом месте реализация Sklearn, потом Catboost, затем LightGBM и XGBoost завершающий;

• После настройки параметров на первое место вышел Catboost, Sklearn переместился на второе, а третье поделили XGBoost и LightGBM с минимальным отставанием от второго места );

• Учитывая близость результатов можно сказать, что современные реализации алгоритма градиентного бустинга достаточно эффективны и все рассмотренные алгоритмы прекрасно справляются со своей задачей.