Источник данных - https://www.kaggle.com/ivanloginov/the-broken-machine
Автор: Потапов Даниил (Slack: @sharthZ23)
В данном проекте ставится задача прогнозирования поломки оборудования при помощи его индикаторов, их здесь почти 60 и все они безымянные. Данная задача - это пример использования алгоритмов машинного обучения в системах обнаружения неполадок. Область применения широка: заводское производство, спутники или другие автономные объекты и так далее. Причем тут надо рассмотреть 2 варианта решений: одно направленное на интерпретируемость, а другое на точность.
Датасет содержит 900000 объектов, каждый из которых содержит 59 признаков, в том числе 1 целевой (он находится в ytrain.csv
, столбец x
).
Как я уже сказал выше, признаки полностью анонимные, единственное, что их объединяет, это то, что они все float64
.
Целевой признак - индикатор того, сломана ли машина (1 - да, 0 - нет).
Импортируем нужные библиотеки
import os
import time
import pandas as pd
import numpy as np
import xgboost as xgb
import lightgbm as lgb
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier, BaggingClassifier, VotingClassifier
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.metrics import roc_auc_score, recall_score
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn_pandas import DataFrameMapper
import category_encoders as ce
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
import statsmodels.discrete.discrete_model as sm
from scipy import stats
stats.chisqprob = lambda chisq, df: stats.chi2.sf(chisq, df)
Укажем пути до данных и загрузим их
DATA_PATH = '/data/kaggle/broken_machine'
X_PATH = os.path.join(DATA_PATH, 'xtrain.csv')
Y_PATH = os.path.join(DATA_PATH, 'ytrain.csv')
X = pd.read_csv(X_PATH)
y = pd.read_csv(Y_PATH)['x']
print(X.shape, y.shape)
Увеличим максимальное кол-во отображаемых столбцов в pandas
до 58
pd.set_option('display.max_columns', 100)
Посмотрим на информацию о нашем датасете
X.info()
Видно, что все параметры - числа с плавающей точкой и что около 1/9 части каждого параметра отсутствует. Посмотрим на сами данные.
X.head()
У же на первых 5 строчках видно, что с NaN'ми придется что-то делать. Попробуем просто их удалить.
print('Initial size: {}'.format(X.shape))
print('After NaN omit size: {}'.format(X.dropna().shape))
Ой, мы потеряли 99% данных. Не получилось и ладно, вспомним о них на стадии feature engineering, а пока посмотрим на первичные числовые признаки каждого параметра.
X.describe()
Сразу бросается в глаза, что все параметры можно разделить на 2 части: те у кого (min
, N%
, max
) не содержат числа после запятой (отличные от 0) и наоборот.
Это говорит о том, что у некоторых параметров мало уникальных значений и они, возможно, больше категориальные, чем числовые. Проверим эту гипотезу.
nunique_x = X.nunique()
axes = nunique_x.plot(kind='barh', figsize=(16, 12))
Гипотеза подтвердилась, наши признаки разделились на 2 группы: в одной у каждого параметра уникальных значений много, около 800000, в другой же их экстремально мало, по сравнению с первой группой. По графику мы можем их четко разделить по середине, поэтому проделаем это.
num_cat_mask = nunique_x > 400000
num_cols = nunique_x[num_cat_mask].keys().tolist()
cat_cols = nunique_x[~num_cat_mask].keys().tolist()
Посмотрим распределения числовых признаков.
nrows, ncols = 8, 4
figsize = (nrows * 3, ncols * 6)
axes = X[num_cols].hist(figsize=figsize, bins=100, color='dodgerblue')
plt.tight_layout()
Красота-то какая, почти все признаки нормально распределены. И даже больше, их распределения очень схожи, пиковых значений примерной одинаковое значение, значит тут можно и нужно(не всегда) применить PCA. Но признаки 12
и 37
выбиваются из общей нормальной массы, рассмотрим их подробнее.
axes = X[['12', '37']].hist(figsize=(12, 8), bins=100, color='dodgerblue')
Видно, что 12
признак распределен более менее равномерно, а вот 37
нет. Более того, у 37
подавляющее большинство значений лежит около 0.
Делаем вывод, что к 12
мы можем применить любой Scaler, а вот к 37
кроме StandarScaler ничего не подойдет.
Теперь рассмотрим категориальные значения.
plt.rcParams['axes.labelsize'] = 20
plt.rcParams['axes.titlesize'] = 20
plt.rcParams['font.size'] = 14
nrows, ncols = 7, 4
figsize = (nrows * 3, ncols * 6)
fig, axes = plt.subplots(ncols=ncols, nrows=nrows, figsize=figsize)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=0.4)
for i in range(len(cat_cols)):
ax = axes[i//4, i%4]
col = cat_cols[i]
X[col].value_counts(normalize=True) \
.plot(kind='bar', label=col, ax=ax, color='dodgerblue')
ax.set_title(col)
plt.tight_layout()
А тут у нас везде что-то похожее на распределение Парето (Брэдфорда). Посмотрим на целевую пременную.
print(y.value_counts(normalize=True))
axes = y.value_counts(normalize=True).plot(kind='bar', color='dodgerblue')
Видим дисбаланс классов, 70:30. Это надо будет учесть, когда будем выбирать, какие алгоритмы машинного обучения применять.
Теперь построим и визуализируем матрицу корреляций. Добавим в X
целевую переменную, чтобы увидеть, если зависимости у остальных признаков с ней.
X['y'] = y
%%time
corr_mat = X.corr()
plt.figure(figsize=(18, 12))
axes = sns.heatmap(corr_mat, cmap="YlGnBu")
Увы, сплошное ничего. Может попробуем другой метод в df.corr
?
%%time
corr_mat = X.corr(method='spearman')
plt.figure(figsize=(18, 12))
axes = sns.heatmap(corr_mat, cmap="YlGnBu")
Не стоило, опять все около 0, но времени ушло гораздо, гораздо больше. С method='kendall'
тоже самое :(
Теперь посмотрим, как распределены признаки относительно целевого. Начнем c численных признаков.
%%time
nrows, ncols = 8, 4
figsize = (nrows * 2, ncols * 6)
fig, axes = plt.subplots(ncols=ncols, nrows=nrows, figsize=figsize)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=0.4)
for i in range(len(num_cols)):
ax = axes[i//4, i%4]
col = num_cols[i]
X.plot(x=col, y='y', kind='scatter', label=col, ax=ax, color='dodgerblue')
ax.set_title(col)
plt.tight_layout()
Никаких зависимостей не видно, зато есть точки-кандидаты на метку "выброс". Но их кол-во незначительно, поэтому пока оставим их в покое.
Теперь тоже самое для категориальных.
%%time
nrows, ncols = 7, 4
figsize = (nrows * 3, ncols * 6)
fig, axes = plt.subplots(ncols=ncols, nrows=nrows, figsize=figsize)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=0.4)
for i in range(len(cat_cols)):
ax = axes[i//4, i%4]
col = cat_cols[i]
sns.countplot(x=col, hue='y', data=X, ax=ax, color='dodgerblue')
ax.set_title(col)
plt.tight_layout()
И тут пустота, на взгляд никаких корреляций. Попробуем посмотреть через долю позитивного класса.
nrows, ncols = 7, 4
figsize = (nrows * 3, ncols * 6)
fig, axes = plt.subplots(ncols=ncols, nrows=nrows, figsize=figsize)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=0.4)
for i in range(len(cat_cols)):
ax = axes[i//4, i%4]
col = cat_cols[i]
X.groupby(col)['y'].mean().plot(kind='bar', ax=ax, color='dodgerblue')
ax.set_title(col)
plt.tight_layout()
А тут уже интересней. Половина признаков почти не несет никакой информации о целевом классе, но у некоторых признаков есть значение, которое прямо сигнализирует о принадлежности объекта к целевому классу. Здесь должны хорошо отработать OHE и TargetEncoding.
Сложно говорить о каких-то инсайтах с анонимными фичами, но все же:
X['48']==13.0
почти у всех объектов с X['y']==1
) , т.е. которые хорошо описывают целевой класс.Основной метрикой для оценки качества модели будет ROC-AUC. Она хорошо справляется с несбалансированными классами. Но также не стоит забывать, что мы делаем алгоритм для систем нахождения неполадок и здесь ошибки первого (FalsePositive) и второго рода (FalseNegative) не раноправны . Ведь затраты на лишнюю проверку обычно ниже, чем убыток от вовремя несработавшей системы. Исходя из этого, надо также внимательно следить за метрикой Recall.
Итого:
В данном датасете есть как числовые, так и категориальные признаки, поэтому градиентный бустинг - наш выбор. Они показывают хорошие результаты на данных такого типа и в тоже время интерпретируемы. Но это тяжеловесные методы, не будем забывать о старой доброй логистической регрессии и случайном лесе, возможно выигрыш по времени будет гораздо значительней, чем по точности.
Итого:
Настало время вспомнить о том, что у нас очень много NaN
. Поступим с ними так: в числовых признаках заменим пустые значения средним, а в категориальных - медианой.
for col in num_cols:
X[col] = X[col].fillna(X[col].mean())
for col in cat_cols:
X[col] = X[col].fillna(X[col].median())
print("Count of NaN's in X - {}".format(X.isnull().sum().sum()))
Прежде чем приступать к обучению моделей, попробуем снизить размерность, датасет не маленький ведь. Для этого обучим логистическую регрессию, но не из scikit-learn
, а из statsmodel
. Это делаем для того, чтобы получить более подробный отчет о значимости признаков. Но больше всего там нас будет интересовать 2 вещи:
model = sm.Logit(y.values, X.drop(columns='y'))
result = model.fit()
print(result.summary())
Видно, что Pseudo R-squ.
очень и очень мал, что означает наша модель не лучше, чем просто предсказывать среднее значение. Попробуем теперь взять только те значение, у которых p-value (P>|z|
) меньше или равна 10%.
sig_columns = [i for i,x in enumerate(result.pvalues.ravel()) if x<=0.1]
model = sm.Logit(y, X.iloc[:,sig_columns])
result = model.fit()
print(result.summary())
Мда, мы сократили кол-во признаков до 11, но получили негативный Pseudo R-squ.
. Это означает, что наша модель теперь хуже, чем просто предсказывать среднее значение. Тоже само будет, если мы ограничем p-value 1 процентом.
Теперь займемся обработкой. Числовые признаки обработаем с помощью StandardScaler
.
scaler = StandardScaler()
scaler.fit(X[num_cols])
#X_nums = scaler.transform(X[num_cols])
X[num_cols] = scaler.transform(X[num_cols])
Категориальные признаки обработаем с помощью OneHotEncoder
и TargetEncoder
.
cat_union = FeatureUnion([
('ohe', ce.OneHotEncoder()),
('target', ce.TargetEncoder())
], n_jobs=-1)
cat_union.fit(X[cat_cols], y)
X_cats = cat_union.transform(X[cat_cols])
print(X_cats.shape)
Теперь объединим обработанные данные с помощью np.hstack
X_train = np.hstack((X_nums, X_cats))
print(X_train.shape)