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, et la comparaion des résultats obtenus avec R. Il propose également de tester l'analyse discriminante pas reprise dans ce calepin car les résultats sont en retrait vis-à-vis des autres méthodes.
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 [ ]:
visaData=read.table("vispremv.dat")
Q Combien d'individus et combien de variables sont concernées?
Vérifier ci-dessous que la plupart des variables ont deux versions, l'une quantitative et l'autre qualitative.
In [ ]:
var=names(visaData)
var
Deux bases sont constituées: une contenant les variables explicatives initiales: qualitatives (csp, sexe, famille) et quantitatives , une autre celles qualitatives obtenues après recodage, les deux intègrent la variable à modéliser CARVP.
In [ ]:
varquant=var[c(1:3,26:54)]
varqual=var[c(1:25,54)]
visaDataR=visaData[,varquant]
visaDataQ=visaData[,varqual]
Comme annoncé ci-dessus, la plupart des variables apparaissent deux fois: sous forme quantitative ou qualitative. Ainsi, la variable AGER est quantitative en années tandis que ageq est qualitative, de type facteur à trois niveaux. Ceci permet de s’adapter facilement à certaines méthodes qui n’acceptent que des variables explicatives quan-
titatives. Néanmoins, le choix de découper en classe une variable quantitative peut avoir un impact important sur la qualité de la prévision de certaines méthodes:
Cette possibilité impacte une méthode comme la régression logistique mais peut nuire à la construction d'un arbre qui recherche par construction des segmentations "optimales" des ensembles de valeurs.
Il reste à comparer les stratégies pour retenir celle conduisant à la meilleure prévision.
In [ ]:
xxx=111 # modifer cette valeur
set.seed(xxx)
# modifier 111
npop=nrow(visaData)
# tirage de 200 indices sans remise
testi=sample(1:npop,200)
#Liste des indices restant qui n’ont pas été tirés
appri=setdiff(1:npop,testi)
# Extraction échantillons d’apprentissage
visApptQ=visaDataQ[appri,]
visApptR=visaDataR[appri,]
# Extraction échantillons de test
visTestQ=visaDataQ[testi,]
visTestR=visaDataR[testi,]
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 chez Criteo ou CDiscount. Elle est étudiée de façon plus détaillée pour comparer les deux stratégies: sur variables quantitatives ou qualitatives.
Pluasieurs algorithmes de sélection de variables peuvent être testés: ceux forward, backward, step-wise en minimisant le critère AIC ou encore les algorithmes avec pénalisation ridge Lasso ou elastic net.
Ces derniers algorihtmes de sélection par pénalisation sont testés de manière plus approfondie dans la version Python de ce scénario. Seule la sélection step-wise est testée pour simplifier les comparaisons. Les interactions ne sont pas prises en compte. Certains services ont pour pratique de mixer arbre de classification et régression logistique pour sélectionner et intégrer les interactions les plus évidentes: Les premières segmentations d'un arbre de décision déterminent une nouvelle variable qui est intégrée à la régression logistique.
In [ ]:
varquant
In [ ]:
# modèle trivial de départ avec le seul terme constant
visaR.logit=glm(CARVP~1,data=visApptR,family=binomial,na.action=na.omit)
# sélection step-wise
# le paramètre trace=0 peut être enlevé pour suivre la progression
visaR.step=step(visaR.logit,direction="both",scope=list(lower=~1, upper=~SEXEQ+FAMIQ+PCSPQ+RELAT+AGER+OPGNBL+
MOYRVL+TAVEPL+ENDETL+GAGETL+GAGECL+GAGEML+KVUNB+QSMOY+QCREDL+DMVTPL+BOPPNL+FACANL+LGAGTL+
VIENB+VIEMTL+UEMNB+XLGNB+XLGMTL+YLVNB+YLVMTL+ROCNB+NPTAG+ITAVCL+HAVEFL+JNBJDL), trace=0)
# observer les p-valeurs
anova(visaR.step,test="Chisq")
Estimation de l'erreur par validation croisée puis sur l'échantillon test.
In [ ]:
cv.glm(visApptR,visaR.step,K=10)$delta[1]
In [ ]:
predTestR=predict(visaR.step,newdata=visTestR)>0.5
table(predTestR,visTestR$CARVP=="Coui")
In [ ]:
varqual
In [ ]:
# modèle trivial de départ avec le seul terme constant
visaQ.logit=glm(CARVP~1,data=visApptQ,family=binomial,na.action=na.omit)
# sélection step-wise
# le paramètre trace=0 peut être enlevé pour suivre la progression
visaQ.step=step(visaQ.logit,direction="both",scope=list(lower=~1, upper=~SEXEQ +FAMIQ +PCSPQ + kvunbq + vienbq+uemnbq +
xlgnbq+ylvnbq+rocnbq+nptagq + endetq + gagetq + facanq +lgagtq+ havefq + ageq+relatq + qsmoyq + opgnbq +
moyrvq + tavepq + dmvtpq + boppnq + jnbjdq + itavcq), trace=0)
# observer les p-valeurs
anova(visaQ.step,test="Chisq")
Q Quelles sont les variables importantes? Comment interpréter?
Estimation de l'erreur de prévision par validation croisée puis sur l'échantillon test.
In [ ]:
cv.glm(visApptQ,visaQ.step,K=10)$delta[1]
In [ ]:
predTestQ=predict(visaQ.step,newdata=visTestQ)>0.5
table(predTestQ,visTestQ$CARVP=="Coui")
Conclusions (qui dépendent du choix de l'initialisation de xxx):
Q quel modèle est le moins complexe (mois de paramètre)?
Q Quel modèle est le plus simple à interpréter?
Q Quel modèle prédit le mieux?
Les arbres binaires de décision concurrencent la régression logistique et gardent une pace de choix dans le sservices 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.
In [ ]:
library(rpart)
visaR.tree=rpart(CARVP~.,data=visApptR,parms=list(split='information'),cp=0.001)
# Choix élémentaire du coefficient de complexité
options(repr.plot.width=8, repr.plot.height=4)
plotcp(visaR.tree)
Q Comment lire le graphique précédent?
In [ ]:
# Tracé de l’arbre élagué
visaR.treeOpt=prune(visaR.tree,cp=0.018)
library(partykit)
options(repr.plot.width=8, repr.plot.height=6)
plot(as.party(visaR.treeOpt), type="simple")
In [ ]:
visaQ.tree=rpart(CARVP~.,data=visApptQ,parms=list(split='information'),cp=0.001)
# Choix élémentaire du coefficient de complexité
options(repr.plot.width=8, repr.plot.height=4)
plotcp(visaQ.tree)
In [ ]:
# Tracé de l’arbre élagué
visaQ.treeOpt=prune(visaQ.tree,cp=0.012)
library(partykit)
options(repr.plot.width=8, repr.plot.height=6)
plot(as.party(visaQ.treeOpt), type="simple")
Une procédure d'élagage par validation croisée serait sans doute plus précise. Elle est intégrée ci-après dan sla comparaison automatique utilisant le package caret.
Erreurs de prévision sur le test.
In [ ]:
predTestR=predict(visaR.treeOpt,newdata=visTestR,type="class")
table(predTestR,visTestR$CARVP=="Coui")
In [ ]:
predTestQ=predict(visaQ.treeOpt,newdata=visTestQ,type="class")
table(predTestQ,visTestQ$CARVP=="Coui")
Comparer les résultats
Q Quel est le modèle le plus simple / facile à interpréter
Q Quel est celui fournissant les meilleures prévisions.
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 [ ]:
library(ROCR)
ROCdislogR=predict(visaR.step,newdata=visTestR)
preddislogR=prediction(ROCdislogR,visTestR$CARVP=="Coui")
perfdislogR=performance(preddislogR,"tpr","fpr")
ROCdislogQ=predict(visaQ.step,newdata=visTestQ)
preddislogQ=prediction(ROCdislogQ,visTestQ$CARVP=="Coui")
perfdislogQ=performance(preddislogQ,"tpr","fpr")
ROCdistreeR=predict(visaR.treeOpt,newdata=visTestR,type="prob")[,2]
preddistreeR=prediction(ROCdistreeR,visTestR$CARVP=="Coui")
perfdistreeR=performance(preddistreeR,"tpr","fpr")
ROCdistreeQ=predict(visaQ.treeOpt,newdata=visTestQ,type="prob")[,2]
preddistreeQ=prediction(ROCdistreeQ,visTestQ$CARVP=="Coui")
perfdistreeQ=performance(preddistreeQ,"tpr","fpr")
# tracer les courbes ROC en les superposant
# pour mieux comparer
plot(perfdislogR,col=1)
plot(perfdislogQ,col=2,add=TRUE)
plot(perfdistreeR,col=3,add=TRUE)
plot(perfdistreeQ,col=4,add=TRUE)
legend("bottomright",legend=c("LogitR","LogitQ","TreeR","TreeQ"),col=c(1:4),pch="_")
Commenter les résultats.
Q Intérêt de la régression logistique par rapport à un arbre.
Q Conséquence du croisemen des courbes ROC sur l'évaluation 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;
caretUn avantage de R est le nombre considérables d'utilisateurs qui participent au développement des librairies. cet avantage a un revers: le manque d'homogénéité de celles-ci. Pour y remédier dans les applications d'apprentissage machine, la (méta)librairie caret de Max Kuhn (2008) intègre dans un même usage, une même syntaxe, l'ensemble des fonctionnalités d'apprentissage et propose une approche unifiée des procédures d'optimisation des paramètres.
De très nombreuses stratégies sont possibles. Il suffit de les mettre en oeuvre afin de comparer les qualités prédictives: considérer les seules variables initiales, celles transformées, la réunion des deux ensembles, les composantes issues de l'AFCM de l'étape précédente ... à croiser avec plu sde 200 options et méthodes d'apprentissage de caret, considérer aussi caretEnsemble qui permet de combiner des modèles entre eux comme cela est courant pour les concours kaggle.
Une sélection drastique est opérée en ne considérant que quelques méthodes parmi les grandes familles et sur les seules variables initiales.
Par ailleurs, même sous windows, caret offre simplement des possibilités de parallèlisation en utilisant la package doParallel. Même si les algorithmes des différentes méthodes d'apprentissage ne sont pas parallélisés, les itérations des calculs de validations croiser pour l'optimisation des paramètres sont effectivement parallélisés avec un gain de temps très appréciable fonciton du nombre de processeurs. Ceci est obtenu en exécutant les commandes suivantes en supposant que 4 processeurs sont disponibles.
In [ ]:
library(doParallel)
cl <- makeCluster(4)
registerDoParallel(cl)
In [ ]:
library(caret)
# extraction des données
# Variable cible
Y=visaData[,"CARVP"]
# Variables explicatives
X=visaData[,varquant[-32]]
# Transformation des facteurs en indicatrices pour utiliser certains algorithmes
# notamment xgboost
library(FactoMineR)
X=data.frame(tab.disjonctif(X[,c(1:3)]),X[,-c(1:3)])
summary(Y);summary(X)
Préparation des échantillons
In [ ]:
xx=11 # Changer cette valeur pour personnaliser l'échantillonnage
set.seed(xx)
inTrain = createDataPartition(X[,1],p = 0.8, list = FALSE)
# Extraction des échantillons
trainDescr=X[inTrain,]
testDescr=X[-inTrain,]
testY=Y[-inTrain]
trainY=Y[inTrain]
Certaines méthodes sont sensibles à des effets de variance ou d'unité des variables. Il est préférable d'introduire une normalisation.
In [ ]:
# Normalisation calculée sur les paramètres de l'échantillon d'apprentissage
xTrans=preProcess(trainDescr)
trainDescr=predict(xTrans,trainDescr)
# Puis appliquée également à l'échantillon test
testDescr=predict(xTrans,testDescr)
# Choix de la validation croisée
cvControl=trainControl(method="cv",number=10)
Q Pour chacun des modèles ou méthodes utilisés ci-dessous, identifier la méthode, préciser les paramètres associés et noter celui ou ceux optimisés par défaut par caret.
In [ ]:
#1 Régression logistique
# Attention, la régression logistique sans interaction (linéaire) est estimée ci-dessous
set.seed(2)
rlogFit = train(trainDescr, trainY,method = "glmStepAIC", trace=FALSE)
rlogFit
In [ ]:
#2 Arbre de décision
set.seed(2)
rpartFit = train(trainDescr, trainY, method = "rpart", tuneLength = 10,
trControl = cvControl)
rpartFit
plot(rpartFit)
In [ ]:
#3 Réseau de neurones
set.seed(2)
nnetFit = train(trainDescr, trainY, method = "nnet", tuneLength = 6,
trControl = cvControl, trace=FALSE)
nnetFit
plot(nnetFit)
In [ ]:
#4 Random forest
set.seed(2)
rfFit = train(trainDescr, trainY,method = "rf", tuneLength = 8,
trControl = cvControl, trace=FALSE)
rfFit
plot(rfFit)
In [ ]:
#5 Boosting
set.seed(2)
gbmFit = train(trainDescr, trainY,method = "gbm", tuneLength = 8,
trControl = cvControl)
gbmFit
plot(gbmFit)
In [ ]:
#6 Extrême boosting
set.seed(2)
xgbFit = train(trainDescr, trainY,method = "xgbTree", tuneLength = 6,
trControl = cvControl, trace=FALSE)
xgbFit
plot(xgbFit)
In [ ]:
#7 SVM
set.seed(2)
svmFit = train(trainDescr, trainY,method = "svmRadial", tuneLength = 6,
trControl = cvControl, trace=FALSE)
svmFit
plot(svmFit)
In [ ]:
models=list(logit=rlogFit,cart=rpartFit,nnet=nnetFit,rf=rfFit,gbm=gbmFit,xgb=xgbFit,svm=svmFit)
testPred=predict(models, newdata = testDescr)
# taux de bien classés
lapply(testPred,function(x)mean(x==testY))
Tracer les courbes ROC pour analyser spécificité et sensibilité des différentes méthodes.
In [ ]:
library(ROCR)
models=list(logit=rlogFit,cart=rpartFit,nnet=nnetFit,rf=rfFit,gbm=gbmFit,xgb=xgbFit)
testProb=predict(models, newdata = testDescr,type="prob")
predroc=lapply(testProb,function(x)prediction(x[,1],testY=="Cnon"))
perfroc=lapply(predroc,
function(x)performance(x, "tpr", "fpr"))
plot(perfroc$logit,col=1)
plot(perfroc$cart,col=2,add=TRUE)
plot(perfroc$nnet,col=3,add=TRUE)
plot(perfroc$rf,col=4,add=TRUE)
plot(perfroc$gbm,col=5,add=TRUE)
plot(perfroc$xgb,col=6,add=TRUE)
legend("bottomright",legend=c("logit","CART","nnet","RF","boost","xgBoost"),col=c(1:6),pch="_")
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 [ ]:
# Choisir la liste des méthodes et l’effort d’optimisation
models=c("glmStepAIC","rpart","nnet","rf","gbm","xgbTree")
noptim=c(6,6,6,6,6,6)
# Initialiser le générateur et fixer le nombre d’itérations
# Changer ces valeurs. Attention au temps de calcul! Être patient!
Niter=3 ; Init=11
# Appel de la fonction définie en annexe
pred.ozone=pred.autom(X,Y,methodes=models,N=Niter,xinit=Init,size=noptim,type="prob")
Puis calculer et représenter les erreurs pour les méthodes considérées
In [ ]:
# Calcul des taux de bien classés
obs=pred.ozone$obs
prev.ozone=pred.ozone$pred
res.ozone=lapply(prev.ozone,function(x)apply((x>0.5)==(obs==1),2,mean))
# Moyennes des taux de bien classés par méthode
lapply(res.ozone,mean)
# distributions des taux de bien classés
boxplot(data.frame(res.ozone))
Les commandes suivandes tracent les courbes ROC moyennes.
In [ ]:
## Comparaison des méthodes par le
# tracer des courbes ROC moyennes
# Problème pas identifié avec rlogit!
predroc.ozone=lapply(prev.ozone,function(x)prediction(x,obs==1))
perfroc.ozone=lapply(predroc.ozone,function(x)performance(x,"tpr","fpr"))
plot(perfroc.ozone$gbm,col=1,lwd=2,avg="vertical")
plot(perfroc.ozone$rf,col=2,add=TRUE,lwd=1.5,avg="vertical")
plot(perfroc.ozone$nnet,add=TRUE,col=3,lwd=1.5,avg="vertical")
plot(perfroc.ozone$xgbTree,add=TRUE,col=4,lwd=2,avg="vertical")
plot(perfroc.ozone$glmStepAIC,add=TRUE,col=5,lwd=1.5,avg="vertical")
legend("bottomright",legend=c("boost","RF", "nnet","xgBoost","logit"),col=c(1:5),pch="_")
Q Quelle méthode retenir, en fonction du taux de faux positif acceptable, pour prévoir le dépassement du seuil? Et si le comanditaire veut une solution explicable?
N.B. L'algorithme xgboost nécessiterait des efforts plus important d'optimisation des paramètres.
N réplications des estimations / prévisions
In [ ]:
pred.autom=function(X,Y,p=1/2,methodes=c("knn",
"rf"),size=c(10,2),xinit=11,N=10,typerr="cv",
number=4,type="raw") {
# Fonction de prévision de N échantillons tests
# par une liste de méthodes de régression
# ou classification (uniquement 2 classes)
# Optimisation des paramètres par validation
# croisée (défaut) ou bootstrap ou... (cf. caret)
# X : matrice ou frame des variables explicatives
# Y : variable cible quantitative ou qualitative
# p : proportion entre apprentissage et test
# methodes : liste des méthodes de rdiscrimination
# size : e grille des paramètres à optimiser
# xinit : générateur de nombres aléatoires
# N : nombre de réplications apprentissage/test
# typerr : "cv" ou "boo" ou "oob"
# number : nombre de répétitions CV ou bootstrap
# pred : liste des matrices de prévision
# type d’erreur
Control=trainControl(method=typerr,number=number)
# initialisation du générateur
set.seed(xinit)
# liste de matrices stockant les prévisions
# une par méthode
inTrain=createDataPartition(Y,p=p,list=FALSE)
ntest=length(Y[-inTrain])
pred=vector("list",length(methodes))
names(pred)=methodes
pred=lapply(pred,function(x)x=matrix(0,
nrow=ntest,ncol=N))
obs=matrix(0,ntest,N)
set.seed(xinit)
for(i in 1:N) {
# N itérations
# indices de l’échantillon d’apprentissage
inTrain=createDataPartition(Y,p=p,list=FALSE)
# Extraction des échantillons
trainDescr=X[inTrain,]
testDescr=X[-inTrain,]
trainY=Y[inTrain]
testY=Y[-inTrain]
# stockage des observés de testY
obs[,i]=testY
# centrage et réduction des variables
xTrans=preProcess(trainDescr)
trainDescr=predict(xTrans,trainDescr)
testDescr=predict(xTrans,testDescr)
# estimation et optimisation des modèles
# pour chaque méthode de la liste
for(j in 1:length(methodes)) {
# modélisation
modFit = train(trainDescr, trainY,method = methodes[j], tuneLength = size[j],
trControl = Control)
# prévisions
if (type=="prob") pred[[j]][,i]=predict(modFit,
newdata = testDescr,type=type)[,1]
else pred[[j]][,i]=predict(modFit,
newdata = testDescr)
}}
list(pred=pred,obs=obs)
# résultats
}