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

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

• Sklearn;

• XGBoost;

• LightGBM;

• Catboost;

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

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

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

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

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

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 с минимальным отставанием от второго места );

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