feature-selection/model-based/sklearn-feature-selection-shap.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scikit-Learn Gradient Boosted Tree Feature Selection with Shapley Importance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook explains how to use Shapley importance from `SHAP` and a `scikit-learn` tree-based model to perform feature selection.\n",
    "\n",
    "This notebook will work with an OpenML dataset to predict who pays for internet with 10108 observations and 69 columns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Packages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial uses:\n",
    "* [pandas](https://pandas.pydata.org/docs/)\n",
    "* [statsmodels](https://www.statsmodels.org/stable/index.html)\n",
    "    * [statsmodels.api](https://www.statsmodels.org/stable/api.html#statsmodels-api)\n",
    "* [numpy](https://numpy.org/doc/stable/)\n",
    "* [scikit-learn](https://scikit-learn.org/stable/)\n",
    "    * [sklearn.metrics](https://scikit-learn.org/stable/modules/model_evaluation.html)\n",
    "    * [sklearn.model_selection](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html)\n",
    "    * [sklearn.ensemble](https://scikit-learn.org/stable/modules/ensemble.html#ensemble)\n",
    "* [shap](https://shap.readthedocs.io/en/latest/index.html)\n",
    "* [category_encoders](https://contrib.scikit-learn.org/category_encoders/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import train_test_split\n",
    "import shap\n",
    "import category_encoders as ce\n",
    "\n",
    "from sklearn.ensemble import GradientBoostingRegressor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The data is from `rdatasets` imported using the Python package `statsmodels`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 336776 entries, 0 to 336775\n",
      "Data columns (total 19 columns):\n",
      " #   Column          Non-Null Count   Dtype  \n",
      "---  ------          --------------   -----  \n",
      " 0   year            336776 non-null  int64  \n",
      " 1   month           336776 non-null  int64  \n",
      " 2   day             336776 non-null  int64  \n",
      " 3   dep_time        328521 non-null  float64\n",
      " 4   sched_dep_time  336776 non-null  int64  \n",
      " 5   dep_delay       328521 non-null  float64\n",
      " 6   arr_time        328063 non-null  float64\n",
      " 7   sched_arr_time  336776 non-null  int64  \n",
      " 8   arr_delay       327346 non-null  float64\n",
      " 9   carrier         336776 non-null  object \n",
      " 10  flight          336776 non-null  int64  \n",
      " 11  tailnum         334264 non-null  object \n",
      " 12  origin          336776 non-null  object \n",
      " 13  dest            336776 non-null  object \n",
      " 14  air_time        327346 non-null  float64\n",
      " 15  distance        336776 non-null  int64  \n",
      " 16  hour            336776 non-null  int64  \n",
      " 17  minute          336776 non-null  int64  \n",
      " 18  time_hour       336776 non-null  object \n",
      "dtypes: float64(5), int64(9), object(5)\n",
      "memory usage: 48.8+ MB\n"
     ]
    }
   ],
   "source": [
    "df = sm.datasets.get_rdataset('flights', 'nycflights13').data\n",
    "df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature Engineering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Handle null values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As this model will predict arrival delay, the `Null` values are caused by flights did were cancelled or diverted. These can be excluded from this analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.dropna(inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert the times from floats or ints to hour and minutes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['arr_hour'] = df.arr_time.apply(lambda x: int(np.floor(x/100)))\n",
    "df['arr_minute'] = df.arr_time.apply(lambda x: int(x - np.floor(x/100)*100))\n",
    "df['sched_arr_hour'] = df.sched_arr_time.apply(lambda x: int(np.floor(x/100)))\n",
    "df['sched_arr_minute'] = df.sched_arr_time.apply(lambda x: int(x - np.floor(x/100)*100))\n",
    "df['sched_dep_hour'] = df.sched_dep_time.apply(lambda x: int(np.floor(x/100)))\n",
    "df['sched_dep_minute'] = df.sched_dep_time.apply(lambda x: int(x - np.floor(x/100)*100))\n",
    "df.rename(columns={'hour': 'dep_hour',\n",
    "                   'minute': 'dep_minute'}, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare data for modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up train-test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "target = 'arr_delay'\n",
    "y = df[target]\n",
    "X = df.drop(columns=[target, 'flight', 'tailnum', 'time_hour', 'year', 'dep_time', 'sched_dep_time', 'arr_time', 'sched_arr_time', 'dep_delay'])\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2, random_state=1066)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Encode categorical variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use a leave-one-out encoder as it creates a single column for each categorical variable instead of creating a column for each level of the categorical variable like one-hot-encoding.  This makes interpreting the impact of categorical variables with feature impact easier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "encoder = ce.LeaveOneOutEncoder(return_df=True)\n",
    "X_train_loo = encoder.fit_transform(X_train, y_train)\n",
    "X_test_loo = encoder.transform(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43.10690848187833"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = GradientBoostingRegressor(learning_rate=0.05, max_depth=5, n_estimators=500, min_samples_split=5, n_iter_no_change=10)\n",
    "model.fit(X_train_loo, y_train)\n",
    "\n",
    "rmse = np.sqrt(mean_squared_error(y_test, model.predict(X_test_loo)))\n",
    "rmse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature selection using mean absolute value of the SHAP values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a data frame to hold the SHAP values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>features</th>\n",
       "      <th>importance</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>arr_hour</td>\n",
       "      <td>2.331481</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>dep_hour</td>\n",
       "      <td>1.112755</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>origin</td>\n",
       "      <td>1.081143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>sched_dep_hour</td>\n",
       "      <td>0.747311</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>sched_arr_hour</td>\n",
       "      <td>0.358844</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>distance</td>\n",
       "      <td>0.169818</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>carrier</td>\n",
       "      <td>0.115962</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>dest</td>\n",
       "      <td>0.086835</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>sched_arr_minute</td>\n",
       "      <td>0.078803</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>arr_minute</td>\n",
       "      <td>0.057935</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>air_time</td>\n",
       "      <td>0.027669</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>dep_minute</td>\n",
       "      <td>0.014226</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>day</td>\n",
       "      <td>0.012283</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>sched_dep_minute</td>\n",
       "      <td>0.007768</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>month</td>\n",
       "      <td>0.001376</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            features  importance\n",
       "9           arr_hour    2.331481\n",
       "7           dep_hour    1.112755\n",
       "3             origin    1.081143\n",
       "13    sched_dep_hour    0.747311\n",
       "11    sched_arr_hour    0.358844\n",
       "6           distance    0.169818\n",
       "2            carrier    0.115962\n",
       "4               dest    0.086835\n",
       "12  sched_arr_minute    0.078803\n",
       "10        arr_minute    0.057935\n",
       "5           air_time    0.027669\n",
       "8         dep_minute    0.014226\n",
       "1                day    0.012283\n",
       "14  sched_dep_minute    0.007768\n",
       "0              month    0.001376"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "explainer = shap.Explainer(model)\n",
    "shap_values = explainer(X_test_loo)\n",
    "shap_importance = shap_values.abs.mean(0).values\n",
    "importance_df = pd.DataFrame({'features': X_train_loo.columns,\n",
    "                              'importance': shap_importance})\n",
    "importance_df.sort_values(by='importance', ascending=False, inplace=True)\n",
    "importance_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a list of the features with Gini importance greater than `0.5` and use that list to retrain the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['arr_hour', 'dep_hour', 'origin', 'sched_dep_hour']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "feature_list = importance_df[importance_df.importance > 0.5]['features'].tolist()\n",
    "feature_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, to keep the top 5 features, use the following instead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['arr_hour', 'dep_hour', 'origin', 'sched_dep_hour', 'sched_arr_hour']"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "feature_list = importance_df['features'].head(5).tolist()\n",
    "feature_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_loo_new = X_train_loo[feature_list]\n",
    "X_test_loo_new = X_test_loo[feature_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43.6599708368277"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reduced_model = GradientBoostingRegressor(learning_rate=0.05, max_depth=5, n_estimators=500, min_samples_split=5, n_iter_no_change=10)\n",
    "reduced_model.fit(X_train_loo_new, y_train)\n",
    "\n",
    "rmse = np.sqrt(mean_squared_error(y_test, reduced_model.predict(X_test_loo_new)))\n",
    "rmse"
   ]
  }
 ],
 "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}