Как правильно провести анализ классификации с помощью sklearn, когда ваш набор данных несбалансирован, и улучшить его результаты.
Представим, что у вас есть набор данных с десятком функций, и вам нужно классифицировать каждое наблюдение. Это может быть проблема с двумя классами (ваш результат - 1 или 0; истина или ложь) или проблема с несколькими классами (возможно более двух альтернатив). Однако в этом случае есть изюминка. Данные несбалансированы. Подумайте о пациентах, которые могут быть или не болеть раком (большинство, вероятно, не болеет), или о решении продлить кредитную линию (большинство клиентов банка получают продление). Ваш алгоритм машинного обучения будет «переэкспонирован» для одного класса и «недоэкспонирован» для другого. В Интернете есть множество статей об этой проблеме с разными подходами. Здесь я объединю несколько из них, чтобы получить надежное решение.
Мы начнем с более простых линейных моделей, затем добавим пару настроек, перейдем к усиленным деревьям и, наконец, к нейронным сетям.
Вы можете скачать набор данных и рабочие тетради здесь.
Мы будем следовать этим шагам:
- Загрузите данные и некоторые из наиболее необходимых зависимостей
- Сделайте минимальную предварительную обработку
- Постройте пару простых моделей, которые послужат основой для дальнейших попыток.
- Потратьте время на разработку функций
- Применять различные методы, которые помогут нам справиться с несбалансированными данными
Загрузка и понимание
Для этого упражнения мы воспользуемся набором данных Framingham Heart Study от Kaggle. Он представляет собой задачу бинарной классификации, в которой нам нужно предсказать значение переменной TenYearCHD (ноль или единица), которая показывает, разовьется ли у пациента болезнь сердца.
import pandas as pd import numpy as np import matplotlib.pyplot as plt import scipy.stats as st import seaborn as sns import pandas_profiling %matplotlib inline df = pd.read_csv(r'path to dataset')
Сделаем предварительное исследование данных немного удобнее. Недавно я наткнулся на статью под названием 10 простых приемов для ускорения анализа данных в Python, написанную Парул Пандей, и установил пакет Профилирование.
С помощью всего одной команды он делает действительно полезные вещи:
# looking at stats pandas_profiling.ProfileReport(df)
Как видите, данные довольно чистые: у нас есть пропущенные значения кое-где, но мы скоро с ними справимся.
Корреляция между функциями в приведенной ниже таблице тоже не говорит о многом:
Два красных прямоугольника вокруг курения и диастолического / систолического артериального давления достаточно очевидны: обычно курение и высокое артериальное давление связаны
Названия большинства функций дают хорошее представление о том, что означает каждая переменная, но некоторые из них менее очевидны:
«Образование» может быть:
1 - для какой-то средней школы
2 - для средней школы GED
3 - какой-нибудь колледж или профессиональное училище
4 - колледж
BPMeds - это двоичная переменная, в которой
0 - означает, что пациент не принимает лекарства от кровяного давления.
1 - означает обратное
Давайте проверим целевую переменную вручную:
df['TenYearCHD'].value_counts(normalize = True)
А если вы визуальный человек:
sns.countplot(x='TenYearCHD',data=df)
Это не сильно несбалансированный набор, но он может исказить как линейные, так и нелинейные алгоритмы.
Предварительная обработка данных
Здесь я могу подумать о нескольких неотложных проблемах:
- отсутствующие значения и NaN
- Колонка «Образование» может рассматриваться или не рассматриваться как порядковая.
- Нормализация или стандартизация данных (обычно это относится к этой части, но я сделаю это позже, когда мы начнем строить модель, требующую нормализованных функций)
В столбце «BPMeds» отсутствуют некоторые значения, но ~ 96% столбцов равны нулю (лекарства от артериального давления не принимаются). Итак, было бы справедливо заполнить NaN нулями.
«Глюкоза» и «totChol» являются непрерывными переменными, и мы будем использовать среднее значение столбца для заполнения пустых ячеек.
«ИМТ» и «частота сердечных сокращений», кажется, тоже хорошо работают со средним значением.
«CigsPerDay» требует двухэтапного подхода. В нем 29 пропущенных значений. Вы можете подумать, что использование среднего - это нормально. Однако cigsPerDay связан с другой двоичной функцией currentSmoker. Итак, вы можете иметь NaN для некурящих и назначить этому человеку некоторое среднее количество сигарет. Мы хотим избежать этого сценария. Посмотрим, как это можно сделать.
df['cigsPerDay'].value_counts(normalize = True).plot(kind="bar")
Большинство людей - некурящие, мы не хотим давать им сигареты.
df['cigsPerDay'][df['currentSmoker']==0].isna().sum()
Эта команда возвращает ноль, поэтому кажется, что у нас нет NaN для некурящих. Но использовать .mean () по-прежнему плохо, так как он будет смещен в сторону нуля. Мы хотим применить команду fillna () только к курильщикам:
# creating a boolean array of smokers smoke = (df['currentSmoker']==1) # applying mean to NaNs in cigsPerDay but using a set of smokers only df.loc[smoke,'cigsPerDay'] = df.loc[smoke,'cigsPerDay'].fillna(df.loc[smoke,'cigsPerDay'].mean())
В результате среднее количество выкуриваемых сигарет сейчас чуть больше 18:
df['cigsPerDay'][df['currentSmoker']==1].mean()
Если бы мы решили пойти более легким путем, не исключая некурящих, это значение было бы около девяти.
Вы можете быть уверены, что некурящим не присвоены никакие значения:
df['cigsPerDay'][df['currentSmoker']==0].mean()
Образование - интересный вопрос. Хотя это категориальная переменная сама по себе, она уже закодирована для нас, и более высокий уровень образования соответствует большему числу. Потенциальная проблема здесь заключается в том, что интервал между колледжем и некоторым колледжем или профессиональным училищем, вероятно, не такой, как между средней школой или GED и некоторым колледжем или профессиональным училищем. Однако числа говорят об обратном: расстояние между четырьмя и тремя точно такое же, как между двумя и тремя. Вы можете прочитать об этой проблеме, набрав в Google Порядковые значения в логистической регрессии или просто просмотрев эту статью. Предположим пока, что расстояния одинаковы. Тем не менее, мы просто избавимся от NaN в образовании, чтобы регрессия могла работать с столбцом.
df['education'].value_counts(normalize = True).plot(kind="bar")
Кажется справедливым заменить здесь NaN на «1».
Внесение сразу всех необходимых изменений:
# Filling out missing values df['BPMeds'].fillna(0, inplace = True) df['glucose'].fillna(df.glucose.mean(), inplace = True) df['totChol'].fillna(df.totChol.mean(), inplace = True) df['education'].fillna(1, inplace = True) df['BMI'].fillna(df.BMI.mean(), inplace = True) df['heartRate'].fillna(df.heartRate.mean(), inplace = True)
Проверим, прошло ли это:
df.isna().sum()
И мы в порядке:
Базовые модели
Мы начнем со случайного леса, а также рассмотрим показатели важности наших функций.
Не забудьте отделить функции от целевой переменной:
features = df.iloc[:,:-1] result = df.iloc[:,-1]
Сама модель ниже.
Сначала зависимости:
import sklearn from sklearn.model_selection import train_test_split from sklearn.linear_model import LogisticRegression from sklearn.metrics import classification_report from sklearn.metrics import confusion_matrix from sklearn.metrics import accuracy_score from sklearn.metrics import roc_auc_score from sklearn.preprocessing import StandardScaler
Разделение на тренировку / тест:
X_train, X_test, y_train, y_test = train_test_split(features, result, test_size = 0.2, random_state = 14)
Установка:
rf = RandomForestClassifier() rf.fit(X_train, y_train) # Making predictions on unseen data predictions_rf = rf.predict(X_test)
Важность функций (иногда ради скорости вы хотите построить модель, используя меньшее количество функций):
# what features are the most important? plt.plot(rf.feature_importances_) plt.xticks(np.arange(X_train.shape[1]), X_train.columns.tolist(), rotation=90)
Если вы хотите, чтобы это было списком:
# View a list of the features and their importance scores list(zip(features, rf.feature_importances_))
Оценка модели:
print(classification_report(y_test, predictions_rf)) print(confusion_matrix(y_test, predictions_rf)) # Under ROC curve prob_rf = rf.predict_proba(X_test) prob_rf = [p[1] for p in prob_rf] print(roc_auc_score(y_test, prob_rf))
Здесь два вывода:
- Показатель точности 84,90% едва ли превосходит случайное предположение. Помните, что 84,81% данных помечены как ноль. Таким образом, мы превзошли прогноз на 0,09%.
- Кроме того, матрица неточностей показывает большое количество ложноотрицательных результатов:
Мы хотим избежать, потому что модель неверно утверждает, что у пациента нет проблемы, в то время как на самом деле она есть! Мы скоро решим эту проблему, а пока давайте перестроим модель, используя только важные функции.
from sklearn.feature_selection import SelectFromModel
from sklearn.ensemble import RandomForestClassifier
# Create a selector object that will use the random forest classifier to identify
# features that have an importance of more than 0.12
sfm = SelectFromModel(clf, threshold=0.12)
# Train the selector
sfm.fit(X_train_std, y_train)
feat_labels = list(features.columns.values) # creating a list with features' names
for feature_list_index in sfm.get_support(indices=True):
print(feat_labels[feature_list_index])
importances = clf.feature_importances_
std = np.std([tree.feature_importances_ for tree in clf.estimators_],
axis=0)
indices = np.argsort(importances)[::-1]
print("Feature ranking:")
for f in range(X_train_std.shape[1]):
print("%d. feature %d (%f)" % (f + 1, indices[f], importances[indices[f]]))
# Plot the feature importances of the forest
plt.figure()
plt.title("Feature importances")
plt.bar(range(X_train_std.shape[1]), importances[indices],
color="r", yerr=std[indices], align="center")
plt.xticks(range(X_train_std.shape[1]), indices)
plt.xlim([-1, X_train_std.shape[1]])
plt.show()
# with only important features. Can check X_important_train.shape[1]
X_important_train = sfm.transform(X_train_std)
X_important_test = sfm.transform(X_test_std)
clf_important = RandomForestClassifier(n_estimators=10000, random_state=0, n_jobs=-1)
clf_important.fit(X_important_train, y_train)
predictions_y_4 = clf_important.predict(X_important_test)
print(classification_report(y_test, predictions_y_4))
print(confusion_matrix(y_test, predictions_y_4))
accuracy_score(y_test, predictions_y_4)
# Under ROC curve
prob_y_4 = clf_important.predict_proba(X_important_test)
prob_y_4 = [p[1] for p in prob_y_4]
print(roc_auc_score(y_test, prob_y_4))
Я выбрал объекты с важностью более 0,12 и перестроил случайный лес, используя только эти столбцы. Результаты были примерно одинаковыми, но нам удалось сэкономить некоторую вычислительную мощность.
Логистическая регрессия
Для регрессий необходимы стандартизованные функции:
scaler = StandardScaler() X_train_std = scaler.fit_transform(X_train) X_test_std = scaler.fit_transform(X_test)
И модель:
logmodel = LogisticRegression(solver='liblinear') logmodel.fit(X_train_std, y_train) predictions_y_2 = logmodel.predict(X_test_std)
И оценка результатов:
print(classification_report(y_test, predictions_y_2)) print(confusion_matrix(y_test, predictions_y_2)) # Under ROC curve prob_y_2 = logmodel.predict_proba(X_test_std) prob_y_2 = [p[1] for p in prob_y_2] print(roc_auc_score(y_test, prob_y_2))
Результаты очень похожи на результаты случайного леса:
Немного лучше справился с ложными негативами. Но этого все еще недостаточно.
LogisticRegression имеет параметр class_weight, который должен помочь точнее. Давайте посмотрим:
logmodel = LogisticRegression(solver='liblinear', class_weight='balanced') logmodel.fit(X_train_std, y_train) predictions_y_3 = logmodel.predict(X_test_std) print(classification_report(y_test, predictions_y_3)) print(confusion_matrix(y_test, predictions_y_3)) accuracy_score(y_test, predictions_y_3) # Under ROC curve prob_y_3 = logmodel.predict_proba(X_test_std) prob_y_3 = [p[1] for p in prob_y_3] print(roc_auc_score(y_test, prob_y_3))
Это бесполезно: мы уменьшаем количество ошибок типа II за счет некоторой потери общей точности. Похоже, нам нужно перейти к чему-то другому.
Мы также можем реализовать class_weight вручную, передав разбивку из нуля и единицы значений, которые вы можете быстро получить с помощью value_counts:
df['TenYearCHD'].value_counts(normalize = True)
weights = {0 : '0.848113', 1 : '0.151887'}
logmodel_auto = LogisticRegression(class_weight = weights, solver = 'liblinear')
logmodel_auto.fit(X_train_std, y_train)
predictions_std_auto = logmodel_auto.predict(X_test_std)
print(classification_report(y_test, predictions_std_auto))
print(confusion_matrix(y_test, predictions_std_auto))
accuracy_score(y_test, predictions_std_auto)
# Under ROC curve
prob_y_4 = logmodel.predict_proba(X_test_std)
prob_y_4 = [p[1] for p in prob_y_4]
print(roc_auc_score(y_test, prob_y_4))
Хорошая работа над ошибками типа I, на самом деле:
Вы можете попробовать разные комбинации. Методом проб и ошибок (или с помощью GridSearch) вы можете найти тот, который соответствует вашим целям, например:
weights = {0 : '0.09042', 1 : '0.90958'}
Если вы решите использовать GridSearch для поиска подходящих весов, это можно реализовать следующим образом:
from sklearn.model_selection import GridSearchCV
weights = np.linspace(0.03, 0.97, 55)
scaler = StandardScaler()
features_std = scaler.fit_transform(features)
gsc = GridSearchCV(
estimator=LogisticRegression(solver='liblinear'),
param_grid={
'class_weight': [{0: x, 1: 1.0-x} for x in weights]
},
scoring='roc_auc',
cv=3
)
grid_result = gsc.fit(features_std, result)
В приведенном выше коде мы тестируем различные комбинации нулей и единиц для переменной TenYearCHD в диапазоне от 3-97% до 97-3%. GridSearchCV может оценивать различные оценщики, и они выбираются и контролируются в поле скоринг. Обычно вы используете что-то вроде точности для классификации и R-квадрат для регрессии, но в нашем случае точность не очень информативна, поэтому я решил выбрать площадь под кривой AUC. Все доступные параметры скоринга описаны здесь.
Затем мы можем просто распечатать лучшие параметры и передать их модели:
print("Best parameters : %s" % grid_result.best_params_)
# passing weights found above
rf_w = RandomForestClassifier(class_weight = {0:0.882962962962963, 1:0.11703703703703705})
rf_w.fit(X_train, y_train)
print(classification_report(y_test, predictions_rf_w))
print(confusion_matrix(y_test, predictions_rf_w))
accuracy_score(y_test, predictions_rf_w)
Или к логистической регрессии:
weights = {0 : '0.882962962962963', 1 : '0.11703703703703705'}
logmodel_auto_gridsearch = LogisticRegression(class_weight = weights, solver = 'liblinear')
logmodel_auto_gridsearch.fit(X_train_std, y_train)
predictions_std_auto_gridsearch = logmodel_auto_gridsearch.predict(X_test_std)
print(classification_report(y_test, predictions_std_auto_gridsearch))
print(confusion_matrix(y_test, predictions_std_auto_gridsearch))
accuracy_score(y_test, predictions_std_auto_gridsearch)
# Under ROC curve
prob_y_3_gridsearch = logmodel_auto_gridsearch.predict_proba(X_test_std)
prob_y_3_gridsearch= [p[1] for p in prob_y_3_gridsearch]
print(roc_auc_score(y_test, prob_y_3_gridsearch))
Специальная работа с несбалансированными данными
Один из самых популярных методов - это восходящая (нисходящая) выборка. Проблема с нашим набором данных, как вы, возможно, помните, заключается в том, что алгоритмы видят значительно больше отрицательных (нулевых) случаев, чем положительных на этапе обучения, и это делает модели менее точными. Следующее, что нужно сделать, это обучить модель одинаковому количеству положительных (1) и отрицательных (0) случаев.
from sklearn.utils import resample df_minority = df[df.TenYearCHD==1] df_majority = df[df.TenYearCHD==0] df['TenYearCHD'].value_counts()
В приведенном выше коде мы разделяем существующий набор данных на два: наблюдения, при которых развилось заболевание сердца, и все остальное.
# sample with replacement to match majority class and get #reproducible results
df_minority_upsampled = resample(df_minority,
replace=True,
n_samples=3596,
random_state=123)
# Display new class counts
df_upsampled.TenYearCHD.value_counts()
Как вы теперь видите, df_upsampled имеет одинаковое количество нулевых и единичных наблюдений.
Теперь нормализуем его и снова добавим в модель:
# Train/test, normalize the new data set features_upsampled = df_upsampled.iloc[:,:-1] result_upsampled = df_upsampled.iloc[:,-1] X_train_upsampled, X_test_upsampled, y_train_upsampled, y_test_upsampled = train_test_split(features_upsampled, result_upsampled, test_size = 0.2, random_state = 14) X_train_std_upsampled = scaler.fit_transform(X_train_upsampled) X_test_std_upsampled = scaler.fit_transform(X_test_upsampled) # new log model for upsampled data logmodel_upsampled = LogisticRegression(solver='liblinear') logmodel_upsampled.fit(X_train_std_upsampled, y_train_upsampled) predictions_y_2_upsampled = logmodel_upsampled.predict(X_test_std_upsampled) print(classification_report(y_test_upsampled, predictions_y_2_upsampled)) print(confusion_matrix(y_test_upsampled, predictions_y_2_upsampled)) accuracy_score(y_test_upsampled, predictions_y_2_upsampled) # Under ROC curve prob_y_2_upsampled = logmodel_upsampled.predict_proba(X_test_std_upsampled) prob_y_2_upsampled = [p[1] for p in prob_y_2_upsampled] print(roc_auc_score(y_test_upsampled, prob_y_2_upsampled))
В нашем случае результаты ухудшились. Почему так? Мы исследуем это в следующий раз. Хотя обычно это помогает.
Мы уже упоминали, что Ошибки II более критичны в этом случае, чем Ошибки I. Другой способ учесть это - сместить порог (который теперь установлен на 0,5).
logmodel_lowering = LogisticRegression(solver='liblinear')
logmodel_lowering.fit(X_train_std, y_train)
from sklearn.preprocessing import binarize
for i in range(1,7):
cm2=0
predictions_y_2_lowering = logmodel_lowering.predict_proba(X_test_std)
y_pred2_lowering=binarize(predictions_y_2_lowering,i/10)[:,1]
cm2=confusion_matrix(y_test,y_pred2_lowering)
print ('With',i/10,'threshold the Confusion Matrix is ','\n',cm2,'\n',
'with',cm2[0,0]+cm2[1,1],'correct predictions and',cm2[1,0],'Type II errors( False Negatives)','\n\n',
'Sensitivity: ',cm2[1,1]/(float(cm2[1,1]+cm2[1,0])),'Specificity: ',cm2[0,0]/(float(cm2[0,0]+cm2[0,1])),'\n\n\n')
Мы написали здесь цикл, чтобы перейти от порогового значения 10% к порогу 70%, а затем показать результаты. Для этого упражнения я использую базовую модель логистической регрессии.
Пока мы достигли лишь ограниченного прогресса. Вероятно, пора обратиться к XGBoost - оружию, которое выбирают для большинства соревнований Kaggle.
import xgboost as xgb
from sklearn.metrics import mean_squared_error
xg_reg = xgb.XGBRegressor(objective ='binary:logistic', colsample_bytree = 0.3, learning_rate = 0.05,
max_depth = 9, alpha = 10, n_estimators = 20)
eval_set = [(X_test_std, y_test)]
xg_reg.fit(X_train_std, y_train, eval_metric="error", eval_set = eval_set, verbose = True)
rmse = np.sqrt(mean_squared_error(y_test, prediction_y_5))
print("RMSE: %f" % (rmse))
В отличие от других методов, XGBoost может сообщать и оценивать производительность набора тестов на этапе подбора. Я создал и eval_set передал его методу fit и установил verbose = True для просмотра подробностей в реальном времени.
Стандартно он возвращает среднеквадратичную ошибку (RMSE) 0,3678. Давай поработаем над его понижением. Для этого нам нужно настроить некоторые из многочисленных параметров, которые передаются в алгоритм. Вам действительно нужно прочитать о них в официальной документации, потому что это ключ к достижению превосходных характеристик вашей модели. Сегодня я остановлюсь на трех доступных параметрах. И снова GridSearchCV сделает это за нас.
Мы собираемся протестировать различные входы этих параметров:
- n_estimators (определяет размер леса для обучения)
- max_depth (максимальная глубина дерева для базовых учащихся)
- learning_rate (ну, это скорость обучения)
n_estimators = [10, 20, 30, 40, 50, 60] max_depth = [2, 4, 5, 6, 7, 8] learning_rate = [0.0001, 0.001, 0.01, 0.1, 0.2, 0.3] param_grid = dict(max_depth = max_depth, n_estimators = n_estimators, learning_rate=learning_rate) kfold = StratifiedKFold(n_splits = 10, shuffle = True, random_state = 10) grid_search_xg = GridSearchCV(xg_reg, param_grid, scoring = 'roc_auc', n_jobs = -1, cv=kfold, verbose = 1) result_gcv_xgb = grid_search_xg.fit(X_train_std, y_train)
Здесь мы выполняем перекрестную проверку, используя StratifiedKFold и GridSearchCV для перебора параметров. Чем больше параметров вы хотите протестировать, тем больше перестановок выполняется на вашем компьютере и тем больше времени на это потребуется. Будьте осторожны, если делаете это на своем ноутбуке.
Как только процесс будет завершен, давайте посмотрим, каковы лучшие параметры:
print("Best: %f using %s" % (result_gcv_xgb.best_score_, result_gcv_xgb.best_params_))
means = result_gcv_xgb.cv_results_['mean_test_score']
stds = result_gcv_xgb.cv_results_['std_test_score']
params = result_gcv_xgb.cv_results_['params']
Затем вы можете повторно запустить модель, используя то, что вы получили выше:
# rebuild using best params
xg_reg = xgb.XGBRegressor(objective ='binary:logistic', colsample_bytree = 0.3, learning_rate = 0.1,
max_depth = 2, alpha = 10, n_estimators = 50)
eval_set = [(X_test_std, y_test)]
xg_reg.fit(X_train_std, y_train, eval_metric="error", eval_set = eval_set, verbose = False)
prediction_y_5 = xg_reg.predict(X_test_std)
rmse = np.sqrt(mean_squared_error(y_test, prediction_y_5))
print("RMSE: %f" % (rmse))
Мы видим улучшение: RMSE снизился до 0,3384.
XGBoost возвращает вероятности, а не фактические прогнозы. Однако нам понадобятся реальные прогнозы для построения матрицы путаницы.
prediction_y_5_01 = prediction_y_5 prediction_y_5_01[prediction_y_5 > 0.5] = 1 prediction_y_5_01[prediction_y_5 <= 0.5] = 0 print(classification_report(y_test, prediction_y_5_01)) print(confusion_matrix(y_test, prediction_y_5_01)) accuracy_score(y_test, prediction_y_5_01)
Мы попытались использовать некоторые традиционные и более конкретные методы для построения точной модели прогнозирования. В зависимости от характера ваших данных, этого может быть уже достаточно для повышения производительности модели. Во второй части статьи мы продолжим наше путешествие, построив кривую обучения, ансамбль из нескольких моделей и, наконец, повторно проведя анализ с помощью Keras.
