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