Shap_XGboost_training_SA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Shap for SA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import xgboost as xgb\n",
    "import shap\n",
    "# Initialize your Jupyter notebook with initjs(), otherwise you will get an error message.\n",
    "shap.initjs()\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.model_selection import TimeSeriesSplit\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define our search space for grid search\n",
    "random_grid = {\n",
    "    'objective' : ['reg:squarederror'],\n",
    "    'n_estimators': [50, 100, 150, 200],\n",
    "    'learning_rate': [0.01, 0.1, 0.2, 0.3],\n",
    "    'max_depth': range(3, 10),\n",
    "    'colsample_bytree': [i/10.0 for i in range(1, 3)],\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Used to fix the shap error in the code\n",
    "def myfun(self=None):\n",
    "        return model_bytearray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def non_zero_variables(df):\n",
    "    cols = df.columns\n",
    "    cols_non_0 = []\n",
    "    for col in cols:\n",
    "        if ((df[col].eq(0).sum(axis=0)) < 0.4 * (len(df[col]))): #if we have more than 60% zeros\n",
    "            cols_non_0.append(col)\n",
    "    return cols_non_0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train the model and get the predictions through the SHAP methodology:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_set = 0.5\n",
    "\n",
    "for country in ['SA']: \n",
    "     \n",
    "    #Read all variables\n",
    "    all_var = pd.read_csv('../../all_variables_and_GPI_monthly_all_countries/all_variables_%s.csv'\n",
    "                          %country, index_col = 0)\n",
    "\n",
    "    #Delete all columns that have more than 60% of their values 0\n",
    "    variables_non_0 = non_zero_variables(all_var) #Filter the variables that have many zeros\n",
    "    df_country = all_var[variables_non_0]\n",
    "    \n",
    "    print(country, len(all_var))\n",
    "\n",
    "    df_country = all_var\n",
    "\n",
    "    #Set the target variable\n",
    "    Y = df_country['GPI']\n",
    "\n",
    "    #Set the independent variables\n",
    "    X = df_country.loc[:, df_country.columns != 'GPI']\n",
    "\n",
    "    #Set the training sets:\n",
    "    Y_train = Y[:int(Y.shape[0]*train_set)]\n",
    "    X_train = X[:int(X.shape[0]*train_set)]\n",
    "\n",
    "    #Set the test sets\n",
    "    X_test = X[int(X.shape[0]*train_set):]\n",
    "    Y_test = Y[int(Y.shape[0]*train_set):]\n",
    "    \n",
    "        \n",
    "    #Create a dataframe to add the predictions\n",
    "    Predictions = pd.DataFrame(columns = ['MonthYear','prediction1','prediction2','prediction3','prediction4',\n",
    "                                              'prediction5','prediction6','prediction7','prediction8','prediction9',\n",
    "                                              'prediction10','prediction11','prediction12'])\n",
    "    \n",
    "    k = 0\n",
    "    \n",
    "    for i in range(0, len(X_test)):\n",
    "        \n",
    "        print(i)\n",
    "    \n",
    "        tscv = TimeSeriesSplit(n_splits=10).split(X_train)\n",
    "        \n",
    "        l_predictions = []\n",
    "\n",
    "        #Train the model\n",
    "\n",
    "        xgb_reg = xgb.XGBRegressor() #model to tune\n",
    "\n",
    "        xgb_reg_random = GridSearchCV(estimator = xgb_reg, param_grid = random_grid,\n",
    "                                       cv = tscv,  n_jobs = -1)\n",
    "\n",
    "        #Best model\n",
    "        bmodel = xgb_reg_random.fit(X_train, Y_train)   \n",
    "\n",
    "        #Pull the best estimated model from the gridsearch and send it to the TreeExplainer\n",
    "        model = bmodel.best_estimator_\n",
    "\n",
    "        #Fix a shap error\n",
    "        mybooster = model.get_booster()\n",
    "        model_bytearray = mybooster.save_raw()[4:]\n",
    "        mybooster.save_raw = myfun\n",
    "        #Finish Fix a shap error\n",
    "        \n",
    "        #If we have less than 12 data points ahead to predict, then increase the value of the k variable\n",
    "        if (len(X_test) - i) < 12: \n",
    "            k = k + 1 \n",
    "            #print('k:',k)\n",
    "\n",
    "        #Make the prediction(s)\n",
    "        for j in range(i, i+12-k):\n",
    "            Y_pred = model.predict(X_test.iloc[[j]])\n",
    "            l_predictions.append(Y_pred[0])\n",
    "\n",
    "        #k variable helps me understand how many months ahead I can predict\n",
    "        if k>0:\n",
    "            l_predictions2 = []\n",
    "            for l in range(0, k):\n",
    "                l_predictions2.append('-')\n",
    "            #Concatenate the predictions list and the '-' list to add them on the dataframe\n",
    "            l_predictions = l_predictions + l_predictions2\n",
    "\n",
    "        #Concatenate the month and Year of the last training with the prediction list to add them on the dataframe\n",
    "        l_predictions = [str(X_train.index[-1])] + l_predictions\n",
    "\n",
    "        #Add monthyear and predictions on the dataframe\n",
    "        Predictions_length = len(Predictions)\n",
    "        Predictions.loc[Predictions_length] = l_predictions\n",
    "\n",
    "         \n",
    "        #Set the new training sets\n",
    "        X_train = X_train.iloc[1:]\n",
    "        X_train = X_train.append(X_test.iloc[i])\n",
    "        Y_train = Y_train[1:]\n",
    "        Y_train = Y_train.append(pd.Series(Y_test.iloc[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the shap predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Predictions.to_csv('../../%s_shap_predictions.csv' %country)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Variables importance analysis through SHAP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for country in ['SA']: \n",
    "     \n",
    "    #Read all variables\n",
    "    all_var = pd.read_csv('../../all_variables_%s.csv' %country, index_col = 0)\n",
    "\n",
    "    #Delete all columns that have more than 60% of their values 0\n",
    "    variables_non_0 = non_zero_variables(all_var) #Filter the variables that have many zeros\n",
    "    df_country = all_var[variables_non_0]\n",
    "    \n",
    "    print(country, len(all_var))\n",
    "\n",
    "    df_country = all_var\n",
    "    X_future = df_country.loc[df_country.index >= 201805] \n",
    "    df_country = df_country.loc[df_country.index <= 201804]\n",
    "    \n",
    "\n",
    "    #Set the target variable\n",
    "    Y = df_country['GPI']\n",
    "\n",
    "    #Set the independent variables\n",
    "    X = df_country.loc[:, df_country.columns != 'GPI']\n",
    "\n",
    "    #Set the training sets:\n",
    "    X_train = X.tail(72) #corresponds to the lenght of the dataset\n",
    "    Y_train = Y.tail(72) #corresponds to the lenght of the dataset\n",
    "    \n",
    "    tscv = TimeSeriesSplit(n_splits=10).split(X_train)\n",
    "\n",
    "    #Train the model\n",
    "\n",
    "    xgb_reg = xgb.XGBRegressor() #model to tune\n",
    "\n",
    "    xgb_reg_random = GridSearchCV(estimator = xgb_reg, param_grid = random_grid,\n",
    "                                   cv = tscv,  n_jobs = -1)\n",
    "\n",
    "    #Best model\n",
    "    bmodel = xgb_reg_random.fit(X_train, Y_train)   \n",
    "    \n",
    "    #Pull the best estimated model from the gridsearch and send it to the TreeExplainer\n",
    "    model = bmodel.best_estimator_\n",
    "    \n",
    "    #Fix a shap error\n",
    "    mybooster = model.get_booster()\n",
    "    model_bytearray = mybooster.save_raw()[4:]\n",
    "    mybooster.save_raw = myfun\n",
    "    #Finish Fix a shap error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Global variable importance plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get the shap values\n",
    "shap_values = shap.TreeExplainer(mybooster).shap_values(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "shap.summary_plot(shap_values, X_train, alpha=1, plot_type=\"bar\", max_display=10, show=False)\n",
    "                            \n",
    "plt.title('Global variable importance Saudi Arabia XGBoost model ', loc='right', fontsize=25, fontweight='bold' )\n",
    "\n",
    "plt.savefig('../../shap_summary_bar_%s.pdf' %country, bbox_inches='tight')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Individual Shap Value plot for October 2018:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shap_plot(j):\n",
    "    explainerModel = shap.TreeExplainer(mybooster)\n",
    "  \n",
    "    shap_values_Model = explainerModel.shap_values(X_future[X_future.columns[1:]])\n",
    "    p = shap.force_plot(explainerModel.expected_value, shap_values_Model[j], X_future[X_future.columns[1:]].iloc[[j]].round(0).astype(object), show=False,matplotlib=True)\n",
    "                                               \n",
    "    p.set_figwidth(40)\n",
    "    p.set_figheight(3)\n",
    "\n",
    "    plt.tick_params(axis='x', labelsize = 12)\n",
    "    plt.suptitle('Saudi Arabia XGBoost model\\n Prediction for October 2018', fontsize=25, fontweight='bold', x=0.5, y=1.27)\n",
    "    p.savefig('../../shap_plot_%s.pdf'%country, bbox_inches='tight')\n",
    "    return(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shap_plot(5)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}