Cajamar UniversityHack 2018 Datathon: Salesforce Predictive Modelling¶
Donkeys | Universidad de Barcelona¶
- Alex Escolá Nixon
- Jaume Puigbó
- Florent Micand
Agenda¶
- Introducción
- Análisis Exploratorio
- Manipulación de variables
- Modelo y ajuste de parámetros
- Evaluación
- Modelo final y generación de predicciones
- Interpretabilidad y conclusiones
1. Introducción¶
Se entiende por poder adquisitivo la capacidad de una persona por comprar bienes y servicios. Es crucial para una entidad financiera poder estimar el poder adquisitivo de sus clientes para ser capaz de ofrecerles los productos financieros más adecuados en su etapa como clientes de la entidad. Además, es una variable de gran utilidad para establecer el riesgo financiero asociado.
En este proyecto presentamos un modelo predictivo que busca estimar de manera precisa el poder adquisitivo de cada cliente, a partir de un dataset con históricos de productos, consumos, saldos, número de operaciones y datos sociodemográficos de más de 300000 clientes. Se provee un fichero Train.txt para entrenar un modelo de aprendizaje automático de 363834 registros y 88 características, más la variable Poder Adquisitivo, que nos permitirá estimar el poder adquisitivo de los 156315 usuarios del fichero Test.txt.
Este problema de Machine Learning posee algunas particularidades, para las cuales avanzamos algunas hipótesis a continuación. En primer lugar, la variable a estimar, el poder adquisitivo, no es una magnitud real que pueda ser observada, por la propia definición o naturaleza de la misma. Los valores de poder adquisitivo que aparecen en el dataset Train.txt han sido ellos mismos calculados o estimados por la entidad financiera, probablemente a partir de reglas, algoritmos y mayor cantidad de datos, en un amplio rango temporal, utilizando su know-how y experiencia en el sector.
En segundo lugar, dado el orden de magnitud del poder adquisitivo y de otras variables numéricas, entendemos las características del dataset como la situación financiera en un preciso momento de los clientes, una foto instantánea, con algunos de los saldos, consumos, y número de operaciones agregados para un cierto periodo, por ejemplo mensualmente, trimestralmente o anualmente.
Si la entidad financiera ya es capaz de establecer el poder adquisitivo con suficiente fiabilidad como para lanzar un reto nacional, ¿en qué medida podemos aportar valor presentando una solución a este reto? Pensamos que el valor radica en obtener el poder adquisitivo con las limitadas características de este dataset, marcadas por su restricción temporal. Ser capaz de establecer el poder adquisitivo de nuevos clientes, de los que se tengan pocos meses de información, para ofrecerles cuanto antes los productos financieros acordes con su poder adquisitivo, en cuanto a interés del cliente y riesgo financiero.
Esta propuesta nos lleva a reflexionar sobre la evaluación más apropiada de la solución presentada a continuación. Aunque un error en términos absolutos o relativos es un buen indicador, y él que utilizaremos, podríamos aportar mayor valor de negocio conociendo umbrales de concesión de productos en término de poder adquisitivo, y así presentar un sistema de evaluación basado en probabilidad de idoneidad de producto según umbrales.
Presentamos en primer lugar un análisis exploratorio de los datos, que justifican la manipulación de variables llevada a cabo, para después probar diferentes modelos predictivos. Tras comentar la bondad del modelo final mediante el sistema de evaluación elegido, generamos el fichero de entrega con las predicciones de poder adquisitivo para los clientes incluidos en Test.txt.
2. Análisis Exploratorio¶
%reset -f
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.sparse as sp
from sklearn.preprocessing import MinMaxScaler, OneHotEncoder
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.linear_model import LinearRegression
import xgboost as xgb
path = "data/"
train_filename = "Dataset_Salesforce_Predictive_Modelling_TRAIN.txt"
df = pd.read_csv(path+train_filename,low_memory=False)
#Remplazamos un valor hexadecimal en la variable categórica Socio_Demo_01 por su correspondiente decimal
#para que no de problemas de tipo de datos luego.
df["Socio_Demo_01"].replace(to_replace='0X301', value=769, inplace=True)
df["Socio_Demo_01"] = df["Socio_Demo_01"].astype(np.float64)
2.1 Poder adquisitivo¶
Empecemos por analizar los parámetros estadísticos y distribución del poder adquisitivo:
df['Poder_Adquisitivo'].describe()
fig, ax =plt.subplots(1,1, figsize=(8,5))
facet =sns.boxplot(y="Poder_Adquisitivo", data=df, palette="RdBu_r", ax = ax)
facet.set(facecolor ="white")
plt.ylim(ymin=0,ymax=50000)
plt.show()
fig, ax =plt.subplots(1,1, figsize=(8,5))
sns.distplot(df['Poder_Adquisitivo'][df['Poder_Adquisitivo'] < 100000], ax = ax, bins = 120)
plt.xlim(xmax=50000)
plt.show()
Lo primero que podemos observar que la distribución tiene una larga cola derecha. Con lo que una buena opción, para quitar parte del peso que los outliers pueden suponer en el modelo, es transformar la variable dependiente mediante el logaritmo. De esta forma tendremos una distribución que se asemeja más a una distribución normal, por lo tanto más balanceada.
Debido a que hay una cierta cantidad de datos mucho más elevados que los demás, hemos considerado que la solución que podria tener más lógica para este problema, desde el punto de vista entidad financiera, es la de tratar estos datos por separado. Además estos valores suponen un importante empeoramiento de los resultados del regresor utilizado.
Teniendo ello en cuenta, hemos decidido prescindir de los clientes cuyo poder adquisitivo sea superior a aproximadamente los 200.000€, lo que supone quitar el percentil 0.1% de la parte superior de la distribución. Esto podría empeorar la predicción del modelo para clientes de poder adquisitivo muy alto, proponiendo por ejemplo valores cercanos a 200000€ y no a 1 millón. Este error no es relevante teniendo en cuenta que la entidad financiera seguramente provee una atención personalizada a clientes de muy alto poder adquisitivo a través de un asesor financiero.
fig, ax =plt.subplots(1,1, figsize=(8,5))
sns.distplot(np.log(df['Poder_Adquisitivo'][df['Poder_Adquisitivo'] < 100000]), ax = ax, bins = 120)
plt.show()
2.2 Tenencia de los productos financieros </h3>¶
Realizamos un analisis exploratorio de las variables categoricas que indican la tenencia de los distintos productos financieros. Observamos las distribuciones de la cantidad de ceros, unos y doses que contienen cada una de la variables de Ind_Prod, las cuales indican respectivamente si no existe tenencia en los distintos productos financieros, si sí existe, o si ha existido en los últimos tres meses.
Ind_Prods = df.loc[:,'Ind_Prod_01':'Ind_Prod_24']
for col in Ind_Prods.columns:
Ind_Prods.loc[:,col] = Ind_Prods.loc[:,col].astype(float)
fig, axes =plt.subplots(6,4, figsize=(13,13), sharex=True)
axes = axes.flatten()
for ax, catplot in zip(axes, Ind_Prods.columns):
sns.countplot(y=catplot, data=Ind_Prods, ax=ax,
order=np.unique(Ind_Prods.values), palette="husl")
ax.set_xlabel('')
ax.set_xticklabels([])
plt.tight_layout()
plt.show()
Por lo general, son pocos los casos en los que exista tenencia de productos financieros. La cantidad de doses en todos los casos es todavía más minoritaria. Teniendo en cuenta esta información, veamos más en detalle que utilidad pueden tener cada una de estas variables de cara a predecir el poder adquisitivo de un cliente.
Para ello en la siguiente figura vemos para cada variable, la distribución del poder adquisitivo correspondiente a los clientes seccionados según los distintos valores que tóma cada Ind_Prod, es decir, 0,1 o 2.
fig, axes = plt.subplots(6,4, figsize=(13,13), sharex=False)
axes = axes.flatten()
for ax, catplot in zip(axes, Ind_Prods.columns):
data_temp = pd.concat([df['Poder_Adquisitivo'], df[catplot]], axis=1)
sns.boxplot(x=catplot, y="Poder_Adquisitivo", data=data_temp, ax=ax, palette="husl")
ax.set_ylim(bottom=0,top=55000)
plt.tight_layout()
plt.show()
Si dejamos de lado por un momento los diagramas de cajas correspondientes a los doses de cada variable, vemos que las correspondientes distribuciones del poder adquisitivo varían bastante dependiendo de la tenencia de los distintos productos financieros. Ello nos está indicando por lo tanto que por lo general (exceptuando algunos casos como por ejemplo 'Ind_Prod_18'), dichas variables pueden resultar útiles para predecir el poder adquisitivo.
También se observa que en la mayoría de casos la distribución del poder adquisitivo es parecida para los casos en los que sí se posee el producto financiero y para aquellos en los que la ha habido los últimos 3 meses. Por lo que, si añadimos el hecho de que la cantidad de doses es muy minoritaria, parece una buena opción tratar éstos dos útlimos como una sola variable.
Por lo que a excepción de Ind_Prod_24, para la que sí se aprecia un importante augmento de media del poder adquisitivo en los clientes con un 2, los doses serán tratados como unos, de esta forma, debido a que las variables categóricas deben ser convertidas a OneHot, pasaríamos de tener tres columnas para cada variable (una para cada posible valor), a una sola.
A continuación analicemos si existen correlaciones entre productos:
corr = Ind_Prods.corr()
mask = np.zeros_like(corr)
mask[np.triu_indices_from(mask)] = True
with sns.axes_style("white"):
fig, ax = plt.subplots(figsize = (15,15))
sns.heatmap(corr, mask=mask, square=True, annot=False, cmap='Blues', ax=ax)
plt.xticks(rotation=90, ha='center');
plt.show()
A primera vista las correlaciones cruzadas entre las variables son bajas. A excepción de alguna como es la correlación de 'Ind_Prod_08' con 'Ind_Prod_24', con una correlación del 89%:
np.corrcoef(Ind_Prods['Ind_Prod_24'],Ind_Prods['Ind_Prod_08'])
No obstante si analizamos un poco más en detalle ambas variables, nos damos cuenta de que la correlacón se debe a que ambas variables contienen los ceros en las mismas posiciones.:
matches_of_zeros = len(np.intersect1d(list(np.where(Ind_Prods['Ind_Prod_24'] == 0)[0]),
list(np.where(Ind_Prods['Ind_Prod_08'] == 0)[0])))
print 'Posiciones en que coinciden los ceros:', matches_of_zeros
print "Cantidad de ceros:", len(Ind_Prods[Ind_Prods['Ind_Prod_24'] == 0])
Por lo que, teniendo en cuenta que 'Ind_Prod_24' contiene una gran cantidad de doses, y que la distribución del poder adquisitivo es más elevada para los usuarios con un dos como se ha argumentado previamente, ello nos indica que dicha información puede resultar valiosa de cara a la predicción. Por lo que podemos suponer que la variable 'Ind_Prod_08' no nos está aportando información disponiendo de 'Ind_Prod_24'.
También vemos una correlación significativa entre las variables 'Ind_Prod_14' y 'Ind_Prod_22'. En este caso se debe a que las posiciones contienen unos en las mismas posiciones. Pero los demás valores que contienen las categorías aportan información distinta por lo que no podemos realizar asunciones como con el caso previo.
Observando la gráfica de diagramas de caja de los valores categóricos vs el poder adquisitivo, vemos que 'Ind_Prod_18', además de contener unos casi en la totalidad de los casos, las distribuciones de los ceros y unos en dicha variable son muy parecidas. Por lo que cabe suponer que esta variable no aporta mucha información de cara a la predicción.
2.3 Importe de los saldos de los distintos productos financieros </h3>¶
Para el análisis de dichas variables, empecemos por ver la cantidad de ceros que tiene cada una:
Imp_Sals = df.loc[:,'Imp_Sal_01':'Imp_Sal_21']
for i in Imp_Sals.columns:
print i, 'contiene', len(df[i][df[i]!=0]), 'valores distintos a 0'
Vemos claramente, que muchas de ellas contienen muy poca información y no aportarían valor de cara a la regresión. Ello ha sido además corroborado con el análisis de la importancia de cada variable obtenido a partir del mismo regresor XGBoostRegressor con el que se ha realizado la predicción. No obstante, como comentamos, ello será corroborado con la importancia de las variables obtenidas con el mismo modelo, ya que, para algunos casos como es el de Imp_Sal_07, aun contener muy poca información sí resulta útil para el modelo. Por lo que dichas variables no serán utilizadas en el modelo.
to_drop = ['Imp_Sal_01', 'Imp_Sal_03', 'Imp_Sal_11', 'Imp_Sal_14', 'Imp_Sal_17']
low_values = ['Imp_Sal_06', 'Imp_Sal_07', 'Imp_Sal_12', 'Imp_Sal_13', 'Imp_Sal_18']
Imp_Sal_hist = pd.concat([Imp_Sals, df.loc[:,'Poder_Adquisitivo']], axis = 1)\
.drop(to_drop, axis = 1).drop(low_values, axis = 1)
A continuación veamos los histogramas de las siguientes variables:
fig, axes = plt.subplots(4,3, figsize=(10,10), sharex=False)
axes = axes.flatten()
for ax, catplot in zip(axes, Imp_Sal_hist.columns[:-1]):
train_subplot = Imp_Sal_hist[Imp_Sal_hist[catplot] != 0]
data_to_plot = train_subplot[catplot][train_subplot[catplot] < np.percentile(train_subplot[catplot], 90)]
data_to_plot.hist(bins = 250, ax=ax)
ax.set_title(catplot)
ax.xaxis.set_ticks(np.arange(min(data_to_plot), max(data_to_plot)+1, (max(data_to_plot) - min(data_to_plot)+1)/5.))
ax.set_xlim(xmin = data_to_plot.min(), xmax=np.percentile(data_to_plot, 90))
plt.tight_layout()
plt.show()
Lo primero que observamos es la evidente asimetría estadística que poseen todas las variables. Vemos que el importe de los saldos de los distintos productos financieros son por lo general bajos, con una frecuencia menguante para importes más altos. Por lo que también podría plantearse la tranformación logarítmica de estas variables, con tal de quitar peso a la mayor variabilidad que tendrán los valores más elevados de éstas, para ver si benefician al modelo.
Log_Imp_Sal = [u'Imp_Sal_02', u'Imp_Sal_08',
u'Imp_Sal_09', u'Imp_Sal_10', u'Imp_Sal_15', u'Imp_Sal_16',
u'Imp_Sal_19', u'Imp_Sal_20', u'Imp_Sal_21']
def log_transform(df,col):
return pd.DataFrame(np.log(df[col] + 1 + abs(min(df[col]))), columns = [col])
fig, axes = plt.subplots(3,3, figsize=(10,10), sharex=False)
axes = axes.flatten()
for ax, catplot in zip(axes, Log_Imp_Sal):
train_subplot = Imp_Sal_hist.loc[Imp_Sal_hist[catplot] != 0]
train_subplot.loc[:,catplot] = log_transform(train_subplot, catplot)
train_subplot.loc[:,catplot].hist(bins = 250, ax=ax)
ax.set_xlim(xmin = np.percentile(train_subplot[catplot], 2), xmax=np.percentile(train_subplot[catplot], 99.9))
ax.set_title(catplot)
plt.tight_layout()
plt.show()
Vemos que de esta forma también resulta más fácil extraer información relevante de las distribibuciones al asemejarse más a distribuciones normales. Vemos por ejemplo que Imp_Sal_20 posee una distribución claramente bimodal:
fig, ax =plt.subplots(1,1, figsize=(8,5))
log_transform(df,'Imp_Sal_20').hist(bins = 400, ax=ax)
plt.xlim(xmin=4, xmax = 13)
plt.ylim(ymax = 3000)
plt.show()
Por lo que, teniendo en cuenta que ésta es una variable importante para el modelo como mostramos más adelante, una opción que hemos analizado es tratar ambas modas de dicha variable por separado, con lo que cada una se asemejaría más a una distribución normal. En este preciso caso no ha mejorado la predicción.
Analicemos ahora algunas de las correlaciones existentes entre las variables de saldo más destacadas:
sns.set()
g = sns.pairplot(df[['Imp_Sal_08', 'Imp_Sal_09', 'Imp_Sal_19', 'Imp_Sal_21']], size = 3)
g.set(xticklabels=[], facecolor ="white")
plt.show();
Observamos que existe cierta relación entre las variables. Por ejemplo, el importe del salario del producto financiero 19 marca un cierto umbral inferior en el del 21. Sin embargo, la correlación entre ambas no es suficientemente alta como para no utilizar alguna de ellas.
No obstante sí existe una alta correlación el 8 i el 9:
np.corrcoef(df.Imp_Sal_08,df.Imp_Sal_09)
Al estar correlacionadas casi al 90%, en este caso sí podriamos considerar el utilizar tan solo una de ellas para evitar que ello afecte al modelo, o dar demasiado peso a esta variable que, en cierto modo podemos considerar que aporta la misma información. Sin embargo la incidencia de tener variables correlacionadas es posteriormente analizada a partir de los resultados debido a que el XGBRegressor es robusto a correlaciones entre variables.
Viendo más en detalle la distribución de ambas, vemos que son casi idénticas:
fig, ax = plt.subplots(figsize = (8,5))
for a in [log_transform(df,'Imp_Sal_08'), log_transform(df,'Imp_Sal_09')]:
g = sns.distplot(a, bins=50, ax=ax, label=a.columns[0])
g.set(facecolor = "white")
ax.set_xlim([2, 12])
ax.set_xlabel('')
plt.legend()
plt.show()
También analizando las correlaciones de las variables con el poder adquisitivo, hemos considerado interesante comentar el comportamiento de Imp_Sal_21. Como podemos observar en el siguiente diagrama de dispersión de dicha variable con el poder adquisitivo, vemos que el importe del salario del producto 21 marca un cierto umbral inferior de poder adquisitivo. Ello nos indica que dicha variable posiblemente sea útil de cara a predecir el poder adquisitivo, ya que para ciertos importes el poder adquisitivo, como se observa, no baja de cierto mínimo.
fig, ax = plt.subplots(figsize = (8,5))
g = df.plot.scatter(x='Imp_Sal_21', y='Poder_Adquisitivo', ax = ax)
g.set(facecolor = "white")
plt.xticks(np.arange(min(df.Imp_Sal_21), max(df.Imp_Sal_21)+1, (max(df.Imp_Sal_21)+1 - min(df.Imp_Sal_21))/5.))
plt.ylim(ymin=-50000)
plt.xlim(xmin=-50000)
plt.show()
2.4 Importe de consumos habituales del cliente en base a sus operaciones con tarjetas y domiciliaciones más comunes </h3>¶
Veamos que distribuciones tienen estas variables:
Imp_Cons = df.loc[:,'Imp_Cons_01':'Imp_Cons_17']
fig, ax = plt.subplots(figsize = (15,8))
g = sns.boxplot(data=np.array(Imp_Cons), ax = ax)
g.set(facecolor = "white")
ax.set(ylim=(-10, 600))
plt.xticks(np.arange(len(Imp_Cons.columns)), (Imp_Cons.columns), rotation = 90)
plt.tick_params(axis='both', which='major', labelsize=10)
plt.show()
Toman rangos de valores bastante distintos, con algunas de ellas muy concentradas en 0, debido a la gran cantidad de ceros que tienen algunas de ellas.
Veamos las distribuciones de estas variables en forma de histogramas también, tomando la transformación logarítmica para una mejor observación de sus respectivas distribuciones:
fig, axes = plt.subplots(5,4, figsize=(13,13), sharex=False)
axes = axes.flatten()
for ax, catplot in zip(axes, Imp_Cons.columns):
train_subplot = Imp_Cons[Imp_Cons[catplot] != 0]
train_subplot.loc[:,catplot] = np.log(train_subplot[catplot] + 1 + abs(min(train_subplot[catplot])))
g = train_subplot[catplot][train_subplot[catplot] < np.percentile(train_subplot[catplot], 99)].hist(bins = 150, ax=ax)
g.set(facecolor = "white")
ax.set_xlim(xmin = train_subplot[catplot].min(), xmax=np.percentile(train_subplot[catplot], 99))
ax.set_title(catplot)
plt.tight_layout()
plt.show()
Veamos ahora si existen correlaciones destacables entre las variables:
corr = Imp_Cons.corr()
mask = np.zeros_like(corr)
mask[np.triu_indices_from(mask)] = True
with sns.axes_style("white"):
fig, ax = plt.subplots(figsize = (15,15))
sns.heatmap(corr, mask=mask, square=True, annot=False, cmap='Blues', ax = ax)
plt.xticks(rotation=90, ha='center');
plt.show()
En este caso vemos que existe muy poca correlación entre las variables.
2.5 Datos socio-demográficos¶
Como viene especificado, Socio_Demo_01 es una variable categórica. Por lo que una opción es generar dummies con cada uno de los posible valores:
len(df.Socio_Demo_01.unique())
Sin embargo ello sigue suponiendo generar más de 500 variables nuevas. Por lo que una opción es quedarse con la información más relevante mediante un bucketing de la variable.
Para explicar mejor el proceso, veamos como depende el poder adquisitivo del cliente en función de dicha variable.
mean_adq_power = df.groupby('Socio_Demo_01').mean()['Poder_Adquisitivo']
fig, ax = plt.subplots(1,1, figsize = (8,5))
ax.set_facecolor('white')
g = plt.scatter(mean_adq_power.index,mean_adq_power.values)
plt.ylim(ymin=0,ymax=60000)
plt.xlim(xmin=0,xmax=10000)
plt.xlabel('Socio Demo 01 (rango numerico)')
plt.ylabel('Poder adquisitivo')
plt.show()
Vemos que existe una gran variabilidad en el poder adquisitivo en función del valor categórico que toma Socio_Demo_01. Por lo que para no desperdiciar información relevante que puede proporcionar esta variable, una opción es seleccionar las categorias de que toma dicha variable para las cuales el poder adquisitivo es más extremo.
Para ello hemos creado dos variables nuevas, una en la cual hemos introducido unos si la categoría forma parte de las 100 categorías con mayor poder adquisitivo de media, y la otra variable equivalentemente pero para el menor poder adquisitivo.
def buckets(df):
sumation_SD1 = df.Poder_Adquisitivo.groupby(df.Socio_Demo_01).sum()
count_SD1 = df.Poder_Adquisitivo.groupby(df.Socio_Demo_01).count()
avg_sum_SD1 = sumation_SD1/count_SD1
temp = avg_sum_SD1.argsort().reset_index()
bucket_1 = temp.loc[temp.Poder_Adquisitivo,:].Socio_Demo_01.values[-30:]
bucket_2 = temp.loc[temp.Poder_Adquisitivo,:].Socio_Demo_01.values[-100:-30]
bucket_3 = temp.loc[temp.Poder_Adquisitivo,:].Socio_Demo_01.values[:100]
return bucket_1,bucket_2,bucket_3
bucket_1,bucket_2,bucket_3 = buckets(df)
mean_adq_power = df.groupby('Socio_Demo_01').mean()['Poder_Adquisitivo']
high = mean_adq_power[bucket_1]
medium_high = mean_adq_power[bucket_2]
low = mean_adq_power[bucket_3]
others = mean_adq_power[(~mean_adq_power.index.isin(bucket_1)) & (~mean_adq_power.index.isin(bucket_2))\
& (~mean_adq_power.index.isin(bucket_3))]
fig, ax = plt.subplots(1,1, figsize = (8,5))
ax.set_facecolor('white')
blue = plt.scatter(others.index,others.values, c='b')
yellow = plt.scatter(medium_high.index,medium_high.values, c='y')
red = plt.scatter(high.index,high.values, c='r')
green = plt.scatter(low.index,low.values, c='g')
plt.legend((red,yellow,blue,green),
('High AP','Medium-High','Interm.', 'Low AP'))
plt.ylim(ymin=0,ymax=60000)
plt.xlim(xmin=0,xmax=10000)
plt.xlabel('Socio Demo 01 (rango numerico)')
plt.ylabel('Poder adquisitivo')
plt.show()
De esta forma evitamos generar una gran cantidad de variables nuevas y mantenemos cierta información que nos puede proporcionar esta variable.
2.6 Número de operaciones a través de los distintos productos financieros </h3>¶
Empecemos por ver primero, la cantidad de entradas distintas a 0 que tienen estas variables:
Num_Opers = df.loc[:,'Num_Oper_01':'Num_Oper_20']
for i in df.loc[:,'Num_Oper_01':'Num_Oper_20'].columns:
print i, 'has', len(df[i][df[i]!=0]), 'entries different than 0'
Vemos que, del mismo modo que con el importe de los saldos de los productos, hay algunas con muy poca información, por lo que podemos prescindir de ellas.
Veamos ahora las distribuciones de estas variables:
fig, axes = plt.subplots(5,4, figsize=(13,13), sharex=False)
axes = axes.flatten()
for ax, catplot in zip(axes, Num_Opers.columns):
train_subplot = Num_Opers[Num_Opers[catplot] != 0]
g = train_subplot[catplot][train_subplot[catplot] < np.percentile(train_subplot[catplot], 99)].hist(bins = 150, ax=ax)
g.set(facecolor = "white")
ax.set_xlim(xmin = train_subplot[catplot].min(), xmax=np.percentile(train_subplot[catplot], 99))
ax.set_title(catplot)
plt.tight_layout()
plt.show()
2.7 Correlaciones con el poder adquisitivo </h3>¶
corr = df.corr()
k = 10 #number of variables for heatmap
cols = corr.nlargest(k, 'Poder_Adquisitivo')['Poder_Adquisitivo'].index
mask = np.zeros_like(df[list(cols.values)].corr())
mask[np.triu_indices_from(mask, k=1)] = True
cm = np.corrcoef(df[cols].values.T)
sns.set(font_scale=1.25)
with sns.axes_style("white"):
fig, ax = plt.subplots(figsize = (13,13))
hm = sns.heatmap(cm, mask=mask, annot=True, square=True, annot_kws={'size': 10}, yticklabels=cols.values, xticklabels=cols.values,cmap='Blues')
plt.show()
Esta matriz de correlación nos está dando una visión de qué variables están más correlacionadas con el poder adquisitivo. Vemos que, las variables que más correlación muestran con la variable dependiente son algunos de los importes de saldos de distintos productos. En especial el 8 y 9, los cuales no es de sorprender que tengan valores similares de correlación debido a que también están altamente correlacionados entre sí. Y vemos también que el número de operaciones 5 y 6 también muestran cierta correlación, por lo que presumiblemente también serán importantes para el modelo.
También hemos probamos de añadir nuevas variables cogiendo la suma, media, desviación estándar, etc de los grupos de las variables numéricas de Imp_Cons e Imp_Sal. Por otra parte, probamos de hacer un k-means de 10 y 11 clústers de dichas variables respectivamente. Elejimos estas k's despúes de una evaluación de la curva Elbow.
Finalmente, no usamos ninguna de estas variables ya que no ayudaban al modelo.
3. Manipulación de variables¶
%reset -f
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.sparse as sp
from sklearn.preprocessing import MinMaxScaler, OneHotEncoder
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.model_selection import TimeSeriesSplit
from sklearn.linear_model import LinearRegression
import pickle
import operator
import xgboost as xgb
path = "data/"
train_filename = "Dataset_Salesforce_Predictive_Modelling_TRAIN.txt"
df = pd.read_csv(path+train_filename,low_memory=False)
#Remplazamos un valor hexadecimal en la variable categórica Socio_Demo_01 por su correspondiente decimal
#para que no de problemas de tipo de datos luego.
df["Socio_Demo_01"].replace(to_replace='0X301', value=769, inplace=True)
df["Socio_Demo_01"] = df["Socio_Demo_01"].astype(np.float64)
A partir de los insights encontrados en el análisis exploratorio, aplicamos las siguientes manipulaciones de variables:
- Drop de variables no relevantes para el modelo
- Drop de usuarios cuyo poder adquisitivo es mayor que el percentil 99.9%
- Conversión de "2" a "1" en productos (excepto producto 24)
- Normalización de variables numéricas
- OneHot Encoding de variables catégoricas no binarias
- Conversión logarítmica del poder adquisitivo
- Buckets para la variable Socio_Demo_01 en función de su relación con el poder adquisitivo
- Creación de variables comentadas en la sección anterior
# NA
print df.shape
df = df.dropna().reset_index(drop = True)
print "Duplicated: "+str(np.any(df.duplicated().values))
print df.shape
# Outliers
print df.shape
df = df[df['Poder_Adquisitivo'] < np.percentile(df['Poder_Adquisitivo'], 99.9)]
print df.shape
En un primer momento, realizamos un split temporal para calcular los buckets de la variable Socio_Demo_01 sobre el futuro conjunto de train unicamente; calculando estos buckets sin usar el futuro conjunto de test.
df_train_temp, df_test_temp, y_train_temp, y_test_temp = train_test_split(df, df["Poder_Adquisitivo"].values, random_state=44)
# Función auxiliar para generar buckets de Socio_Demo_01 usuada durante el análisis exploratorio
def buckets(df):
sumation_SD1 = df.Poder_Adquisitivo.groupby(df.Socio_Demo_01).sum()
count_SD1 = df.Poder_Adquisitivo.groupby(df.Socio_Demo_01).count()
avg_sum_SD1 = sumation_SD1/count_SD1
temp = avg_sum_SD1.argsort().reset_index()
bucket_1 = temp.loc[temp.Poder_Adquisitivo,:].Socio_Demo_01.values[-30:]
bucket_2 = temp.loc[temp.Poder_Adquisitivo,:].Socio_Demo_01.values[-100:-30]
bucket_3 = temp.loc[temp.Poder_Adquisitivo,:].Socio_Demo_01.values[:100]
bucket_others = np.setdiff1d(np.unique(df_train_temp["Socio_Demo_01"].values), np.append(np.append(bucket_1,bucket_2),bucket_3))
return bucket_1,bucket_2,bucket_3,bucket_others
bucket_1,bucket_2,bucket_3,bucket_others = buckets(df_train_temp)
Reproducimos el split anterior con 2/3 para el conjunto de train y 1/3 para el conjunto de test, ya de manera definitiva, con el mismo random_state anterior.
#Labels con transformacion logaritmica
y = np.log(df["Poder_Adquisitivo"].values)
print(y.shape)
#Dataset principal sin labels
X = df.drop("Poder_Adquisitivo", axis = 1)
print(X.shape)
#Split 2/3 y 1/3
df_train, df_test, y_train, y_test = train_test_split(X, y, random_state=44)
Detallamos los grupos de variables que usaremos en el modelo.
categorical = np.array([u'Ind_Prod_01', u'Ind_Prod_02',
u'Ind_Prod_03', u'Ind_Prod_04', u'Ind_Prod_05', u'Ind_Prod_06',
u'Ind_Prod_07', u'Ind_Prod_08', u'Ind_Prod_09', u'Ind_Prod_10',
u'Ind_Prod_11', u'Ind_Prod_12', u'Ind_Prod_13', u'Ind_Prod_14',
u'Ind_Prod_15', u'Ind_Prod_16', u'Ind_Prod_17', u'Ind_Prod_18',
u'Ind_Prod_19', u'Ind_Prod_20', u'Ind_Prod_21', u'Ind_Prod_22',
u'Ind_Prod_23',"Socio_Demo_02"])
categorical_to_onehot = np.array(['Ind_Prod_24',"Socio_Demo_01"])
numerical = np.array([u'Imp_Cons_01', u'Imp_Cons_02', u'Imp_Cons_03',
u'Imp_Cons_04', u'Imp_Cons_05', u'Imp_Cons_06', u'Imp_Cons_07',
u'Imp_Cons_08', u'Imp_Cons_09', u'Imp_Cons_10', u'Imp_Cons_11',
u'Imp_Cons_12', u'Imp_Cons_13', u'Imp_Cons_14', u'Imp_Cons_15',
u'Imp_Cons_16', u'Imp_Cons_17', u'Imp_Sal_02',
u'Imp_Sal_04', u'Imp_Sal_05', u'Imp_Sal_06',
u'Imp_Sal_07', u'Imp_Sal_08', u'Imp_Sal_09', u'Imp_Sal_10',
u'Imp_Sal_12', u'Imp_Sal_13', u'Imp_Sal_15',
u'Imp_Sal_16', u'Imp_Sal_18', u'Imp_Sal_19',
u'Imp_Sal_21', u'Num_Oper_01', u'Num_Oper_02',
u'Num_Oper_03', u'Num_Oper_05', u'Num_Oper_06',
u'Num_Oper_07', u'Num_Oper_08', u'Num_Oper_09', u'Num_Oper_10',
u'Num_Oper_13', u'Num_Oper_14',
u'Num_Oper_15', u'Num_Oper_17', u'Num_Oper_18',
u'Num_Oper_19', u'Socio_Demo_03',
u'Socio_Demo_04', u'Socio_Demo_05'])
Aplicamos las transformaciones mencionadas, sobre ambos conjuntos. La función siguiente devuelve matrices numpy en formato sparse, lo que incrementará la velocidad del modelo.
#Transformaciones
def transform(df,numerical,categorical,minmax,onehot,bucket_1,bucket_2,bucket_3,bucket_others, train = True):
# Restamos 1 a la categórica socio_demo_02 para que tenga la forma 0,1 y no 1,2
df.loc[:,"Socio_Demo_02"] = df.loc[:,"Socio_Demo_02"]-1
# Remplazamos 2 por 1 en los productos con pocos 2
for col in df.loc[:,'Ind_Prod_01':'Ind_Prod_23'].columns:
df.loc[:,col] = df.loc[:,col].astype(int)
df.loc[:,col].replace(to_replace=2, value=1, inplace=True)
# Bucketing de Socio_Demo_01
df["Socio_Demo_01"].replace(to_replace = bucket_1, value = 1, inplace = True)
df["Socio_Demo_01"].replace(to_replace = bucket_2, value = 2, inplace = True)
df["Socio_Demo_01"].replace(to_replace = bucket_3, value = 3, inplace = True)
df["Socio_Demo_01"].replace(to_replace = bucket_others, value = 0, inplace = True)
if train == True:
#MinMax Scaling Fit de las features numéricas y OneHot Encoding Fit de las features
#a las que aplica, entrenadas en el conjunto de train unicamente
minmax = MinMaxScaler()
minmax.fit(df[numerical])
onehot = OneHotEncoder()
onehot.fit(df[categorical_to_onehot])
#MinMax Scaling de las features numéricas y en sparse
df_num_scaled = minmax.transform(df[numerical])
df_num_sparse = sp.csr_matrix(df_num_scaled)
print "Features Numéricas: "+str(df_num_sparse.shape[1])
# OneHot de features categoricas
f_onehot_names = []
for name in categorical_to_onehot:
for i in range(len(df[name].unique())):
aux = name + '_' + str(i)
f_onehot_names.append(aux)
df_cat_encoded = onehot.transform(df[categorical_to_onehot])
print "Features Categóricas sin OneHot encoding: "+str(df[categorical].shape[1])
print "Features Categóricas con OneHot encoding: "+str(df_cat_encoded.shape[1])
# Concatenamos las matrices en formato sparse
X = sp.hstack([df_num_sparse,sp.csr_matrix(df[categorical]),df_cat_encoded])
f_names = np.hstack([numerical,categorical,f_onehot_names])
print "Dataset: "+str(X.shape)
return X, f_names, minmax, onehot
minmax = []
onehot = []
x_train, f_names, minmax1, onehot1 = transform(df_train,numerical,categorical,minmax,onehot,bucket_1,bucket_2,bucket_3,bucket_others)
x_test, trash1, trash2, trash3 = transform(df_test,numerical,categorical,minmax1,onehot1,bucket_1,bucket_2,bucket_3,bucket_others, train = False)
4. Modelo y ajuste de parámetros¶
Este problema estándar de aprendize supervisado para obtener una regresión puede plantearse con distintos modelos muy conocidos, desde una simple regresión linear, que planteamos a continuación a modo de baseline, hasta aplicar Random Forest, XGBoost, o Redes Neuronales, entre otras opciones. Tras probar un modelo Random Forest y una red neuronal de varias capas densas, optamos por el Extreme Gradient Boosting Regressor o XGBoost, un modelo que ha destacado durante los últimos anos en competiciones Kaggle, y con él que hemos obtenido mejores resultados. XGBoost clasifica ejemplos mediante el uso de un conjunto de árboles de decisión. Los árboles se construyen secuencialmente, añadiendo en cada iteración el árbol que mejor compense por los errores de los árboles ya existentes. Se le llama método de gradiente porque el modelo evoluciona en dirección al menor error, árbol a árbol.
Aplicamos un sencillo modelo de regresión linear como baseline.
LR = LinearRegression()
LR.fit(x_train,y_train)
y_pred = LR.predict(x_test)
# Evaluacion: Mean Absolute Error
def mae(true, pred):
y_diff = np.absolute(true - pred)
print "MAE: " + str(np.mean(y_diff))
mae(np.exp(y_test), np.exp(y_pred))
Se obtiene un error absoluto medio cercano a los 5600€, sabiendo que el promedio de poder adquisitivo es de 16000€ y la mediana de 13000€. Utilizamos esta regresión linear como característica de entrada al modelo XGBoost ya que mejorará su capacidad predictora.
LR = LinearRegression()
LR.fit(x_train,y_train)
lr_train = pd.DataFrame(data=LR.predict(x_train),columns = ["lr"])
lr_test = pd.DataFrame(data=LR.predict(x_test),columns = ["lr"])
print lr_train.shape
print lr_test.shape
x_train = sp.hstack([x_train,lr_train])
x_test = sp.hstack([x_test,lr_test])
XGBoost posee varios hyperparámetros que deben ser tuneados para encontrar la combinación de aquellos que conducen al menor error absoluto medio: este proceso se denomina Grid Search. Lo ejecutamos a continuación, además de una validación cruzada (cross validation) para asegurarnos una correcta generalización del resultado. Además, utilizamos Early Stopping de Gradient Descent para evitar overfitting.
En la celda siguiente lanzamos un Grid Search de ejemplo, ya con los mejores parámetros obtenidos en una ejecucción anterior, para que la presente ejecucción no tarde demasiado.
#Grid Search original
paramGrid = {'learning_rate': [0.01, 0.05, 0.1], 'max_depth': [15,20,25], 'min_child_weight':[3,5],
'gamma': [0.3,1,3], 'subsample':[0.75], 'colsample_bytree':[0.75],
'n_estimators': [2000]}
#Grid Search Best Parameters (Comentar la linea siguiente en caso de querer
#ejecutar el Grid Search original)
paramGrid = {'learning_rate': [0.01], 'max_depth': [25], 'min_child_weight':[3],
'gamma': [0.3], 'subsample':[0.75], 'colsample_bytree':[0.75],
'n_estimators': [1000]}
fit_params={"early_stopping_rounds": 20,
"eval_metric" : "mae",
"eval_set" : [[x_test, y_test]]}
cv = 5
model = xgb.XGBRegressor(nthread=-1)
grid = GridSearchCV(model, paramGrid, scoring = 'neg_mean_absolute_error', verbose=3,
fit_params=fit_params,
cv=TimeSeriesSplit(n_splits=cv).get_n_splits([x_train,y_train]))
grid.fit(x_train,y_train)
print grid.best_params_
Sería conveniente explorar más combinaciones de parámetros para conseguir un ajuste aún más óptimo, pero tenemos ciertas limitaciones computacionales. Es un punto a considerar para productizar la solución.
El Grid Search genera los siguientes mejores parámetros:
- colsample_bytree: 0.75
- learning_rate: 0.01
- min_child_weight: 3
- n_estimators: 1000
- subsample: 0.75
- max_depth: 25
- gamma: 0.3
Guardamos el objeto Grid para posibles futuras operaciones.
pickle.dump(grid, open("grid.dat", "wb"))
#grid = pickle.load(open("grid.dat", "rb"))
5. Evaluación¶
Proponemos dos métricas para evaluar el error del modelo :
Mean Absolute Error: $$MAE = \frac{1}{n} \sum_{i=1}^{n} \big|y_{real} - y_{pred}\big|$$
Mean Absolute Percentage Error: $$MAPE = \frac{1}{n} \sum_{i=1}^{n} \frac{\big|y_{real} - y_{pred}\big|}{\big|y_{real}\big|}$$
donde $y_{real} \in \mathbb{R^n}$ son los valores reales, $y_{pred} \in \mathbb{R^n}$ son los valores de la predicción.
Nos hemos decido por estas métricas ya que son las más naturales para medir el error. El MAE está expresado en euros, mientras que el MAPE es un porcentaje de error. El MAPE puede ser más representativo al tener distintos órdenes de magnitud en la variable de poder adquisitivo.
Generamos las predicciones para su evaluación.
y_pred = grid.best_estimator_.predict(x_test)
# MAE
def eval_mae(true, pred, q=1.0) :
y_diff = np.absolute(true - pred)
y_diff_q=pd.Series(y_diff)
y_diff_q = y_diff_q[y_diff_q<y_diff_q.quantile(q)]
mae = np.mean(y_diff_q)
print "MAE (" +str(int(q*100)) + "%): " + str(mae)
plt.hist(y_diff_q, bins = 100)
plt.title('MAE Distribution')
plt.ylabel('Frequency')
plt.xlabel('MAE')
plt.show()
return mae
mae_100 = eval_mae(np.exp(y_pred), np.exp(y_test))
Obtenemos un error medio de 3700 euros aproximadamente. Para ver mejor la distribución del error nos quedamos con aquellos errores que son inferior al percentil 0.95. De esta manera obtenemos un error medio de 2506 euros aproximadamente. Como podemos observar en la siguiente gráfica el 95% de los casos el error es menor a 12000 euros.
# MAE percentile 0.95
mae_95 = eval_mae(np.exp(y_pred), np.exp(y_test), q=0.95)
Veamos ahora los resultados del MAPE y veremos como obtenemos otra interpretación de los resultados.
# MAPE
def eval_mape(true, pred, q=1.0) :
mape = np.absolute(true - pred)/np.abs(true)
y_diff_q=pd.Series(mape)
y_diff_q = y_diff_q[y_diff_q<y_diff_q.quantile(q)]
mape = np.mean(y_diff_q)
print "MAPE (" +str(int(q*100)) + "%): " + str(mape)
plt.hist(y_diff_q, bins = 100)
plt.title('MAPE Distribution')
plt.ylabel('Frequency')
plt.xlabel('MAPE')
plt.show()
return mape
mape_100 = eval_mape(np.exp(y_pred), np.exp(y_test))
Obtenemos un error medio de aproximadamente un 23%. Es decir, la predicción es un 23% mayor o menor respecto al valor real de media. Si nos fijamos en el 95% de los mejores casos obtenemos un MAPE de 18% aproximadamente y podemos ver que el 95% el error del 12000 euros se corresponde como máximo a menos del 70% de error.
# MAPE percentile 0.95
mape_95 = eval_mape(np.exp(y_pred), np.exp(y_test), q=0.95)
En la siguiente gráfica podemos comprobar si nos equivocamos más a la alta o la baja. Realmente no se aprecia que el modelo se decante más hacia un lado que al otro.
#Pred minus real
y_diff = np.exp(y_pred)- np.exp(y_test)
plt.hist(y_diff, bins=1000)
plt.xlim(xmin=-10000, xmax=10000)
plt.title('Error Distribution')
plt.ylabel('Frequency')
plt.xlabel('MAE')
plt.show()
Analizemos ahora donde cometemos los errores, es decir relacionaremos el error con el poder adquisitivo.
# MAE vs Poder Adquisitivo
y_diff = np.abs(np.exp(y_pred)- np.exp(y_test))
plt.scatter(y_diff, np.exp(y_test))
plt.title('MAE vs Poder Adquisitivo')
plt.ylabel('Poder Adquisitivo')
plt.xlabel('MAE')
plt.show()
Como podemos observar a mayor error mayor es el poder adquisitivo como cabe esperar. Veamos ahora el mismo gráfico usando el MAPE.
# MAPE vs Poder Adquisitivo
mape = np.absolute(np.exp(y_test) - np.exp(y_pred))/np.abs(np.exp(y_test))
plt.scatter(mape, np.exp(y_test))
plt.title('MAPE vs Poder Adquisitivo')
plt.ylabel('Poder Adquisitivo')
plt.xlabel('MAPE')
plt.show()
De esta manera podemos ver que cometemos errores más importantes cuándo el valor del poder adquisitivo es pequeño. En cambio, a medida que el valor adquisitivo crece no obtenemos tantos errores grandes proporcionalmente hablando.
Esta situación no nos preocupa si suponemos que la entidad bancaria usa esta predicción en función de umbrales: el modelo puede predicir un poder adquisitivo de 8000€ para un usuario que realmente tiene un poder adquisitivo de 4000€ y cometer un error del 100%, pero no es un problema si el banco propone productos financieros para clientes con un poder adquisitivo por debajo del umbral de 10000€. Esta necesidad de negocio basada en umbrales motiva la creación de un modelo de adicional de evaluación. En la función siguiente podemos establecer el umbral mínimo, umbral máximo, y el "step" o intervalo para crear umbrales en el recorrido entre mínimo y máximo. Para cada intervalo, se obtienen los clientes del conjunto de test cuyo poder adquisitivo esté en el rango, y de éstos se obtienen el subgrupo de clientes cuyo poder adquisitivo predicho también se situa en el mismo rango.
def thresholds_evaluation(y_test, y_pred, minimo, maximo, step):
thresholds = np.arange(minimo,maximo,step)
df_index = []
df_acc = []
df_test = []
df_pred = []
for limit in thresholds:
# límite inferior
if limit == thresholds[0]:
test_indices = np.where(np.exp(y_test)<limit)[0]
test_num = test_indices.shape[0]
pred_num = np.where(np.exp(y_pred[test_indices])<limit)[0].shape[0]
acc = np.round(100.*pred_num/test_num,2)
df_index.append("< "+str(limit)+"€")
df_acc.append(acc)
df_test.append(test_num)
df_pred.append(pred_num)
# límite superior
elif limit == thresholds[-1]:
test_indices = np.where(np.exp(y_test)>limit)[0]
test_num = test_indices.shape[0]
pred_num = np.where(np.exp(y_pred[test_indices])>limit)[0].shape[0]
acc = np.round(100.*pred_num/test_num,2)
##print np.exp(y_test[test_indices])
#print np.exp(y_pred[test_indices])
df_index.append("> "+str(limit)+"€")
df_acc.append(acc)
df_test.append(test_num)
df_pred.append(pred_num)
# límites intermedios
else:
test_indices = np.where((np.exp(y_test)<limit) & (np.exp(y_test)>(limit-step)))[0]
test_num = test_indices.shape[0]
pred_num = np.where(np.exp(y_pred[test_indices])<limit)[0].shape[0]
acc = np.round(100.*pred_num/test_num,2)
df_index.append(str(limit-step)+"€ - "+str(limit)+"€")
df_acc.append(acc)
df_test.append(test_num)
df_pred.append(pred_num)
df = pd.DataFrame(data=df_test,index=df_index, columns=["Muestras de Test"])
df["Muestras de Pred correspondientes"] = df_pred
df["Precisión por Umbral (%)"] = df_acc
return df
Para tener una visión general, mostramos la distribución con un umbral mínimo de 5000€, máximo de 60000€, y los umbrales intermediarios obtenidos cada 5000€. Observamos como el modelo tiene más dificultad con los extremos, sobretodo en el extremo inferior.
thresholds_evaluation(y_test, y_pred, 5000, 60000, 5000)
La distribución del poder adquisitivo centrada entorno a los 13000€ justifica mirar con más detalle esa zona. Reducimos el intervalo a 2000€ por lo que la precisión disminuye, y conseguimos conocer el rendimiento del modelo cuando variaciones de 2000€ no son relevantes desde el punto de visto de negocio.
Podría agregarse un parámetro de tolerancia para tratar mejor los valores que se encuentran muy cerca de los umbrales.
thresholds_evaluation(y_test, y_pred, 5000, 30000, 2000)
6. Modelo final y generación de predicciones¶
Con el objetivo de generar las predicciones definitivas de los clientes incluidos en el fichero Test.txt, entrenamos el modelo con los parámetros más adecuados, obtenidos del Grid Search anterior, sobre el conjunto completo de los datos de Train.txt.
train_filename = "Dataset_Salesforce_Predictive_Modelling_TRAIN.txt"
df = pd.read_csv(path+train_filename,low_memory=False)
#Remplazamos un valor hexadecimal en la variable categórica Socio_Demo_01 por su correspondiente decimal
#para que no de problemas de tipo de datos luego.
df["Socio_Demo_01"].replace(to_replace='0X301', value=769, inplace=True)
df["Socio_Demo_01"] = df["Socio_Demo_01"].astype(np.float64)
# NA
print df.shape
df = df.dropna().reset_index(drop = True)
print "Duplicated: "+str(np.any(df.duplicated().values))
print df.shape
# Outliers
print df.shape
df = df[df['Poder_Adquisitivo'] < np.percentile(df['Poder_Adquisitivo'], 99.9)]
print df.shape
#Creación de buckets para Socio_Demo_01
bucket_1,bucket_2,bucket_3,bucket_others = buckets(df)
#Labels con transformacion logaritmica
y = np.log(df["Poder_Adquisitivo"].values)
#Dataset principal sin labels
df = df.drop("Poder_Adquisitivo", axis = 1)
#Aplicamos las transformaciones
minmax = []
onehot = []
X, f_names, minmax1, onehot1 = transform(df,numerical,categorical,minmax,onehot,bucket_1,bucket_2,bucket_3,bucket_others)
LR = LinearRegression()
LR.fit(X, y)
lr = pd.DataFrame(data=LR.predict(X),columns = ["lr"])
X = sp.hstack([X,lr])
print X.shape
Entrenamos el modelo final.
model = xgb.train({"max_depth": 25,
"learning_rate": 0.01,
"min_child_weight":3,
"gamma":0.3,
"subsample":0.75,
"colsample_bytree":0.75,
"seed":11},
xgb.DMatrix(X, label=y),1000)
# save model to file
pickle.dump(model, open("model.dat", "wb"))
#model = pickle.load(open("model.dat", "rb"))
Cargamos el fichero Test.txt y aplicamos las mismas transformaciones de variables aplicadas sobre el conjunto de Train, teniendo especial cuidado de utilizar la misma función de MinMaxScaler, ya entrenada y empleada en el paso justamente anterior.
test_filename = "Dataset_Salesforce_Predictive_Modelling_TEST.txt"
df_target = pd.read_csv(path+test_filename,low_memory=False)
#Remplazamos un valor hexadecimal en la variable categórica Socio_Demo_01 por su correspondiente decimal
#para que no de problemas de tipo de datos luego.
df_target["Socio_Demo_01"].replace(to_replace='0X301', value=769, inplace=True)
df_target["Socio_Demo_01"] = df_target["Socio_Demo_01"].astype(np.float64)
# Quitamos NA si aplica
print df_target.shape
df_target = df_target.dropna().reset_index(drop = True)
print "Duplicated: "+str(np.any(df_target.duplicated().values))
print df_target.shape
#Aplicamos las transformaciones
X_to_predict, trash1, trash2, trash3 = transform(df_target,numerical,categorical,minmax1,onehot1,bucket_1,bucket_2,bucket_3,bucket_others, train = False)
Obtenemos las predicciones de la regresión linear, lo añadimos al dataset, obtenemos las predicciones del modelo final, y aplicamos una transformación exponencial.
lr_to_predict = pd.DataFrame(data=LR.predict(X_to_predict),columns = ["lr"])
X_to_predict = sp.hstack([X_to_predict,lr_to_predict])
y_pred = model.predict(xgb.DMatrix(X_to_predict))
print y_pred
y_pred = np.exp(y_pred)
print y_pred
Exportamos el resultado a un fichero Test_Mission.txt en el formato solicitado.
result_filename = "Test_Mission.txt"
y_pred_data=pd.DataFrame(data=y_pred, columns=["PA_Est"])
customerid=pd.DataFrame(data=df_target["ID_Customer"])
output = customerid.assign(PA_Est=y_pred_data)
output.to_csv(path+result_filename,sep=",", encoding="utf-8", index=False,quoting=None)
print output.shape
7. Interpretabilidad y conclusiones¶
import shap
shap.initjs()
La explicabilidad de un modelo predictivo es un aspecto clave para que una entidad financiera pueda utilizarlo y ser transparente con sus clientes, más teniendo en cuenta la última reglamentación europea al respecto (GDPR).
El proyecto SHAP (SHapley Additive exPlanations) plantea explicar la salida de cualquier modelo de Machine Learning basado en las contribuciones de las características empleadas en el modelo, basado en teoría de juegos. Puede encontrarse más información en:
https://github.com/slundberg/shap
Obtenemos las contribuciones o "shap values" de cada característica a las predicciones.
shap_values = model.predict(xgb.DMatrix(X_to_predict.todense()[0:1000]), pred_contribs=True)
X_to_predict_df = pd.DataFrame(data=X_to_predict.todense()[0:1000],columns = np.append([f_names],'lr'))
Para una predicción individual, observamos las características que más pesan en incrementar esta predicción por encima del base value (promedio de todas las predicciones) en rojo, y las que más contribuyen a dismimuir el valor de la predicción en azul. Es interesante ver como, aunque algunas características se repiten para diferentes muestras, cada una aporta diferentes contribuciones, de forma positiva o negativa. Debe tomarse el logaritmo de los valores "base value" y "output value" en estas figuras para conseguir el poder adquisitivo en euro.
n=0
shap.force_plot(shap_values[n,:], X_to_predict_df.iloc[n,:])
n=6
shap.force_plot(shap_values[n,:], X_to_predict_df.iloc[n,:])
Para tener una idea genérica de la importancia de las características, mostramos la distribución del impacto que aporta cada caracterítica en la salida del modelo para cada feature, ordenadas por importancia. A mayor importes de los saldos (por ejemplo "Imp_Sal_21"), mayor es la contribución a aumentar el poder adquisitivo. La variable creada de la regresión lineal está entre las más relevantes, al igual que algunos saldos, consumos, y la variable Socio_Demo_03 que muy probablemente sea la edad (por su distribución).
shap.summary_plot(shap_values, X_to_predict_df)
A lo largo de este proyecto hemos realizado una exploración exhaustiva de los datos, que ha justificado una serie de transformaciones y creaciones de nuevas variables, como por ejemplo aplicar el logaritmo al Poder Adquisitivo. Tras probar varios modelos predictivos, aplicamos un modelo de ensamblado de árboles de decisión para regresión, XGBoost, que evaluamos mediante Mean Absolute Error y Mean Absolute Percentage Error principalmente. Entrenamos un modelo final sobre el conjunto completo de train y generamos las predicciones de Poder Adquisitivo sobre el conjunto de test, adjuntas a este documento. Mostramos finalmente las métricas obtenidas.
results = pd.DataFrame(data = [mae_100, mae_95, mape_100, mape_95], index = ["MAE (100%)","MAE (95%)","MAPE (100%)","MAPE (95%)"], columns = ["Evaluation Table"])
results