Regresión logística para predicciones bancarias¶
Este estudio trata de construir un modelo predictivo que responda a siguiente pregunta:
Para ello hemos utilizado un modelo de regresión logística. Tras el análisis exploratorio de las variables, hemos seleccionado las 12 siguientes, las más significativas, como variables predictoras para la construcción del modelo:
Variables predictoras:
- previous: número de veces que hemos contactado con el cliente antes de esta campaña
- euribor3m: tipos de interés Euribor con un plazo de vigencia de 3 meses
Las otras 10 son variables dummies categóricas creadas para este estudio. El valor 1 indica "sí", el valor 0 indica "no":
- job_entrepreneur: si el cliente es emprendedor
- job_management: si el cliente es administrador
- job_self-employed: si el cliente es autónomo
- month_apr: si el último contacto con el cliente fue en el mes de abril
- month_dec: si el último contacto con el cliente fue en el mes de diciembre
- month_mar: si el último contacto con el cliente fue en el mes de marzo
- month_may: si el último contacto con el cliente fue en el mes de mayo
- month_nov: si el último contacto con el cliente fue en el mes de noviembre
- poutcome_failure: si no compró en la campaña de marketing previa (failure in previous outcome)
- poutcome_success: si sí compró en la campaña de marketing previa (success in previous outcome)
Variable de salida:
- y: si el cliente ha comprado o no (1: sí, 0: no)
1. Análisis exploratorio de los datos¶
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
data = pd.read_csv("datasets/bank/bank.csv", sep=";")
data.head()
data.shape
data.columns.values
Limpieza de datos:
Vamos a poner los datos en un formato uniforme y a darle un valor numérico a la variable de salida
data["y"] = (data["y"]=="yes").astype(int)
data.tail()
data["education"].unique()
#where(haya_esto, "me_pones_esto", y en caso contrario este 3er argumento)
data["education"] = np.where(data["education"]=="basic.4y", "Basic", data["education"])
data["education"] = np.where(data["education"]=="basic.6y", "Basic", data["education"])
data["education"] = np.where(data["education"]=="basic.9y", "Basic", data["education"])
data["education"] = np.where(data["education"]=="high.school", "High School", data["education"])
data["education"] = np.where(data["education"]=="professional.course", "Professional Course", data["education"])
data["education"] = np.where(data["education"]=="university.degree", "University Degree", data["education"])
data["education"] = np.where(data["education"]=="illiterate", "Illiterate", data["education"])
data["education"] = np.where(data["education"]=="unknown", "Unknown", data["education"])
data["education"].unique()
1.1 Información estadística por categorías:¶
Comprobamos los valores medios de las medidas entre los clientes que compran (1) y los que no (0) y el número total de compras, así como la frecuencia de compras entre distintos grupos según su nivel educativo, estado civil, día de la semana del último contacto con el cliente y mes del último contacto con el cliente.
#Medias
data.groupby("y").mean()
#Total de compras (1) y no compras (0)
data["y"].value_counts()
Por nivel educativo:
#Medias por niveles educativos
data.groupby("education").mean()
#Total de compras (1) y no compras (0) por grupos de nivel educativo
pd.crosstab(data["education"], data["y"])
%matplotlib inline
pd.crosstab(data.education, data.y).plot(kind="bar")
plt.title("Frecuencia de compra en función del nivel de educación")
plt.xlabel("Nivel de educación")
plt.ylabel("Frecuencia de compra del producto")
Por estado civil:
table=pd.crosstab(data.marital, data.y)
table.div(table.sum(1).astype(float), axis=0).plot(kind="bar", stacked=True) #stacked = True es para que salga apilado
#Con table.div(table.sum(1).astype(float), axis=0) es como arrastrar en excel la columna (axis = 0)
#para que cada celda quede dividida por la suma de su fila
plt.title("Diagrama apilado de estado civil contra el nivel de compras")
plt.xlabel("Estado civil")
plt.ylabel("Proporción de clientes")
Por día de la semana del último contacto con el cliente:
%matplotlib inline
table= pd.crosstab(data.day_of_week, data.y)
table.div(table.sum(1).astype(float), axis=0).plot(kind="bar", stacked=True)
plt.title("Frecuencia de compra en función del día de la semana")
plt.xlabel("Día de la semana")
plt.ylabel("Frecuencia de compra del producto")
Por mes del último contacto con el cliente:
%matplotlib inline
table= pd.crosstab(data.month, data.y)
table.div(table.sum(1).astype(float), axis=0).plot(kind="bar", stacked=True)
plt.title("Frecuencia de compra en función del mes")
plt.xlabel("Mes del año")
plt.ylabel("Frecuencia de compra del producto")
Tal vez es que en marzo y diciembre haya menos datos y por eso la proporción de compra / no compra es mayor, sin que eso signifique realmente un volumen de compras mayor
%matplotlib inline
table.plot(kind="bar", stacked=False)
plt.title("Frecuencia de compra en función del mes")
plt.xlabel("Mes del año")
plt.ylabel("Frecuencia de compra del producto")
Efectivamente vemos que en diciembre y marzo apenas hay datos, y las mayor proporción de compras frente a rechazos no significaba un mayor volumen de compras en términos absolutos. Vemos que los meses en los que se hizo campaña fueron en primavera y verano (mayo, junio, julio, agosto) y hubo mucho rechazos
Por edad del cliente:
%matplotlib inline
data.age.hist()
plt.title("Histograma de la Edad")
plt.xlabel("Edad")
plt.ylabel("Número de Clientes")
Los clientes tienen edades entre 25 y 60 años mayoritariamente.
pd.crosstab(data.age, data.y).plot(kind="bar", figsize=(16,6), stacked=True)
compras_abs = pd.crosstab(data.age, data.y)
frec_abs_edades = data.age.value_counts().sort_index()
compras_rel = compras_abs.div(frec_abs_edades, axis=0)
compras_rel.plot(kind="bar", stacked=True, figsize=(16,6))
Los mayores de 60, aunque son pocos, tienen en general mayor probabilidad de compra
Por frecuencia de compras anteriores:
#poutcome es la frecuencia de compras anteriores
pd.crosstab(data.poutcome, data.y).plot(kind="bar")
poutcome es buena variable para predecir si alguien va a llevar a cabo una inversión. Vemos que la gente que ha tenido malas experiencias (failure) pasadas, no repite, pero los que tuvieron éxito sí que repiten. Los que no han invertido, por defecto no van a invertir, con algunas excepciones
2. Preprocesado de datos
2.1. Conversión de variables categóricas a variables dummies
categories = ["job", "marital", "education", "default", "housing", "loan", "contact",
"month", "day_of_week", "poutcome"]
for category in categories:
cat_dummies = pd.get_dummies(data[category], prefix=category)
data = data.join(cat_dummies)
data.columns.values
Eliminación de las categorías originales:
data_vars = data.columns.values.tolist()
to_keep = [v for v in data_vars if v not in categories] #Me quiero quedar con estos
bank_data = data[to_keep]
2.2. Selección de variables para la creación del modelo
¿Qué variables son las más significativas para la predicción del modelo?
n = 12
from sklearn import datasets
from sklearn.feature_selection import RFE #Recursive Feature Elimination
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(solver="lbfgs", max_iter = 10000) #Me daba warnings si no metía esos argumentos
bank_data_list = bank_data.columns.values.tolist()
Y_list = ['y']
X_list = [v for v in bank_data_list if v not in Y_list]
rfe = RFE(lr, n) #Con el método RFE elegimos, con el modelo logístico, el número de variables con el que queremos quedarnos
rfe = rfe.fit(bank_data[X_list], bank_data[Y_list].values.ravel())
rfe.support_ ¿Cuales son las 12 variables que se van a quedar en el modelo?
rfe.support_
rfe.support_.sum() #Hay 12 Trues, ya que especifiqué que quería quedarme con 12 variables
rfe.ranking_ ¿Cuales tengo que ir añadiendo si quiero sumar una 13ª, 14ª, etc.? Con rfe.ranking_ obtenemos las 12 que me entran con un 1, y las siguientes que me entrarían con un 2,3,... por orden de prioridad
print(rfe.ranking_)
zip Para ver juntos los True o False junto con los nombres de las variables correspondientes que entran en el modelo o se quedan fuera (respectivamente):
z=zip(X_list,rfe.support_, rfe.ranking_)
z_list = list(z)
z_list
Creo una lista con las 12 variables que han sido seleccionadas:
support = rfe.support_ #Array de True y False
#X_list era una lista con los nombres de las columnas de las variables predictoras (excluida la "y")
columnas = [nombre for nombre in X_list if support[X_list.index(nombre)]==np.bool_("True")]
columnas
Nota: Habíamos hecho una regresión logística con todas las variables, para después ver qué 12 son las mejores predictoras. Ahora finalmente construimos el modelo final sólo con estas 12
3. Implementación del modelo de regresión logística con statsmodel.api
import statsmodels.api as sm
X = bank_data[columnas]
Y = bank_data["y"]
logit_model = sm.Logit(Y, X)
result = logit_model.fit() #Método de Newton-Raphson
result.summary2()
En regresiones logísticas, el estadístico z nos indica si los coeficientes para cada variable predictora son significativamente distintos de 0 y por tanto están haciendo una contribución significativa al modelo para predecir y. "Z" es el número de veces la desviación estándar que tenemos que movernos desde la media (el valor dado para "coef") hasta el 0. Esto coincide con el valor de z en la transformación de una distribución normal con media != 0 y std != 0 a una distribución normal, lo cual se hace, por ejemplo, en el contraste de hipótesis.
$$z = \frac{x - \mu}{\sigma} = \frac{0 - coef}{str.err} \leftarrow \rightarrow \sigma · z = -coef $$El z en este caso, sin embargo, es simplemente esta transformación, pero no tiene nada que ver con un contraste de hipótesis que esté haciendo la librería ni nada por el estilo.
P>|z| es el área bajo la curva desde el 0 en adelante (si coef es negativo) o desde el 0 hacia atrás (si coef es positivo).
Cuanto mayor sea z (en valor absoluto), menos probable es que el coeficiente sea 0. A continuación comprobamos que dividiendo los coef entre las std.err obtenemos los valores de z:
std_err = np.array([0.1410, 0.0204, 0.3812, 0.2249, 0.3258, 0.1932, 0.4654, 0.3154, 0.1241, 0.1935, 0.2501, 0.2927])
result.params.div(std_err) #result.params da los coef.
Las pequeñas discrepancias se deben a los valores redondeados que hemos tomado de std_err, tal y como aparecen en la tabla. Por ejemplo, en el caso del euribor3m, ese str.err = 0.2024 es un redondeo del valor exacto 0.02035.... Vemos que dividiendo entre un valor así sí se obtiene el z de la tabla:
result.params['euribor3m']/0.02035045
4. Implementación del modelo de regresión logística con scikit-learn
from sklearn import linear_model
logit_model = linear_model.LogisticRegression(solver="lbfgs")
logit_model.fit(X,Y)
score: Factor $R^{2}$. Nos dice la calidad del ajuste:
logit_model.score(X,Y)
1-Y.mean()
$1-Y.mean()$ me da el tanto por 1 de gente que no compra (y=0). Si hicieramos las predicciones diciendo "no va a comprar", acierto en el 89,05% de los casos. Este sería el máximo % de valor que puedo acertar sin un modelo y nuestro valor de partida a mejorar.
Con el modelo que hemos creado acertamos el 90,22%
pd.DataFrame(list(zip(X.columns, np.transpose(logit_model.coef_))))
Vemos que los coeficientes que obtenemos no son exactamente los mismos con los modulos de scikit learn que con los de statsmodel.api, aunque los resultados de la "y" que queremos inferir no diferirán demasiado. Pero tenemos que hacer un análisis exhaustivo de las variables, no podemos dejar que el modelo decida por nosotros mismos
5. Validación del modelo logístico
Vamos a hacer una division entre conjunto de entrenamiento y de testing y a validar el modelo con el conjunto de testing
from sklearn.model_selection import train_test_split
X_train, X_test, Y_train, Y_test = train_test_split(X,Y, test_size = 0.3, random_state=0)
#random_state=0 es para fijar la semilla de la división aleatoria de los conjuntos
lm = linear_model.LogisticRegression(solver = "lbfgs")
lm.fit(X_train, Y_train)
El modelo predice probabilidades. Luego nosotros establecemos un umbral (por ejemplo, 0,5), a partir del cual se considerará como 1, y si la probabilidad es más pequeña que el umbral se considerará como 0
from IPython.display import display, Math, Latex
display(Math(r'Y_p=\begin{cases}0& si\ p\leq0.5\\1&si\ p >0.5\end{cases}'))
probs = lm.predict_proba(X_test)
El módulo lm.predictproba devuelve las probabilidades de cada clase (en este caso hay 2 clases, y = 0, y = 1 (no compra, compra)), según el orden dado en lm.classes
probs #1ª col: probabilidades de que no compra. 2ª col: probabilidades de que sí compre
prediction = lm.predict(X_test)
Porcentaje de datos capaz de predecir correctamente
from sklearn import metrics
metrics.accuracy_score(Y_test, prediction, normalize=False)
$Número\ de\ veces\ que\ y_{pred} = y_{i}$
metrics.accuracy_score(Y_test, prediction) #Por defecto normalize=True
$\large \frac{Número\ de\ veces\ que\ y_{pred} = y_{i}}{Total\ de\ puntos\ (len(Y_{test}))}$
Ha mejorado un poco respecto al anterior modelo en el que no dividimos los datos en los conjuntos de entrenamiento y testing. Con este modelo predecimos un 90,37% del valor
Compras con probabilidad p > 50%
El umbral que hemos establecido es el que python toma por defecto con el módulo predict: p = 0,5
prediction
np.sum(prediction)/prediction.shape[0]
Tenemos un 2,7% de compras considerando como "compra" una probabilidad superior al 50% de que el cliente compre
Compras con probabilidad p > 10%
Podemos elegir nosotros otro umbral. Como tenemos sólo un 10% de clientes que compran, puede que un umbral de 0,10 sea un buen umbral de decisión. Así, si las probabilidades de que compre son superiores al 10%, lo catalogaremos como venta:
display(Math(r'\varepsilon\in (0,1), Y_p=\begin{cases}0& si\ p\leq \varepsilon\\1&si\ p >\varepsilon\end{cases}'))
prob_yes = probs[:,1] #Todas las filas, segunda columna (probabilidades de que sí compre, y = 1)
prob_df = pd.DataFrame(prob_yes)
threshold = 0.1
prob_df["prediction"] = np.where(prob_df[0]>threshold, 1, 0)
prob_df.tail()
Si la probabilidad de que compre es superior al 10%, lo catalogaremos como "compra"
pd.crosstab(prob_df.prediction, columns="contador")
336/len(prob_df)*100
27% de compras considerando "compra" una probabilidad de más del 10%
Compras con probabilidad p > 15%
threshold = 0.15
prob_df["prediction"] = np.where(prob_df[0]>threshold, 1, 0)
pd.crosstab(prob_df.prediction, columns="count")
199/len(prob_df)*100
16% de compras considerando "compra" una probabilidad de más del 15%
Compras con probabilidad p > 5%
Nos puede interesar establecer un umbral más bajo y captar al mayor número de clientes posible
threshold = 0.05
prob_df["prediction"] = np.where(prob_df[0]>threshold, 1, 0)
pd.crosstab(prob_df.prediction, columns="count")
920/len(prob_df)*100
74% de compras considerando "compra" una probabilidad de más del 5%
5.1. Validación cruzada
from sklearn.model_selection import cross_val_score
scores = cross_val_score(linear_model.LogisticRegression(solver = "lbfgs"), X, Y, scoring="accuracy", cv=10)
- scoring = "accuracy" --> el objetivo es calcular la mejor eficacia
- cv: número de partes (k) en el que se divide el dataset para hacer la cross validation
scores
scores.mean()
No es muy diferente al valor que obtuvimos, luego podemos decir que el modelo es bueno, no había problemas de overfitting ni nada
5.2. Matrices de confusión y curvas ROC
# probs = lm.predict_proba(X_test)
probs
prob_yes=probs[:,1]
prob_df = pd.DataFrame(probs)
threshold = 0.1
prob_df["prediction"] = np.where(prob_df[1]>=threshold, 1, 0)
prob_df["actual"] = list(Y_test)
prob_df.head()
confusion_matrix = pd.crosstab(prob_df.prediction, prob_df.actual)
confusion_matrix
#OJO: matrix[columna][fila]
TN=confusion_matrix[0][0]
TP=confusion_matrix[1][1]
FP=confusion_matrix[0][1]
FN=confusion_matrix[1][0]
sens = TP/(TP+FN) #Sensibilidad (% de valores actuales positivos predichos como tal por el modelo)
sens
espc_1 = 1-TN/(TN+FP) #El % de todos los valores actuales negativos que son Falsos Positivos --> Also: FP/(TN+FP)
espc_1
thresholds = [0.04, 0.05, 0.07, 0.10, 0.12, 0.15, 0.18, 0.20, 0.25, 0.3, 0.4, 0.5]
sensitivities = [1]
especifities_1 = [1]
for t in thresholds:
prob_df["prediction"] = np.where(prob_df[1]>=t, 1, 0)
prob_df["actual"] = list(Y_test)
prob_df.head()
confusion_matrix = pd.crosstab(prob_df.prediction, prob_df.actual)
TN=confusion_matrix[0][0]
TP=confusion_matrix[1][1]
FP=confusion_matrix[0][1]
FN=confusion_matrix[1][0]
sens = TP/(TP+FN)
sensitivities.append(sens)
espc_1 = 1-TN/(TN+FP)
especifities_1.append(espc_1)
sensitivities.append(0)
especifities_1.append(0)
#Añadimos el 1 y el 0 para luego la representación.
sensitivities #(True Positives)/(total valores positivos actuales)
especifities_1 #(False Positives)/(total valores negativos actuales)
import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(especifities_1, sensitivities, marker="o", linestyle="--", color="r")
x=[i*0.01 for i in range(100)]
y=[i*0.01 for i in range(100)]
plt.plot(x,y)
plt.xlabel("1-Especifidad")
plt.ylabel("Sensibilidad")
plt.title("Curva ROC")
- Thresholds altos corresponden a puntos abajo a la izquierda, y threshold bajos corresponden a puntos arriba a la derecha. Puntos muy aglutinados indican probabilidades muy extremas (cercanas a 0/1), y puntos agrupados en 2 grupos muy separados corresponderían a probabilidades muy intermedias.
- La diagonal corresponde a un modelo que, haya el número de predicciones iguales a 0/1 que haya (lo cual depende del threshold), estos no siguen a los datos actuales y se distribuyen aleatoriamente. Un threshold alto correspondería a puntos de la recta por abajo: (0.1, 0.1), (0.2, 0.2), etc., y viceversa. Que los puntos de la recta estuvieran distribuidos más uniformemente o menos (es decir, más equidistantes o más aglutinados), depende de lo extremas que sean las probabilidades (predict_proba).
- El peor de los modelos (peor que el aleatorio, uno que yerre a drede) sería uno que fuera por y = 0 y luego subiera en vertical en x = 1. En general, por debajo de la diagonal el modelo "falla a propósito", siendo mejor cualquier modelo aleatorio.
- En el mejor de los modelos, en x=0 línea vertical hasta y=1, y el último punto en (x,y)=(1,1)
- El modelo será mejor cuanto mayor sea el área bajo la curva, que puede ir desde 0 en el peor de los casos hasta 1.
A la vista de la curva ROC, un buen threshold podría ser el 4º o hasta el 7º, es decir threshold = 0.07 / thresold = 0.15
No obstante, nos puede interesar acertar un gran número de compras a expensas de un gran porcentaje de falsos positivos, por lo que un threshold = 0.05 sería interesante
from sklearn import metrics
from ggplot import *
prob_yes
espc_1, sensit, _ = metrics.roc_curve(Y_test, prob_yes)
metrics.roc_curve(Y_test, prob_yes) nos devuelve 3 parámetros:
- Array con 1 - especifidad
- Array con Sensibilidad
- Otro que no vamos a utilizar, por lo que ponemos una "_"
La división de thresholds que hace es de 0,01 en 0,01
df = pd.DataFrame({
"esp":espc_1,
"sens":sensit
})
df.head()
ggplot:
Es una librería que se utiliza tanto en Python como en R y que vale la pena echarle un vistazo. Está diseñado para el lenguaje de programación R, que utiliza el signo "+" para añadir características a la gráfica
ggplot(df, aes(x="esp", y="sens")) + geom_line() + geom_abline(linetype="dashed")+xlim(-0.01,1.01)+ylim(-0.01,1.01)+xlab("1-Especifidad")+ylab("Sensibilidad")
Área bajo la curva:
auc = metrics.auc(espc_1, sensit)
auc #Area Under the Curve
ggplot(df, aes(x="esp", y="sens")) + geom_area(alpha=0.25)+geom_line(aes(y="sens"))+ggtitle("Curva ROC y AUC=%s"%str(auc))
