Estadistica-Orange.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Librerías manejo de datos\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# Librerías de gráficos\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Carga de datos y ETL\n",
    "\n",
    "Vamos a cargar un dataset sobre el peso de los niños al nacer.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#birth_weight = pd.read_csv('birthweight.csv')\n",
    "birth_weight = pd.read_csv('https://github.com/calabozo/orange-stats/raw/master/birthweight.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "birth_weight.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "birth_weight[\"Birthweight_kg\"]=birth_weight[\"Birthweight\"]*0.453592"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Media\n",
    "\n",
    "\\\\[ \\mu_{x} = E[x] =  \\frac{1}{N} \\sum _{i=1}^N x_i\\\\]\n",
    "\n",
    "Rellenar con formula y texto"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean(vector):\n",
    "    pass\n",
    "mean(birth_weight[\"Birthweight_kg\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(birth_weight[\"Birthweight_kg\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Varianza\n",
    "\n",
    "Es la esperanza del cuadrado de la desviación de dicha variable respecto a su media. Otra definición es la media de los residuos al cuadrado.\n",
    "\n",
    "Su formula es:\n",
    "\n",
    "\\\\[ Var[x] =  \\frac{1}{N} \\sum _{i=1}^N \\left( x_i-\\bar{x} \\right)^2 \\\\]\n",
    "\n",
    "Pero como no sabremos $\\bar{x}$ sino que siempre tendremoso que estimar la media E[x] la formula del estimador insesgado de la varianza para una muestra de la población es:\n",
    "\n",
    "\\\\[ Var[x] =  \\frac{1}{N-1} \\sum _{i=1}^N \\left( x_i-E[x] \\right)^2 \\\\]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def var(vector):\n",
    "    pass\n",
    "var(birth_weight[\"Birthweight_kg\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.var(birth_weight[\"Birthweight_kg\"],ddof=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mediana\n",
    "\n",
    "La mediana es el valor de una distribución tal que la mitad de los elementos son superiores y la mitad son inferiores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def median(vector: pd.Series) -> float:\n",
    "    pass\n",
    "    \n",
    "print(median(pd.Series([1,2,3,3.5,4,6])))\n",
    "print(median(pd.Series([1,2,3,4,6])))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "median(birth_weight[\"Birthweight_kg\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.median(birth_weight[\"Birthweight_kg\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cuantiles\n",
    "\n",
    "Los cuantiles suelen usarse por grupos que dividen la distribución en partes iguales; entendidas estas como intervalos que comprenden la misma proporción de valores. Los más usados son:\n",
    "\n",
    "* Los cuartiles, que dividen a la distribución en cuatro partes (corresponden a los cuantiles 0,25; 0,50 y 0,75);\n",
    "* Los quintiles, que dividen a la distribución en cinco partes (corresponden a los cuantiles 0,20; 0,40; 0,60 y 0,80);\n",
    "* Los deciles, que dividen a la distribución en diez partes;\n",
    "* Los percentiles, que dividen a la distribución en cien partes.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.quantile(birth_weight[\"Birthweight_kg\"],[0.25,0.5,0.75])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boxplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.boxplot(x=birth_weight[\"Birthweight_kg\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Histograma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(birth_weight[\"Birthweight_kg\"], binwidth=0.5,label=\"0.5 binwidth\")\n",
    "sns.histplot(birth_weight[\"Birthweight_kg\"], bins=20,color=\"red\",label=\"20 bins\")\n",
    "plt.xlabel(\"Weight [kg]\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Curva de densidad (kernel density estimation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.kdeplot(birth_weight[\"Birthweight_kg\"], bw_adjust=0.5,label=\"bw_adjust=0.5\")\n",
    "sns.kdeplot(birth_weight[\"Birthweight_kg\"], bw_adjust=2,color=\"red\",label=\"bw_adjust=1\")\n",
    "plt.xlabel(\"Weight [kg]\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Violin plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.violinplot(data=birth_weight,x=\"Birthweight_kg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tests estadísticos\n",
    "\n",
    "¿Es la media de estas dos poblaciones iguales?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weight_NO_smoke = birth_weight.loc[birth_weight[\"smoker\"]==0,\"Birthweight_kg\"]\n",
    "weight_YES_smoke = birth_weight.loc[birth_weight[\"smoker\"]==1,\"Birthweight_kg\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weight_NO_smoke.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weight_YES_smoke.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.kdeplot(weight_NO_smoke,label=\"No fumadora\")\n",
    "sns.kdeplot(weight_YES_smoke,color=\"red\",label=\"Fumadora\")\n",
    "plt.xlabel(\"Weight [kg]\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_NO_YES=np.mean(weight_NO_smoke.values)-np.mean(weight_YES_smoke.values)\n",
    "diff_YES_NO=np.mean(weight_YES_smoke.values)-np.mean(weight_NO_smoke.values)\n",
    "print(f\"La diferencia de medias NO_smoker-YES_smoker es:{np.round(diff_NO_YES,2)}\")\n",
    "print(f\"La diferencia de medias YES_smoker-NO_smoker es:{np.round(diff_YES_NO,2)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "¿Cómo de probable es que estos cambios hayan ocurrido al azar?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "X1=weight_NO_smoke\n",
    "X2=weight_YES_smoke\n",
    "\n",
    "def perm_func(X1: pd.Series, X2:pd.Series) -> (pd.Series,pd.Series): \n",
    "    X = X1.append(X2).reset_index(drop=True)\n",
    "\n",
    "    X1_reorder=X.sample(frac=0.5)\n",
    "    X2_reorder=X.drop(X1_reorder.index)\n",
    "    return (X1_reorder.reset_index(drop=True),X2_reorder.reset_index(drop=True))\n",
    "\n",
    "X1_reorder, X2_reorder= perm_func(weight_NO_smoke,weight_YES_smoke)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_values = []\n",
    "num_mean_X1_larger_X2 = 0\n",
    "num_mean_X2_larger_X1 = 0\n",
    "\n",
    "total_trials = 10000\n",
    "for i in range(total_trials):\n",
    "    X1_reorder, X2_reorder= perm_func(weight_NO_smoke,weight_YES_smoke)\n",
    "    \n",
    "    mean_diff = np.mean(X1_reorder.values) - np.mean(X2_reorder.values)\n",
    "    mean_values.append(mean_diff)\n",
    "    \n",
    "    if mean_diff>diff_NO_YES:\n",
    "        num_mean_X1_larger_X2 += 1\n",
    "    \n",
    "    if mean_diff<diff_YES_NO:\n",
    "        num_mean_X2_larger_X1 += 1\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"\"\"La probabilidad por puro azar la media de X1 menos X2 sea mayor que la diferencia original: \n",
    "    {np.round(num_mean_X1_larger_X2/total_trials*100,2)}%\"\"\")\n",
    "\n",
    "print(f\"\"\"La probabilidad por puro azar la media de X1 menos X2 sea menor que la diferencia original: \n",
    "    {np.round(num_mean_X2_larger_X1/total_trials*100,2)}%\"\"\")\n",
    "\n",
    "print(f\"\"\"La probabilidad de que por puro azar los valores de NO_smoker y YES_smoker se hayan reordenado\n",
    "para dar diferencias tan grandes es:{np.round((num_mean_X2_larger_X1+num_mean_X1_larger_X2)/total_trials*100,2)}%\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "#sns.displot(mean_values)\n",
    "sns.kdeplot(mean_values)\n",
    "plt.axvline(diff_NO_YES, 0,2,color=\"#FF0000\")\n",
    "plt.axvline(diff_YES_NO, 0,2,color=\"#FF0000\")\n",
    "x=np.linspace(-1, 1, 100)\n",
    "plt.plot( x, norm.pdf(x, 0, np.std(mean_values,ddof=1)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.var(mean_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vemos que la probabilidad de que esos datos se hayan conseguido por puro hazar es muy baja, inferior al 5%-"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### T-student"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El t-test es usado cuando tienes dos grupos de datos y quieres saber si la media de ambos grupos es igual o no. Tenemos dos hipótesis:\n",
    "\\\\[H_0 : \\mu(X) = \\mu(Y) \\\\]\n",
    "\\\\[H_1 : \\mu(X) \\ne \\mu(Y) \\\\]\n",
    "\n",
    "Empezamos con las siguientes suposiciones:\n",
    "1. Las muestras han sido seleccionas por muestreo aleatorio simple dentro de la población.\n",
    "2. Las poblaciones siguen una distribución normal. (sino usar test Wilcoxon-Mann-Whitney)\n",
    "3. Ambas poblaciones tienen la misma varianza. (sino usar test Welch)\n",
    "4. Las poblaciones no están correladas. (sino usar t-test emparejado)\n",
    "\n",
    "Primero calculamos la media y varianza de ambos grupos. El test para verificar si la hipótesis nula $H_0$ es cierta puede ser calculado como sigue:\n",
    "\n",
    "\\\\[ t=\\frac{E[X]-E[Y]}{s_p·\\sqrt{\\frac{1}{n_x}+\\frac{1}{n_y}}} \\\\]\n",
    "\n",
    "Donde $s_p$ es la desviación estandar compuesta, calculada como:\n",
    "\n",
    "\\\\[ s^2_p=\\frac{(n_x-1)Var[X]+(n_y-1)Var[Y]}{n_x+n_y-2} \\\\]\n",
    "\n",
    "Donde $n_x$ y $n_y$ son los números de muestras en cada uno de los grupos muestreados. El número de grados de libertad es \n",
    "$d.f.=n_x+n_y-2$. Podemos asumir que las dos varianzas son iguales si ambas pasan el test de varianza de Fisher(F-test).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "t, p = stats.ttest_ind(weight_NO_smoke,weight_YES_smoke)\n",
    "p\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ejemplo practico\n",
    "\n",
    "¿Cuantos caramelos hay en estos grupos?\n",
    "\n",
    "![](https://github.com/calabozo/orange-stats/raw/master/caramelos.jpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regresiones lineales\n",
    "\n",
    "Modelo matemático usado para aproximar la relación de dependencia entre una variable dependiente $Y$, diferentes variables independientes definidas en la matriz $X$ y un término aleatorio $\\varepsilon$. Este modelo puede ser expresado como:\n",
    "\n",
    "\\\\[\n",
    "\\hat{Y}=\\beta_1 X_1+\\beta_2 X_2+\\cdots +\\beta_p X_p = \\sum \\beta_k X_k\n",
    "\\\\]\n",
    "\n",
    "\n",
    "Tratamos de encontrar la predicción $\\hat{Y}$ tal que minimize el error cuadrático medio:\n",
    "\\\\[\n",
    "MSE = {1 \\over n} \\sum_{i=0}^n{(Y-\\hat{Y})^2}\n",
    "\\\\]\n",
    "\n",
    "\n",
    "\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/thumb/3/3a/Linear_regression.svg/400px-Linear_regression.svg.png)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "train,test = train_test_split(birth_weight,test_size=0.33, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.formula.api import ols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = ols(formula='Birthweight_kg ~ mheight + motherage + fheight + fage + mnocig', data=birth_weight)\n",
    "res = model.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = ols(formula='Birthweight_kg ~ mheight+ mnocig', data=birth_weight)\n",
    "res = model.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clasificación: Regresión logística"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En lugar de realizar una predicción de un valor queremos hacer un clasificador.\n",
    "Podemos tratar de asignar una probabilidad. Pero hay un problema porque la regresión lineal va entre 0 y 1.\n",
    "\n",
    "Para ello transformarmos la regresión lineal mediante la función logit:\n",
    "\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/thumb/c/c8/Logit.svg/350px-Logit.svg.png)\n",
    "\n",
    "#### función de enlace (link function)\n",
    "\n",
    "Para pasar del dominio de números reales $(-\\infty,\\infty)$ al de probabilidades $[0,1]$ a vamos a utilizar la **función logística**:\n",
    "\\\\[\n",
    "p = h(x)=  \\frac{1}{1+e^{-x}}\n",
    "\\\\]\n",
    "Su inversa se conoce como la función **logit**:\n",
    "\\\\[\n",
    "h^{-1}(p) = log \\left( \\frac{p}{1-p} \\right)\n",
    "\\\\]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Es decir, cuando estemos trabajando con una **distribución binomial** un modelo lineal del tipo:\n",
    "\\\\[\n",
    "y = \\beta \\vec{x}+\\beta_0\n",
    "\\\\]\n",
    "lo podemos trasnformar en:\n",
    "\\\\[\n",
    "y = p(x) = \\frac{1}{1+e^{-\\beta \\vec{x}-\\beta_0}} \n",
    "\\\\]\n",
    "Ahora $p(x)$ es una función que muestra valores en el rango $[0,1]$, puede ser considerada como una probabilidad.\n",
    "\n",
    "Y definiendo un umbral podríamos tener el siguiente clasificador:\n",
    "* Seleccionamos clase 1 si p(x)>=0.5\n",
    "* Seleccionamos clase 0 si p(x)< 0.5\n",
    "\n",
    "\n",
    "\n",
    "Es decir, tenemos una probabilidad, su valor está en el rango $[0,1]$:\n",
    "\\\\[\n",
    "    p = \\frac{1}{1-e^{-\\hat{Y}}}= \\frac{1}{1-e^{-(\\beta_1 X_1+\\beta_2 X_2+\\cdots +\\beta_p X_p)}}   \n",
    "\\\\]\n",
    "\n",
    "Definimos la razón de monomios (Odds ratio) como el cociente entre dos probabilidades, su valor está en el rango $[0,\\infty]$:\n",
    "\n",
    "\\\\[\n",
    " Odds = \\frac{p}{1-p}=\\frac{\\frac{1}{1-e^{-(\\beta_1 X_1+\\beta_2 X_2+\\cdots +\\beta_p X_p)}}}{\\frac{e^{-(\\beta_1 X_1+\\beta_2 X_2+\\cdots +\\beta_p X_p)}}{1-e^{-(\\beta_1 X_1+\\beta_2 X_2+\\cdots +\\beta_p X_p)}}}=e^{(\\beta_1 X_1+\\beta_2 X_2+\\cdots +\\beta_p X_p)}\n",
    "\\\\]\n",
    "\n",
    "Si aplicamos el logaritmo a la razón de monomios tenemos un valor que está en el rango $[-\\infty,\\infty]$:\n",
    "\\\\[\n",
    " log(Odds)= log \\left(\\frac{p}{1-p} \\right) = \\beta_1 X_1+\\beta_2 X_2+\\cdots +\\beta_p X_p\n",
    "\\\\]\n",
    "\n",
    "La función de coste que vamos a tratar de minimizar será:\n",
    "\\\\[\n",
    "\\begin{split}\n",
    "Cost(p(x),y) &= {1 \\over n} \\sum_{i=0}^n{(y-\\hat{y})^2}\\\\\n",
    "Cost(p(x),y) &= {1 \\over n} \\sum_{i=0}^n{(y-p(x))^2}\n",
    "\\end{split}\n",
    "\\\\]\n",
    "Que transformaremos en:\n",
    "\\\\[\n",
    "Cost(p(x),y) = -y ·log(p(x))-(1-y)·log(1-p(x))\n",
    "\\\\]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo:Churn rate\n",
    "\n",
    "Vamos a utilizar un dataset publicado por IBM en [kaggle](https://www.kaggle.com/blastchar/telco-customer-churn).\n",
    "\n",
    "\n",
    "En este ejemplo vamos a cargar el dataset proporcionado y ver si somos capaces de ver qué usuarios son los que corren más riesgo de irse.\n",
    "\n",
    "El conjunto de datos incluye información sobre:\n",
    "\n",
    "* Clientes que se fueron en el último mes: la columna se llama Churn\n",
    "* Servicios para los que se ha registrado cada cliente: teléfono, líneas múltiples, Internet, seguridad en línea, copia de seguridad en línea, protección de dispositivos, soporte técnico y transmisión de TV y películas\n",
    "* Información de la cuenta del cliente: cuánto tiempo han sido cliente (columna tenure), contrato, método de pago, facturación electrónica, cargos mensuales y cargos totales\n",
    "* Información demográfica sobre los clientes: sexo, rango de edad y si tienen socios y dependientes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_churn=pd.read_csv(\"WA_Fn-UseC_-Telco-Customer-Churn.csv\")\n",
    "df_churn=pd.read_csv(\"https://github.com/calabozo/orange-stats/raw/master/WA_Fn-UseC_-Telco-Customer-Churn.csv\")\n",
    "df_churn.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_churn.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_churn=df_churn.drop(\n",
    "    columns=[\"customerID\",\"OnlineSecurity\",\"OnlineBackup\",\n",
    "             \"DeviceProtection\",\"TechSupport\",\"StreamingTV\",\"StreamingMovies\",\"TotalCharges\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_churn.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.formula.api import logit\n",
    "# https://www.statsmodels.org/stable/generated/statsmodels.formula.api.logit.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_churn[\"Churn\"].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_churn[\"Churn\"]=df_churn[\"Churn\"].astype('category')\n",
    "df_churn[\"Churn_id\"]=df_churn[\"Churn\"].cat.codes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "df_churn_train,df_churn_test = train_test_split(df_churn,test_size=0.33, random_state=42)\n",
    "df_churn_test=df_churn_test.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "churn_model=logit(formula=\"Churn_id ~ tenure+C(InternetService)+MonthlyCharges+C(PaymentMethod)\",\n",
    "                  data=df_churn_train).fit()\n",
    "churn_model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Salida de un modelo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La salida de esta predicción es el logaritmo de la razón de monomios:\n",
    "\\\\[\n",
    " log(Odds)= log \\left(\\frac{p}{1-p} \\right) = \\beta_1 X_1+\\beta_2 X_2+\\cdots +\\beta_p X_p\n",
    "\\\\]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "user_test = df_churn_test.iloc[5:6,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_odds=churn_model.predict(user_test, linear=True).values[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Es {np.round(np.exp(log_odds),2)} veces más probable que user_test haga Churn que permanezca en la compañia\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob_est = churn_model.predict(user_test, linear=False).values[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"\"\"Aunque no se corresponde con la probabilidad teórica, \n",
    "una aproximación a la probabilidad de churn para el usuario test rondaría {np.round(prob_est,2)}\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    return 1/(1+np.exp(-x))\n",
    "\n",
    "prob_est == sigmoid(log_odds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Medidas de calidad de un modelo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnf_mat=pd.crosstab(index=df_churn_test[\"Churn\"],columns=churn_model.predict(df_churn_test)>0.5 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_churn_test[\"pred\"]=churn_model.predict(df_churn_test)>0.5\n",
    "df_churn_test[\"pred\"] = df_churn_test[\"pred\"].replace({True: 'Yes', False: 'No'})\n",
    "df_churn_test=df_churn_test.dropna(subset=[\"Churn\"])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "cnf_mat=confusion_matrix(df_churn_test[\"Churn\"], df_churn_test[\"pred\"], labels=[\"No\",\"Yes\"] )\n",
    "cnf_mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FP = cnf_mat.sum(axis=0) - np.diag(cnf_mat) \n",
    "FN = cnf_mat.sum(axis=1) - np.diag(cnf_mat)\n",
    "TP = np.diag(cnf_mat)\n",
    "TN = cnf_mat.sum() - (FP + FN + TP)\n",
    "\n",
    "# Sensitivity, hit rate, recall, or true positive rate\n",
    "TPR = TP/(TP+FN)\n",
    "print(f\"Recall o True Positive Rate: {TPR}\")\n",
    "# Specificity or true negative rate\n",
    "TNR = TN/(TN+FP) \n",
    "print(f\"True Negative rate: {TPR}\")\n",
    "# Precision or positive predictive value\n",
    "PPV = TP/(TP+FP)\n",
    "print(f\"Precision: {PPV}\")\n",
    "# Negative predictive value\n",
    "NPV = TN/(TN+FN)\n",
    "# Fall out or false positive rate\n",
    "FPR = FP/(FP+TN)\n",
    "# False negative rate\n",
    "FNR = FN/(TP+FN)\n",
    "# False discovery rate\n",
    "FDR = FP/(TP+FP)\n",
    "\n",
    "# Overall accuracy\n",
    "ACC = (TP+TN)/(TP+FP+FN+TN)\n",
    "print(f\"Accuracy: {ACC}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import roc_curve\n",
    "fpr, tpr, thresholds = roc_curve(df_churn_test[\"Churn\"]==\"Yes\", churn_model.predict(df_churn_test))\n",
    "plt.plot(fpr, tpr)\n",
    "plt.title(\"ROC Curve\")\n",
    "plt.xlabel(\"False Positive Rate\")\n",
    "plt.ylabel(\"True Positive Rate\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.metrics import roc_auc_score\n",
    "\n",
    "roc_auc_score(df_churn_test[\"Churn\"], churn_model.predict(df_churn_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mejorando la calidad del modelo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_payment=pd.crosstab(index=df_churn_train[\"Churn\"],columns=df_churn_train[\"PaymentMethod\"])\n",
    "df_payment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "from statsmodels.stats.proportion import proportion_confint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proportion_confint(3, nobs=15, method='binom_test')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_payment.apply(lambda x:proportion_confint(x[0],nobs=x[0]+x[1],method='binom_test'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_churn[\"ElectronicCheck\"]=df_churn.PaymentMethod==\"Electronic check\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_churn.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_churn_train,df_churn_test = train_test_split(df_churn,test_size=0.33, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "churn_new_model=logit(formula=\"Churn_id ~ tenure+C(InternetService)+C(PaperlessBilling)+C(Contract)+MonthlyCharges+C(ElectronicCheck)+C(SeniorCitizen)\",\n",
    "                  data=df_churn_train).fit()\n",
    "churn_new_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fpr_new, tpr_new, thresholds_new = roc_curve(df_churn_test[\"Churn\"]==\"Yes\", churn_new_model.predict(df_churn_test))\n",
    "plt.plot(fpr, tpr)\n",
    "plt.plot(fpr_new, tpr_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}