38-logistic_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Logistic regression\r\n",
    "\r\n",
    "## Introduction\r\n",
    "\r\n",
    "Logistic regression is used when the outcome (dependent variable) has two possible values or categories. A simple linear regression is incompatible with a binary outcome as the distribution of the $\\epsilon$ random term follows a Gaussian distrubtion and therefore the whole term lies between $(-\\infty, +\\infty)$.\r\n",
    "\r\n",
    "Because it is almost always used with two or more independent variables, is should be called _multiple logistic regression_."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "| Type of regression   | Dependent of dep variable (Y) | Examples                                                      |\r\n",
    "| :------------------- | :---------------------------- | :------------------------------------------------------------ | \r\n",
    "| Linear               | Continuous (interval, ratio)  | Enzyme activity, renal function, weight etc.                  |\r\n",
    "| Logistic             | Binary (dichotomous)          | Death during surgery, graduation, recurrence of cancer etc.   |\r\n",
    "| Proportional hazards | Elapsed time to event         | Months until death, quarters in school before graduation etc. |"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The odds of an event occuring equals the probability that that event will occur divided by the probability that it will not occur. Every probability can be expressed as odds, every odds can be expressed as a probability.\r\n",
    "\r\n",
    "The logistic regression model computes the odds from baseline odds and from odds ratios computed for each independent variable:\r\n",
    "$$ \\text{Odds} = (\\text{Baseline Odds}).\\text{OR}_1.\\text{OR}_2.(...)\\text{OR}_n $$\r\n",
    "The baseline odds answers this question:\r\n",
    ">if every single independent X variables equaled 0 what are te odds of a particular category?\r\n",
    "\r\n",
    "To make the baseline odds meaningful, we can encode that variable as `Age - 20` so that X=0 would encode people who are 20.\r\n",
    "\r\n",
    "For each continuous variable such as the age, the corresponding odds ratio answers the question:\r\n",
    "> for each additional year of age, by how much do the odds increase or decrease?\r\n",
    "\r\n",
    "* If the OR associated with age equals 1.0, then the age is not related\r\n",
    "* If OR > 1 then the odds increase by a set percentage for each additional year of age\r\n",
    "\r\n",
    "For instance, an OR of 1.03 would mean that the odds of Y increase by 3% as a person grows older by one year. In another example with OR = 0.99, every year, the odds ratio of `suicide.hr` declines by ca. 1%. The OR for a 40-y compared to a 25-y is $1 \\times 0.99 \\times 0.99 \\times ... 0.99 = 0.99^{15} \\approxeq 0.86$, with the exponent equals $40-25=15$. That means a 40-y persone has about 14% lower odds of being `suicide.hr` than a 25-y.\r\n",
    "\r\n",
    "The same assumptions and attentions as for multiple regression applied to logistic regression (e.g. at least 5-10 events per variable, not too many independent variables in the model, etc.)\r\n",
    "\r\n",
    "## The equation\r\n",
    "\r\n",
    "$$ Y_i = \\beta_0 + \\beta_1 X_{i,1} + \\beta_2 X_{i,2} + ... + \\beta_n X_{i,n} $$\r\n",
    "\r\n",
    "with $Y_i$ the **natural log (ln)**, not the common Log ($\\log_{10}$), of the odds for a particular participant, $\\beta_0$ the natural log of the baseline odds, $\\beta_2$ the natural log of the odds ratio for the first independent variable etc. For example:\r\n",
    "\r\n",
    "$$ \\log \\frac{\\text{p}_\\text{suicide.hr=1}}{1-\\text{p}_\\text{suicide.hr=1}} = \\beta_0 + \\beta_1 \\times \\text{duree} + \\beta_2 \\times \\text{abus} $$\r\n",
    "\r\n",
    "In fact, the log-odds are given by the logit function, which maps a probability $p$ of the response variable being `1` from $ [0,1) $ to $(-\\infty, +\\infty)$, with $ \\text{logit}(p) = \\ln \\frac{p}{1-p} = \\beta_0 + \\beta X $.\r\n",
    "\r\n",
    "The $Y$ value is the natural log of odds, which can be transformed to a probability. Since it implicitly embodies uncertainty, there is no need to explicitly add a random term to the model. Because there is no Gaussian distribution, the method of least squares is not used, instead logistic regression finds the values of the odds ratios using what is called a _maximum likelihood method_ (MLE).\r\n",
    "\r\n",
    "The odds of the response variable being `1` can be obtained by exponentiating the log-odds $ \\frac{p}{1-p} = e^{\\beta_0 + \\beta X}$, and the probability of the response variable being `1` is given by the logistic function $ p = \\frac{1}{1 + e^{-(\\beta_0 + \\beta X)}} $. The first coefficient $\\beta_0$ is always the constant term (intercept) of the model."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Fitting a simple model\r\n",
    "\r\n",
    "We use the dataset from the MOOC 'Introducton à la statistique avec R' which deals with multiple regression."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "source": [
    "import pandas as pd\r\n",
    "import numpy as np\r\n",
    "import statsmodels.api as sm"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "source": [
    "data = pd.read_csv(\r\n",
    "    'https://raw.githubusercontent.com/sbwiecko/intro_statistique_R/master/data/smp2.tsv',\r\n",
    "    delimiter='\\t',\r\n",
    ")\r\n",
    "data.head()"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "    age                prof  duree  discip  n.enfant  n.fratrie  ecole  \\\n",
       "0  31.0               autre    4.0     0.0       2.0          4    1.0   \n",
       "1  49.0                 NaN    NaN     0.0       7.0          3    2.0   \n",
       "2  50.0  prof.intermediaire    5.0     0.0       2.0          2    2.0   \n",
       "3  47.0             ouvrier    NaN     0.0       0.0          6    1.0   \n",
       "4  23.0         sans emploi    4.0     1.0       1.0          6    1.0   \n",
       "\n",
       "   separation  juge.enfant  place  ...  subst.cons  scz.cons  char   rs   ed  \\\n",
       "0         0.0          0.0    0.0  ...           0         0   1.0  2.0  1.0   \n",
       "1         1.0          0.0    0.0  ...           0         0   1.0  2.0  2.0   \n",
       "2         0.0          0.0    0.0  ...           0         0   1.0  2.0  3.0   \n",
       "3         1.0          0.0    1.0  ...           0         0   1.0  2.0  2.0   \n",
       "4         1.0          NaN    1.0  ...           0         0   1.0  2.0  2.0   \n",
       "\n",
       "    dr  suicide.s  suicide.hr  suicide.past  dur.interv  \n",
       "0  1.0        0.0         0.0           0.0         NaN  \n",
       "1  1.0        0.0         0.0           0.0        70.0  \n",
       "2  2.0        0.0         0.0           0.0         NaN  \n",
       "3  2.0        1.0         0.0           0.0       105.0  \n",
       "4  2.0        0.0         0.0           1.0         NaN  \n",
       "\n",
       "[5 rows x 26 columns]"
      ],
      "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>prof</th>\n",
       "      <th>duree</th>\n",
       "      <th>discip</th>\n",
       "      <th>n.enfant</th>\n",
       "      <th>n.fratrie</th>\n",
       "      <th>ecole</th>\n",
       "      <th>separation</th>\n",
       "      <th>juge.enfant</th>\n",
       "      <th>place</th>\n",
       "      <th>...</th>\n",
       "      <th>subst.cons</th>\n",
       "      <th>scz.cons</th>\n",
       "      <th>char</th>\n",
       "      <th>rs</th>\n",
       "      <th>ed</th>\n",
       "      <th>dr</th>\n",
       "      <th>suicide.s</th>\n",
       "      <th>suicide.hr</th>\n",
       "      <th>suicide.past</th>\n",
       "      <th>dur.interv</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>31.0</td>\n",
       "      <td>autre</td>\n",
       "      <td>4.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>4</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>49.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>3</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>70.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>50.0</td>\n",
       "      <td>prof.intermediaire</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>47.0</td>\n",
       "      <td>ouvrier</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>6</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>105.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>23.0</td>\n",
       "      <td>sans emploi</td>\n",
       "      <td>4.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>6</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 26 columns</p>\n",
       "</div>"
      ]
     },
     "metadata": {},
     "execution_count": 102
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "source": [
    "data_ = data.dropna(subset=['suicide.hr', 'abus'], how='any')\r\n",
    "y = data_['suicide.hr']\r\n",
    "X = data_['abus']\r\n",
    "X = sm.add_constant(X)\r\n",
    "model = sm.Logit(y, X) # same as sm.GLM(y, X, family=sm.families.Binomial())\r\n",
    "result = model.fit()"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.494196\n",
      "         Iterations 5\n"
     ]
    },
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "C:\\Users\\WIECKOWS\\AppData\\Local\\Programs\\Python\\Python39\\lib\\site-packages\\statsmodels\\tsa\\tsatools.py:142: FutureWarning: In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only\n",
      "  x = pd.concat(x[::order], 1)\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "source": [
    "print(result.summary())"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:             suicide.hr   No. Observations:                  753\n",
      "Model:                          Logit   Df Residuals:                      751\n",
      "Method:                           MLE   Df Model:                            1\n",
      "Date:                Thu, 12 Aug 2021   Pseudo R-squ.:                 0.02098\n",
      "Time:                        14:58:16   Log-Likelihood:                -372.13\n",
      "converged:                       True   LL-Null:                       -380.11\n",
      "Covariance Type:            nonrobust   LLR p-value:                 6.494e-05\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -1.6161      0.115    -14.003      0.000      -1.842      -1.390\n",
      "abus           0.7688      0.190      4.052      0.000       0.397       1.141\n",
      "==============================================================================\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "source": [
    "result.params"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "const   -1.616082\n",
       "abus     0.768785\n",
       "dtype: float64"
      ]
     },
     "metadata": {},
     "execution_count": 105
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "source": [
    "result.conf_int()"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "              0         1\n",
       "const -1.842275 -1.389890\n",
       "abus   0.396937  1.140632"
      ],
      "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>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>const</th>\n",
       "      <td>-1.842275</td>\n",
       "      <td>-1.389890</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>abus</th>\n",
       "      <td>0.396937</td>\n",
       "      <td>1.140632</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ]
     },
     "metadata": {},
     "execution_count": 106
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "source": [
    "OR1 = np.exp(result.params[1])\r\n",
    "CI_OR1 = np.round(np.exp(result.conf_int().loc['abus']), 2)\r\n",
    "print(f\"Odds ratio for 'abus': {OR1:3.3f}\")\r\n",
    "print(f\"CI of the OR for 'abus': {CI_OR1[0]} - {CI_OR1[1]}\")"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Odds ratio for 'abus': 2.157\n",
      "CI of the OR for 'abus': 1.49 - 3.13\n"
     ]
    }
   ],
   "metadata": {
    "scrolled": true
   }
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "source": [
    "data_['abus'].unique()"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "array([0., 1.])"
      ]
     },
     "metadata": {},
     "execution_count": 108
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The independent variable has two possible values. The corresponding odds ratio is 1.0 for `abus=0` and 2.16 for `abus=1`, meaning that participant with `abus=1`, **but sharing the other attributes**, has a bit more than twice the odds of being `suicide.hr`, with CI ranging from 1.49 to 3.13. If the H0 were true, the OR would equal 1.0, but the CI above doesn't include 1.0, so the corresponding P value must be less than 0.05.\r\n",
    "\r\n",
    "In the summary table P < 0.0001, therefore the association between `abus` and `suicide.hr` is statistically significant."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Analysis of the corresponding contingency table"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "source": [
    "table = pd.crosstab(data_['abus'], data_['suicide.hr'])\r\n",
    "sm.stats.Table2x2(table).summary()"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ],
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>        <th>Estimate</th>  <th>SE</th>    <th>LCB</th>   <th>UCB</th>  <th>p-value</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Odds ratio</th>        <td>2.157</td>      <td></td> <td>1.487</td> <td>3.129</td>   <td>0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Log odds ratio</th>    <td>0.769</td> <td>0.190</td> <td>0.397</td> <td>1.141</td>   <td>0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Risk ratio</th>        <td>1.192</td>      <td></td> <td>1.083</td> <td>1.312</td>   <td>0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Log risk ratio</th>    <td>0.175</td> <td>0.049</td> <td>0.079</td> <td>0.272</td>   <td>0.000</td>\n",
       "</tr>\n",
       "</table>"
      ]
     },
     "metadata": {},
     "execution_count": 109
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The Odds Ratio obtained with the analysis of the contingency table equals the Odds Ratio computed using the logistic regression."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Analysis using pingouin"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "source": [
    "import pingouin as pg\r\n",
    "\r\n",
    "lom = pg.logistic_regression(\r\n",
    "    X=X,\r\n",
    "    y=y\r\n",
    ")\r\n",
    "\r\n",
    "lom.round(3)"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "       names   coef     se       z  pval  CI[2.5%]  CI[97.5%]\n",
       "0  Intercept -1.616  0.115 -14.003   0.0    -1.842     -1.390\n",
       "1       abus  0.769  0.190   4.052   0.0     0.397      1.141"
      ],
      "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>names</th>\n",
       "      <th>coef</th>\n",
       "      <th>se</th>\n",
       "      <th>z</th>\n",
       "      <th>pval</th>\n",
       "      <th>CI[2.5%]</th>\n",
       "      <th>CI[97.5%]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Intercept</td>\n",
       "      <td>-1.616</td>\n",
       "      <td>0.115</td>\n",
       "      <td>-14.003</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-1.842</td>\n",
       "      <td>-1.390</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>abus</td>\n",
       "      <td>0.769</td>\n",
       "      <td>0.190</td>\n",
       "      <td>4.052</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.397</td>\n",
       "      <td>1.141</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ]
     },
     "metadata": {},
     "execution_count": 110
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Interpretation of the intercept coefficient"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "source": [
    "intercept = lom.loc[0, 'coef']\r\n",
    "\r\n",
    "print(\"Coefficients for intercept (`abus=0`):\")\r\n",
    "print(\"--------------------------------------\")\r\n",
    "print(f\"Log-odds\\t {intercept:.3f}\")\r\n",
    "print(f\"Odds    \\t {np.exp(intercept):.2f}\")\r\n",
    "print(f\"Ratio   \\t 1:{1/np.exp(intercept):.0f}\")\r\n",
    "print(f\"Proba   \\t {1/(1+np.exp(-(intercept))):.3f}\")"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Coefficients for intercept (`abus=0`):\n",
      "--------------------------------------\n",
      "Log-odds\t -1.616\n",
      "Odds    \t 0.20\n",
      "Ratio   \t 1:5\n",
      "Proba   \t 0.166\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "source": [
    "coeff_abus = lom.loc[1, 'coef']\r\n",
    "np.exp(coeff_abus)"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "2.1571428116306226"
      ]
     },
     "metadata": {},
     "execution_count": 112
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "source": [
    "print(\"Coefficients for `abus=1`):\")\r\n",
    "print(\"---------------------------\")\r\n",
    "print(f\"Log-odds \\t {intercept + coeff_abus:.3f}\")\r\n",
    "print(f\"Odds     \\t {np.exp(intercept + coeff_abus):.2f}\")\r\n",
    "print(f\"Ratio ca.\\t 1:{1/np.exp(intercept + coeff_abus):.0f}\")\r\n",
    "print(f\"Proba    \\t {1/(1+np.exp(-(intercept + coeff_abus))):.3f}\")"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Coefficients for `abus=1`):\n",
      "---------------------------\n",
      "Log-odds \t -0.847\n",
      "Odds     \t 0.43\n",
      "Ratio ca.\t 1:2\n",
      "Proba    \t 0.300\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In the previous example, the _intercept_ or _constant_ coefficient (-1.616) is the log-odds of `suicide.hr=1` when `abus=0`. Thus, the odds ratio of _suicide when no abus_ is `0.20:1`, i.e. a probability of ca. 16.6%.\r\n",
    "\r\n",
    "The _abus_ coefficient (0.769) means that for each additional _abus_ (altough for this binary parameter there are only 2 possible values), the log-odds of _suicide.hr_ increases by 0.769, therefore  the odds are multiplied by $e^{0.769}$, i.e. the odds are multiplied by ca. 2.16. The table below shows the probability of `suicide.hr` for both values of `abus`:\r\n",
    "\r\n",
    "| abus | Log-odds | Odds | Ratio | Proba |\r\n",
    "| ---- | -------- | ---- | ----- | ----- |\r\n",
    "|  0   |  -1.616  | 0.20 |  1:5  | 0.166 |\r\n",
    "|  1   |  -0.847  | 0.43 |  1:2  | 0.300 |"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Fitting a more complex model"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "source": [
    "data2_ = data.dropna(\r\n",
    "    subset=['suicide.hr', 'abus', 'discip', 'duree'], how='any')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "source": [
    "y2 = data2_['suicide.hr']\r\n",
    "X2 = data2_[['abus', 'discip', 'duree']]\r\n",
    "# duree is graduated from 1 to 5\r\n",
    "X2 = sm.add_constant(X2)\r\n",
    "model2 = sm.Logit(y2, X2)\r\n",
    "result2 = model2.fit()\r\n",
    "print(result2.summary2())"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.484782\n",
      "         Iterations 6\n",
      "                         Results: Logit\n",
      "=================================================================\n",
      "Model:              Logit            Pseudo R-squared: 0.041     \n",
      "Dependent Variable: suicide.hr       AIC:              541.2607  \n",
      "Date:               2021-08-12 15:01 BIC:              558.5004  \n",
      "No. Observations:   550              Log-Likelihood:   -266.63   \n",
      "Df Model:           3                LL-Null:          -277.97   \n",
      "Df Residuals:       546              LLR p-value:      4.7043e-05\n",
      "Converged:          1.0000           Scale:            1.0000    \n",
      "No. Iterations:     6.0000                                       \n",
      "-------------------------------------------------------------------\n",
      "           Coef.    Std.Err.      z      P>|z|     [0.025    0.975]\n",
      "-------------------------------------------------------------------\n",
      "const     -0.0246     0.4963   -0.0496   0.9604   -0.9974    0.9482\n",
      "abus       0.6229     0.2276    2.7363   0.0062    0.1767    1.0691\n",
      "discip     0.5281     0.2377    2.2219   0.0263    0.0623    0.9939\n",
      "duree     -0.3986     0.1172   -3.4004   0.0007   -0.6284   -0.1689\n",
      "=================================================================\n",
      "\n"
     ]
    },
    {
     "output_type": "stream",
     "name": "stderr",
     "text": [
      "C:\\Users\\WIECKOWS\\AppData\\Local\\Programs\\Python\\Python39\\lib\\site-packages\\statsmodels\\tsa\\tsatools.py:142: FutureWarning: In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only\n",
      "  x = pd.concat(x[::order], 1)\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can already interpret the sign of the Log odds coefficients, with negative coefficient associated with a decrease in risk, as for example with the explanatory variable `duree`. For measuring the amplitude of the effect, further manipulation is required as shown below."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "source": [
    "OR2 = np.exp(result2.params)\r\n",
    "OR2"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "const     0.975680\n",
       "abus      1.864315\n",
       "discip    1.695687\n",
       "duree     0.671249\n",
       "dtype: float64"
      ]
     },
     "metadata": {},
     "execution_count": 119
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "source": [
    ".99**15"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "0.8600583546412884"
      ]
     },
     "metadata": {},
     "execution_count": 63
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Table of the probability of `suicide.hr=1` for the different `duree` values"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "source": [
    "lom2 = pg.logistic_regression(\r\n",
    "    X=X2,\r\n",
    "    y=y2\r\n",
    ")\r\n",
    "\r\n",
    "lom2.round(4)"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "       names    coef      se       z    pval  CI[2.5%]  CI[97.5%]\n",
       "0  Intercept -0.0246  0.4963 -0.0496  0.9604   -0.9974     0.9482\n",
       "1       abus  0.6229  0.2276  2.7363  0.0062    0.1767     1.0691\n",
       "2     discip  0.5281  0.2377  2.2219  0.0263    0.0623     0.9939\n",
       "3      duree -0.3986  0.1172 -3.4004  0.0007   -0.6284    -0.1689"
      ],
      "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>names</th>\n",
       "      <th>coef</th>\n",
       "      <th>se</th>\n",
       "      <th>z</th>\n",
       "      <th>pval</th>\n",
       "      <th>CI[2.5%]</th>\n",
       "      <th>CI[97.5%]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Intercept</td>\n",
       "      <td>-0.0246</td>\n",
       "      <td>0.4963</td>\n",
       "      <td>-0.0496</td>\n",
       "      <td>0.9604</td>\n",
       "      <td>-0.9974</td>\n",
       "      <td>0.9482</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>abus</td>\n",
       "      <td>0.6229</td>\n",
       "      <td>0.2276</td>\n",
       "      <td>2.7363</td>\n",
       "      <td>0.0062</td>\n",
       "      <td>0.1767</td>\n",
       "      <td>1.0691</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>discip</td>\n",
       "      <td>0.5281</td>\n",
       "      <td>0.2377</td>\n",
       "      <td>2.2219</td>\n",
       "      <td>0.0263</td>\n",
       "      <td>0.0623</td>\n",
       "      <td>0.9939</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>duree</td>\n",
       "      <td>-0.3986</td>\n",
       "      <td>0.1172</td>\n",
       "      <td>-3.4004</td>\n",
       "      <td>0.0007</td>\n",
       "      <td>-0.6284</td>\n",
       "      <td>-0.1689</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ]
     },
     "metadata": {},
     "execution_count": 125
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "source": [
    "intercept2  = lom2.loc[0, 'coef']\r\n",
    "coeff_duree = lom2.loc[3, 'coef']\r\n",
    "\r\n",
    "print(\"`duree`\\tLogodds\\tOdds\\tRatio\\tProba\")\r\n",
    "print(\"-------------------------------------\")\r\n",
    "\r\n",
    "for duree in range(6):  # for values 0 to 5\r\n",
    "    print(duree, end='\\t')\r\n",
    "    print(f\"{intercept + duree*coeff_duree:.3f}\", end='\\t')\r\n",
    "    print(f\"{np.exp(intercept + duree*coeff_duree):.3f}\", end='\\t')\r\n",
    "    print(f\"1:{1/np.exp(intercept + duree*coeff_duree):.0f}\", end='\\t')\r\n",
    "    print(f\"{1/(1+np.exp(-(intercept + duree*coeff_duree))):.3f}\")"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "`duree`\tLogodds\tOdds\tRatio\tProba\n",
      "-------------------------------------\n",
      "0\t-1.616\t0.199\t1:5\t0.166\n",
      "1\t-2.015\t0.133\t1:7\t0.118\n",
      "2\t-2.413\t0.090\t1:11\t0.082\n",
      "3\t-2.812\t0.060\t1:17\t0.057\n",
      "4\t-3.211\t0.040\t1:25\t0.039\n",
      "5\t-3.609\t0.027\t1:37\t0.026\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Interactions (synergies) - R-style formula"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "source": [
    "import statsmodels.formula.api as smf"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "source": [
    "model3 = smf.logit(\r\n",
    "    formula=\"Q('suicide.hr') ~ abus + discip*duree\", \r\n",
    "    data=data_\r\n",
    ")\r\n",
    "# Q() is a way to 'quote' variable names, especially ones that do not otherwise \r\n",
    "# meet Python's variable name rules, such as with a dot in the variable name\r\n",
    "\r\n",
    "# \":\" adds a new column to the design matrix with the product of the other two columns\r\n",
    "# \"*\" will also include the individual columns that were multiplied together\r\n",
    "# thus \"Q('suicide.hr') ~ abus + discip*duree\" is eq \"Q('suicide.hr') ~ abus + discip + dure + discip:duree\"\r\n",
    "\r\n",
    "print(model3.fit().summary2())"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.484782\n",
      "         Iterations 6\n",
      "                         Results: Logit\n",
      "=================================================================\n",
      "Model:              Logit            Pseudo R-squared: 0.041     \n",
      "Dependent Variable: Q('suicide.hr')  AIC:              543.2599  \n",
      "Date:               2021-08-12 15:03 BIC:              564.8095  \n",
      "No. Observations:   550              Log-Likelihood:   -266.63   \n",
      "Df Model:           4                LL-Null:          -277.97   \n",
      "Df Residuals:       545              LLR p-value:      0.00014651\n",
      "Converged:          1.0000           Scale:            1.0000    \n",
      "No. Iterations:     6.0000                                       \n",
      "------------------------------------------------------------------\n",
      "               Coef.   Std.Err.     z     P>|z|    [0.025   0.975]\n",
      "------------------------------------------------------------------\n",
      "Intercept     -0.0315    0.5512  -0.0571  0.9544  -1.1119   1.0489\n",
      "abus           0.6229    0.2276   2.7362  0.0062   0.1767   1.0691\n",
      "discip         0.5639    1.2686   0.4445  0.6567  -1.9226   3.0503\n",
      "duree         -0.3969    0.1311  -3.0286  0.0025  -0.6538  -0.1401\n",
      "discip:duree  -0.0084    0.2917  -0.0287  0.9771  -0.5802   0.5634\n",
      "=================================================================\n",
      "\n"
     ]
    }
   ],
   "metadata": {}
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}