{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 10: Logistic Regression\n",
    "\n",
    "> (c) 2019 Galit Shmueli, Peter C. Bruce, Peter Gedeck \n",
    ">\n",
    "> Code included in\n",
    ">\n",
    "> _Data Mining for Business Analytics: Concepts, Techniques, and Applications in Python_ (First Edition) \n",
    "> Galit Shmueli, Peter C. Bruce, Peter Gedeck, and Nitin R. Patel. 2019.\n",
    "\n",
    "## Import required packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:05.117357Z",
     "iopub.status.busy": "2022-02-05T23:04:05.111412Z",
     "iopub.status.idle": "2022-02-05T23:04:06.295120Z",
     "shell.execute_reply": "2022-02-05T23:04:06.295431Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no display found. Using non-interactive Agg backend\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.9/site-packages/statsmodels/compat/pandas.py:65: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.\n",
      "  from pandas import Int64Index as NumericIndex\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.linear_model import LogisticRegression, LogisticRegressionCV\n",
    "from sklearn.model_selection import train_test_split\n",
    "import statsmodels.api as sm\n",
    "from mord import LogisticIT\n",
    "import matplotlib.pylab as plt\n",
    "import seaborn as sns\n",
    "import dmba\n",
    "from dmba import classificationSummary, gainsChart, liftChart\n",
    "from dmba.metric import AIC_score\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 10.2\n",
    "Load the `UniversalBank.csv` dataset. The columns `ID` and `ZIP Code` are not relevant for model building and therefore removed. Treat Education as categorical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:06.302724Z",
     "iopub.status.busy": "2022-02-05T23:04:06.302336Z",
     "iopub.status.idle": "2022-02-05T23:04:06.450640Z",
     "shell.execute_reply": "2022-02-05T23:04:06.450214Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intercept  -12.493436061176814\n",
      "            Age  Experience    Income    Family     CCAvg  Mortgage  \\\n",
      "coeff -0.037685    0.039202  0.058844  0.612251  0.240489  0.001012   \n",
      "\n",
      "       Securities_Account  CD_Account    Online  CreditCard  \\\n",
      "coeff            -1.01428    3.649097 -0.678306   -0.958283   \n",
      "\n",
      "       Education_Graduate  Education_Advanced/Professional  \n",
      "coeff            4.202148                         4.355761  \n",
      "\n",
      "AIC -709.1524769205962\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_181/1278374867.py:8: FutureWarning: The `inplace` parameter in pandas.Categorical.rename_categories is deprecated and will be removed in a future version. Removing unused categories will always return a new Categorical object.\n",
      "  bank_df.Education.cat.rename_categories(new_categories, inplace=True)\n"
     ]
    }
   ],
   "source": [
    "bank_df = dmba.load_data('UniversalBank.csv')\n",
    "bank_df.drop(columns=['ID', 'ZIP Code'], inplace=True)\n",
    "bank_df.columns = [c.replace(' ', '_') for c in bank_df.columns]\n",
    "\n",
    "# Treat education as categorical, convert to dummy variables\n",
    "bank_df['Education'] = bank_df['Education'].astype('category')\n",
    "new_categories = {1: 'Undergrad', 2: 'Graduate', 3: 'Advanced/Professional'}\n",
    "bank_df.Education.cat.rename_categories(new_categories, inplace=True)\n",
    "bank_df = pd.get_dummies(bank_df, prefix_sep='_', drop_first=True)\n",
    "\n",
    "y = bank_df['Personal_Loan']\n",
    "X = bank_df.drop(columns=['Personal_Loan'])\n",
    "\n",
    "# partition data\n",
    "train_X, valid_X, train_y, valid_y = train_test_split(X, y, test_size=0.4, random_state=1)\n",
    "\n",
    "# fit a logistic regression (set penalty=l2 and C=1e42 to avoid regularization)\n",
    "logit_reg = LogisticRegression(penalty=\"l2\", C=1e42, solver='liblinear')\n",
    "logit_reg.fit(train_X, train_y)\n",
    "\n",
    "print('intercept ', logit_reg.intercept_[0])\n",
    "print(pd.DataFrame({'coeff': logit_reg.coef_[0]}, index=X.columns).transpose())\n",
    "print()\n",
    "print('AIC', AIC_score(valid_y, logit_reg.predict(valid_X), df = len(train_X.columns) + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:06.455834Z",
     "iopub.status.busy": "2022-02-05T23:04:06.455369Z",
     "iopub.status.idle": "2022-02-05T23:04:06.457277Z",
     "shell.execute_reply": "2022-02-05T23:04:06.457633Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set to True to calculate logistic regression using Income only\n",
    "if False:\n",
    "    predictors = ['Income']\n",
    "    outcome = 'Personal_Loan'\n",
    "\n",
    "    y = bank_df[outcome]\n",
    "    X = bank_df[predictors]\n",
    "\n",
    "    # partition data\n",
    "    train_X, valid_X, train_y, valid_y = train_test_split(X, y, test_size=0.4, random_state=1)\n",
    "\n",
    "    # fit a logistic regression (set penalty=l2 and C=1e42 to avoid regularization)\n",
    "    logit_reg_income = LogisticRegression(penalty=\"l2\", C=1e42, solver='liblinear')\n",
    "    logit_reg_income.fit(train_X, train_y)\n",
    "\n",
    "    print('intercept ', logit_reg_income.intercept_[0])\n",
    "    print(pd.DataFrame({'coefficient': logit_reg_income.coef_[0]}, index=X.columns).transpose())\n",
    "    print()\n",
    "    print('AIC', AIC_score(valid_y, logit_reg_income.predict(valid_X), df = len(train_X.columns) + 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 10.3\n",
    "Predict to get the probabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:06.462136Z",
     "iopub.status.busy": "2022-02-05T23:04:06.459579Z",
     "iopub.status.idle": "2022-02-05T23:04:06.470937Z",
     "shell.execute_reply": "2022-02-05T23:04:06.470621Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      actual      p(0)      p(1)  predicted\n",
      "2764       0  0.976092  0.023908          0\n",
      "932        0  0.331000  0.669000          1\n",
      "2721       1  0.031430  0.968570          1\n",
      "702        1  0.985893  0.014107          0\n"
     ]
    }
   ],
   "source": [
    "logit_reg_pred = logit_reg.predict(valid_X)\n",
    "logit_reg_proba = logit_reg.predict_proba(valid_X)\n",
    "logit_result = pd.DataFrame({'actual': valid_y, \n",
    "                             'p(0)': [p[0] for p in logit_reg_proba],\n",
    "                             'p(1)': [p[1] for p in logit_reg_proba],\n",
    "                             'predicted': logit_reg_pred })\n",
    "\n",
    "# display four different cases\n",
    "interestingCases = [2764, 932, 2721, 702]\n",
    "print(logit_result.loc[interestingCases])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 10.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:06.475005Z",
     "iopub.status.busy": "2022-02-05T23:04:06.472946Z",
     "iopub.status.idle": "2022-02-05T23:04:06.484664Z",
     "shell.execute_reply": "2022-02-05T23:04:06.484228Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Confusion Matrix (Accuracy 0.9600)\n",
      "\n",
      "       Prediction\n",
      "Actual    0    1\n",
      "     0 2683   30\n",
      "     1   90  197\n",
      "Confusion Matrix (Accuracy 0.9595)\n",
      "\n",
      "       Prediction\n",
      "Actual    0    1\n",
      "     0 1791   16\n",
      "     1   65  128\n"
     ]
    }
   ],
   "source": [
    "classificationSummary(train_y, logit_reg.predict(train_X))\n",
    "classificationSummary(valid_y, logit_reg.predict(valid_X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 10.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:06.488611Z",
     "iopub.status.busy": "2022-02-05T23:04:06.486527Z",
     "iopub.status.idle": "2022-02-05T23:04:06.816580Z",
     "shell.execute_reply": "2022-02-05T23:04:06.816931Z"
    }
   },
   "outputs": [],
   "source": [
    "df = logit_result.sort_values(by=['p(1)'], ascending=False)\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 4))\n",
    "\n",
    "gainsChart(df.actual, ax=axes[0])\n",
    "liftChart(df.actual, title=False, ax=axes[1])\n",
    "    \n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 10.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:06.819070Z",
     "iopub.status.busy": "2022-02-05T23:04:06.818583Z",
     "iopub.status.idle": "2022-02-05T23:04:06.874140Z",
     "shell.execute_reply": "2022-02-05T23:04:06.874394Z"
    }
   },
   "outputs": [],
   "source": [
    "delays_df = dmba.load_data('FlightDelays.csv')\n",
    "# Create an indicator variable\n",
    "delays_df['isDelayed'] = [1 if status == 'delayed' else 0 for status in delays_df['Flight Status']]\n",
    "\n",
    "# group information by day of week and determine the average delay\n",
    "averageDelay = delays_df.groupby(['DAY_WEEK'])['isDelayed'].mean()\n",
    "\n",
    "# create a bar chart\n",
    "ax = averageDelay.plot.bar(color='C0')\n",
    "ax.set_xlabel('Day of week')\n",
    "ax.set_ylabel('Average Delay')\n",
    "_ = ax.set_xticklabels(['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:06.876695Z",
     "iopub.status.busy": "2022-02-05T23:04:06.876345Z",
     "iopub.status.idle": "2022-02-05T23:04:07.319498Z",
     "shell.execute_reply": "2022-02-05T23:04:07.319788Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "delays_df = dmba.load_data('FlightDelays.csv')\n",
    "# Create an indicator variable\n",
    "delays_df['isDelayed'] = [1 if status == 'delayed' else 0 \n",
    "                          for status in delays_df['Flight Status']]\n",
    "\n",
    "def createGraph(group, xlabel, axis):\n",
    "    groupAverage = delays_df.groupby([group])['isDelayed'].mean()\n",
    "    if group == 'DAY_WEEK': # rotate so that display starts on Sunday\n",
    "        groupAverage = groupAverage.reindex(index=np.roll(groupAverage.index,1))\n",
    "        groupAverage.index = ['Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat']\n",
    "    ax = groupAverage.plot.bar(color='C0', ax=axis)\n",
    "    ax.set_ylabel('Average Delay')\n",
    "    ax.set_xlabel(xlabel)\n",
    "    return ax\n",
    "\n",
    "def graphDepartureTime(xlabel, axis):\n",
    "    temp_df = pd.DataFrame({'CRS_DEP_TIME': delays_df['CRS_DEP_TIME'] // 100, \n",
    "                            'isDelayed': delays_df['isDelayed']})\n",
    "    groupAverage = temp_df.groupby(['CRS_DEP_TIME'])['isDelayed'].mean()\n",
    "    ax = groupAverage.plot.bar(color='C0', ax=axis)\n",
    "    ax.set_xlabel(xlabel); ax.set_ylabel('Average Delay')\n",
    "    \n",
    "fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(10, 9))\n",
    "\n",
    "createGraph('DAY_WEEK', 'Day of week', axis=axes[0][0])\n",
    "createGraph('DEST', 'Destination', axis=axes[0][1])\n",
    "graphDepartureTime('Departure time', axis=axes[1][0])\n",
    "createGraph('CARRIER', 'Carrier', axis=axes[1][1])\n",
    "createGraph('ORIGIN', 'Origin', axis=axes[2][0])\n",
    "createGraph('Weather', 'Weather', axis=axes[2][1])\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 10.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:07.321959Z",
     "iopub.status.busy": "2022-02-05T23:04:07.321497Z",
     "iopub.status.idle": "2022-02-05T23:04:07.630749Z",
     "shell.execute_reply": "2022-02-05T23:04:07.631889Z"
    }
   },
   "outputs": [],
   "source": [
    "agg = delays_df.groupby(['ORIGIN', 'DAY_WEEK', 'CARRIER']).isDelayed.mean()\n",
    "agg = agg.reset_index()\n",
    "\n",
    "# Define the layout of the graph\n",
    "height_ratios = []\n",
    "for i, origin in enumerate(sorted(delays_df.ORIGIN.unique())):\n",
    "    height_ratios.append(len(agg[agg.ORIGIN == origin].CARRIER.unique()))\n",
    "gridspec_kw = {'height_ratios': height_ratios, 'width_ratios': [15, 1]}\n",
    "fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(10, 6), \n",
    "                         gridspec_kw = gridspec_kw)\n",
    "axes[0, 1].axis('off')\n",
    "axes[2, 1].axis('off')\n",
    "\n",
    "maxIsDelay = agg.isDelayed.max()\n",
    "for i, origin in enumerate(sorted(delays_df.ORIGIN.unique())):\n",
    "    data = pd.pivot_table(agg[agg.ORIGIN == origin], values='isDelayed', aggfunc=np.sum, \n",
    "                          index=['CARRIER'], columns=['DAY_WEEK'])\n",
    "    data = data[[7, 1, 2, 3, 4, 5, 6]]  # Shift last columns to first\n",
    "    ax = sns.heatmap(data, ax=axes[i][0], vmin=0, vmax=maxIsDelay, \n",
    "                     cbar_ax=axes[1][1], cmap=sns.light_palette(\"navy\"))\n",
    "    ax.set_xticklabels(['Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat'])\n",
    "    if i != 2: \n",
    "        ax.get_xaxis().set_visible(False)\n",
    "    ax.set_ylabel('Airport ' + origin)\n",
    "    \n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 10.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:07.634827Z",
     "iopub.status.busy": "2022-02-05T23:04:07.634494Z",
     "iopub.status.idle": "2022-02-05T23:04:07.668342Z",
     "shell.execute_reply": "2022-02-05T23:04:07.668037Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intercept  -1.2191013737953593\n",
      "       Weather  DAY_WEEK_2  DAY_WEEK_3  DAY_WEEK_4  DAY_WEEK_5  DAY_WEEK_6  DAY_WEEK_7  \\\n",
      "coeff    9.325      -0.598      -0.705      -0.799      -0.296      -1.129      -0.135   \n",
      "\n",
      "       CRS_DEP_TIME_7  CRS_DEP_TIME_8  CRS_DEP_TIME_9  CRS_DEP_TIME_10  CRS_DEP_TIME_11  \\\n",
      "coeff           0.631           0.382          -0.365            0.337            0.078   \n",
      "\n",
      "       CRS_DEP_TIME_12  CRS_DEP_TIME_13  CRS_DEP_TIME_14  CRS_DEP_TIME_15  CRS_DEP_TIME_16  \\\n",
      "coeff            0.399            0.175            0.202            1.265            0.628   \n",
      "\n",
      "       CRS_DEP_TIME_17  CRS_DEP_TIME_18  CRS_DEP_TIME_19  CRS_DEP_TIME_20  CRS_DEP_TIME_21  \\\n",
      "coeff            1.093            0.285            1.655            1.023            1.077   \n",
      "\n",
      "       ORIGIN_DCA  ORIGIN_IAD  DEST_JFK  DEST_LGA  CARRIER_DH  CARRIER_DL  CARRIER_MQ  \\\n",
      "coeff       -0.01      -0.134    -0.524    -0.546       0.352      -0.685       0.743   \n",
      "\n",
      "       CARRIER_OH  CARRIER_RU  CARRIER_UA  CARRIER_US  \n",
      "coeff      -0.711      -0.194       0.315      -0.971  \n",
      "\n",
      "AIC 1004.5346225948085\n"
     ]
    }
   ],
   "source": [
    "delays_df = dmba.load_data('FlightDelays.csv')\n",
    "# Create an indicator variable\n",
    "delays_df['isDelayed'] = [1 if status == 'delayed' else 0 \n",
    "                          for status in delays_df['Flight Status']]\n",
    "\n",
    "# convert to categorical\n",
    "delays_df.DAY_WEEK = delays_df.DAY_WEEK.astype('category')\n",
    "\n",
    "# create hourly bins departure time \n",
    "delays_df.CRS_DEP_TIME = [round(t / 100) for t in delays_df.CRS_DEP_TIME]\n",
    "delays_df.CRS_DEP_TIME = delays_df.CRS_DEP_TIME.astype('category')\n",
    "\n",
    "predictors = ['DAY_WEEK', 'CRS_DEP_TIME', 'ORIGIN', 'DEST', 'CARRIER', 'Weather']\n",
    "outcome = 'isDelayed'\n",
    "\n",
    "X = pd.get_dummies(delays_df[predictors], drop_first=True)\n",
    "y = delays_df[outcome]\n",
    "classes = ['ontime', 'delayed']\n",
    "\n",
    "# split into training and validation\n",
    "train_X, valid_X, train_y, valid_y = train_test_split(X, y, test_size=0.4, \n",
    "                                                      random_state=1)\n",
    "\n",
    "logit_full = LogisticRegression(penalty=\"l2\", C=1e42, solver='liblinear')\n",
    "logit_full.fit(train_X, train_y)\n",
    "\n",
    "pd.set_option('display.width', 95)\n",
    "pd.set_option('display.precision',3)\n",
    "pd.set_option('display.max_columns', 33)\n",
    "print('intercept ', logit_full.intercept_[0])\n",
    "\n",
    "print(pd.DataFrame({'coeff': logit_full.coef_[0]}, index=X.columns).transpose())\n",
    "print()\n",
    "print('AIC', AIC_score(valid_y, logit_full.predict(valid_X), df=len(train_X.columns) + 1))\n",
    "pd.reset_option('display.width')\n",
    "pd.reset_option('display.precision')\n",
    "pd.reset_option('display.max_columns')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 10.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:07.673316Z",
     "iopub.status.busy": "2022-02-05T23:04:07.672952Z",
     "iopub.status.idle": "2022-02-05T23:04:07.771417Z",
     "shell.execute_reply": "2022-02-05T23:04:07.771958Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Confusion Matrix (Accuracy 0.8309)\n",
      "\n",
      "        Prediction\n",
      " Actual  ontime delayed\n",
      " ontime     705       9\n",
      "delayed     140      27\n"
     ]
    }
   ],
   "source": [
    "logit_reg_pred = logit_full.predict_proba(valid_X)\n",
    "full_result = pd.DataFrame({'actual': valid_y, \n",
    "                            'p(0)': [p[0] for p in logit_reg_pred],\n",
    "                            'p(1)': [p[1] for p in logit_reg_pred],\n",
    "                            'predicted': logit_full.predict(valid_X)})\n",
    "full_result = full_result.sort_values(by=['p(1)'], ascending=False)\n",
    "\n",
    "# confusion matrix\n",
    "classificationSummary(full_result.actual, full_result.predicted, class_names=classes)\n",
    "\n",
    "gainsChart(full_result.actual, figsize=[5, 5])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 10.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:07.779588Z",
     "iopub.status.busy": "2022-02-05T23:04:07.779210Z",
     "iopub.status.idle": "2022-02-05T23:04:07.872247Z",
     "shell.execute_reply": "2022-02-05T23:04:07.871984Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "regularization [2.7825594]\n",
      "intercept  -2.287363183029368\n",
      "        Sun_Mon   Weather  CARRIER_CO_MQ_DH_RU   MORNING      NOON   AFTER2P   EVENING\n",
      "coeff  0.577972  4.977926             1.298819 -0.583306 -0.665845 -0.055188  0.560893\n",
      "AIC 934.6153607819033\n",
      "Confusion Matrix (Accuracy 0.8343)\n",
      "\n",
      "        Prediction\n",
      " Actual  ontime delayed\n",
      " ontime     711       3\n",
      "delayed     143      24\n"
     ]
    }
   ],
   "source": [
    "delays_df = dmba.load_data('FlightDelays.csv')\n",
    "delays_df['isDelayed'] = [1 if status == 'delayed' else 0 \n",
    "                          for status in delays_df['Flight Status']]\n",
    "delays_df['CRS_DEP_TIME'] = [round(t / 100) for t in delays_df['CRS_DEP_TIME']]\n",
    "delays_red_df = pd.DataFrame({\n",
    "    'Sun_Mon' : [1 if d in (1, 7) else 0 for d in delays_df.DAY_WEEK],\n",
    "    'Weather' : delays_df.Weather,\n",
    "    'CARRIER_CO_MQ_DH_RU' : [1 if d in (\"CO\", \"MQ\", \"DH\", \"RU\") else 0 \n",
    "                             for d in delays_df.CARRIER],\n",
    "    'MORNING' : [1 if d in (6, 7, 8, 9) else 0 for d in delays_df.CRS_DEP_TIME],\n",
    "    'NOON' : [1 if d in (10, 11, 12, 13) else 0 for d in delays_df.CRS_DEP_TIME],\n",
    "    'AFTER2P' : [1 if d in (14, 15, 16, 17, 18) else 0 for d in delays_df.CRS_DEP_TIME],\n",
    "    'EVENING' : [1 if d in (19, 20) else 0 for d in delays_df.CRS_DEP_TIME],\n",
    "    'isDelayed' : [1 if status == 'delayed' else 0 for status in delays_df['Flight Status']],\n",
    "})\n",
    "\n",
    "X = delays_red_df.drop(columns=['isDelayed'])\n",
    "y = delays_red_df['isDelayed']\n",
    "classes = ['ontime', 'delayed']\n",
    "\n",
    "# split into training and validation\n",
    "train_X, valid_X, train_y, valid_y = train_test_split(X, y, test_size=0.4, \n",
    "                                                      random_state=1)\n",
    "\n",
    "logit_red = LogisticRegressionCV(penalty=\"l1\", solver='liblinear', cv=5)\n",
    "logit_red.fit(train_X, train_y)\n",
    "\n",
    "pd.set_option('display.width', 100)\n",
    "print('regularization', logit_red.C_)\n",
    "print('intercept ', logit_red.intercept_[0])\n",
    "print(pd.DataFrame({'coeff': logit_red.coef_[0]}, index=X.columns).transpose())\n",
    "pd.reset_option('display.width')\n",
    "print('AIC', AIC_score(valid_y, logit_red.predict(valid_X), df=len(train_X.columns) + 1))\n",
    "\n",
    "\n",
    "# confusion matrix\n",
    "classificationSummary(valid_y, logit_red.predict(valid_X), class_names=classes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:07.876908Z",
     "iopub.status.busy": "2022-02-05T23:04:07.876584Z",
     "iopub.status.idle": "2022-02-05T23:04:07.966855Z",
     "shell.execute_reply": "2022-02-05T23:04:07.966549Z"
    }
   },
   "outputs": [],
   "source": [
    "logit_reg_proba = logit_red.predict_proba(valid_X)\n",
    "red_result = pd.DataFrame({'actual': valid_y, \n",
    "                            'p(0)': [p[0] for p in logit_reg_proba],\n",
    "                            'p(1)': [p[1] for p in logit_reg_proba],\n",
    "                            'predicted': logit_red.predict(valid_X),\n",
    "                          })\n",
    "red_result = red_result.sort_values(by=['p(1)'], ascending=False)\n",
    "\n",
    "ax = gainsChart(full_result.actual, label='Full model', color='C1', figsize=[5, 5])\n",
    "ax = gainsChart(red_result.actual, label='Reduced model', color='C0', ax=ax)\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 10.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:07.970421Z",
     "iopub.status.busy": "2022-02-05T23:04:07.970060Z",
     "iopub.status.idle": "2022-02-05T23:04:07.997448Z",
     "shell.execute_reply": "2022-02-05T23:04:07.997788Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:          Personal_Loan   No. Observations:                 3000\n",
      "Model:                            GLM   Df Residuals:                     2987\n",
      "Model Family:                Binomial   Df Model:                           12\n",
      "Link Function:                  Logit   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -340.15\n",
      "Date:                Sat, 05 Feb 2022   Deviance:                       680.30\n",
      "Time:                        23:04:07   Pearson chi2:                 8.10e+03\n",
      "No. Iterations:                     8   Pseudo R-squ. (CS):             0.3325\n",
      "Covariance Type:            nonrobust                                         \n",
      "===================================================================================================\n",
      "                                      coef    std err          z      P>|z|      [0.025      0.975]\n",
      "---------------------------------------------------------------------------------------------------\n",
      "const                             -12.5634      2.336     -5.377      0.000     -17.143      -7.984\n",
      "Age                                -0.0354      0.086     -0.412      0.680      -0.204       0.133\n",
      "Experience                          0.0369      0.086      0.431      0.666      -0.131       0.205\n",
      "Income                              0.0589      0.004     15.044      0.000       0.051       0.067\n",
      "Family                              0.6128      0.103      5.931      0.000       0.410       0.815\n",
      "CCAvg                               0.2408      0.060      4.032      0.000       0.124       0.358\n",
      "Mortgage                            0.0010      0.001      1.301      0.193      -0.001       0.003\n",
      "Securities_Account                 -1.0305      0.422     -2.443      0.015      -1.857      -0.204\n",
      "CD_Account                          3.6628      0.460      7.961      0.000       2.761       4.565\n",
      "Online                             -0.6794      0.216     -3.140      0.002      -1.103      -0.255\n",
      "CreditCard                         -0.9609      0.274     -3.507      0.000      -1.498      -0.424\n",
      "Education_Graduate                  4.2075      0.364     11.573      0.000       3.495       4.920\n",
      "Education_Advanced/Professional     4.3580      0.365     11.937      0.000       3.642       5.074\n",
      "===================================================================================================\n"
     ]
    }
   ],
   "source": [
    "# same initial preprocessing and creating dummies\n",
    "\n",
    "# add constant column\n",
    "bank_df = sm.add_constant(bank_df, prepend=True)\n",
    "\n",
    "y = bank_df['Personal_Loan']\n",
    "X = bank_df.drop(columns=['Personal_Loan'])\n",
    "\n",
    "# partition data\n",
    "train_X, valid_X, train_y, valid_y = train_test_split(X, y, test_size=0.4, random_state=1)\n",
    "\n",
    "# use GLM (general linear model) with the binomial family to fit a logistic regression\n",
    "logit_reg = sm.GLM(train_y, train_X, family=sm.families.Binomial())\n",
    "logit_result = logit_reg.fit()\n",
    "print(logit_result.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:08.003530Z",
     "iopub.status.busy": "2022-02-05T23:04:08.003135Z",
     "iopub.status.idle": "2022-02-05T23:04:08.479953Z",
     "shell.execute_reply": "2022-02-05T23:04:08.479662Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Nominal logistic regression\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  intercept [-0.09100315  0.9036454  -0.81264225]\n",
      "  coefficients [[ 0.51606685  0.3391015 ]\n",
      " [ 0.14900396  0.09543369]\n",
      " [-0.66507082 -0.43453518]]\n",
      "\n",
      "   actual  predicted      P(0)      P(1)      P(2)\n",
      "0       1          1  0.490649  0.498989  0.010362\n",
      "1       0          0  0.553461  0.441147  0.005392\n",
      "2       0          0  0.553461  0.441147  0.005392\n",
      "3       0          1  0.490649  0.498989  0.010362\n",
      "4       0          1  0.394192  0.578684  0.027124\n",
      "\n",
      "Ordinal logistic regression\n",
      "  theta [-1.06916285  2.77444326]\n",
      "  coefficients [-0.40112008 -0.25174207]\n",
      "\n",
      "   actual  predicted      P(0)      P(1)      P(2)\n",
      "0       1          1  0.496205  0.482514  0.021281\n",
      "1       0          0  0.558866  0.424510  0.016624\n",
      "2       0          0  0.558866  0.424510  0.016624\n",
      "3       0          1  0.496205  0.482514  0.021281\n",
      "4       0          1  0.397402  0.571145  0.031453\n"
     ]
    }
   ],
   "source": [
    "data = dmba.load_data('accidentsFull.csv')\n",
    "outcome = 'MAX_SEV_IR'\n",
    "predictors = ['ALCHL_I', 'WEATHER_R']\n",
    "\n",
    "y = data[outcome]\n",
    "X = data[predictors]\n",
    "train_X, train_y = X, y\n",
    "classes = sorted(y.unique())\n",
    "\n",
    "print('Nominal logistic regression')\n",
    "logit = LogisticRegression(penalty=\"l2\", solver='lbfgs', C=1e24, multi_class='multinomial')\n",
    "logit.fit(X, y)\n",
    "print('  intercept', logit.intercept_)\n",
    "print('  coefficients', logit.coef_)\n",
    "print()\n",
    "probs = logit.predict_proba(X)\n",
    "results = pd.DataFrame({\n",
    "    'actual': y, 'predicted': logit.predict(X),\n",
    "    'P(0)': [p[0] for p in probs],\n",
    "    'P(1)': [p[1] for p in probs],\n",
    "    'P(2)': [p[2] for p in probs],\n",
    "})\n",
    "print(results.head())\n",
    "print()\n",
    "\n",
    "\n",
    "print('Ordinal logistic regression')\n",
    "logit = LogisticIT(alpha=0)\n",
    "logit.fit(X, y)\n",
    "print('  theta', logit.theta_)\n",
    "print('  coefficients', logit.coef_)\n",
    "print()\n",
    "probs = logit.predict_proba(X)\n",
    "results = pd.DataFrame({\n",
    "    'actual': y, 'predicted': logit.predict(X),\n",
    "    'P(0)': [p[0] for p in probs],\n",
    "    'P(1)': [p[1] for p in probs],\n",
    "    'P(2)': [p[2] for p in probs],\n",
    "})\n",
    "print(results.head())\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nominal logistic regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:08.484641Z",
     "iopub.status.busy": "2022-02-05T23:04:08.484279Z",
     "iopub.status.idle": "2022-02-05T23:04:08.806404Z",
     "shell.execute_reply": "2022-02-05T23:04:08.806042Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "predictors ['ALCHL_I', 'WEATHER_R']\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intercept [-0.09100315  0.9036454  -0.81264225]\n",
      "coef [[ 0.51606685  0.3391015 ]\n",
      " [ 0.14900396  0.09543369]\n",
      " [-0.66507082 -0.43453518]]\n",
      "classes [0 1 2]\n",
      "   actual  predicted      P(a)      P(b)      P(c)\n",
      "0       1          1  0.490649  0.498989  0.010362\n",
      "1       0          0  0.553461  0.441147  0.005392\n",
      "2       0          0  0.553461  0.441147  0.005392\n",
      "3       0          1  0.490649  0.498989  0.010362\n",
      "4       0          1  0.394192  0.578684  0.027124\n"
     ]
    }
   ],
   "source": [
    "data = dmba.load_data('accidentsFull.csv')\n",
    "outcome = 'MAX_SEV_IR'\n",
    "predictors = ['ALCHL_I', 'WEATHER_R']\n",
    "print('predictors', predictors)\n",
    "print()\n",
    "y = data[outcome]\n",
    "X = data[predictors]\n",
    "train_X, train_y = X, y\n",
    "classes = sorted(y.unique())\n",
    "\n",
    "logit = LogisticRegression(penalty=\"l2\", solver='lbfgs', C=1e24, multi_class='multinomial')\n",
    "logit.fit(X, y)\n",
    "\n",
    "print('intercept', logit.intercept_)\n",
    "print('coef', logit.coef_)\n",
    "print('classes', logit.classes_)\n",
    "\n",
    "probs = logit.predict_proba(X)\n",
    "results = pd.DataFrame({\n",
    "    'actual': y,\n",
    "    'predicted': logit.predict(X),\n",
    "    'P(a)': [p[0] for p in probs],\n",
    "    'P(b)': [p[1] for p in probs],\n",
    "    'P(c)': [p[2] for p in probs],\n",
    "})\n",
    "print(results.head())\n",
    "# classificationSummary(y, results.predicted, class_names=classes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:08.814066Z",
     "iopub.status.busy": "2022-02-05T23:04:08.813675Z",
     "iopub.status.idle": "2022-02-05T23:04:09.282956Z",
     "shell.execute_reply": "2022-02-05T23:04:09.282652Z"
    }
   },
   "outputs": [],
   "source": [
    "def plotProbabilities(model):\n",
    "    n = 100\n",
    "    Xtest = pd.DataFrame({\n",
    "        'ALCHL_I': [(i % n) * 0.1 + (j // n) * 0.1 - 5 for i in range(n) for j in range(n)],\n",
    "        'WEATHER_R': [(i // n) * 0.1 + (j % n) * 0.1 - 5 for i in range(n) for j in range(n)],\n",
    "    })\n",
    "    probs = logit.predict_proba(Xtest[predictors])\n",
    "\n",
    "    results = pd.DataFrame({\n",
    "        'predicted': logit.predict(Xtest[predictors]),\n",
    "        'P(a)': [round(10 * p[0]) / 10 for p in probs],\n",
    "        'P(b)': [round(10 * p[1]) / 10 for p in probs],\n",
    "        'P(c)': [round(10 * p[2]) / 10 for p in probs],\n",
    "    })\n",
    "    mapped = pd.concat([results, Xtest], axis=1, sort=False)\n",
    "\n",
    "    fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 3))\n",
    "    for i, c in enumerate(['P(a)', 'P(b)', 'P(c)']):\n",
    "        ax = mapped.plot.scatter(x='ALCHL_I', y='WEATHER_R', c=c, title=c, colormap='cividis', ax=axes[i])\n",
    "        if i > 0: ax.set_ylabel('')\n",
    "plotProbabilities(logit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ordinal logistic regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:09.288367Z",
     "iopub.status.busy": "2022-02-05T23:04:09.287987Z",
     "iopub.status.idle": "2022-02-05T23:04:09.598942Z",
     "shell.execute_reply": "2022-02-05T23:04:09.598665Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theta [-1.06916285  2.77444326]\n",
      "coef [-0.40112008 -0.25174207]\n",
      "classes [0 1 2]\n",
      "\n",
      "     actual  predicted      P(a)      P(b)      P(c)\n",
      "0  1.000046          1  0.495504  0.482286  0.022849\n",
      "1  0.000133          0  0.559977  0.422456  0.015655\n",
      "2  0.000810          0  0.559380  0.425080  0.014003\n",
      "3  0.000306          1  0.493719  0.482300  0.021894\n",
      "4 -0.000340          1  0.397216  0.571625  0.031833\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "from mord import LogisticIT\n",
    "\n",
    "data = dmba.load_data('accidentsFull.csv')\n",
    "outcome = 'MAX_SEV_IR'\n",
    "predictors = ['ALCHL_I', 'WEATHER_R']\n",
    "y = data[outcome]\n",
    "X = data[predictors]\n",
    "X['ALCHL_I']\n",
    "train_X, train_y = X, y\n",
    "classes = sorted(y.unique())\n",
    "\n",
    "logit = LogisticIT(alpha=0)\n",
    "logit.fit(X, y)\n",
    "\n",
    "print('theta', logit.theta_)\n",
    "print('coef', logit.coef_)\n",
    "print('classes', logit.classes_)\n",
    "print()\n",
    "\n",
    "\n",
    "probs = logit.predict_proba(X)\n",
    "\n",
    "results = pd.DataFrame({\n",
    "    'actual': [yi + random.gauss(0, 0.001) for yi in y],\n",
    "    'predicted': logit.predict(X),\n",
    "    'P(a)': [p[0] + random.gauss(0, 0.001) for p in probs],\n",
    "    'P(b)': [p[1] + random.gauss(0, 0.001) for p in probs],\n",
    "    'P(c)': [p[2] + random.gauss(0, 0.001) for p in probs],\n",
    "})\n",
    "print(results.head())\n",
    "\n",
    "# classificationSummary(y, results.predicted, class_names=classes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-05T23:04:09.606123Z",
     "iopub.status.busy": "2022-02-05T23:04:09.605708Z",
     "iopub.status.idle": "2022-02-05T23:04:10.121869Z",
     "shell.execute_reply": "2022-02-05T23:04:10.121586Z"
    }
   },
   "outputs": [],
   "source": [
    "plotProbabilities(logit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
