{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# импортируем необходимые библиотеки и классы\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from statsmodels.formula.api import logit\n",
    "np.set_printoptions(precision=None, suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>age</th>\n",
       "      <th>job</th>\n",
       "      <th>debtinc</th>\n",
       "      <th>creddebt</th>\n",
       "      <th>othdebt</th>\n",
       "      <th>default</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>28</td>\n",
       "      <td>working - other</td>\n",
       "      <td>17.7</td>\n",
       "      <td>2.990592</td>\n",
       "      <td>4.797408</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>64</td>\n",
       "      <td>working - production</td>\n",
       "      <td>14.7</td>\n",
       "      <td>5.047392</td>\n",
       "      <td>12.004608</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>40</td>\n",
       "      <td>working - IT</td>\n",
       "      <td>4.8</td>\n",
       "      <td>1.042368</td>\n",
       "      <td>1.885632</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>30</td>\n",
       "      <td>working - IT</td>\n",
       "      <td>34.5</td>\n",
       "      <td>1.751220</td>\n",
       "      <td>7.563780</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>25</td>\n",
       "      <td>working - IT</td>\n",
       "      <td>22.4</td>\n",
       "      <td>0.759360</td>\n",
       "      <td>5.960640</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   age                   job  debtinc  creddebt    othdebt  default\n",
       "0   28       working - other     17.7  2.990592   4.797408        0\n",
       "1   64  working - production     14.7  5.047392  12.004608        0\n",
       "2   40          working - IT      4.8  1.042368   1.885632        0\n",
       "3   30          working - IT     34.5  1.751220   7.563780        0\n",
       "4   25          working - IT     22.4  0.759360   5.960640        1"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# загружаем и смотрим данные \n",
    "data = pd.read_csv('Data/Bankloan.csv', sep=';', decimal=',')\n",
    "data.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# формируем массив меток и массив предикторов\n",
    "y = data.pop('default')\n",
    "# создаем список количественных предикторов\n",
    "num_columns = data.select_dtypes(exclude='object').columns\n",
    "# создаем экземпляр класса StandardScaler\n",
    "scal = StandardScaler()\n",
    "# выполняем стандартизацию\n",
    "data[num_columns] = scal.fit_transform(data[num_columns])\n",
    "# выполняем дамми-кодирование\n",
    "X = pd.get_dummies(data, drop_first=True)\n",
    "# записываем названия предикторов\n",
    "feat_labels = X.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# создаем экземпляр класса LogisticRegression, \n",
    "# отключаем регуляризации\n",
    "logreg_sklearn = LogisticRegression(solver='lbfgs', \n",
    "                                    penalty='none')\n",
    "# обучаем модель\n",
    "logreg_sklearn.fit(X, y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Найдем стандартные ошибки коэффициентов модели, которые вычисляются как    \n",
    "квадратный корень из диагональных элементов матрицы ковариации предикторов. <br>\n",
    "\n",
    "Матрица ковариации может быть получена следующим образом:\n",
    "(X.T * V * X) ** (-1), где  \n",
    "> X - матрица признаков с добавленным в начало столбцом из единиц (для константы)  \n",
    "> V - диагональная матрица, c диагональными элементами v[i, i] =  pi * (1 - pi),  \n",
    "> где pi - предсказанная вероятность положительного класса для i-го наблюдения "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1500, 2)\n",
      "\n",
      "[[0.289 0.711]\n",
      " [0.914 0.086]\n",
      " [0.921 0.079]\n",
      " ...\n",
      " [0.722 0.278]\n",
      " [0.687 0.313]\n",
      " [0.773 0.227]]\n"
     ]
    }
   ],
   "source": [
    "# вычисляем матрицу спрогнозированных вероятностей классов\n",
    "predProbs = logreg_sklearn.predict_proba(X)\n",
    "print(predProbs.shape)\n",
    "print('')\n",
    "print(np.round(predProbs, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1500, 9)\n",
      "\n",
      "[[ 1.    -0.47   1.165 ...  0.     1.     0.   ]\n",
      " [ 1.     2.27   0.715 ...  0.     0.     1.   ]\n",
      " [ 1.     0.443 -0.769 ...  1.     0.     0.   ]\n",
      " ...\n",
      " [ 1.    -0.698 -1.009 ...  0.     0.     1.   ]\n",
      " [ 1.    -0.013  0.685 ...  1.     0.     0.   ]\n",
      " [ 1.    -0.546 -0.619 ...  1.     0.     0.   ]]\n"
     ]
    }
   ],
   "source": [
    "# создаем матрицу плана, добавляем \n",
    "# в начало столбец с единицами \n",
    "X_ = np.hstack([np.ones((X.shape[0], 1)), X])\n",
    "print(X_.shape)\n",
    "print('')\n",
    "print(np.round(X_, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1500, 1500)\n",
      "\n",
      "[[0.20530004 0.         0.         ... 0.         0.         0.        ]\n",
      " [0.         0.07837879 0.         ... 0.         0.         0.        ]\n",
      " [0.         0.         0.07273397 ... 0.         0.         0.        ]\n",
      " ...\n",
      " [0.         0.         0.         ... 0.20090302 0.         0.        ]\n",
      " [0.         0.         0.         ... 0.         0.21486664 0.        ]\n",
      " [0.         0.         0.         ... 0.         0.         0.17549169]]\n"
     ]
    }
   ],
   "source": [
    "# создаем диагональную матрицу, помещаем на главную \n",
    "# диагональ произведение вероятностей отрицательного и\n",
    "# положительного классов для каждого наблюдения\n",
    "V = np.diagflat(np.product(predProbs, axis=1))\n",
    "print(V.shape)\n",
    "print(\"\")\n",
    "print(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ковариационная матрица:\n",
      " [[ 0.0185  0.0022  0.0008 -0.     -0.0018 -0.018  -0.0186 -0.0177 -0.0185]\n",
      " [ 0.0022  0.0104  0.0001 -0.0041 -0.0018  0.0008 -0.0011  0.0019 -0.0022]\n",
      " [ 0.0008  0.0001  0.0069 -0.0006 -0.0033 -0.0016 -0.0029 -0.0014 -0.0006]\n",
      " [-0.     -0.0041 -0.0006  0.0133 -0.0065 -0.001  -0.     -0.0013 -0.001 ]\n",
      " [-0.0018 -0.0018 -0.0033 -0.0065  0.0114  0.0021  0.003   0.0018  0.0017]\n",
      " [-0.018   0.0008 -0.0016 -0.001   0.0021  0.0364  0.0187  0.0186  0.0179]\n",
      " [-0.0186 -0.0011 -0.0029 -0.      0.003   0.0187  0.048   0.0185  0.0184]\n",
      " [-0.0177  0.0019 -0.0014 -0.0013  0.0018  0.0186  0.0185  0.0299  0.0176]\n",
      " [-0.0185 -0.0022 -0.0006 -0.001   0.0017  0.0179  0.0184  0.0176  0.0912]]\n"
     ]
    }
   ],
   "source": [
    "# теперь вычисляем и выводим ковариационную матрицу\n",
    "\n",
    "'''\n",
    "Обратите внимание, что в Python 3.5+ оператор @ выполняет \n",
    "умножение матриц, таким образом, если вы используете \n",
    "Python 3.5+, вы можете заменить строчку кода ниже\n",
    "на более читаемый вариант:\n",
    "cov_matrix = np.linalg.pinv(X_.T @ V @ X_)\n",
    "'''\n",
    "cov_matrix = np.linalg.inv(np.dot(np.dot(X_.T, V), X_))\n",
    "print(\"Ковариационная матрица:\\n\", np.round(cov_matrix, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# вычисляем коэффициенты, z-значения, статистики Вальда, p-значения\n",
    "coefs = np.round(np.insert(\n",
    "    logreg_sklearn.coef_, 0, logreg_sklearn.intercept_), 4)\n",
    "se = np.round(np.sqrt(np.diag(cov_matrix)), 4)\n",
    "z_values = np.round(coefs / se, 4)\n",
    "wald_values = np.round((coefs / se) ** 2, 4)\n",
    "p_values = np.round((1 - stats.norm.cdf(abs(z_values))) * 2, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "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>coef</th>\n",
       "      <th>se</th>\n",
       "      <th>z_value</th>\n",
       "      <th>Wald</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>-0.3205</td>\n",
       "      <td>0.1362</td>\n",
       "      <td>-2.3532</td>\n",
       "      <td>5.5373</td>\n",
       "      <td>0.0186</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>-1.1907</td>\n",
       "      <td>0.1019</td>\n",
       "      <td>-11.6850</td>\n",
       "      <td>136.5389</td>\n",
       "      <td>0.0000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>debtinc</th>\n",
       "      <td>0.8541</td>\n",
       "      <td>0.0833</td>\n",
       "      <td>10.2533</td>\n",
       "      <td>105.1302</td>\n",
       "      <td>0.0000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>creddebt</th>\n",
       "      <td>0.8140</td>\n",
       "      <td>0.1155</td>\n",
       "      <td>7.0476</td>\n",
       "      <td>49.6689</td>\n",
       "      <td>0.0000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>othdebt</th>\n",
       "      <td>-0.2967</td>\n",
       "      <td>0.1070</td>\n",
       "      <td>-2.7729</td>\n",
       "      <td>7.6890</td>\n",
       "      <td>0.0056</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>job_own business</th>\n",
       "      <td>-0.5412</td>\n",
       "      <td>0.1908</td>\n",
       "      <td>-2.8365</td>\n",
       "      <td>8.0456</td>\n",
       "      <td>0.0046</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>job_working - IT</th>\n",
       "      <td>-0.8158</td>\n",
       "      <td>0.2191</td>\n",
       "      <td>-3.7234</td>\n",
       "      <td>13.8638</td>\n",
       "      <td>0.0002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>job_working - other</th>\n",
       "      <td>-0.5678</td>\n",
       "      <td>0.1729</td>\n",
       "      <td>-3.2840</td>\n",
       "      <td>10.7845</td>\n",
       "      <td>0.0010</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>job_working - production</th>\n",
       "      <td>-0.3522</td>\n",
       "      <td>0.3020</td>\n",
       "      <td>-1.1662</td>\n",
       "      <td>1.3601</td>\n",
       "      <td>0.2435</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            coef      se  z_value      Wald  p_value\n",
       "Intercept                -0.3205  0.1362  -2.3532    5.5373   0.0186\n",
       "age                      -1.1907  0.1019 -11.6850  136.5389   0.0000\n",
       "debtinc                   0.8541  0.0833  10.2533  105.1302   0.0000\n",
       "creddebt                  0.8140  0.1155   7.0476   49.6689   0.0000\n",
       "othdebt                  -0.2967  0.1070  -2.7729    7.6890   0.0056\n",
       "job_own business         -0.5412  0.1908  -2.8365    8.0456   0.0046\n",
       "job_working - IT         -0.8158  0.2191  -3.7234   13.8638   0.0002\n",
       "job_working - other      -0.5678  0.1729  -3.2840   10.7845   0.0010\n",
       "job_working - production -0.3522  0.3020  -1.1662    1.3601   0.2435"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# сохраняем результаты вычислений в датафрейм\n",
    "feat_labels = np.insert(feat_labels, 0, 'Intercept')\n",
    "summary = np.array([coefs, se, z_values, wald_values, p_values])\n",
    "columns = ['coef', 'se', 'z_value', 'Wald', 'p_value']\n",
    "logit_skl_summary = pd.DataFrame(summary.T, \n",
    "                                 index=feat_labels,\n",
    "                                 columns=columns)\n",
    "logit_skl_summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'default ~ age + debtinc + creddebt + othdebt + job_own_business + job_working___IT + job_working___other + job_working___production + '"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# конкатенируем массив признаков и массив меток\n",
    "data = pd.concat([X, y], axis=1)\n",
    "# преобразовываем имена столбцов в формат, понятный для statmodels\n",
    "data.rename(columns=lambda x: x.replace(' ', '_').replace('-', '_'), \n",
    "            inplace=True)\n",
    "# создаем и выводим формулу\n",
    "columns = data.columns.drop('default')\n",
    "formula = 'default ~ '\n",
    "for col in columns:\n",
    "    formula = formula + str(col) + ' + '\n",
    "formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'default ~ age + debtinc + creddebt + othdebt + job_own_business + job_working___IT + job_working___other + job_working___production'"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# редактируем формулу, удалив три лишних символа с конца\n",
    "formula = formula[0: -3]\n",
    "formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.500446\n",
      "         Iterations 6\n"
     ]
    }
   ],
   "source": [
    "# создаем экземпляр класса logit и обучаем \n",
    "logreg_statsmodels = logit(formula, data).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                default   No. Observations:                 1500\n",
      "Model:                          Logit   Df Residuals:                     1491\n",
      "Method:                           MLE   Df Model:                            8\n",
      "Date:                Sat, 11 Feb 2023   Pseudo R-squ.:                  0.2376\n",
      "Time:                        15:48:08   Log-Likelihood:                -750.67\n",
      "converged:                       True   LL-Null:                       -984.64\n",
      "Covariance Type:            nonrobust   LLR p-value:                 5.294e-96\n",
      "============================================================================================\n",
      "                               coef    std err          z      P>|z|      [0.025      0.975]\n",
      "--------------------------------------------------------------------------------------------\n",
      "Intercept                   -0.3205      0.136     -2.354      0.019      -0.587      -0.054\n",
      "age                         -1.1907      0.102    -11.680      0.000      -1.391      -0.991\n",
      "debtinc                      0.8541      0.083     10.256      0.000       0.691       1.017\n",
      "creddebt                     0.8140      0.116      7.045      0.000       0.588       1.040\n",
      "othdebt                     -0.2967      0.107     -2.774      0.006      -0.506      -0.087\n",
      "job_own_business            -0.5412      0.191     -2.836      0.005      -0.915      -0.167\n",
      "job_working___IT            -0.8158      0.219     -3.723      0.000      -1.245      -0.386\n",
      "job_working___other         -0.5678      0.173     -3.285      0.001      -0.907      -0.229\n",
      "job_working___production    -0.3522      0.302     -1.167      0.243      -0.944       0.240\n",
      "============================================================================================\n"
     ]
    }
   ],
   "source": [
    "# смотрим полученные результаты\n",
    "print(logreg_statsmodels.summary())"
   ]
  },
  {
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
