Les données sont composées de 825 clients d'une banque décrits par 32 variables concernant leurs avoirs, et utilisations de leurs comptes. Après avoir réalisé, avec R ou Python, le premier objectif de segmentation ou profilage des types de comportement des clients, le 2ème consiste à estimer puis prévoir un score d'appétence pour un produit bancaie, ici la carte visa premier. Comparaison des différentes méthodes et algorihtmes d'apprentissage pour atteindre cet objectif de la régression logistique au boosting (extrem gradient) en passant par les arbres, les SVM ou random forest. Une procédure de validation croisée généralisée est itérée sur une selection de ces méthodes. Celles d'agrégation de modèles conduisent aux meilleurs résultats.
Un calepin), qu'il est préférable d'exécuter préalablement, décrit le premier objectif d'exploration puis segmentation ou profilage des types de comportement des clients d'une banque.
Ce deuxième calepin propose de construire un score d'appétence pour la carte Visa Premier. Ce deuxième objectif est intégré à la saison 3 (Apprentissage Statistique). Il s'agit d'un score d'appétence mais ce pourrait être le score d'attrition (churn) d'un opérateur téléphonique ou encore un score de défaillance d'un emprunteur ou de faillite d'une entreprise; les outils de modélisation sont les mêmes et sont très largement utilisés dans tout le secteur tertiaire pour l'aide à la décision.
Remarque: un scénario plus ancien propose une exécution avec SAS, encore utilisé dans de nombreuses grandes entreprises.
La liste des variables est issue d'une base de données retraçant l'historique mensuel bancaire et les caractéristiques de tous les clients. Un sondage a été réalisé afin d'alléger les traitements ainsi qu'une première sélection de variables. Les variables contenues dans le fichier initial sont décrites dans le tableau ci-dessous. Elles sont observées sur 1425 clients.
Tableau: Liste des variables initiales et de leur libellé Attention, certains sont écrits en majuscules dans les programmes puis en minuscules après transfomation des données (logarithme, recodage) au cours d ela phase d'exploration. Les noms des variables logarithmes des variables quantitatives se terminent par Lles variables qualitatives se terminent par Qou q.
| Identifiant | Libellé |
|---|---|
sexeq |
Sexe (qualitatif) |
ager |
Age en années |
famiq |
Situation familiale: Fmar Fcel Fdiv Fuli Fsep Fveu |
relat |
Ancienneté de relation en mois |
pcspq |
Catégorie socio-professionnelle (code num) |
opgnb |
Nombre d'opérations par guichet dans le mois |
moyrv |
Moyenne des mouvements nets créditeurs des 3 mois en Kf |
tavep |
Total des avoirs épargne monétaire en francs |
endet |
Taux d'endettement |
gaget |
Total des engagements en francs |
gagec |
Total des engagements court terme en francs |
gagem |
Total des engagements moyen terme en francs |
kvunb |
Nombre de comptes à vue |
qsmoy |
Moyenne des soldes moyens sur 3 mois |
qcred |
Moyenne des mouvements créditeurs en Kf |
dmvtp |
Age du dernier mouvement (en jours)\hline |
boppn |
Nombre d'opérations à M-1 |
facan |
Montant facturé dans l'année en francs |
lgagt |
Engagement long terme |
vienb |
Nombre de produits contrats vie |
viemt |
Montant des produits contrats vie en francs |
uemnb |
Nombre de produits épargne monétaire |
xlgnb |
Nombre de produits d'épargne logement |
xlgmt |
Montant des produits d'épargne logement en francs |
ylvnb |
Nombre de comptes sur livret |
ylvmt |
Montant des comptes sur livret en francs |
rocnb |
Nombre de paiements par carte bancaire à M-1 |
nptag |
Nombre de cartes point argent |
itavc |
Total des avoirs sur tous les comptes |
havef |
Total des avoirs épargne financière en francs |
| `jnbjd | Nombre de jours à débit à M |
carvp |
Possession de la carte VISA Premier |
Réponde aux questions en s'aidant des résultats des exécutions
Les données sont disponibles dans le répertoire de ce calepin et chargées en même temps. Elles sont issues de la première phase de prétraitement ou data munging pour détecter, corriger les erreurs et incohérences, éliminer des redondances, traiter les données manquantes, transformer certaines variables.
In [ ]:
# Importation des librairies.
import numpy as np
import pandas as pd
import random as rd
import matplotlib.pyplot as plt
%matplotlib inline
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV
In [ ]:
# Lecture d'un data frame
vispremv = pd.read_table('vispremv.dat', delimiter=' ')
vispremv.shape
In [ ]:
vispremv.head()
In [ ]:
# Variables quantitatives
vispremv.describe()
Vérifier ci-dessous que la plupart des variables ont deux versions, l'une quantitative et l'autre qualitative. La version en R de ce calepin compare deux stratégies: l'une basée sur l'utilisation des variables explicatives initiales quantitatives, l'autre sur celles qualitatives après découpage en classes. Compte tenu des résultats et des contraintes de python, ou plus précisement de scikit-learn , il est plus adapté de ne considérer que les variables quantitatives.
Les variables qualitatives (sexe, csp, famille) sont transformées en indicatrices à l'exception de la cible CARVP.
In [ ]:
vispremv.dtypes
In [ ]:
# Transformation en indicatrices
vispremDum=pd.get_dummies(vispremv[["SEXEQ","FAMIQ","PCSPQ"]])
# Une seule est conservée pour les variables binaires
vispremDum.drop(["SEXEQ_Sfem","FAMIQ_Fseu"], axis = 1, inplace = True)
# Sélection des variables numériques
vispremNum = vispremv.select_dtypes(exclude=['object'])
# Concaténation des variables retenues
vispremR=pd.concat([vispremDum,vispremNum],axis=1)
vispremR.columns
Q Combien d'individus et combien de variables sont finalement concernés?
In [ ]:
vispremR.shape
In [ ]:
# La variable à expliquer est recodée
y=vispremv["CARVP"].map(lambda x: 0 if x=="Cnon" else 1)
In [ ]:
rd_seed=111 # Modifier cette valeur d'initialisation
npop=len(vispremv)
xApp,xTest,yApp,yTest=train_test_split(vispremR,y,test_size=200,random_state=rd_seed)
xApp.shape
Cette ancienne méthode reste toujours très utilisée. D'abord par habitude mais aussi par efficacité pour le traitement de données très volumineuses lors de l'estimation de très gros modèles (beaucoup de variables) notamment par exemple chez Criteo ou CDiscount.
La procédure de sélection de modèle est celle par pénalisation: ridge, Lasso ou une combinaison (elastic net). Contrairement aux procédures disponibles en R (stepwise, backward, forward) optimisant un critère comme AIC, l'algorithme proposé dans scikit-learn nepermet pas une prise en compte simple des interactions. D'autre part les compléments usuels (test de Wald ou du rapport de vraisemblance) ne sont pas directement fournis.
In [ ]:
from sklearn.linear_model import LogisticRegression
# Grille de valeurs du paramètre de pénalisaiton
param=[{"C":[0.5,1,5,10,12,15,30]}]
logitL = GridSearchCV(LogisticRegression(penalty="l1"), param,cv=5,n_jobs=-1)
logitLasso=logitL.fit(xApp, yApp)
# Sélection du paramètre optimal
logitLasso.best_params_["C"]
print("Meilleur score (apprentissage) = %f, Meilleur paramètre = %s" %
(1.-logitLasso.best_score_,logitLasso.best_params_))
Erreur de prévision
In [ ]:
# Prévision
yChap = logitLasso.predict(xTest)
# matrice de confusion
table=pd.crosstab(yChap,yTest)
print(table)
# Erreur sur l'échantillon test
print("Erreur de test régression Lasso = %f" % (1-logitLasso.score(xTest, yTest)))
In [ ]:
# Grilles de valeurs du paramètre de pénalisation
param=[{"C":[0.5,1,5,10,12,15,30]}]
logitR = GridSearchCV(LogisticRegression(penalty="l2"), param,cv=5,n_jobs=-1)
logitRidge=logitR.fit(xApp, yApp)
# Sélection du paramètre optimal
logitRidge.best_params_["C"]
print("Meilleur score = %f, Meilleur paramètre = %s" % (1.-logitRidge.best_score_,logitRidge.best_params_))
In [ ]:
# Prévision
yChap = logitRidge.predict(xTest)
# matrice de confusion
table=pd.crosstab(yChap,yTest)
print(table)
# Erreur sur l'échantillon test
print("Erreur de test régression Ridge = %f" % (1-logitRidge.score(xTest, yTest)))
Q Noter l'erreur de prévision; Comparer avec celle estimée par validation croisée.
L'objet LassoOpt issu de GridSearchCV ne retient pas les paramètres estimés dans le modèle. Il faut donc ré-estimer avec la valeur optimale du paramètre de pénalisation si l'on souhaite afficher ces coefficients.
In [ ]:
LassoOpt=LogisticRegression(penalty="l1",C=12)
LassoOpt=LassoOpt.fit(xApp, yApp)
# Récupération des coefficients
vect_coef=np.matrix.transpose(LassoOpt.coef_)
vect_coef=vect_coef.ravel()
#Affichage des 25 plus importants
coef=pd.Series(abs(vect_coef),index=xApp.columns).sort_values(ascending=False)
print(coef)
In [ ]:
plt.figure(figsize=(7,4))
coef.plot(kind='bar')
plt.title('Coeffients')
plt.tight_layout()
plt.show()
In [ ]:
from sklearn.metrics import roc_curve
listMethod=[["Lasso",logitLasso],["Ridge",logitRidge]]
for method in enumerate(listMethod):
probas_ = method[1][1].predict_proba(xTest)
fpr, tpr, thresholds = roc_curve(yTest, probas_[:,1])
plt.plot(fpr, tpr, lw=1,label="%s"%method[1][0])
plt.xlabel('Taux de faux positifs')
plt.ylabel('Taux de vrais positifs')
plt.legend(loc="best")
plt.show()
In [ ]:
from sklearn import discriminant_analysis
from sklearn.neighbors import KNeighborsClassifier
In [ ]:
lda = discriminant_analysis.LinearDiscriminantAnalysis()
disLin=lda.fit(xApp, yApp)
# Prévision de l'échantillon test
yChap = disLin.predict(xTest)
# matrice de confusion
table=pd.crosstab(yChap,yTest)
print(table)
# Erreur de prévision sur le test
print("Erreur de test lda = %f" % (1-disLin.score(xTest,yTest)))
In [ ]:
qda = discriminant_analysis.QuadraticDiscriminantAnalysis()
disQua=qda.fit(xApp, yApp)
In [ ]:
# Prévision de l'échantillon test
yChap = disQua.predict(xTest)
# matrice de confusion
table=pd.crosstab(yChap,yTest)
print(table)
# Erreur de prévision sur le test
print("Erreur de test qda = %f" % (1-disQua.score(xTest,yTest)))
In [ ]:
knn=KNeighborsClassifier(n_neighbors=10)
# Définition du modèle
disKnn=knn.fit(xApp, yApp)
# Prévision de l'échantillon test
yChap = disKnn.predict(xTest)
# matrice de confusion
table=pd.crosstab(yChap,yTest)
print(table)
# Erreur de prévision sur le test
print("Erreur de test knn = %f" % (1-disKnn.score(xTest,yTest)))
In [ ]:
yChap
In [ ]:
#Optimisation du paramètre de complexité k
#Grille de valeurs
param_grid=[{"n_neighbors":list(range(1,15))}]
disKnn=GridSearchCV(KNeighborsClassifier(),param_grid,cv=5,n_jobs=-1)
disKnnOpt=disKnn.fit(xApp, yApp) # GridSearchCV est lui même un estimateur
# paramètre optimal
disKnnOpt.best_params_["n_neighbors"]
print("Meilleur score = %f, Meilleur paramètre = %s" % (1.-disKnnOpt.best_score_,disKnnOpt.best_params_))
In [ ]:
# Prévision de l'échantillon test
yChap = disKnnOpt.predict(xTest)
# matrice de confusion
table=pd.crosstab(yChap,yTest)
print(table)
# Estimation de l'erreur de prévision sur l'échantillon test
print("Erreur de test knn_opt = %f" % (1-disKnnOpt.score(xTest,yTest)))
Courbes ROC
In [ ]:
from sklearn.metrics import roc_curve
# Liste des méthodes
listMethod=[["lda",disLin],["qda",disQua],["knn",disKnnOpt]]
# Tracé des courbes
for method in enumerate(listMethod):
probas_ = method[1][1].predict_proba(xTest)
fpr, tpr, thresholds = roc_curve(yTest, probas_[:,1])
plt.plot(fpr, tpr, lw=1,label="%s"%method[1][0])
plt.xlabel('Taux de faux positifs')
plt.ylabel('Taux de vrais positifs')
plt.legend(loc="best")
plt.show()
Les arbres binaires de décision concurrencent la régression logistique et gardent une place de choix dans les services de Gestion de la Relation Client, maintenant de Science des Données, par la facilité d'interprétation des modèles qui en découlent. L'optimisation de la complexité d'un artbre peut être délicate à opérer cr très sensible aux fluctuations de l'échantillon.
In [ ]:
from sklearn.tree import DecisionTreeClassifier
In [ ]:
# définition du modèle
tree= DecisionTreeClassifier()
treeC=tree.fit(xApp, yApp)
Q Quel est le critère d'homogénéité des noeuds utilisé par défaut?
Q Quel est le problème concernant l'élagage de l'arbre dans Scikkit-learn vis à vis des possibliités de la librairie rpart de R?
In [ ]:
# Optimisation de la profondeur de l'arbre
param=[{"max_depth":list(range(2,10))}]
tree= GridSearchCV(DecisionTreeClassifier(),param,cv=10,n_jobs=-1)
treeOpt=tree.fit(xApp, yApp)
# paramètre optimal
print("Meilleur score = %f, Meilleur paramètre = %s" % (1. - treeOpt.best_score_,treeOpt.best_params_))
In [ ]:
# Prévision de l'échantillon test
yChap = treeOpt.predict(xTest)
# matrice de confusion
table=pd.crosstab(yChap,yTest)
print(table)# Erreur de prévision sur le test
print("Erreur de test tree qualitatif = %f" % (1-treeOpt.score(xTest,yTest)))
In [ ]:
from sklearn.tree import export_graphviz
from sklearn.externals.six import StringIO
import pydotplus
treeG=DecisionTreeClassifier(max_depth=treeOpt.best_params_['max_depth'])
treeG.fit(xApp,yApp)
dot_data = StringIO()
export_graphviz(treeG, out_file=dot_data)
graph=pydotplus.graph_from_dot_data(dot_data.getvalue())
graph.write_png("treeOpt.png")
In [ ]:
from IPython.display import Image
Image(filename='treeOpt.png')
Comparaison des méthodes précédentes.
La valeur de seuil par défaut (0.5) n'étant pas nécessairement celle "optimale", il est important de comparer les courbes ROC.
In [ ]:
# Liste des méthodes
listMethod=[["Logit",logitLasso],["lda",disLin],["Arbre",treeOpt]]
# Tracé des courbes
for method in enumerate(listMethod):
probas_ = method[1][1].predict_proba(xTest)
fpr, tpr, thresholds = roc_curve(yTest, probas_[:,1])
plt.plot(fpr, tpr, lw=1,label="%s"%method[1][0])
plt.xlabel('Taux de faux positifs')
plt.ylabel('Taux de vrais positifs')
plt.legend(loc="best")
plt.show()
Commenter les résultats.
Q Intérêt de la régression logistique par rapport à l'analyse discriminante linéaire?
Q Conséquence du croisement des courbes ROC sur l'évaluation de l'AUC.
L'échantillon test reste de taille modeste (200). une étude plus systématique est nécessaire ainsi que la prise en compte des autres méthodes.
Il s'agit de comparer les principaux algorithmes issus de l'apprentissage machine: bagging, random forest, boosting.
Q Quel est par défaut l'estimateur qui est agrégé?
Q Quel est le nombre d'estimateurs par défaut? Est-il nécessaire de l'optimiser?
In [ ]:
from sklearn.ensemble import BaggingClassifier
bag= BaggingClassifier(n_estimators=100,oob_score=False)
bagC=bag.fit(xApp, yApp)
# Prévision de l'échantillon test
yChap = bagC.predict(xTest)
# matrice de confusion
table=pd.crosstab(yChap,yTest)
print(table)
# Erreur de prévision sur le test
print("Erreur de test avec le bagging = %f" % (1-bagC.score(xTest,yTest)))
In [ ]:
from sklearn.ensemble import RandomForestClassifier
In [ ]:
# Optimisation de max_features
param=[{"max_features":list(range(2,10,1))}]
rf= GridSearchCV(RandomForestClassifier(n_estimators=100),param,cv=5,n_jobs=-1)
rfOpt=rf.fit(xApp, yApp)
# paramètre optimal
print("Meilleur score = %f, Meilleur paramètre = %s" % (1. - rfOpt.best_score_,rfOpt.best_params_))
In [ ]:
# Prévision de l'échantillon test
yChap = rfOpt.predict(xTest)
# matrice de confusion
table=pd.crosstab(yChap,yTest)
print(table)
# Erreur de prévision sur le test
print("Erreur de test random forest opt -quantitatif = %f" % (1-rfOpt.score(xTest,yTest)))
In [ ]:
from sklearn.ensemble import GradientBoostingClassifier
# Optimisation de deux paramètres
paramGrid = [
{'n_estimators': list(range(100,601,50)), 'learning_rate': [0.1,0.2,0.3,0.4]}
]
gbmC= GridSearchCV(GradientBoostingClassifier(),paramGrid,cv=5,n_jobs=-1)
gbmOpt=gbmC.fit(xApp, yApp)
# paramètre optimal
print("Meilleur score = %f, Meilleur paramètre = %s" % (1. - gbmOpt.best_score_,gbmOpt.best_params_))
In [ ]:
# Prévision de l'échantillon test
yChap = gbmOpt.predict(xTest)
# matrice de confusion
table=pd.crosstab(yChap,yTest)
print(table)
# Erreur de prévision sur le test
print("Erreur de test gbm opt = %f" % (1-gbmOpt.score(xTest,yTest)))
In [ ]:
# Liste des méthodes
listMethod=[["Logit",logitLasso],["lda",disLin],["Arbre",treeOpt],["RF",rfOpt],["GBM",gbmOpt]]
# Tracé des courbes
for method in enumerate(listMethod):
probas_ = method[1][1].predict_proba(xTest)
fpr, tpr, thresholds = roc_curve(yTest, probas_[:,1])
plt.plot(fpr, tpr, lw=1,label="%s"%method[1][0])
plt.xlabel('Taux de faux positifs')
plt.ylabel('Taux de vrais positifs')
plt.legend(loc="best")
plt.show()
Q Quelles meilleure méthode interprétable? Quelle meilleure méthode?
Q Que dire de l'extrem gradient boosting ? Du nombre de paramètres à optimiser? De son implémentation en Python par rapport à R? De sa disponibilité sous Windows?
Exercice Ajouter les réseaux de neurones et les SVM dans la comparaison.
L'échantillon est de faible taille (#200), et les estimations des taux de bien classés comme le tracé des courbes ROC sont très dépendants de l’échantillon test; on peut s’interroger sur l’identité du modèle le plus performant ainsi que sur la significativité des différences observées entre les méthodes. Il est donc important d’itérer le processus (validation croisée Monte Carlo) sur plusieurs échantillons tests. Exécuter la fonction en annexe en choisissant les méthodes semblant les plus performantes.
In [ ]:
from sklearn.utils import check_random_state
import time
check_random_state(13)
tps0=time.clock()
# définition des estimateurs
logit = LogisticRegression(penalty="l1")
lda = discriminant_analysis.LinearDiscriminantAnalysis()
arbre = DecisionTreeClassifier()
rf = RandomForestClassifier(n_estimators=200)
gbm = GradientBoostingClassifier()
# Nombre d'itérations
B=3 # pour utiliser le programme, mettre plutôt B=30
# définition des grilles de paramètres
listMethGrid=[
[logit,{"C":[0.5,1,5,10,12,15,30]}],
[lda,{}],
[arbre,{"max_depth":[2,3,4,5,6,7,8,9,10]}],
[rf,{"max_features":[2,3,4,5,6]}],
[gbm,{"n_estimators": list(range(100,601,50)),"learning_rate": [0.1,0.2,0.3,0.4]}]
]
# Initialisation à 0 des erreurs pour chaque méthode (colonne) et chaque itération (ligne)
arrayErreur=np.empty((B,5))
for i in range(B): # itérations sur B échantillons test
# extraction apprentissage et test
xApp,xTest,yApp,yTest=train_test_split(vispremR,y,test_size=200)
# optimisation de chaque méthode et calcul de l'erreur sur le test
for j,(method, grid_list) in enumerate(listMethGrid):
methodGrid=GridSearchCV(method,grid_list,cv=5,n_jobs=-1).fit(xApp, yApp)
methodOpt = methodGrid.best_estimator_
methFit=methodOpt.fit(xApp, yApp)
arrayErreur[i,j]=1-methFit.score(xTest,yTest)
tps1=time.clock()
print("Temps execution en mn :",(tps1 - tps0)/60)
In [ ]:
dataframeErreur=pd.DataFrame(arrayErreur,columns=["Logit","LDA","Arbre","RF","GBM"])
In [ ]:
# Distribution des erreurs
dataframeErreur[["Logit","LDA","Arbre","RF","GBM"]].boxplot(return_type='dict')
plt.show()
Q Finalement, quelle meilleure méthode? Quelle meilleure méthode interprétable?
Exercice: tester XGBoost.
In [ ]:
# Moyennes
dataframeErreur.mean()