capitulo2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f141f0d3-89c3-4988-91ea-6ac0edb48220",
   "metadata": {},
   "source": [
    "### Librerias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "63e687c8-683d-411d-b675-62dc4281c3fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2536b74-ae81-445d-88ef-bce950f99f21",
   "metadata": {},
   "source": [
    "### Librerias para interactuar con R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7d526e6a-91e7-4bb6-b0fb-deffc3c091d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rpy2.robjects as ro\n",
    "from rpy2.robjects.packages import importr\n",
    "from rpy2.robjects import pandas2ri\n",
    "from rpy2.robjects import r\n",
    "\n",
    "\n",
    "from rpy2.robjects.conversion import localconverter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4c4cb78-ed48-431f-9ad4-fc75468de82a",
   "metadata": {},
   "source": [
    "### Convertidor DataFrame R a Pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "34b1700b-a2c8-44bc-be68-111d5fe9eb70",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dset(datasetr):\n",
    "    x=r(datasetr)\n",
    "    with localconverter(ro.default_converter + pandas2ri.converter):\n",
    "        pd_from_r_df = ro.conversion.rpy2py(x)\n",
    "        return pd_from_r_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa38a502-b778-490e-b222-1b1659578026",
   "metadata": {},
   "source": [
    "### Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "516c0bfe-528c-4edc-abfd-34c8677a2997",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <span>StrVector with 9 elements.</span>\n",
       "        <table>\n",
       "        <tbody>\n",
       "          <tr>\n",
       "          \n",
       "            <td>\n",
       "            'aplore3'\n",
       "            </td>\n",
       "          \n",
       "            <td>\n",
       "            'tools'\n",
       "            </td>\n",
       "          \n",
       "            <td>\n",
       "            'stats'\n",
       "            </td>\n",
       "          \n",
       "            <td>\n",
       "            ...\n",
       "            </td>\n",
       "          \n",
       "            <td>\n",
       "            'datasets'\n",
       "            </td>\n",
       "          \n",
       "            <td>\n",
       "            'methods'\n",
       "            </td>\n",
       "          \n",
       "            <td>\n",
       "            'base'\n",
       "            </td>\n",
       "          \n",
       "          </tr>\n",
       "        </tbody>\n",
       "        </table>\n",
       "        "
      ],
      "text/plain": [
       "<rpy2.robjects.vectors.StrVector object at 0x00000228DDB1F280> [RTYPES.STRSXP]\n",
       "R classes: ('character',)\n",
       "['aplore3', 'tools', 'stats', 'graphics', ..., 'utils', 'datasets', 'methods', 'base']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r('library(aplore3)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17cc2309-0f33-4da6-b468-37586d0dabb2",
   "metadata": {},
   "source": [
    "## 2. Regresión Logística Multiple\n",
    "\n",
    "Consideramos una coleción de $p$ variables denotada $X^T = (X_1,X_2,\\cdots, X_p)$. Denotamos la probabilidad condicional $\\mathbb{P}(Y=1|X) = \\pi(X)$\n",
    "\n",
    "$$\\pi(X) = \\frac{e^{\\beta_0+\\beta_1x_1 + \\beta_2x_2+\\cdots+\\beta_px_p}}{1+e^{\\beta_0+\\beta_1x_1 + \\beta_2x_2+\\cdots+\\beta_px_p}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4dee0e4-3459-42cd-b9d5-fa9744683fc1",
   "metadata": {},
   "source": [
    "Si tenemos una muestra de n observaciones independientes $(x_i,y_i)$ con $i=1,2,\\cdots,n$. $x_i$ un vector con $p$ valores e $y_i$ toma el valor 0 o 1.\n",
    "\n",
    "Las ecuaciones del likelihood son las mismas que para el caso univariado. \n",
    "\n",
    "$$\\sum_{i=1}^n y_i - \\pi(x_i) = 0$$\n",
    "\n",
    "$$\\sum_{i=1}^n x_{ij}(y_i - \\pi(x_i)) = 0$$\n",
    "\n",
    "para $j=1,2,\\cdots,p$\n",
    "\n",
    "Para la estimación de las varianzas y covarianzas de los coeficientes necesitamos las segundas derivadas del loglikelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa23cd51-ebaf-48bd-bf17-bcffb15751fd",
   "metadata": {},
   "source": [
    "$$\\frac{\\partial^2 L}{\\partial \\beta_j^2} = -\\sum_{i=1}^nx_{ij}^2\\pi(x_i)(1-\\pi(x_i))$$\n",
    "\n",
    "$$\\frac{\\partial^2 L}{\\partial \\beta_j\\beta_k} = -\\sum_{i=1}^nx_{ij}x_{ik}\\pi(x_i)(1-\\pi(x_i))$$\n",
    "\n",
    "$j,k = 0,1,2,\\cdots,p$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa12c49f-a7be-4060-8878-8319c6b4c331",
   "metadata": {},
   "source": [
    "El estimador de la varianza está dado por la inversa de la matriz de información observada $I$\n",
    "\n",
    "$$\\hat{Var(\\beta)} = \\begin{pmatrix} -\\frac{\\partial^2L}{\\partial \\beta_1^2} & -\\frac{\\partial^2L}{\\partial \\beta_1\\beta_2} & \\cdots & -\\frac{\\partial^2L}{\\partial \\beta_1\\beta_{p+1}}\\\\\n",
    "-\\frac{\\partial^2L}{\\partial \\beta_2\\beta_{1}} & -\\frac{\\partial^2L}{\\partial \\beta_2^2} & \\cdots & -\\frac{\\partial^2L}{\\partial \\beta_2\\beta_{p+1}}\\\\\n",
    "\\vdots & \\ddots & \\ddots & \\vdots\\\\\n",
    "-\\frac{\\partial^2L}{\\partial \\beta_{p+1}\\beta_{1}} &\\cdots & \\cdots & -\\frac{\\partial^2L}{\\partial \\beta_{p+1}^2}\n",
    "\\end{pmatrix}^{-1}$$\n",
    "\n",
    "Podemos definir \n",
    "\n",
    "$$X=\\begin{pmatrix} 1 & x_{11} & x_{12}&\\cdots&x_{1p}\\\\\n",
    "1& x_{21} & x_{22} &\\cdots &x_{2p}\\\\\n",
    "\\vdots & \\vdots &\\vdots&\\ddots&\\vdots\\\\\n",
    "1 & x_{n1}&x_{n2} & \\cdots & x_{np}\\end{pmatrix}$$\n",
    "\n",
    "y la matriz \n",
    "\n",
    "$$\\hat{V} = \\begin{pmatrix} \\hat{\\pi}_1(1-\\hat{\\pi}_1) & 0 &\\cdots &0\\\\\n",
    "0 &  \\hat{\\pi}_2(1-\\hat{\\pi}_2)&\\cdots & 0\\\\\n",
    "\\vdots  &0 &\\ddots & \\vdots\\\\\n",
    "0 &\\cdots & 0 & \\hat{\\pi}_n(1-\\hat{\\pi}_n)\\end{pmatrix}$$\n",
    "\n",
    "donde $\\pi_i = \\pi(x_i)$\n",
    "\n",
    "Luego podemos estimar la matriz de información observada como \n",
    "\n",
    "$$\\hat{I}(\\hat{\\beta}) = X^T\\hat{V}X$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17575f0a-1900-407c-8982-dfe327db3b50",
   "metadata": {},
   "source": [
    "## Ejemplo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e719616e-9aec-4900-8446-2a9c8eb2e153",
   "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>sub_id</th>\n",
       "      <th>site_id</th>\n",
       "      <th>phy_id</th>\n",
       "      <th>priorfrac</th>\n",
       "      <th>age</th>\n",
       "      <th>weight</th>\n",
       "      <th>height</th>\n",
       "      <th>bmi</th>\n",
       "      <th>premeno</th>\n",
       "      <th>momfrac</th>\n",
       "      <th>armassist</th>\n",
       "      <th>smoke</th>\n",
       "      <th>raterisk</th>\n",
       "      <th>fracscore</th>\n",
       "      <th>fracture</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>14</td>\n",
       "      <td>No</td>\n",
       "      <td>62</td>\n",
       "      <td>70.3</td>\n",
       "      <td>158</td>\n",
       "      <td>28.16055</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>Same</td>\n",
       "      <td>1</td>\n",
       "      <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>284</td>\n",
       "      <td>No</td>\n",
       "      <td>65</td>\n",
       "      <td>87.1</td>\n",
       "      <td>160</td>\n",
       "      <td>34.02344</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>Same</td>\n",
       "      <td>2</td>\n",
       "      <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>6</td>\n",
       "      <td>305</td>\n",
       "      <td>Yes</td>\n",
       "      <td>88</td>\n",
       "      <td>50.8</td>\n",
       "      <td>157</td>\n",
       "      <td>20.60936</td>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>Yes</td>\n",
       "      <td>No</td>\n",
       "      <td>Less</td>\n",
       "      <td>11</td>\n",
       "      <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>6</td>\n",
       "      <td>309</td>\n",
       "      <td>No</td>\n",
       "      <td>82</td>\n",
       "      <td>62.1</td>\n",
       "      <td>160</td>\n",
       "      <td>24.25781</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>Less</td>\n",
       "      <td>5</td>\n",
       "      <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>37</td>\n",
       "      <td>No</td>\n",
       "      <td>61</td>\n",
       "      <td>68.0</td>\n",
       "      <td>152</td>\n",
       "      <td>29.43213</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>Same</td>\n",
       "      <td>1</td>\n",
       "      <td>No</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   sub_id  site_id  phy_id priorfrac  age  weight  height       bmi premeno  \\\n",
       "1       1        1      14        No   62    70.3     158  28.16055      No   \n",
       "2       2        4     284        No   65    87.1     160  34.02344      No   \n",
       "3       3        6     305       Yes   88    50.8     157  20.60936      No   \n",
       "4       4        6     309        No   82    62.1     160  24.25781      No   \n",
       "5       5        1      37        No   61    68.0     152  29.43213      No   \n",
       "\n",
       "  momfrac armassist smoke raterisk  fracscore fracture  \n",
       "1      No        No    No     Same          1       No  \n",
       "2      No        No    No     Same          2       No  \n",
       "3     Yes       Yes    No     Less         11       No  \n",
       "4      No        No    No     Less          5       No  \n",
       "5      No        No    No     Same          1       No  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "GLOW = dset('glow500')\n",
    "GLOW.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8cb1dcd1-1fd8-46a9-993d-a7756ebc3055",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import LabelEncoder\n",
    "from sklearn.preprocessing import OneHotEncoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "bbcda603-6d18-4129-9dcd-1258c65dc43d",
   "metadata": {},
   "outputs": [],
   "source": [
    "labelencoder = LabelEncoder()\n",
    "GLOW[\"priorfrac\"]=labelencoder.fit_transform(GLOW[\"priorfrac\"])\n",
    "GLOW[\"premeno\"]=labelencoder.fit_transform(GLOW[\"premeno\"])\n",
    "GLOW[\"momfrac\"]=labelencoder.fit_transform(GLOW[\"momfrac\"])\n",
    "GLOW[\"armassist\"]=labelencoder.fit_transform(GLOW[\"armassist\"])   \n",
    "GLOW[\"smoke\"]=labelencoder.fit_transform(GLOW[\"smoke\"])\n",
    "GLOW[\"fracture\"]=labelencoder.fit_transform(GLOW[\"fracture\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a4c24f5c-2aa1-4cd5-88b4-753489fd1350",
   "metadata": {},
   "outputs": [],
   "source": [
    "onehotencoder = OneHotEncoder()\n",
    "ratk = onehotencoder.fit_transform(GLOW[\"raterisk\"].to_numpy().reshape(-1,1)).toarray()\n",
    "\n",
    "GLOW[\"raterisk_1\"] = ratk[:,0]\n",
    "GLOW[\"raterisk_2\"] = ratk[:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "06bb4f04-6a2c-4b98-b8f8-b8587bbd227d",
   "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>sub_id</th>\n",
       "      <th>site_id</th>\n",
       "      <th>phy_id</th>\n",
       "      <th>priorfrac</th>\n",
       "      <th>age</th>\n",
       "      <th>weight</th>\n",
       "      <th>height</th>\n",
       "      <th>bmi</th>\n",
       "      <th>premeno</th>\n",
       "      <th>momfrac</th>\n",
       "      <th>armassist</th>\n",
       "      <th>smoke</th>\n",
       "      <th>raterisk</th>\n",
       "      <th>fracscore</th>\n",
       "      <th>fracture</th>\n",
       "      <th>raterisk_1</th>\n",
       "      <th>raterisk_2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>14</td>\n",
       "      <td>0</td>\n",
       "      <td>62</td>\n",
       "      <td>70.3</td>\n",
       "      <td>158</td>\n",
       "      <td>28.16055</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>Same</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>284</td>\n",
       "      <td>0</td>\n",
       "      <td>65</td>\n",
       "      <td>87.1</td>\n",
       "      <td>160</td>\n",
       "      <td>34.02344</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>Same</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>6</td>\n",
       "      <td>305</td>\n",
       "      <td>1</td>\n",
       "      <td>88</td>\n",
       "      <td>50.8</td>\n",
       "      <td>157</td>\n",
       "      <td>20.60936</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>Less</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>6</td>\n",
       "      <td>309</td>\n",
       "      <td>0</td>\n",
       "      <td>82</td>\n",
       "      <td>62.1</td>\n",
       "      <td>160</td>\n",
       "      <td>24.25781</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>Less</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>37</td>\n",
       "      <td>0</td>\n",
       "      <td>61</td>\n",
       "      <td>68.0</td>\n",
       "      <td>152</td>\n",
       "      <td>29.43213</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>Same</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   sub_id  site_id  phy_id  priorfrac  age  weight  height       bmi  premeno  \\\n",
       "1       1        1      14          0   62    70.3     158  28.16055        0   \n",
       "2       2        4     284          0   65    87.1     160  34.02344        0   \n",
       "3       3        6     305          1   88    50.8     157  20.60936        0   \n",
       "4       4        6     309          0   82    62.1     160  24.25781        0   \n",
       "5       5        1      37          0   61    68.0     152  29.43213        0   \n",
       "\n",
       "   momfrac  armassist  smoke raterisk  fracscore  fracture  raterisk_1  \\\n",
       "1        0          0      0     Same          1         0         0.0   \n",
       "2        0          0      0     Same          2         0         0.0   \n",
       "3        1          1      0     Less         11         0         0.0   \n",
       "4        0          0      0     Less          5         0         0.0   \n",
       "5        0          0      0     Same          1         0         0.0   \n",
       "\n",
       "   raterisk_2  \n",
       "1         1.0  \n",
       "2         1.0  \n",
       "3         0.0  \n",
       "4         0.0  \n",
       "5         1.0  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "GLOW.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e469c8a6-0dd8-4bf7-8d97-9f3d48431bd0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.518075\n",
      "         Iterations 6\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:               fracture   No. Observations:                  500\n",
      "Model:                          Logit   Df Residuals:                      493\n",
      "Method:                           MLE   Df Model:                            6\n",
      "Date:                Thu, 08 Jul 2021   Pseudo R-squ.:                 0.07871\n",
      "Time:                        14:21:37   Log-Likelihood:                -259.04\n",
      "converged:                       True   LL-Null:                       -281.17\n",
      "Covariance Type:            nonrobust   LLR p-value:                 6.565e-08\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -5.6058      1.221     -4.592      0.000      -7.998      -3.213\n",
      "age            0.0501      0.013      3.737      0.000       0.024       0.076\n",
      "weight         0.0041      0.007      0.589      0.556      -0.009       0.018\n",
      "priorfrac      0.6795      0.242      2.803      0.005       0.204       1.155\n",
      "premeno        0.1870      0.277      0.676      0.499      -0.355       0.729\n",
      "raterisk_1     0.8741      0.289      3.023      0.003       0.307       1.441\n",
      "raterisk_2     0.5345      0.276      1.938      0.053      -0.006       1.075\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "from statsmodels.tools.tools import add_constant\n",
    "\n",
    "log_reg = sm.Logit(GLOW[\"fracture\"],add_constant(GLOW[[\"age\",\"weight\",\"priorfrac\",\"premeno\",\"raterisk_1\",\"raterisk_2\"]])).fit()\n",
    "print(log_reg.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f072b1f-794b-494d-a400-a0080e121c5d",
   "metadata": {},
   "source": [
    "## Testeo de significancia de las variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8959df07-d299-48e2-bfa6-0e0b9b480bb7",
   "metadata": {},
   "source": [
    "El test de significancia sigue la misma logica del caso univariado. Solo que ahora la $G$ tiene distribución chi cuadrado con $p$ grados de libertad. El test ratio test ahora compara el modelo con todos los parametros con el modelo con solo la constante.\n",
    "\n",
    "## Continuación ejemplo \n",
    "\n",
    "Mediante la tabla anterior obtenemos\n",
    "\n",
    "$$LikelihoodRatioTest= -2[-281.17-(-259.04)] = 44.2598$$\n",
    "\n",
    "con un p-value $6.565\\cdot 10^{-8}$, luego las variables son significantes. Lo que nos dice que al menos una de las variables no es cero. \n",
    "\n",
    "Si bien esto no es muy util, ahora podemos usar el wald test parámetro por parámetro\n",
    "\n",
    "$$W_j = \\frac{\\hat{\\beta}_j}{\\hat{SE}(\\hat{\\beta})_j}$$\n",
    "\n",
    "Estos valores pueden ser obtenidos directo desde la tabla que retorna statsmodel junto a los p-values.\n",
    "\n",
    "Siguiendo con el ejemplo podemos notar que $WEIGHT$ y $PREMENO$ no son significantes (tomando un nivel de significancia de 0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16a76fef-5e6b-4537-ba20-448cfac13694",
   "metadata": {},
   "source": [
    "Ahora ajustaremos el modelo sin las variables no significantes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ad0dc658-e0da-4236-bae3-162cf3f55e5e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.518899\n",
      "         Iterations 6\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:               fracture   No. Observations:                  500\n",
      "Model:                          Logit   Df Residuals:                      495\n",
      "Method:                           MLE   Df Model:                            4\n",
      "Date:                Thu, 08 Jul 2021   Pseudo R-squ.:                 0.07724\n",
      "Time:                        14:22:07   Log-Likelihood:                -259.45\n",
      "converged:                       True   LL-Null:                       -281.17\n",
      "Covariance Type:            nonrobust   LLR p-value:                 8.400e-09\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -4.9906      0.903     -5.529      0.000      -6.760      -3.221\n",
      "age            0.0459      0.012      3.690      0.000       0.022       0.070\n",
      "priorfrac      0.7002      0.241      2.904      0.004       0.228       1.173\n",
      "raterisk_1     0.8658      0.286      3.025      0.002       0.305       1.427\n",
      "raterisk_2     0.5486      0.275      1.995      0.046       0.010       1.088\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "log_reg_2 = sm.Logit(GLOW[\"fracture\"],add_constant(GLOW[[\"age\",\"priorfrac\",\"raterisk_1\",\"raterisk_2\"]])).fit()\n",
    "print(log_reg_2.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "078bd92c-d317-4c6a-ac96-43f6cd770929",
   "metadata": {},
   "source": [
    "Si hacemos un likelihood ratio test entre los dos modelos (chi cuadrado con 2 grados de libertad) obtenemos\n",
    "\n",
    "$$G = -2[-259.45 - (-259.04)] = 0.82$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "fc8a0128-4f61-40ee-8709-60548919d141",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6636502501363195"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import chi2\n",
    "\n",
    "1-chi2.cdf(0.82,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f73895d-050b-4627-8ba8-21d409b59974",
   "metadata": {},
   "source": [
    "Ya que este valor es mayor a 0.05 concluimos que el modelo con todos los parametros no es mejor que el modelo sin los parametros"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd376e53-04fe-4d24-9783-ec35914ac8ff",
   "metadata": {},
   "source": [
    "## Wald test multivariado\n",
    "\n",
    "El Wald-test multivariado está dado por \n",
    "\n",
    "$$W = \\hat{\\beta}^T[\\hat{\\mathbb{V}}(\\hat{\\beta})]^{-1}\\hat{\\beta}$$\n",
    "\n",
    "este test distribuye chi cuadrado con $p+1$ grados de libertad, bajo la hipotesis de que los $p+1$ coeficientes son ceros. \n",
    "\n",
    "Es equivalente al likelihood ratio test si eliminamos la constante, pero es mucho mas costoso ya que hay que invertir una matriz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09f5533c-4b40-43a1-b723-1b03a2579106",
   "metadata": {},
   "source": [
    "## Intervalos de confianza\n",
    "\n",
    "El intervalo de confianza basado en el Wald test es el mismo visto anteriormente."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22e50a2f-a2ea-4ee3-b0ee-852a555d8971",
   "metadata": {},
   "source": [
    "## Intervalo de confianza para el logit\n",
    "\n",
    "El estimador del logit es \n",
    "\n",
    "$$\\hat{g}(x) = \\hat{\\beta}_0 + \\hat{\\beta}_1 x_1+\\cdots+\\hat{\\beta}_px_p$$\n",
    "\n",
    "Reescribimos \n",
    "\n",
    "$$\\hat{g}(x) = x^T\\hat{\\beta}$$\n",
    "\n",
    "Luego el estimador de la varianza es \n",
    "\n",
    "$$\\hat{\\mathbb{V}}(\\hat{g}(x))=\\sum_{i=0}^p x_i^2\\hat{\\mathbb{V}}(\\hat{\\beta}_i)+\\sum_{j=0}^p\\sum_{i=j+1}^p 2x_{i}x_{j}\\hat{cov}(\\hat{\\beta}_j,\\hat{\\beta}_i)$$\n",
    "\n",
    "En forma matricial queda \n",
    "\n",
    "$$\\hat{\\mathbb{V}}(\\hat{g}(x)) = x^T\\hat{\\mathbb{V}}(\\hat{\\beta})x = x^T(X^T\\hat{V}X)^{-1}x$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc798752-6c68-4789-9c37-a85aa1b4bb0f",
   "metadata": {},
   "source": [
    "Mediante statsmodels podemos obtener la matriz de covarianza "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "7e7da3f0-980b-41d0-9e39-ba95170a92ba",
   "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>const</th>\n",
       "      <th>age</th>\n",
       "      <th>priorfrac</th>\n",
       "      <th>raterisk_1</th>\n",
       "      <th>raterisk_2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>const</th>\n",
       "      <td>0.81487</td>\n",
       "      <td>-0.01089</td>\n",
       "      <td>0.04450</td>\n",
       "      <td>-0.08055</td>\n",
       "      <td>-0.06039</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>-0.01089</td>\n",
       "      <td>0.00015</td>\n",
       "      <td>-0.00083</td>\n",
       "      <td>0.00054</td>\n",
       "      <td>0.00022</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>priorfrac</th>\n",
       "      <td>0.04450</td>\n",
       "      <td>-0.00083</td>\n",
       "      <td>0.05816</td>\n",
       "      <td>-0.01184</td>\n",
       "      <td>-0.00313</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>raterisk_1</th>\n",
       "      <td>-0.08055</td>\n",
       "      <td>0.00054</td>\n",
       "      <td>-0.01184</td>\n",
       "      <td>0.08191</td>\n",
       "      <td>0.04624</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>raterisk_2</th>\n",
       "      <td>-0.06039</td>\n",
       "      <td>0.00022</td>\n",
       "      <td>-0.00313</td>\n",
       "      <td>0.04624</td>\n",
       "      <td>0.07563</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              const      age  priorfrac  raterisk_1  raterisk_2\n",
       "const       0.81487 -0.01089    0.04450    -0.08055    -0.06039\n",
       "age        -0.01089  0.00015   -0.00083     0.00054     0.00022\n",
       "priorfrac   0.04450 -0.00083    0.05816    -0.01184    -0.00313\n",
       "raterisk_1 -0.08055  0.00054   -0.01184     0.08191     0.04624\n",
       "raterisk_2 -0.06039  0.00022   -0.00313     0.04624     0.07563"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "round(log_reg_2.cov_params(),5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85fe1fd9-996d-4518-883c-f88171769b8b",
   "metadata": {},
   "source": [
    "## Ejemplo\n",
    "\n",
    "intervalo de confianza de tener Edad = 65, Priorfrac = 1, Raterisk_1=1, Raterisk_2=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "f4a8ff88-3b49-4f76-87c3-cdd43d189d49",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.039509999999999976"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([1,65,1,1,0]).T.dot(round(log_reg_2.cov_params(),5)).dot(np.array([1,65,1,1,0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "5727c890-af6b-4528-a315-b770d309fc43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.04937000000000011"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.81487+(65)**2*0.00015+0.05816 + 0.07563 + 2*65*(-0.01089) + 2*0.04450 + 2*(-0.06039) + 2*65*(-0.00083) + 2*65*0.00022 + 2*(-0.00313)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f4e65c5-1a8c-4f83-82bc-9c873babdc85",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}