notebooks/02_residual_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Residual Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. ML Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import pandas as pd\n",
    "from loguru import logger\n",
    "from sklearn.model_selection import train_test_split\n",
    "import lightgbm as lgb\n",
    "import matplotlib.pyplot as plt\n",
    "from blog.model.regularisation import select_features, get_monotone_constraints\n",
    "from blog.model.evaluation import show_model_results\n",
    "from probatus.interpret import ShapModelInterpreter\n",
    "from blog.data.lendingclub_dataset import LendingClubDataset\n",
    "from blog.model.sampling import oot_split\n",
    "from blog.model.regularisation import get_monotone_constraints\n",
    "from blog.model.feature_transformations import get_simple_feature_transformation\n",
    "from blog.model.optimise import getObjective, tune_model\n",
    "from blog.model.evaluation import get_segment_rocauc\n",
    "from blog.model.target_transformations import get_sample_weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_name = \"lendingclub\"\n",
    "target = \"loan_status\"\n",
    "path = \"../data/lc_accepted.csv\"\n",
    "\n",
    "initial_features = [\n",
    "    \"verification_status\",\n",
    "    \"emp_title\",\n",
    "    \"int_rate\",\n",
    "    \"loan_amnt\",\n",
    "    \"total_rec_int\",\n",
    "    \"total_acc\",\n",
    "    \"tot_cur_bal\",\n",
    "    \"fico_range_low\",\n",
    "    \"fico_range_high\",\n",
    "    \"grade\",\n",
    "    \"total_rev_hi_lim\",\n",
    "    \"sub_grade\",\n",
    "    \"initial_list_status\",\n",
    "    \"purpose\",\n",
    "    \"issue_d\",\n",
    "    \"emp_length\",\n",
    "    \"pub_rec_bankruptcies\",\n",
    "    \"last_pymnt_amnt\",\n",
    "    \"num_actv_bc_tl\",\n",
    "    \"total_pymnt\",\n",
    "    \"loan_status\",\n",
    "    \"term\",\n",
    "    \"home_ownership\",\n",
    "    \"revol_util\",\n",
    "    \"application_type\",\n",
    "    \"addr_state\",\n",
    "    \"inq_last_6mths\",\n",
    "    \"pub_rec\",\n",
    "    \"dti\",\n",
    "    \"mort_acc\",\n",
    "    \"revol_bal\",\n",
    "    \"title\",\n",
    "    \"annual_inc\",\n",
    "    \"out_prncp\",\n",
    "    \"open_acc\",\n",
    "    \"installment\",\n",
    "    \"home_ownership\",\n",
    "    \"avg_cur_bal\",\n",
    "    \"annual_inc\",\n",
    "    \"num_tl_90g_dpd_24m\"\n",
    "]\n",
    "\n",
    "# combine two lists without duplicates\n",
    "initial_features = list(set(initial_features))\n",
    "print(initial_features)\n",
    "# ================== Read and clean the dataset.===========================================\n",
    "logger.info(f\"1.Read data\")\n",
    "lcd = LendingClubDataset()\n",
    "X, y = lcd.get_data(path=path, use_cols=initial_features, dropna=False, sample=-1)\n",
    "logger.info(f\"\\nTarget distribution:\\n{y.value_counts(normalize=True)}\")\n",
    "\n",
    "# Split data into IT & OOT datasets.\n",
    "data = X.copy()\n",
    "data[target] = y\n",
    "\n",
    "X_it, y_it, X_oot, y_oot = oot_split(\n",
    "    df=data, split_date=\"Jun-2016\", split_col=\"issue_d\", target_col=target\n",
    ")\n",
    "# Original data dictionary\n",
    "data_dict = {\"xtrain\": X_it, \"ytrain\": y_it, \"xtest\": X_oot, \"ytest\": y_oot}\n",
    "\n",
    "# =================== Apply feature transformation =========================================\n",
    "logger.info(f\"2.Transform features\")\n",
    "#data_dict = get_simple_feature_transformation(data_dict_org)\n",
    "\n",
    "# =================== Feature Selection =========================================\n",
    "logger.info(f\"3.Starting feature selection.\")\n",
    "# Select the best features based on SHAPRFE CV\n",
    "# selected_features, fs_plot = select_features(data=data_dict, verbose=100)\n",
    "\n",
    "# selected_features=['loan_amnt','int_rate','installment','grade','annual_inc','home_ownership',\n",
    "#                    'emp_length','term','addr_state', \n",
    "#                    'fico_range_low', 'fico_range_high']\n",
    "\n",
    "selected_features = [\n",
    "    \"emp_length\",\n",
    "    \"addr_state\",\n",
    "    \"revol_util\",\n",
    "    \"revol_bal\",\n",
    "    \"term\",\n",
    "    \"num_actv_bc_tl\",\n",
    "    \"total_acc\",\n",
    "    \"application_type\",\n",
    "    \"purpose\",\n",
    "    \"grade\",\n",
    "    'dti',\n",
    "    \"home_ownership\",\n",
    "    \"annual_inc\",\n",
    "    \"num_tl_90g_dpd_24m\"\n",
    "]\n",
    "logger.info(f\"Final features :  {selected_features}\")\n",
    "# Update the data dictionary with the selected features.\n",
    "data_dict[\"xtrain\"] = data_dict[\"xtrain\"][selected_features]\n",
    "data_dict[\"xtest\"] = data_dict[\"xtest\"][selected_features]\n",
    "\n",
    "model_params = {\n",
    "    \"objective\": \"binary\",\n",
    "    \"metric\": \"binary_logloss\",\n",
    "    \"verbosity\": -1,\n",
    "    \"boosting_type\": \"gbdt\",\n",
    "    \"feature_pre_filter\": False,\n",
    "    \"lambda_l1\": 3.221919143923753e-08,\n",
    "    \"lambda_l2\": 0.0016740960905730054,\n",
    "    \"num_leaves\": 66,\n",
    "    \"feature_fraction\": 0.6,\n",
    "    \"bagging_fraction\": 1.0,\n",
    "    \"bagging_freq\": 0,\n",
    "    \"min_child_samples\": 20,\n",
    "}\n",
    "\n",
    "# model_params = tune_model(data_dict[\"xtrain\"], data_dict[\"ytrain\"], \"rocauc\")\n",
    "\n",
    "# =================== Train Model =========================================\n",
    "logger.info(f\"5.Train Model \")\n",
    "logger.info(f\"Creating model with params : {model_params}\")\n",
    "mono_lgb = lgb.LGBMClassifier(**model_params)\n",
    "sample_weight = get_sample_weights(data_dict)\n",
    "mono_model = show_model_results(data=data_dict, model=mono_lgb,sample_weight=sample_weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a validation set, to do the error analysis.\n",
    "X_train, X_val, y_train, y_val = train_test_split(data_dict[\"xtrain\"], data_dict[\"ytrain\"], test_size=0.2,stratify=data_dict[\"ytrain\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train the model\n",
    "lgb_model = lgb.LGBMClassifier(**model_params)\n",
    "sample_weight=get_sample_weights(data_dict={'xtrain':X_train,'ytrain':y_train})\n",
    "lgb_model.fit(X_train,y_train,sample_weight=sample_weight)\n",
    "\n",
    "yhat_name = 'p_loan_status'\n",
    "preds1 = lgb_model.predict_proba(X_val)[:,1]\n",
    "\n",
    "X_val_results = X_val.copy().reset_index(drop=True)\n",
    "X_val_results[yhat_name] = preds1\n",
    "X_val_results[target] = y_val.values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the log loss error.\n",
    "from sklearn.metrics import log_loss\n",
    "import numpy as np\n",
    "resid_name = f'r_{target}' \n",
    "# calculate logloss residuals\n",
    "X_val_results[resid_name] = -X_val_results[target]*np.log(X_val_results[yhat_name]) -(1 - X_val_results[target])*np.log(1 - X_val_results[yhat_name])   \n",
    "# check that logloss is calculated correctly\n",
    "print('Mean logloss residual: %.6f' % X_val_results[resid_name].mean())\n",
    "print('Logloss from sklearn %.6f' % log_loss(y_val, preds1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize figure\n",
    "# Matplotlib figure inline\n",
    "%matplotlib inline   \n",
    "fig, ax_ = plt.subplots(figsize=(8, 8))         \n",
    "valid_yhat_df = X_val_results.copy()\n",
    "# plot groups with appropriate color\n",
    "color_list = ['#000066', '#FF6200'] \n",
    "c_idx = 0\n",
    "groups = valid_yhat_df.groupby(target) # define groups for levels of PAY_0\n",
    "for name, group in groups:\n",
    "    ax_.plot(group.p_loan_status, group.r_loan_status, \n",
    "             label=' '.join([target, str(name)]),\n",
    "             marker='o', linestyle='', color=color_list[c_idx], alpha=0.3)\n",
    "    c_idx += 1\n",
    "    \n",
    "# annotate plot\n",
    "plt.xlabel(yhat_name)\n",
    "plt.ylabel(resid_name)\n",
    "ax_.legend(loc=1)\n",
    "plt.title('Global Logloss Residuals')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Some high-magnitude outlying residuals are visible. \n",
    "* Who are these customers? \n",
    "* Why is the model so wrong about them? \n",
    "* And are they somehow exerting undue influence on other predictions? \n",
    "\n",
    "The model could be retrained without these individuals and retested as a potentially remediation strategy.\n",
    "\n",
    "`loan_status = 1 Residuals`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = [0, 30, 60, 100]\n",
    "labels = ['low','medium','high']\n",
    "valid_yhat_df['dti_binned'] = pd.cut(valid_yhat_df['dti'], bins=bins, labels=labels)\n",
    "valid_yhat_df1 = valid_yhat_df[valid_yhat_df[target] == 1]\n",
    "valid_yhat_df1 = valid_yhat_df1.sort_values(by=f'r_{target}', ascending=False).reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use Seaborn FacetGrid for convenience\n",
    "# some seaborn configs\n",
    "import seaborn as sns\n",
    "sns.set(font_scale=0.9)                                         # legible font size\n",
    "sns.set_style('whitegrid', {'axes.grid': False})                # white background, no grid in plots\n",
    "sns.set_palette(sns.color_palette(['#000066', '#FF6200']))            # consistent colors\n",
    "\n",
    "# facet grid of residuals by PAY_0 \n",
    "sorted_ = valid_yhat_df.sort_values(by='grade')                 # sort for better layout of by-groups\n",
    "g = sns.FacetGrid(sorted_, col='grade', hue=target, col_wrap=4) # init grid\n",
    "_ = g.map(plt.scatter, yhat_name, resid_name, alpha=0.4)        # plot points\n",
    "_ = g.add_legend(bbox_to_anchor=(0.82, 0.2))                    # legend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# facet grid of residuals by PAY_0 \n",
    "sorted_ = valid_yhat_df.sort_values(by='dti_binned')                 # sort for better layout of by-groups\n",
    "g = sns.FacetGrid(sorted_, col='dti_binned', hue=target, col_wrap=4) # init grid\n",
    "_ = g.map(plt.scatter, yhat_name, resid_name, alpha=0.4)        # plot points\n",
    "_ = g.add_legend(bbox_to_anchor=(0.82, 0.2))                    # legend"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Residual Decision Tree\n",
    "\n",
    "A simple decision tree can be built, to highlight the failure modes for the model.\n",
    "We can futher simplfy it by choosing only the highly important features and creating a high cutoff for the residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top5_feats = ['grade','dti','term','home_ownership','num_actv_bc_tl']\n",
    "# Drop missing values\n",
    "valid_yhat_df1 = valid_yhat_df1.dropna()\n",
    "X_err = valid_yhat_df1[top5_feats]\n",
    "y_err = valid_yhat_df1['r_loan_status'].apply(lambda x: 1 if x >= 1.5 else 0)\n",
    "print(X_err.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "valid_yhat_df1['r_loan_status'].hist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.tree import DecisionTreeClassifier,DecisionTreeRegressor\n",
    "from sklearn.model_selection import StratifiedKFold\n",
    "\n",
    "X_ohe = pd.get_dummies(X_err)\n",
    "clf =  DecisionTreeClassifier(max_depth=2,random_state=0)\n",
    "mean_auc = cross_val_score(clf, X_ohe, y_err, \n",
    "                           cv=StratifiedKFold(n_splits=5, shuffle=True, random_state=0),\n",
    "                           scoring='roc_auc').mean()\n",
    "print(f\"AUC on residuals : {mean_auc}\")\n",
    "clf.fit(X_ohe,y_err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dtreeviz.trees import dtreeviz # remember to load the package\n",
    "viz = dtreeviz(clf, X_ohe, y,\n",
    "                target_name=\"target\",\n",
    "                feature_names=X_ohe.columns,\n",
    "                class_names=[\"low residual\", \"high residual\"],)\n",
    "\n",
    "viz"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "9ee18155c7f08bca26f68a2e7822e37aab5c5efb1655e9dcf78a129488236b07"
  },
  "kernelspec": {
   "display_name": "Python 3.6.10 ('gdd3.6')",
   "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.6.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}