Какие из ваших функций переобучения?

"Секреты и уловки"

Какие из ваших функций переобучения?

Откройте для себя «ParShap»: расширенный метод определения того, какие столбцы делают вашу модель неэффективной на новых данных.

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

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

Кто виноват в перетренированности? Другими словами,

Какие функции (столбцы набора данных) мешают модели хорошо обобщать новые данные?

В этой статье — с помощью реального набора данных — мы увидим продвинутый метод ответа на этот вопрос.

На этот раз важность функции не спасет вас!

Если вы ответили на вопрос выше: «Я бы обратил внимание на важность функций», попробуйте усерднее.

Важность функций ничего не говорит о том, как функции будут работать с новыми данными.

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

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

«Здесь когда-то был немецкий госпиталь…»

Чтобы объяснить подход, я буду использовать набор данных, содержащий записи из реестра здоровья Германии с 1984 по 1988 год (набор данных доступен из библиотеки Pydataset по лицензии MIT).

Загрузка данных так же проста, как:

import pydataset
X = pydataset.data('rwm5yr')
y = (X['hospvis'] > 1).astype(int)
X = X.drop('hospvis', axis = 1)

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

Целевая переменная:

`hospvis`: находился ли пациент в больнице более 1 дня в течение соответствующего года.

Утилизируем 16 столбцов:

`id`: идентификатор пациента (1-7028);
`docvis`: количество посещений врача в течение года (0-121);
`год`: год (1984-1988);
`edlevel`: уровень образования (1-4);
`возраст`: возраст (25–64);
`outwork`: 1, если нет работы, 0 в противном случае;
`женщина`: 1, если женщина, 0 иначе;
`женат`: 1, если женат, 0 иначе;
`дети`: 1, если есть дети, 0 иначе;
`hhninc`: годовой доход семьи в марках (в марках);
`educ`: годы формального образования (7–18 лет);
`self`: 1, если работает не по найму, 0 иначе;
`edlevel1`: 1, если вы не закончили среднюю школу, иначе 0;
`edlevel2`: 1, если выпускник средней школы, 0 иначе;
`edlevel3`: 1, если университет/колледж, 0 иначе;
`edlevel4`: 1, если аспирантура, 0 иначе.

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

from sklearn.model_selection import train_test_split
from catboost import CatBoostClassifier
X_train, X_test, y_train, y_test = train_test_split(
  X, y, test_size = .2, stratify = y)
cat = CatBoostClassifier(silent = True).fit(X_train, y_train)

После того, как модель обучена, давайте посмотрим на важность функции:

import pandas as pd
fimpo = pd.Series(cat.feature_importances_, index = X_train.columns)

Неудивительно, что «docvis» — количество посещений врача — важен для прогнозирования того, провел ли пациент в больнице более 1 дня. Также несколько очевидны «возраст» и «hhninc» (доход). Но тот факт, что идентификатор пациента занимает второе место по важности, должен вызывать у нас подозрения, особенно потому, что мы рассматривали его как числовой признак.

Теперь посчитаем производительность модели (площадь под ROC) как на обучающей, так и на тестовой выборке.

from sklearn.metrics import roc_auc_score
roc_train = roc_auc_score(y_train, cat.predict_proba(X_train)[:, 1])
roc_test = roc_auc_score(y_test, cat.predict_proba(X_test)[:, 1])

Это огромная разница! Это свидетельствует о сильном переоснащении. Но какие функции «ответственны» за это?

Вы когда-нибудь слышали о значениях SHAP?

У нас есть много показателей для измерения того, как модель работает с некоторыми данными. Но как мы можем измерить, как функция работает с некоторыми данными?

Самый мощный инструмент для этого называется «Значения SHAP».

В общем, чтобы эффективно вычислять значения SHAP для любой прогностической модели, вы можете использовать выделенную библиотеку Python. Однако в этом примере мы воспользуемся нативным методом Catboost:

from catboost import Pool
shap_train = pd.DataFrame(
  data = cat.get_feature_importance(
    data = Pool(X_train), 
    type = 'ShapValues')[:, :-1], 
  index = X_train.index, 
  columns = X_train.columns
)
shap_test = pd.DataFrame(
  data = cat.get_feature_importance(
    data = Pool(X_test), 
    type = 'ShapValues')[:, :-1], 
  index = X_test.index, 
  columns = X_test.columns
)

Если вы посмотрите на shap_train и shap_test, вы заметите, что они имеют одинаковую форму соответствующих наборов данных.

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

Давайте посмотрим на пару примеров:

У пациента в строке 12071 было 0 посещений. Соответствующее значение SHAP (-0,753) говорит нам, что эта часть информации снижает на -0,753 вероятность (фактически логарифмическую вероятность) того, что он находится в больнице более 1 дня. Напротив, пациентка в ряду 18650 была у врача 4 раза, что увеличивает на 0,918 логарифмическую вероятность того, что она находится в больнице более 1 дня.

Встреча ПарШап

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

Например, если мы хотим рассчитать корреляцию между функцией «docvis» и целевой переменной в наблюдениях, содержащихся в тестовом наборе:

np.corrcoef(shap_test['docvis'], y_test)

Однако значения SHAP являются аддитивными, что означает, что окончательный прогноз представляет собой сумму значений SHAP всех функций. Поэтому имеет больше смысла, если мы удалим влияние других функций перед вычислением корреляции. Это точное определение частичной корреляции. Удобная реализация частичной корреляции содержится в библиотеке Python Pingouin:

import pingouin
pingouin.partial_corr(
  data = pd.concat([shap_test, y_test], axis = 1).astype(float), 
  x = 'docvis', 
  y = y_test.name,
  x_covar = [feature for feature in shap_test.columns if feature != 'docvis'] 
)

Этот фрагмент кода означает «вычислить корреляцию между значениями SHAP функции `docvis` и целевой переменной в наблюдениях тестового набора после удаления влияния всех других функций».

Поскольку всем нужно имя, я назову эту формулу "ParShap" (от "Частичная корреляция значений Shap").

Мы можем повторить эту процедуру для каждой функции, как в наборе поездов, так и в наборе тестов:

parshap_train = partial_correlation(shap_train, y_train)
parshap_test = partial_correlation(shap_test, y_test)

Примечание: вы можете найти определение функции partial_correlation в конце этой статьи.

Теперь построим parshap_train по оси x и parshap_test по оси y.

plt.scatter(parshap_train, parshap_test)

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

Визуально сразу видно, какие функции неэффективны: я обвел их синим цветом.

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

parshap_diff = parshap_test — parshap_train

Как мы должны интерпретировать этот вывод? Основываясь на том, что мы сказали выше, чем отрицательнее этот показатель, тем больше переобучения дает функция.

«Я тебе не верю»

Хорошо: не надо! Можете ли вы придумать способ проверить, верна ли интуиция, стоящая за этой статьей?

По логике, если мы удалим «функции переобучения» из набора данных, мы сможем уменьшить переоснащение (а именно, расстояние между roc_train и roc_test).

Итак, давайте попробуем отбросить одну функцию за раз и посмотреть, как изменится область под ROC.

Слева слева мы удаляем по одной функции за раз, отсортированные по важности функции. Итак, сначала удаляется наименее важный (`edlevel4`), затем удаляются два наименее важных (`edlevel4` и `edlevel1`) и так далее.

Справа справа мы следуем тому же процессу, но порядок удаления определяется ParShap. Итак, сначала удаляется самый отрицательный ParShap ("id"), затем удаляются два самых отрицательных ParShap ("id" и "year") и так далее.

Как мы и надеялись, удаление признаков с наиболее отрицательным значением ParShap привело к сильному уменьшению переобучения. На самом деле roc_train все ближе и ближе к roc_test.

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

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

Полный код, использованный в этой статье (вы можете получить немного другие результаты из-за случайного начального числа):

# Import libraries
import pandas as pd
import pydataset
from sklearn.model_selection import train_test_split
from catboost import CatBoostClassifier, Pool
from sklearn.metrics import roc_auc_score
from pingouin import partial_corr
import matplotlib.pyplot as plt
# Print documentation and read data
print('################# Print docs')
pydataset.data('rwm5yr', show_doc = True)
X = pydataset.data('rwm5yr')
y = (X['hospvis'] > 1).astype(int)
X = X.drop('hospvis', axis = 1)
# Split data in train + test
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = .2, stratify = y)
# Fit model
cat = CatBoostClassifier(silent = True).fit(X_train, y_train)
# Show feature importance
fimpo = pd.Series(cat.feature_importances_, index = X_train.columns)
fig, ax = plt.subplots()
fimpo.sort_values().plot.barh(ax = ax)
fig.savefig('fimpo.png', dpi = 200, bbox_inches="tight")
fig.show()
# Compute metrics
roc_train = roc_auc_score(y_train, cat.predict_proba(X_train)[:, 1])
roc_test = roc_auc_score(y_test, cat.predict_proba(X_test)[:, 1])
print('\n################# Print roc')
print('roc_auc train: {:.2f}'.format(roc_train))
print('roc_auc  test: {:.2f}'.format(roc_test))
# Compute SHAP values  
shap_train = pd.DataFrame(
  data = cat.get_feature_importance(data = Pool(X_train), type = 'ShapValues')[:, :-1],
  index = X_train.index, 
  columns = X_train.columns
)
shap_test = pd.DataFrame(
  data = cat.get_feature_importance(data = Pool(X_test), type = 'ShapValues')[:, :-1],
  index = X_test.index, 
  columns = X_test.columns
)
print('\n################# Print df shapes')
print(f'X_train.shape:    {X_train.shape}')
print(f'X_test.shape:     {X_test.shape}\n')
print(f'shap_train.shape: {shap_train.shape}')
print(f'shap_test.shape:  {shap_test.shape}')
print('\n################# Print data and SHAP')
print('Original data:')
display(X_test.head(3))
print('\nCorresponding SHAP values:')
display(shap_test.head(3).round(3))
# Define function for partial correlation
def partial_correlation(X, y):
  out = pd.Series(index = X.columns, dtype = float)
  for feature_name in X.columns:
    out[feature_name] = partial_corr(
      data = pd.concat([X, y], axis = 1).astype(float), 
      x = feature_name, 
      y = y.name,
      x_covar = [f for f in X.columns if f != feature_name] 
    ).loc['pearson', 'r']
  return out
# Compute ParShap
parshap_train = partial_correlation(shap_train, y_train)
parshap_test = partial_correlation(shap_test, y_test)
parshap_diff = pd.Series(parshap_test - parshap_train, name = 'parshap_diff')
print('\n################# Print parshap_diff')
print(parshap_diff.sort_values())
                         
# Plot parshap
plotmin, plotmax = min(parshap_train.min(), parshap_test.min()), max(parshap_train.max(), parshap_test.max())
plotbuffer = .05 * (plotmax - plotmin)
fig, ax = plt.subplots()
if plotmin < 0:
    ax.vlines(0, plotmin - plotbuffer, plotmax + plotbuffer, color = 'darkgrey', zorder = 0)
    ax.hlines(0, plotmin - plotbuffer, plotmax + plotbuffer, color = 'darkgrey', zorder = 0)
ax.plot(
    [plotmin - plotbuffer, plotmax + plotbuffer], [plotmin - plotbuffer, plotmax + plotbuffer], 
    color = 'darkgrey', zorder = 0
)
sc = ax.scatter(
    parshap_train, parshap_test, 
    edgecolor = 'grey', c = fimpo, s = 50, cmap = plt.cm.get_cmap('Reds'), vmin = 0, vmax = fimpo.max())
ax.set(title = 'Partial correlation bw SHAP and target...', xlabel = '... on Train data', ylabel = '... on Test data')
cbar = fig.colorbar(sc)
cbar.set_ticks([])
for txt in parshap_train.index:
    ax.annotate(txt, (parshap_train[txt], parshap_test[txt] + plotbuffer / 2), ha = 'center', va = 'bottom')
fig.savefig('parshap.png', dpi = 300, bbox_inches="tight")
fig.show()
# Feature selection
n_drop_max = 5
iterations = 4
features = {'parshap': parshap_diff, 'fimpo': fimpo}
features_dropped = {}
roc_auc_scores = {
  'fimpo': {'train': pd.DataFrame(), 'test': pd.DataFrame()},
  'parshap': {'train': pd.DataFrame(), 'test': pd.DataFrame()}
}
for type_ in ['parshap', 'fimpo']:
  for n_drop in range(n_drop_max + 1):
    features_drop = features[type_].sort_values().head(n_drop).index.to_list()
    features_dropped[type_] = features_drop
    X_drop = X.drop(features_drop, axis = 1)
    for i in range(iterations):
      X_train, X_test, y_train, y_test = train_test_split(X_drop, y, test_size = .2, stratify = y)
      cat = CatBoostClassifier(silent = True).fit(X_train, y_train)
      roc_auc_scores[type_]['train'].loc[n_drop, i] = roc_auc_score(y_train, cat.predict_proba(X_train)[:, 1])
      roc_auc_scores[type_]['test'].loc[n_drop, i] = roc_auc_score(y_test, cat.predict_proba(X_test)[:, 1])
        
# Plot feature selection
fig, axs = plt.subplots(1, 2, sharey = True, figsize = (8, 3))
plt.subplots_adjust(wspace = .1)
axs[0].plot(roc_auc_scores['fimpo']['train'].index, roc_auc_scores['fimpo']['train'].mean(axis = 1), lw = 3, label = 'Train')
axs[0].plot(roc_auc_scores['fimpo']['test'].index, roc_auc_scores['fimpo']['test'].mean(axis = 1), lw = 3, label = 'Test')
axs[0].set_xticks(roc_auc_scores['fimpo']['train'].index)
axs[0].set_xticklabels([''] + features_dropped['fimpo'], rotation = 90)
axs[0].set_title('Feature Importance')
axs[0].set_xlabel('Feature dropped')
axs[0].grid()
axs[0].legend(loc = 'center left')
axs[0].set(ylabel = 'ROC-AUC score')
axs[1].plot(roc_auc_scores['parshap']['train'].index, roc_auc_scores['parshap']['train'].mean(axis = 1), lw = 3, label = 'Train')
axs[1].plot(roc_auc_scores['parshap']['test'].index, roc_auc_scores['parshap']['test'].mean(axis = 1), lw = 3, label = 'Test')
axs[1].set_xticks(roc_auc_scores['parshap']['train'].index)
axs[1].set_xticklabels([''] + features_dropped['parshap'], rotation = 90)
axs[1].set_title('ParShap')
axs[1].set_xlabel('Feature dropped')
axs[1].grid()
axs[1].legend(loc = 'center left')
fig.savefig('feature_selection.png', dpi = 300, bbox_inches="tight")
fig.show()

Спасибо за чтение! Надеюсь, вам понравилась эта статья. Если хотите, добавьте меня в Linkedin!

Какие из ваших функций переобучения?

"Секреты и уловки"