Regresión logística con Python
Introducción¶
La regresión logística es un método estadístico que trata de modelar la probabilidad de una variable cualitativa binaria (dos posibles valores) en función de una o más variables independientes. La principal aplicación de la regresión logística es la creación de modelos de clasificación binaria.
Se llama regresión logística simple cuando solo hay una variable independiente y regresión logística múltiple cuando hay más de una. Dependiendo del contexto, a la variable modelada se le conoce como variable dependiente o variable respuesta, y a las variables independientes como regresores, predictores o features.
A lo largo de este documento, se describen de forma progresiva los fundamentos teóricos de la regresión logística, los principales aspectos prácticos a tener en cuenta y ejemplos de cómo crear este tipo de modelos en Python.
Modelos de regresión logística en Python
Dos de las implementaciones de modelos de regresión logística más utilizadas en Python son: scikit-learn y statsmodels. Aunque ambas están muy optimizadas, Scikit-learn está orientada principalmente a la predicción, por lo que no dispone de apenas funcionalidades que muestren las muchas características del modelo que se deben analizar para hacer inferencia. Statsmodels es más completa en este sentido.
Definición matemática¶
El modelo de regresión lineal (Legendre, Gauss, Galton y Pearson) considera que, dado un conjunto de observaciones $\{y_i, x_{i1},...,x_{np}\}^{n}_{i=1}$, la media $\mu$ de la variable respuesta $y$ se relaciona de forma lineal con la o las variables regresoras $x_1$ ... $x_p$ acorde a la ecuación:
$$\mu_y = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_p x_{p}$$El resultado de esta ecuación se conoce como la línea de regresión poblacional, y recoge la relación entre los predictores y la media de la variable respuesta.
Otra definición que se encuentra con frecuencia en los libros de estadística es:
$$y_i= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + ... + \beta_p x_{ip} +\epsilon_i$$En este caso, se está haciendo referencia al valor de $y$ para una observación $i$ concreta. El valor de una observación puntual nunca va a ser exactamente igual al promedio, de ahí que se añada el término de error $\epsilon$.
En ambos casos, la interpretación de los elementos del modelo es la misma:
$\beta_{0}$: es la ordenada en el origen, se corresponde con el valor promedio de la variable respuesta $y$ cuando todos los predictores son cero.
$\beta_{p}$: es el efecto promedio que tiene sobre la variable respuesta el incremento en una unidad de la variable predictora $x_{p}$, manteniéndose constantes el resto de variables. Se conocen como coeficientes parciales de regresión.
$e$: es el residuo o error, la diferencia entre el valor observado y el estimado por el modelo. Recoge el efecto de todas aquellas variables que influyen en $y$ pero que no se incluyen en el modelo como predictores.
En la gran mayoría de casos, los valores $\beta_0$ y $\beta_p$ poblacionales se desconocen, por lo que, a partir de una muestra, se obtienen sus estimaciones $\hat{\beta}_0$ y $\hat{\beta}_p$. Ajustar el modelo consiste en estimar, a partir de los datos disponibles, los valores de los coeficientes de regresión que maximizan la verosimilitud (likelihood), es decir, los que dan lugar al modelo que con mayor probabilidad puede haber generado los datos observados.
El método empleado con más frecuencia es el ajuste por mínimos cuadrados ordinarios (OLS), que identifica como mejor modelo la recta (o plano si es regresión múltiple) que minimiza la suma de las desviaciones verticales entre cada dato de entrenamiento y la recta, elevadas al cuadrado.
El término "lineal" en los modelos de regresión hace referencia al hecho de que los parámetros se incorporan en la ecuación de forma lineal, no a que necesariamente la relación entre cada predictor y la variable respuesta tenga que seguir un patrón de recta.
¿Qué ocurre cuando la variable respuesta es binaria, solo dos posibles valores?
Si una variable cualitativa con dos niveles (binaría) se codifica como 0 y 1, matemáticamente es posible ajustar un modelo de regresión lineal por mínimos cuadrados. Sin embargo, esta aproximación tiene dos problemas:
Al generar una recta (hiperplano si hay múltiples variables), se pueden obtener valores predichos distintos de 0 y 1, lo que entra en contradicción con la definición de la variable respuesta binaria.
Si se quiere interpretar las predicciones del modelo como probabilidades de pertenencia a cada clase, no se cumpliría la condición de que toda probabilidad ha de estar dentro del intervalo [0,1], ya que podrían obtenerse valores fuera de este.
Para evitar estos problemas, la regresión logística (David Cox 1958) transforma el valor devuelto por la regresión lineal con una función cuyo resultado siempre está comprendido entre 0 y 1. Existen varias funciones que cumplen esta descripción, una de las más utilizadas es la función logística (también conocida como función sigmoide):
$$\text{sigmoide} = \sigma(y) = \dfrac{1}{1 + e^{-y}} $$Para valores de $y$ positivos muy grandes, $e^{−y}$ es aproximadamente 0, por lo que el valor de la función sigmoide es 1. Para valores de $y$ negativos muy grandes, $e^{−y}$ tiende a infinito, por lo que el valor de la función sigmoide es 0.
Sustituyendo la $y$ de la ecuación anterior por la función de un modelo lineal $ y = \beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}$ se obtiene que:
$$P(y=1|X=x)= \frac{e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}}{1+e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}} \tag{anexo 1}$$donde $P(y=1|X=x)$ puede interpretarse como la probabilidad de que la variable respuesta $y$ adquiera el valor 1 (clase de referencia), dado los predictores $x_1, ..., x_p$.
El modelo resultante tiene los coeficientes de regresión en los exponentes, por lo que no es un modelo lineal y no puede ajustarse con la estrategia descrita al inicio del documento. ¿Cómo evitar este inconveniente?
La expresión obtenida tiene la característica de ser siempre positiva, ya que la función exponencial solo toma valores positivos y, el cociente de valores positivos, es siempre positivo. Esto hace posible aplicarle el logaritmo:
$$ln(\frac{p(y=1|X=x)}{p(y=0|X=x)}) = \beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p} \tag{anexo 2}$$Al realizar la transformación, en el lado derecho se obtiene la ecuación de un modelo lineal. El término de la izquierda resulta ser el logaritmo de un cociente de probabilidades, lo que se conoce como razón de probabilidad (log of odds).
Como resultado de este proceso se consigue convertir un problema de clasificación no lineal, en un problema de regresión lineal que sí puede ajustarse mediante los métodos convencionales.
Una vez obtenidos los coeficientes del modelo ($\beta_0, \beta_1, \cdots, \beta_p$) se puede obtener la probabilidad de que una nueva observación pertenezca a la clase $y=1$ con la ecuación:
$$p(y=1|X=x) = \frac{e^{(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}}{1+e^{(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}}$$Modelar el log of odds es la estrategia matemática que permite encontrar los valores de los coeficientes (ajustar el modelo).
En el apartado anterior, se han obviado dos demostraciones importantes con el objetivo de facilitar la narrativa. Sin embargo, es importante tenerlas en cuenta para entender bien el funcionamiento de este tipo de modelos.
Anexo 1
Combinación de la función sigmoide y la ecuación de un modelo lineal para obtener un modelo que relacione los predictores con la probabilidad de que $y=1$.
$$P(y=1|X=x)=$$$$= \frac{1}{1+e^{-(\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p})}} =$$$$= \frac{1}{\frac{e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}}{e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}} + \frac{1}{e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}}}=$$$$= \frac{1}{\frac{1 + e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}}{e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}}}=$$$$= \frac{e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}}{1+e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}}$$Anexo 2
Transformación logarítmica para linealizar el modelo.
La probabilidad de que $y=1$ dado los predictores $X=x_1, ..., x_p$ es:
$$P(y=1|X=x) = \frac{e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}}{1+e^{\beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}}} = \frac{e^{\sum \beta_i x_{i}}}{1+e^{\sum \beta_i x_{i}}}$$Al ser un caso de clasificación binaria, la probabilidad de $y=0$ es la probabilidad complementaria:
$$P(y=0|X=x) = 1- P(y=0|X=x) = 1 - \frac{e^{\sum \beta_i x_{i}}}{1+e^{\sum \beta_i x_{i}}} = \frac{1}{1+e^{\sum \beta_i x_{i}}}$$Si se divide una probabilidad entre otra:
$$\frac{P(y=1|X=x)}{P(y=0|X=x)} = \frac{\frac{e^{\sum \beta_i x_{i}}}{1+e^{\sum \beta_i x_{i}}}}{\frac{1}{1+e^{\sum \beta_i x_{i}}}} = e^{\sum \beta_i x_{i}}$$$$ln(\frac{P(y=1|X=x)}{P(y=0|X=x)}) = ln(e^{\sum \beta_i x_{i}}) = \sum \beta_i x_{i}$$$$ln(\frac{P(y=1|X=x)}{P(y=0|X=x)}) = \sum \beta_i x_{i}$$$$ln(\frac{P(y=1|X=x)}{P(y=0|X=x)}) = \beta_0 + \beta_1 x_{1} + ... + \beta_p x_{p}$$Odds y logaritmo de odds¶
Debido a las transformaciones necesarias, la regresión logística modela la probabilidad de que la variable respuesta pertenezca al grupo de referencia de forma indirecta mediante el uso del logaritmo de odds. Aunque la mayoría de implementaciones convierten el logaritmo de odds a probabilidades al mostrar las predicciones, conviene conocer su significado.
Supóngase un evento cuya probabilidad de ser verdadero es de 0.8 ($p=0.8$). Dado que es una variable binaria, la probabilidad ser falso es de $1 - 0.8 = 0.2$, ($q=0.2$). Los odds o razón de probabilidad se definen como el ratio entre la probabilidad de evento verdadero y la probabilidad de evento falso $\frac{p}{q}$. En este caso los odds son 0.8 / 0.2 = 4, lo que equivale a decir que se esperan 4 eventos verdaderos por cada evento falso.
Dado que las probabilidades están acotadas siempre al rango [0,1], los odds están acotados entre [0,$\infty$]. Al aplicarles el logaritmo, el rango de valores pasa a ser de [-$\infty$, +$\infty$].
Interpretación del modelo¶
Los principales elementos que hay que interpretar en un modelo de regresión logística son los coeficientes de los predictores:
$\beta_{0}$ es la ordenada en el origen o intercept. Se corresponde con el valor esperado del logaritmo de odds cuando todos los predictores son cero. Puede transformarse a probabilidad con la ecuación $e^{\beta_0} / (1 + e^{\beta_0})$. Tras la transformación, su valor se corresponde con la probabilidad esperada de pertenecer a la clase 1 cuando todos los predictores son cero.
$\beta_{p}$ los coeficientes de regresión parcial de cada predictor indican el cambio promedio del logaritmo de odds al incrementar en una unidad la variable predictora $x_{p}$, manteniéndose constantes el resto de variables. Esto equivale a decir que, por cada unidad que se incrementa $x_{p}$, se multiplican los odds por $e^{\beta_p}$.
Dado que la relación entre $p(y=1)$ y $x$ no es lineal, los coeficientes de regresión $\beta_{p}$ no se corresponden con el cambio en la probabilidad de $y$ asociada con el incremento en una unidad de $x$. Cuánto se incremente la probabilidad de $y$ por unidad de $x$ depende del valor de $x$, es decir, de la posición en la curva logística en la que se encuentre. Esta es una diferencia muy importante respecto a la interpretación de los coeficientes de un modelo de regresión lineal.
La magnitud de cada coeficiente parcial de regresión depende de las unidades en las que se mida la variable predictora a la que corresponde, por lo que su magnitud no está asociada con la importancia de cada predictor. Para poder determinar qué impacto tienen en el modelo cada una de las variables, se emplean los coeficientes parciales estandarizados, que se obtienen al estandarizar (sustraer la media y dividir entre la desviación estándar) las variables predictoras previo ajuste del modelo.
Si bien los coeficientes de regresión suelen ser el primer objetivo de la interpretación de un modelo lineal, existen muchos otros aspectos (significancia del modelo en su conjunto, significancia de los predictores...). Estos últimos suelen ser tratados con poca detalle cuando el único objetivo del modelo es realizar predicciones, sin embargo, son muy relevantes si se quiere realizar inferencia, es decir, explicar las relaciones entre los predictores y la variable respuesta. A lo largo de este documento se irán introduciendo cada uno de ellos.
Significancia del modelo¶
Likelihood Ratio¶
Uno de los primeros resultados que hay que evaluar al ajustar un modelo de regresión logística es el resultado del test de significancia likelihood ratio (LLR). Este contraste responde a la pregunta de si el modelo en su conjunto es capaz de predecir la variable respuesta mejor de lo esperado por azar, o lo que es equivalente, si al menos uno de los predictores que forman el modelo contribuye de forma significativa. Para realizar este análisis se compara la probabilidad de obtener los valores observados (log likelihood) con el modelo de interés ($m_1$) frente a las de un modelo sin predictores (modelo nulo $m_0$).
$$LLR = -2log(\mathcal{L}(m_0)/\mathcal{L}(m_1)) = 2(log \mathcal{L}(m_1) - log \mathcal{L}(m_0))$$El estadístico LLR tiene una distribución chi-cuadrado con grados de libertad equivalentes a la diferencia de grados de libertad de los dos modelos comparados. Si se compara respecto al modelo nulo, los grados de libertad equivalen al número de predictores.
Con frecuencia, la hipótesis nula y alternativa de este test se describen como:
$H_0$: $\beta_1$ = ... = $\beta_{p} = 0$
$H_a$: al menos un $\beta_i \neq 0$
Si el test resulta significativo, implica que el modelo es útil, pero no que sea el mejor. Podría ocurrir que alguno de sus predictores no fuese necesario.
Pseudo R2¶
A diferencia de los modelos de regresión lineal, en los modelos logísticos no existe un equivalente a $R^2$ que determine exactamente la varianza explicada por el modelo. Se han desarrollado diferentes métodos conocidos como $pseudoR^2$ que intentan aproximarse al concepto de $R^2$ pero que, aunque su rango oscila entre 0 y 1, no se pueden considerar equivalentes. El más conocido es el propuesto por McFadden’s:
$$R^2_{McF} = 1-\frac{log \hat{L}(\text{modelo)}}{log \hat{L}(\text{modelo nulo})}$$siendo $\hat{L}$ el valor de likelihood de cada modelo.
$R^2_{McF}$ tiene valor 0 si el modelo candidato no mejora al modelo nulo, y valor de 1 si se ajusta perfectamente a los datos.
Significancia de los predictores¶
En la mayoría de casos, aunque el estudio de regresión logística se aplica a una muestra, el objetivo último es obtener un modelo que explique la relación entre las variables en toda la población. Esto significa que, el modelo generado, es una estimación de la relación poblacional a partir de la relación que se observa en la muestra y, por lo tanto, está sujeta a variaciones. Para cada uno de los coeficientes de la ecuación de regresión logística ($\beta_p$) se puede calcular su significancia (p-value) y su intervalo de confianza. El test estadístico más empleado es el Wald chi-test.
El test de significancia para los coeficientes ($\beta_p$) del modelo logístico considera como hipótesis:
$H_0$: el predictor $x_p$ no contribuye al modelo ($\beta_p = 0$), en presencia del resto de predictores.
$H_a$: el predictor $x_p$ sí contribuye al modelo ($\beta_p \neq 0$), en presencia del resto de predictores.
En los modelos generados con statsmodels se devuelve, junto con el valor de los coeficientes de regresión, el valor del estadístico $z$ obtenido para cada uno, los p-value y los intervalos de confianza correspondientes. Esto permite saber, además de las estimaciones, si son significativamente distintos de 0, es decir, si contribuyen al modelo.
Variables categóricas como predictores¶
Cuando se introduce una variable categórica como predictor, un nivel se considera el de referencia (normalmente codificado como 0) y el resto de niveles se comparan con él. En el caso de que el predictor categórico tenga más de dos niveles, se generan lo que se conoce como variables dummy o one-hot-encodding, que son variables creadas para cada uno de los niveles del predictor categórico y que pueden tomar el valor de 0 o 1. Cada vez que se emplee el modelo para predecir un valor, solo una variable dummy por predictor adquiere el valor 1 (la que coincida con el valor que adquiere el predictor en ese caso) mientras que el resto se consideran 0. El valor del coeficiente parcial de regresión $\beta_p$ de cada variable dummy indica el porcentaje promedio en el que influye dicho nivel sobre el logaritmo de odds de la variable dependiente $y$ en comparación con el nivel de referencia de dicho predictor.
La idea de variables dummy se entiende mejor con un ejemplo. Supóngase un modelo en el que la variable respuesta ojos verdes (si/no) se predice en función de nacionalidad del sujeto. La variable nacionalidad es cualitativa con 3 niveles (americana, europea y asiática). A pesar de que el predictor inicial es nacionalidad, se crea una variable nueva por cada nivel, cuyo valor puede ser 0 o 1. De tal forma que la ecuación del modelo completo es:
$$p(\text{ojos verdes = si}) = \frac{e^{(\beta_0 + \beta_1americana + \beta_2europea + \beta_3asiatica)}}{1+e^{(\beta_0 + \beta_1americana + \beta_2europea + \beta_3asiatica)}}$$Si el sujeto es europeo, las variables dummy americana y asiática se consideran 0, de forma que el modelo para este caso se queda en:
$$p(\text{ojos verdes = si}) = \frac{e^{(\beta_0 + \beta_2europea)}}{1+e^{(\beta_0 + \beta_2europea)}}$$Predicción¶
Una vez generado un modelo válido, es posible predecir la probabilidad de la variable respuesta $y$ para nuevos valores de las variables predictoras $x$. Es importante tener en cuenta que las predicciones deben, a priori, limitarse al rango de valores dentro del que se encuentran las observaciones con las que se ha entrenado el modelo. Esto es importante puesto que, aunque los modelos de regresión logística pueden extrapolarse, solo en esta región se tiene certeza de que se cumplen las condiciones para que el modelo sea válido. Para calcular las predicciones, se emplea la ecuación:
$$\hat{p}(y=1|X) = \frac{e^{(\hat{\beta}_0 + \hat{\beta}_1 x_1 + \cdots + \hat{\beta}_p x_p)}}{1+e^{(\hat{\beta}_0 + \hat{\beta}_1 x_1 + \cdots + \hat{\beta}_p x_p)}}$$Convertir probabilidad en clasificación¶
Una de las principales aplicaciones de un modelo de regresión logística es clasificar la variable cualitativa en función del valor que tomen los predictores. Dado que la salida de un modelo logístico es una probabilidad, para conseguir la clasificación, es necesario establecer un límite (threshold) a partir del cual se considera que la variable pertenece a uno de los niveles. Por ejemplo, se puede asignar una observación al grupo 1 si la probabilidad estimada es mayor de 0.5 y al grupo 0 de lo contrario.
Validación del modelo¶
Una vez seleccionado el mejor modelo que se puede crear con los datos disponibles, se tiene que comprobar su capacidad prediciendo nuevas observaciones que no se hayan empleado para entrenarlo, de este modo se verifica si el modelo se puede generalizar. Una estrategia comúnmente empleada es dividir aleatoriamente los datos en dos grupos, ajustar el modelo con el primer grupo y estimar la precisión de las predicciones con el segundo.
El tamaño adecuado de las particiones depende en gran medida de la cantidad de datos disponibles y la seguridad que se necesite en la estimación del error, 80% 20% suele dar buenos resultados.
Condiciones para la regresión logística¶
Para que un modelo de regresión logística, y las conclusiones derivadas de él, sean completamente válidas, se deben verificar que se cumplen las asunciones sobre las que se basa su desarrollo matemático. En la práctica, rara vez se cumplen, o se puede demostrar que se cumplen todas, sin embargo esto no significa que el modelo no sea útil. Lo importante es ser consciente de ellas y del impacto que esto tiene en las conclusiones que se extraen del modelo.
No colinealidad o multicolinealidad:
En los modelos de regresión logística múltiple, los predictores deben ser independientes, no debe de haber colinealidad entre ellos. La colinealidad ocurre cuando un predictor está linealmente relacionado con uno o varios de los otros predictores del modelo. Como consecuencia de la colinealidad, no se puede identificar de forma precisa el efecto individual que tiene cada predictor sobre la variable respuesta, lo que se traduce en un incremento de la varianza de los coeficientes de regresión estimados hasta el punto de que resulta imposible establecer su significancia estadística. Además, pequeños cambios en los datos, provocan grandes cambios en las estimaciones de los coeficientes. Si bien la colinealidad propiamente dicha existe solo si el coeficiente de correlación simple o múltiple entre predictores es 1, cosa que raramente ocurre en la realidad, es frecuente encontrar la llamada casi-colinealidad o multicolinealidad no perfecta.
En caso de encontrar colinealidad entre predictores, hay dos posibles soluciones. La primera es excluir uno de los predictores problemáticos intentando conservar el que, a juicio del investigador, está influyendo realmente en la variable respuesta. Esta medida no suele tener mucho impacto en el modelo en cuanto a su capacidad predictiva ya que, al existir colinealidad, la información que aporta uno de los predictores es redundante en presencia del otro. La segunda opción consiste en combinar las variables colineales en un único predictor, aunque con el riesgo de perder su interpretación.
Cuando se intenta establecer relaciones causa-efecto, la colinealidad puede llevar a conclusiones muy erróneas, haciendo creer que una variable es la causa cuando, en realidad, es otra la que está influenciando sobre ese predictor.
Relación lineal entre los predictores numéricos y el logaritmo de odds de la variable respuesta
Cada predictor numérico tiene que estar linealmente relacionado con el logaritmo de odds de la variable respuesta $y$ mientras los demás predictores se mantienen constantes, de lo contrario no se deben introducir en el modelo.
No autocorrelación (Independencia)
Los valores de cada observación son independientes de los otros. Esto es especialmente importante de comprobar cuando se trabaja con mediciones temporales. Se recomienda representar los residuos ordenados acorde al tiempo de registro de las observaciones, si existe un cierto patrón hay indicios de autocorrelación. También se puede emplear el test de hipótesis de Durbin-Watson.
Valores atípicos, con alto leverage o influyentes
Es importante identificar observaciones que sean atípicas o que puedan estar influenciando al modelo. La forma más fácil de detectarlas es a través de los residuos.
Tamaño de la muestra
No se trata de una condición de por sí pero, si no se dispone de suficientes observaciones, predictores que no son realmente influyentes podrían parecerlo. Un recomendación frecuente es que el número de observaciones sea como mínimo entre 10 y 20 veces el número de predictores del modelo.
Parsimonia
Este término hace referencia a que, el mejor modelo, es aquel capaz de explicar con mayor precisión la variabilidad observada en la variable respuesta empleando el menor número de predictores, por lo tanto, con menos asunciones.
La gran mayoría de condiciones se verifican utilizando los residuos, por lo tanto, se suele generar primero el modelo y posteriormente validar las condiciones. De hecho, el ajuste de un modelo debe verse como un proceso iterativo en el que se ajusta el modelo, se evalúan sus residuos y se mejora. Así hasta llegar a un modelo óptimo.
Ejemplo regresión logística simple¶
Un estudio quiere establecer un modelo que permita calcular la probabilidad de obtener matrícula de honor al final del bachillerato en función de la nota que se ha obtenido en matemáticas. La variable matrícula está codificada como 0 si no se tiene matrícula y 1 si se tiene.
# Tratamiento de datos
# ==============================================================================
import pandas as pd
import numpy as np
# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from matplotlib import style
import seaborn as sns
# Preprocesado y modelado
# ==============================================================================
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.weightstats import ttest_ind
# Configuración matplotlib
# ==============================================================================
plt.rcParams['image.cmap'] = "bwr"
#plt.rcParams['figure.dpi'] = "100"
plt.rcParams['savefig.bbox'] = "tight"
style.use('ggplot') or plt.style.use('ggplot')
# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')
Datos¶
# Datos
# ==============================================================================
matricula = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1,
0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0,
1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1,
1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0,
0, 0, 0, 0, 1, 0, 0, 0, 1, 1])
matematicas = np.array([
41, 53, 54, 47, 57, 51, 42, 45, 54, 52, 51, 51, 71, 57, 50, 43,
51, 60, 62, 57, 35, 75, 45, 57, 45, 46, 66, 57, 49, 49, 57, 64,
63, 57, 50, 58, 75, 68, 44, 40, 41, 62, 57, 43, 48, 63, 39, 70,
63, 59, 61, 38, 61, 49, 73, 44, 42, 39, 55, 52, 45, 61, 39, 41,
50, 40, 60, 47, 59, 49, 46, 58, 71, 58, 46, 43, 54, 56, 46, 54,
57, 54, 71, 48, 40, 64, 51, 39, 40, 61, 66, 49, 65, 52, 46, 61,
72, 71, 40, 69, 64, 56, 49, 54, 53, 66, 67, 40, 46, 69, 40, 41,
57, 58, 57, 37, 55, 62, 64, 40, 50, 46, 53, 52, 45, 56, 45, 54,
56, 41, 54, 72, 56, 47, 49, 60, 54, 55, 33, 49, 43, 50, 52, 48,
58, 43, 41, 43, 46, 44, 43, 61, 40, 49, 56, 61, 50, 51, 42, 67,
53, 50, 51, 72, 48, 40, 53, 39, 63, 51, 45, 39, 42, 62, 44, 65,
63, 54, 45, 60, 49, 48, 57, 55, 66, 64, 55, 42, 56, 53, 41, 42,
53, 42, 60, 52, 38, 57, 58, 65])
datos = pd.DataFrame({'matricula': matricula, 'matematicas': matematicas})
datos.head(3)
Exploración gráfica¶
El primer paso antes de generar un modelo de regresión logística simple es representar los datos para poder intuir si existe una relación entre la variable independiente y la variable respuesta.
# Número de obsercaciones por clase
# ==============================================================================
datos.matricula.value_counts().sort_index()
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 3.84))
sns.violinplot(
x = 'matricula',
y = 'matematicas',
data = datos,
#color = "white",
ax = ax
)
ax.set_title('Distribución notas de matemáticas por clase');
# T-test entre clases
# ==============================================================================
res_ttest = ttest_ind(
x1 = matematicas[matricula == 0],
x2 = matematicas[matricula == 1],
alternative='two-sided'
)
print(f"t={res_ttest[0]}, p-value={res_ttest[1]}")
Tanto el gráfico como el t-test muestran evidencias de que existe una diferencia entre la nota de las personas con matrícula y sin matrícula. Esta información es útil para considerar la nota de matemáticas como un buen predictor para el modelo.
Ajuste del modelo¶
Se ajusta un modelo empleando como variable respuesta matricula y como predictor matematicas. Como en todo estudio predictivo, no solo es importante ajustar el modelo, sino también cuantificar su capacidad para predecir nuevas observaciones. Para poder hacer esta evaluación, se dividen los datos en dos grupos, uno de entrenamiento y otro de test.
Scikit-learn¶
# División de los datos en train y test
# ==============================================================================
X = datos[['matematicas']]
y = datos['matricula']
X_train, X_test, y_train, y_test = train_test_split(
X.values.reshape(-1,1),
y.values.reshape(-1,1),
train_size = 0.8,
random_state = 1234,
shuffle = True
)
# Creación del modelo
# ==============================================================================
# Para no incluir ningún tipo de regularización en el modelo se indica
# penalty='none'
modelo = LogisticRegression(penalty='none')
modelo.fit(X = X_train.reshape(-1, 1), y = y_train)
# Información del modelo
# ==============================================================================
print("Intercept:", modelo.intercept_)
print("Coeficiente:", list(zip(X.columns, modelo.coef_.flatten(), )))
print("Accuracy de entrenamiento:", modelo.score(X, y))
Una vez entrenado el modelo, se pueden predecir nuevas observaciones.
# Predicciones probabilísticas
# ==============================================================================
# Con .predict_proba() se obtiene, para cada observación, la probabilidad predicha
# de pertenecer a cada una de las dos clases.
predicciones = modelo.predict_proba(X = X_test)
predicciones = pd.DataFrame(predicciones, columns = modelo.classes_)
predicciones.head(3)
# Predicciones con clasificación final
# ==============================================================================
# Con .predict() se obtiene, para cada observación, la clasificación predicha por
# el modelo. Esta clasificación se corresponde con la clase con mayor probabilidad.
predicciones = modelo.predict(X = X_test)
predicciones
Statsmodels¶
La implementación de regresión logística de Statsmodels, es más completa que la de Scikitlearn ya que, además de ajustar el modelo, permite calcular los test estadísticos y análisis necesarios para verificar que se cumplen las condiciones sobre las que se basa este tipo de modelos. Statsmodels tiene dos formas de entrenar el modelo:
Indicando la fórmula del modelo y pasando los datos de entrenamiento como un dataframe que incluye la variable respuesta y los predictores. Esta forma es similar a la utilizada en R.
Pasar dos matrices, una con los predictores y otra con la variable respuesta. Esta es igual a la empleada por Scikitlearn con la diferencia de que a la matriz de predictores hay que añadirle una primera columna de 1s.
# División de los datos en train y test
# ==============================================================================
X = datos[['matematicas']]
y = datos['matricula']
X_train, X_test, y_train, y_test = train_test_split(
X.values.reshape(-1,1),
y.values.reshape(-1,1),
train_size = 0.8,
random_state = 1234,
shuffle = True
)
# Creación del modelo utilizando el modo fórmula (similar a R)
# ==============================================================================
# datos_train = pd.DataFrame(np.hstack((X_train, y_train)),
# columns=['matematicas', 'matricula'])
# modelo = smf.logit(formula = 'matricula ~matematicas', data = datos_train)
# modelo = modelo.fit()
# print(modelo.summary())
# Creación del modelo utilizando matrices como en scikitlearn
# ==============================================================================
# A la matriz de predictores se le tiene que añadir una columna de 1s para el intercept del modelo
X_train = sm.add_constant(X_train, prepend=True)
modelo = sm.Logit(endog=y_train, exog=X_train,)
modelo = modelo.fit()
print(modelo.summary())
El coeficiente estimado para la intersección (Intercept o const) es el valor esperado del logaritmo de odds de que un estudiante obtenga matrícula teniendo un 0 en matemáticas. Como cabría imaginar, los odds son muy bajos $e^-8.9848 = 0.0001254$, lo que se corresponde con una probabilidad de obtener matrícula de $p=\frac{e^{0.0001254}}{1+e^{0.0001254}}=0.000125$.
Acorde al modelo, el logaritmo de odds de que un estudiante tenga matrícula está positivamente relacionado con la puntuación obtenida en matemáticas (coeficiente de regresión = 0.1439). Esto significa que, por cada unidad que se incrementa la variable matemáticas, se espera que el logaritmo de odds de la variable matrícula se incremente en promedio 0.1439 unidades. Aplicando la inversa del logaritmo natural ($e^{0.1439}=1.155$) se obtiene que, por cada unidad que se incrementa la variable matemáticas, los odds de obtener matrícula se incremente en promedio 1.169 unidades. No hay que confundir esto último con que la probabilidad de matrícula se incremente un 1.169 %.
A diferencia de la regresión lineal en la que $\beta_1$ se corresponde con el cambio promedio en la variable dependiente $y$ debido al incremento en una unidad del predictor $x_1$, en regresión logística, $\beta_1$ indica el cambio en el logaritmo de odds debido al incremento en una unidad de $x_1$, o lo que es lo mismo, multiplica los odds por $e^{\beta_1}$. Dado que la relación entre $p(y=1)$ y $x$ no es lineal, los coeficientes de regresión $\beta_{p}$ no se corresponden con el cambio en la probabilidad de $y$ asociada con el incremento en una unidad de $x$. Cuánto se incremente la probabilidad de $y$ por unidad de $x$ depende del valor de $x$, es decir, de la posición en la curva logística en la que se encuentre.
Además del valor de las estimaciones de los coeficientes parciales de correlación del modelo, es conveniente calcular sus correspondientes intervalos de confianza.
Intervalos de confianza de los coeficientes¶
# Intervalos de confianza para los coeficientes del modelo
# ==============================================================================
intervalos_ci = modelo.conf_int(alpha=0.05)
intervalos_ci = pd.DataFrame(intervalos_ci)
intervalos_ci.columns = ['2.5%', '97.5%']
intervalos_ci
Predicciones¶
Una vez entrenado el modelo, se pueden obtener predicciones para nuevos datos. Los modelos de regresión logística de statsmodels devuelven la probabilidad de pertenecer a la clase de referencia.
# Predicción de probabilidades
# ==============================================================================
predicciones = modelo.predict(exog = X_train)
predicciones[:4]
Para obtener la clasificación final, se convierten los valores de probabilidad mayores de 0.5 a 1 y los mejores a 0.
# Clasificación predicha
# ==============================================================================
clasificacion = np.where(predicciones<0.5, 0, 1)
clasificacion
Representación gráfica del modelo¶
Además de la línea de mínimos cuadrados, es recomendable incluir los límites superior e inferior del intervalo de confianza. Esto permite identificar la región en la que, según el modelo generado y para un determinado nivel de confianza, se encuentra el valor promedio de la variable respuesta.
# Predicciones en todo el rango de X
# ==============================================================================
# Se crea un vector con nuevos valores interpolados en el rango de observaciones.
grid_X = np.linspace(
start = min(datos.matematicas),
stop = max(datos.matematicas),
num = 200
).reshape(-1,1)
grid_X = sm.add_constant(grid_X, prepend=True)
predicciones = modelo.predict(exog = grid_X)
# Gráfico del modelo
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 3.84))
ax.scatter(
X_train[(y_train == 1).flatten(), 1],
y_train[(y_train == 1).flatten()].flatten()
)
ax.scatter(
X_train[(y_train == 0).flatten(), 1],
y_train[(y_train == 0).flatten()].flatten()
)
ax.plot(grid_X[:, 1], predicciones, color = "gray")
ax.set_title("Modelo regresión logística")
ax.set_ylabel("P(matrícula = 1 | matemáticas)")
ax.set_xlabel("Nota matemáticas");
Accuracy de test¶
Se calcula el porcentaje de aciertos que tiene el modelo al predecir las observaciones de test (accuracy).
# Accuracy de test del modelo
# ==============================================================================
X_test = sm.add_constant(X_test, prepend=True)
predicciones = modelo.predict(exog = X_test)
clasificacion = np.where(predicciones<0.5, 0, 1)
accuracy = accuracy_score(
y_true = y_test,
y_pred = clasificacion,
normalize = True
)
print("")
print(f"El accuracy de test es: {100*accuracy}%")
Conclusión¶
El modelo logístico creado para predecir la probabilidad de que un alumno obtenga matrícula de honor a partir de la nota de matemáticas es en conjunto significativo (Likelihood ratio p-value = 9.831e-11). El p-value del predictor matemáticas es significativo (p-value = 7.17e-08).
$$\text{P(matricula)} = \frac{e^{-8.9848 + 0.1439 * \text{nota matematicas}}}{1+e^{-8.9848 + 0.1439 * \text{nota matematicas}}}$$Los resultados obtenidos con el conjunto de test indican que el modelo es capaz de clasificar correctamente el 87.5% de las observaciones.
Ejemplo regresión logística múltiple¶
El set de datos spam, obtenido de UCI Repository Of Machine Learning Databases contiene información sobre 4601 correos electrónicos clasificados como spam y no spam.
Para cada correo electrónico se dispone de 58 variables: las 48 primeras contienen la frecuencia con la que aparecen en el texto del email determinadas palabras. Las variables 49-54 indican la frecuencia con la que aparecen los caracteres ;’, ‘(’, ‘[’, ‘!’, ‘\$’, ‘#’. Las variables 55-57 contienen la media, a longitud máxima y el número total de letras mayúsculas.
Se crear un modelo de regresión logística con el objetivo de clasificar si un correo es spam o no.
# Tratamiento de datos
# ==============================================================================
import pandas as pd
import numpy as np
# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from matplotlib import style
import seaborn as sns
# Preprocesado y modelado
# ==============================================================================
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.metrics import plot_confusion_matrix
import statsmodels.api as sm
import statsmodels.formula.api as smf
# Configuración matplotlib
# ==============================================================================
plt.rcParams['image.cmap'] = "bwr"
#plt.rcParams['figure.dpi'] = "100"
plt.rcParams['savefig.bbox'] = "tight"
style.use('ggplot') or plt.style.use('ggplot')
# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')
Datos¶
# Datos
# ==============================================================================
url = 'https://raw.githubusercontent.com/JoaquinAmatRodrigo/' \
+ 'Estadistica-machine-learning-python/master/data/spam.csv'
datos = pd.read_csv(url)
datos.head(3)
Se codifica la variable respuesta como 1 si es spam y 0 si no lo es, y se identifica cuantas observaciones hay de cada clase.
datos['type'] = np.where(datos['type'] == 'spam', 1, 0)
print("Número de observaciones por clase")
print(datos['type'].value_counts())
print("")
print("Porcentaje de observaciones por clase")
print(100 * datos['type'].value_counts(normalize=True))
El 66.6% de los correos no son spam y el 39.4% sí lo son. Un modelo de clasificación que sea útil debe de ser capaz de predecir correctamente un porcentaje de observaciones por encima del porcentaje de la clase mayoritaria. En este caso, el umbral de referencia que se tiene que superar es del 66.6%.
Ajuste del modelo¶
Se ajusta un modelo de regresión logística múltiple con el objetivo de predecir si un correo es spam en función de todas las variables disponibles.
# División de los datos en train y test
# ==============================================================================
X = datos.drop(columns = 'type')
y = datos['type']
X_train, X_test, y_train, y_test = train_test_split(
X,
y.values.reshape(-1,1),
train_size = 0.8,
random_state = 1234,
shuffle = True
)
# Creación del modelo utilizando matrices como en scikitlearn
# ==============================================================================
# A la matriz de predictores se le tiene que añadir una columna de 1s para el intercept del modelo
X_train = sm.add_constant(X_train, prepend=True)
modelo = sm.Logit(endog=y_train, exog=X_train,)
modelo = modelo.fit()
print(modelo.summary())
Predicciones¶
Una vez entrenado el modelo, se pueden obtener predicciones para nuevos datos. Los modelos de statsmodels permiten calcular los intervalos de confianza asociados a cada predicción.
# Predicciones con intervalo de confianza
# ==============================================================================
predicciones = modelo.predict(exog = X_train)
# Clasificación predicha
# ==============================================================================
clasificacion = np.where(predicciones<0.5, 0, 1)
clasificacion
Accuracy de test¶
Se calcula el porcentaje de aciertos que tiene el modelo al predecir las observaciones de test (accuracy).
# Accuracy de test del modelo
# ==============================================================================
X_test = sm.add_constant(X_test, prepend=True)
predicciones = modelo.predict(exog = X_test)
clasificacion = np.where(predicciones<0.5, 0, 1)
accuracy = accuracy_score(
y_true = y_test,
y_pred = clasificacion,
normalize = True
)
print("")
print(f"El accuracy de test es: {100*accuracy}%")
# Matriz de confusión de las predicciones de test
# ==============================================================================
confusion_matrix = pd.crosstab(
y_test.ravel(),
clasificacion,
rownames=['Real'],
colnames=['Predicción']
)
confusion_matrix
Conclusión¶
El modelo logístico creado para predecir la probabilidad de que un correo sea spam es en conjunto significativo (Likelihood ratio p-value = 0). El porcentaje de clasificación correcta en el conjunto del test es del 92.8%, un valor muy por encima del umbral de 66.6% esperado por azar.
Acorde al p-value individual de cada predictor, solo algunos de ellos aportan información al modelo. Es conveniente identificar cuáles son y excluirlos para simplificar el modelo y evitar la introducción de ruido. Una forma de conseguir reducir la influencia de predictores que no aportan a un modelo de regresión logística es incorporar regularización en su ajuste.
Información de sesión¶
from sinfo import sinfo
sinfo()
Bibliografía¶
Linear Models with R by Julian J.Faraway
An Introduction to Statistical Learning: with Applications in R (Springer Texts in Statistics)
OpenIntro Statistics: Fourth Edition by David Diez, Mine Çetinkaya-Rundel, Christopher Barr
Extending the Linear Model with R: Generalized Linear, Mixed Effects and Nonparametric Regression Models by Julian J.Faraway
Introduction to Machine Learning with Python: A Guide for Data Scientists
Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011.
¿Cómo citar este documento?
Regresión logística con Python por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/py17-regresion-logistica-python.html 
¿Te ha gustado el artículo? Tu ayuda es importante
Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
Este contenido, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.