05-logistic-regression/3-assessing-coefficients-matrix-multiplication.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5e9e49bd",
   "metadata": {},
   "source": [
    "# Evaluate Coefficients of Logistic Regression\n",
    "\n",
    "[Machine Learning Interpretability course](https://www.trainindata.com/p/machine-learning-interpretability)\n",
    "\n",
    "In this notebook, we will evaluate the coefficients of the logistic regression. We will first use `statsmodels` and then `sklearn`.\n",
    "\n",
    "`statsmodels` returns the statistics and p-values out of the box, utilizing MLE. For `sklearn`, we will obtain these parameters manually, also using the derivations from MLE.\n",
    "\n",
    "If the coefficients are significant, then the contribution of the feature towards the probability is meaningful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b619bd11",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "750e70b0",
   "metadata": {},
   "source": [
    "### Load data\n",
    "\n",
    "To obtain the data, check the folder `prepare-data` in this repo or section 2 in the course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "050acc05",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load titanic dataset\n",
    "\n",
    "df = pd.read_csv('../titanic.csv')\n",
    "\n",
    "# split data\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    df.drop(\"survived\", axis=1), \n",
    "    df[\"survived\"],\n",
    "    test_size=0.15,\n",
    "    random_state=1,\n",
    ")\n",
    "\n",
    "# scale the variables\n",
    "scaler = StandardScaler().set_output(transform=\"pandas\")\n",
    "\n",
    "X_train = scaler.fit_transform(X_train)\n",
    "X_test = scaler.transform(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ada99f7",
   "metadata": {},
   "source": [
    "## Statsmodels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ec8ae23c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Our model needs an intercept so we add a column of 1s:\n",
    "\n",
    "X_train_sm = sm.add_constant(X_train)\n",
    "X_test_sm = sm.add_constant(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dfac6e9b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.469479\n",
      "         Iterations 6\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:               survived   No. Observations:                 1112\n",
      "Model:                          Logit   Df Residuals:                     1098\n",
      "Method:                           MLE   Df Model:                           13\n",
      "Date:                Tue, 14 Nov 2023   Pseudo R-squ.:                  0.2923\n",
      "Time:                        10:56:39   Log-Likelihood:                -522.06\n",
      "converged:                       True   LL-Null:                       -737.67\n",
      "Covariance Type:            nonrobust   LLR p-value:                 5.606e-84\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -0.6505      0.079     -8.267      0.000      -0.805      -0.496\n",
      "pclass        -0.5174      0.118     -4.389      0.000      -0.748      -0.286\n",
      "sibsp         -0.2268      0.095     -2.383      0.017      -0.413      -0.040\n",
      "parch          0.0542      0.082      0.660      0.509      -0.107       0.215\n",
      "sex_female     1.1999      0.082     14.620      0.000       1.039       1.361\n",
      "embarked_S    -0.0085      0.125     -0.068      0.946      -0.254       0.237\n",
      "embarked_C     0.2340      0.127      1.836      0.066      -0.016       0.484\n",
      "cabin_B       -0.2198      0.149     -1.479      0.139      -0.511       0.072\n",
      "cabin_C       -0.3765      0.170     -2.210      0.027      -0.710      -0.043\n",
      "cabin_E       -0.0252      0.129     -0.195      0.846      -0.279       0.228\n",
      "cabin_D       -0.1982      0.130     -1.519      0.129      -0.454       0.058\n",
      "cabin_A       -0.1051      0.100     -1.055      0.291      -0.300       0.090\n",
      "cabin_M       -0.6804      0.236     -2.878      0.004      -1.144      -0.217\n",
      "cabin_G       -0.1436      0.073     -1.958      0.050      -0.287       0.000\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "# Fit logistic regression model\n",
    "\n",
    "logit_mod = sm.Logit(y_train, X_train_sm)\n",
    "\n",
    "logit_res = logit_mod.fit()\n",
    "\n",
    "print(logit_res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d666ad47",
   "metadata": {},
   "source": [
    "With statsmodels, the statistical testing of the coefficients comes out of the box.\n",
    "\n",
    "We see that some coefficients are not statistically significant, so we could remove those features from the model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adfae3e9",
   "metadata": {},
   "source": [
    "## Scikit-learn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "89b7867b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {color: black;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>LogisticRegression(penalty=None, random_state=1)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" checked><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">LogisticRegression</label><div class=\"sk-toggleable__content\"><pre>LogisticRegression(penalty=None, random_state=1)</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "LogisticRegression(penalty=None, random_state=1)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# fit model\n",
    "\n",
    "logit = LogisticRegression(penalty=None, random_state=1)\n",
    "\n",
    "logit.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8231597a",
   "metadata": {},
   "source": [
    "### Empirically\n",
    "\n",
    "We follow the derivations shown [here](https://web.stanford.edu/class/archive/stats/stats200/stats200.1172/Lecture26.pdf) and the Python implementation shown [here](https://stats.stackexchange.com/questions/89484/how-to-compute-the-standard-errors-of-a-logistic-regressions-coefficients)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "489ee265",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1], dtype=int64)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the target classes\n",
    "\n",
    "logit.classes_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b0b16b67",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.89651014, 0.10348986],\n",
       "       [0.76111387, 0.23888613],\n",
       "       [0.55621831, 0.44378169],\n",
       "       ...,\n",
       "       [0.93911299, 0.06088701],\n",
       "       [0.71414115, 0.28585885],\n",
       "       [0.41312738, 0.58687262]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# matrix of predicted class probabilities\n",
    "\n",
    "pred_train = logit.predict_proba(X_train)\n",
    "\n",
    "pred_train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "32111fe9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.83427143, -0.47413915, ..., -0.12827074,\n",
       "         0.5399273 , -0.06720649],\n",
       "       [ 1.        , -0.36742524,  0.46905   , ..., -0.12827074,\n",
       "         0.5399273 , -0.06720649],\n",
       "       [ 1.        ,  0.83427143,  2.35542828, ..., -0.12827074,\n",
       "         0.5399273 , -0.06720649],\n",
       "       ...,\n",
       "       [ 1.        ,  0.83427143,  2.35542828, ..., -0.12827074,\n",
       "         0.5399273 , -0.06720649],\n",
       "       [ 1.        , -1.56912191, -0.47413915, ..., -0.12827074,\n",
       "         0.5399273 , -0.06720649],\n",
       "       [ 1.        ,  0.83427143, -0.47413915, ..., -0.12827074,\n",
       "         0.5399273 , -0.06720649]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# add column of 1's at the beginning of X_train\n",
    "\n",
    "X_train_c = np.hstack([np.ones((X_train.shape[0], 1)), X_train])\n",
    "\n",
    "X_train_c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "89268191",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.09277971, 0.        , 0.        , ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       [0.        , 0.18181955, 0.        , ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       [0.        , 0.        , 0.2468395 , ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       ...,\n",
       "       [0.        , 0.        , 0.        , ..., 0.05717978, 0.        ,\n",
       "        0.        ],\n",
       "       [0.        , 0.        , 0.        , ..., 0.        , 0.20414357,\n",
       "        0.        ],\n",
       "       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,\n",
       "        0.24245315]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create matrix of 0's, where diagonal contains each \n",
    "# predicted observation's variance.\n",
    "\n",
    "W = np.diagflat(np.product(pred_train, axis=1))\n",
    "\n",
    "W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f2c97b23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0927797088771804"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The variance is the product of the probabilities\n",
    "\n",
    "# For the first observation:\n",
    "0.89651014 * 0.10348986"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "992915de",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 6.19148576e-03,  9.50979865e-04,  4.04960853e-04,\n",
       "        -1.28885132e-04, -9.73906357e-04,  3.31124990e-04,\n",
       "         1.25713235e-04,  7.24502525e-04,  7.11303716e-04,\n",
       "         5.67242165e-04,  4.37754332e-04,  1.75871235e-04,\n",
       "         8.61118068e-04,  3.04792119e-05],\n",
       "       [ 9.50979865e-04,  1.38932240e-02, -7.85293102e-04,\n",
       "        -5.80503656e-04, -5.68088469e-04,  3.37760169e-03,\n",
       "         3.12774306e-03,  5.16264693e-03,  6.33396731e-03,\n",
       "         3.41794219e-03,  3.88078452e-03,  2.94650855e-03,\n",
       "        -3.87665730e-04, -2.43197064e-04],\n",
       "       [ 4.04960853e-04, -7.85293102e-04,  9.05920372e-03,\n",
       "        -2.29035784e-03, -1.17425375e-03, -1.26791128e-03,\n",
       "        -7.73044204e-04, -3.44877875e-04, -9.14102587e-04,\n",
       "        -1.31800356e-04, -2.24792355e-04, -5.39609035e-05,\n",
       "        -1.18947160e-04,  2.45175335e-04],\n",
       "       [-1.28885132e-04, -5.80503656e-04, -2.29035784e-03,\n",
       "         6.72920647e-03, -1.26879216e-03, -1.42280469e-03,\n",
       "        -1.32916565e-03, -1.39077749e-04,  2.96075240e-04,\n",
       "         1.11698860e-04,  2.22609036e-04,  9.36236925e-05,\n",
       "         7.17364723e-04, -1.58282602e-04],\n",
       "       [-9.73906357e-04, -5.68088469e-04, -1.17425375e-03,\n",
       "        -1.26879216e-03,  6.73616917e-03,  2.03881432e-03,\n",
       "         2.20928078e-03, -8.48356610e-04, -1.21724108e-03,\n",
       "        -3.82684234e-04, -6.67545187e-04, -2.02431098e-04,\n",
       "        -2.04216967e-03, -6.52061191e-04],\n",
       "       [ 3.31124990e-04,  3.37760169e-03, -1.26791128e-03,\n",
       "        -1.42280469e-03,  2.03881432e-03,  1.57473480e-02,\n",
       "         1.25575576e-02,  8.89841420e-04,  1.25642459e-03,\n",
       "         4.57232115e-04,  5.63365097e-04,  4.84621277e-04,\n",
       "         1.68037849e-04, -4.03468109e-04],\n",
       "       [ 1.25713235e-04,  3.12774306e-03, -7.73044204e-04,\n",
       "        -1.32916565e-03,  2.20928078e-03,  1.25575576e-02,\n",
       "         1.62462734e-02, -3.64311995e-04, -1.78656788e-04,\n",
       "         2.23563510e-04, -4.46161755e-04, -2.51488897e-05,\n",
       "        -5.57396840e-04, -2.82215465e-04],\n",
       "       [ 7.24502525e-04,  5.16264693e-03, -3.44877875e-04,\n",
       "        -1.39077749e-04, -8.48356610e-04,  8.89841420e-04,\n",
       "        -3.64311995e-04,  2.20984087e-02,  1.98719938e-02,\n",
       "         1.30344179e-02,  1.39830020e-02,  9.53565539e-03,\n",
       "         2.77620091e-02,  4.37936549e-03],\n",
       "       [ 7.11303716e-04,  6.33396731e-03, -9.14102587e-04,\n",
       "         2.96075240e-04, -1.21724108e-03,  1.25642459e-03,\n",
       "        -1.78656788e-04,  1.98719938e-02,  2.90270997e-02,\n",
       "         1.57875514e-02,  1.69365557e-02,  1.15406090e-02,\n",
       "         3.36614046e-02,  5.29180290e-03],\n",
       "       [ 5.67242165e-04,  3.41794219e-03, -1.31800356e-04,\n",
       "         1.11698860e-04, -3.82684234e-04,  4.57232115e-04,\n",
       "         2.23563510e-04,  1.30344179e-02,  1.57875514e-02,\n",
       "         1.67431229e-02,  1.11563524e-02,  7.62586347e-03,\n",
       "         2.27343334e-02,  3.59461102e-03],\n",
       "       [ 4.37754332e-04,  3.88078452e-03, -2.24792355e-04,\n",
       "         2.22609036e-04, -6.67545187e-04,  5.63365097e-04,\n",
       "        -4.46161755e-04,  1.39830020e-02,  1.69365557e-02,\n",
       "         1.11563524e-02,  1.70221679e-02,  8.14573452e-03,\n",
       "         2.40570826e-02,  3.78588950e-03],\n",
       "       [ 1.75871235e-04,  2.94650855e-03, -5.39609035e-05,\n",
       "         9.36236925e-05, -2.02431098e-04,  4.84621277e-04,\n",
       "        -2.51488897e-05,  9.53565539e-03,  1.15406090e-02,\n",
       "         7.62586347e-03,  8.14573452e-03,  9.92010997e-03,\n",
       "         1.61914449e-02,  2.53667193e-03],\n",
       "       [ 8.61118068e-04, -3.87665730e-04, -1.18947160e-04,\n",
       "         7.17364723e-04, -2.04216967e-03,  1.68037849e-04,\n",
       "        -5.57396840e-04,  2.77620091e-02,  3.36614046e-02,\n",
       "         2.27343334e-02,  2.40570826e-02,  1.61914449e-02,\n",
       "         5.58851472e-02,  8.77693257e-03],\n",
       "       [ 3.04792119e-05, -2.43197064e-04,  2.45175335e-04,\n",
       "        -1.58282602e-04, -6.52061191e-04, -4.03468109e-04,\n",
       "        -2.82215465e-04,  4.37936549e-03,  5.29180290e-03,\n",
       "         3.59461102e-03,  3.78588950e-03,  2.53667193e-03,\n",
       "         8.77693257e-03,  5.38177465e-03]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Covariance matrix, it is derived empirically.\n",
    "\n",
    "# (XT W X)−1\n",
    "\n",
    "cov_m = np.linalg.inv(np.dot(np.dot(X_train_c.T, W), X_train_c))\n",
    "\n",
    "cov_m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "94a88ef2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.07868599, 0.11786952, 0.09517985, 0.08203174, 0.08207417,\n",
       "       0.12548844, 0.12746087, 0.14865534, 0.17037341, 0.12939522,\n",
       "       0.13046903, 0.09959975, 0.2364004 , 0.07336058])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Standard errors\n",
    "\n",
    "se = np.sqrt(np.diag(cov_m))\n",
    "\n",
    "se"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "048f360d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['intercept', 'pclass', 'sibsp', 'parch', 'sex_female',\n",
       "       'embarked_S', 'embarked_C', 'cabin_B', 'cabin_C', 'cabin_E',\n",
       "       'cabin_D', 'cabin_A', 'cabin_M', 'cabin_G'], dtype=object)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# capture name of features\n",
    "\n",
    "features = np.insert(logit.feature_names_in_, 0, \"intercept\")\n",
    "\n",
    "features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8c9e5d36",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.65051925, -0.51737413, -0.22689809,  0.05420292,  1.19989943,\n",
       "       -0.00845913,  0.23404885, -0.21982215, -0.37650021, -0.02520391,\n",
       "       -0.19821006, -0.1050838 , -0.68045168, -0.14362905])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# capture coefficients of the logit including the intercept\n",
    "\n",
    "coefficients = np.insert(logit.coef_[0], 0, logit.intercept_[0])\n",
    "\n",
    "coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "105e944e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-8.26728127, -4.38938008, -2.38388788,  0.66075542, 14.61969636,\n",
       "       -0.06740967,  1.83624077, -1.47873705, -2.20985305, -0.19478241,\n",
       "       -1.5192116 , -1.05506088, -2.8783864 , -1.9578506 ])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Wald's statistic (coefficient / s.e.)\n",
    "\n",
    "z_scores = (coefficients / np.sqrt(np.diag(cov_m)))\n",
    "\n",
    "z_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "4cd51db8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.37048362e-16, 1.13674245e-05, 1.71308208e-02, 5.08769178e-01,\n",
       "       2.10336291e-48, 9.46255570e-01, 6.63220556e-02, 1.39210605e-01,\n",
       "       2.71153623e-02, 8.45563294e-01, 1.28709242e-01, 2.91397480e-01,\n",
       "       3.99715215e-03, 5.02475461e-02])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# p-values\n",
    "\n",
    "p_values = stats.norm.sf(abs(z_scores))*2  # two sided\n",
    "\n",
    "p_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "cc62a1d6",
   "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>feature</th>\n",
       "      <th>coeff</th>\n",
       "      <th>std</th>\n",
       "      <th>z</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>intercept</td>\n",
       "      <td>-0.650519</td>\n",
       "      <td>0.078686</td>\n",
       "      <td>-8.267281</td>\n",
       "      <td>1.370484e-16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>pclass</td>\n",
       "      <td>-0.517374</td>\n",
       "      <td>0.117870</td>\n",
       "      <td>-4.389380</td>\n",
       "      <td>1.136742e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>sibsp</td>\n",
       "      <td>-0.226898</td>\n",
       "      <td>0.095180</td>\n",
       "      <td>-2.383888</td>\n",
       "      <td>1.713082e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>parch</td>\n",
       "      <td>0.054203</td>\n",
       "      <td>0.082032</td>\n",
       "      <td>0.660755</td>\n",
       "      <td>5.087692e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>sex_female</td>\n",
       "      <td>1.199899</td>\n",
       "      <td>0.082074</td>\n",
       "      <td>14.619696</td>\n",
       "      <td>2.103363e-48</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>embarked_S</td>\n",
       "      <td>-0.008459</td>\n",
       "      <td>0.125488</td>\n",
       "      <td>-0.067410</td>\n",
       "      <td>9.462556e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>embarked_C</td>\n",
       "      <td>0.234049</td>\n",
       "      <td>0.127461</td>\n",
       "      <td>1.836241</td>\n",
       "      <td>6.632206e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>cabin_B</td>\n",
       "      <td>-0.219822</td>\n",
       "      <td>0.148655</td>\n",
       "      <td>-1.478737</td>\n",
       "      <td>1.392106e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>cabin_C</td>\n",
       "      <td>-0.376500</td>\n",
       "      <td>0.170373</td>\n",
       "      <td>-2.209853</td>\n",
       "      <td>2.711536e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>cabin_E</td>\n",
       "      <td>-0.025204</td>\n",
       "      <td>0.129395</td>\n",
       "      <td>-0.194782</td>\n",
       "      <td>8.455633e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>cabin_D</td>\n",
       "      <td>-0.198210</td>\n",
       "      <td>0.130469</td>\n",
       "      <td>-1.519212</td>\n",
       "      <td>1.287092e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>cabin_A</td>\n",
       "      <td>-0.105084</td>\n",
       "      <td>0.099600</td>\n",
       "      <td>-1.055061</td>\n",
       "      <td>2.913975e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>cabin_M</td>\n",
       "      <td>-0.680452</td>\n",
       "      <td>0.236400</td>\n",
       "      <td>-2.878386</td>\n",
       "      <td>3.997152e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>cabin_G</td>\n",
       "      <td>-0.143629</td>\n",
       "      <td>0.073361</td>\n",
       "      <td>-1.957851</td>\n",
       "      <td>5.024755e-02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       feature     coeff       std          z       p_value\n",
       "0    intercept -0.650519  0.078686  -8.267281  1.370484e-16\n",
       "1       pclass -0.517374  0.117870  -4.389380  1.136742e-05\n",
       "2        sibsp -0.226898  0.095180  -2.383888  1.713082e-02\n",
       "3        parch  0.054203  0.082032   0.660755  5.087692e-01\n",
       "4   sex_female  1.199899  0.082074  14.619696  2.103363e-48\n",
       "5   embarked_S -0.008459  0.125488  -0.067410  9.462556e-01\n",
       "6   embarked_C  0.234049  0.127461   1.836241  6.632206e-02\n",
       "7      cabin_B -0.219822  0.148655  -1.478737  1.392106e-01\n",
       "8      cabin_C -0.376500  0.170373  -2.209853  2.711536e-02\n",
       "9      cabin_E -0.025204  0.129395  -0.194782  8.455633e-01\n",
       "10     cabin_D -0.198210  0.130469  -1.519212  1.287092e-01\n",
       "11     cabin_A -0.105084  0.099600  -1.055061  2.913975e-01\n",
       "12     cabin_M -0.680452  0.236400  -2.878386  3.997152e-03\n",
       "13     cabin_G -0.143629  0.073361  -1.957851  5.024755e-02"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# put everything together\n",
    "\n",
    "coeffs = pd.DataFrame({\n",
    "    \"feature\": features,\n",
    "    \"coeff\": coefficients,\n",
    "    \"std\": se,\n",
    "    \"z\": z_scores,\n",
    "    \"p_value\": p_values,\n",
    "})\n",
    "\n",
    "coeffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "2293dd2f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:               survived   No. Observations:                 1112\n",
      "Model:                          Logit   Df Residuals:                     1098\n",
      "Method:                           MLE   Df Model:                           13\n",
      "Date:                Tue, 14 Nov 2023   Pseudo R-squ.:                  0.2923\n",
      "Time:                        10:56:40   Log-Likelihood:                -522.06\n",
      "converged:                       True   LL-Null:                       -737.67\n",
      "Covariance Type:            nonrobust   LLR p-value:                 5.606e-84\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -0.6505      0.079     -8.267      0.000      -0.805      -0.496\n",
      "pclass        -0.5174      0.118     -4.389      0.000      -0.748      -0.286\n",
      "sibsp         -0.2268      0.095     -2.383      0.017      -0.413      -0.040\n",
      "parch          0.0542      0.082      0.660      0.509      -0.107       0.215\n",
      "sex_female     1.1999      0.082     14.620      0.000       1.039       1.361\n",
      "embarked_S    -0.0085      0.125     -0.068      0.946      -0.254       0.237\n",
      "embarked_C     0.2340      0.127      1.836      0.066      -0.016       0.484\n",
      "cabin_B       -0.2198      0.149     -1.479      0.139      -0.511       0.072\n",
      "cabin_C       -0.3765      0.170     -2.210      0.027      -0.710      -0.043\n",
      "cabin_E       -0.0252      0.129     -0.195      0.846      -0.279       0.228\n",
      "cabin_D       -0.1982      0.130     -1.519      0.129      -0.454       0.058\n",
      "cabin_A       -0.1051      0.100     -1.055      0.291      -0.300       0.090\n",
      "cabin_M       -0.6804      0.236     -2.878      0.004      -1.144      -0.217\n",
      "cabin_G       -0.1436      0.073     -1.958      0.050      -0.287       0.000\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "# Compare with statsmodels\n",
    "\n",
    "print(logit_res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d003f7e9",
   "metadata": {},
   "source": [
    "We got the same values. Nicely done. Now you know how the values are calculated under the hood in `statsmodels`, and how to determine them if you are using `sklearn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af8250bc",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "fsml",
   "language": "python",
   "name": "fsml"
  },
  "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.10.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}