Après l'exploration et la classification des contenus d'une base de messages électroniques afin de catégoriser les pourriels, il s'agit de construire des scores ou modèles prévoyant la nature d'un message: pourriel ou non, en fonction de son contenu ou plutôt de la présence ou fréquence de certains mots et caractères. Les principales méthodes de modélisation et apprentissage statistique sont testées et comparées, la librairie caret est utilisée pour automatiser les traitements.
Cette étude est un exemple d'analyse textuelle d'un corpus de documents, ici des courriels. Une telle analyse est classiquement basée sur la fréquence d'une sélection de mots. Après l'analyse exploratoire, l'objectif général est de pouvoir discriminer les courriels pertinents ou encore de définir un modèle personnalisé de détection des pourriels (spams), c'est-à-dire adapté au contenu spécifique de la boîte d'un internaute. Il s'agit donc d'un modèle susceptible de prévoir la qualité d'un message reçu en fonction de son contenu. Le déroulement de cette étude est évidemment marqué par le type particulier des données mais celle-ci peut facilement se transposer à d'autres types de données textuelles ou analyse du contenu de livres, pages web, discours politiques, réponses ouvertes à des questionnaires... les exemples sont nombreux en sciences humaines, marketing lorsqu'il s'agit d'estimer / prévoir des scores, par exemple, de satisfaction de clientèle.
L'objectif du présent travail est d'évaluer, comparer différentes stratégies et méthodes pour aboutir à celle la plus efficace de filtrage des courriers indésirables. Il s'agit de déterminer quelle méthode de modélisation parmi celles concurrentes: logit, arbre, réseau de neurones, random forest, boosting, SVM, se montre la plus efficace.
George, ingénieur chez HP dans le département Computer Science a recueilli un échantillon de messages électroniques dans chacun desquels il a évalué le nombre d'occurrences d'une sélection de mots et caractères. Les variables considérées sont, dans une première partie, des rapports: nombre d'occurrences d'un mot spécifique sur le nombre total de mots ou nombre d'occurrences d'un caractère sur le nombre de caractères du message avant d'être, dans une deuxième partie, des indicatrices ou facteurs: présence / absence de mots ou ensemble de caractères. Il a également considéré trois variables prenant en compte la casse (majuscule / minuscule) des caractères et une dernière qualitative binaire indiquant le classement qu'il a fait de chaque message : spam ou Nsp. Les variables d'occurrences sont décrites dans le tableau 1, celles associées à la casse dans le tableau 2. Ces données sont publiques, elles servent régulièrement de benchmark pour la comparaison de méthodes d'apprentissage machine:
Frank A., Asuncion A. (2010). UCI Machine Learning Repository. Irvine, CA: University of California, School of Information and Computer Science.
Ce sont donc finalement 58 variables qui sont observées sur 4601 messages dont 1813 pourriels (spams). La variable binaire Spam est présente à titre illustratif, elle est toujours considérée en supplémentaire dans ce travail exploratoire préliminaire.
Le tableau ci-dessous liste 54 variables exprimant soit:
Tableau 1: Les colonnes contiennent successivement le libellé de la variable, le mot ou ensemble de caractères concernés, le libellé des modalités Présence / Absence utilisées après recodage.
| Variable | Mot ou Carac | Modalités P/A | Variable | Mot ou Carac. | Modalités |
|---|---|---|---|---|---|
| make | make | make / Nmk | X650 | 650 | 650 / N65 |
| address | address | addr / Nad | lab | lab | lab / Nlb |
| all | all | all / Nal | labs | labs | labs / Nls |
| X3d | 3d | 3d / N3d | telnet | telnet | teln / Ntl |
| our | our | our / Nou | X857 | 857 | 857 / N87 |
| over | over | over / Nov | data | data | data / Nda |
| remove | remove | remo / Nrm | X415 | 415 | 415 / N41 |
| internet | internet | inte / Nin | X85 | 85 | 85 / N85 |
| order | order | orde / Nor | technology | technology | tech / Ntc |
| mail / Nma | X1999 | 1999 | 1999/ N19 | ||
| receive | receive | rece / Nrc | parts | parts | part / Npr |
| will | will | will / Nwi | pm | pm | pm / Npm |
| people | people | peop / Npp | direct | direct | dire / Ndr |
| report | report | repo / Nrp | cs | cs | cs / Ncs |
| addresses | addresses | adds / Nas | meeting | meeting | meet/Nmt |
| free | free | free / Nfr | original | original | orig / or |
| business | business | busi / Nbs | project | project | proj / Npj |
| emai / Nem | re | re | re / Nre | ||
| you | you | you / Nyo | edu | edu | edu / Ned |
| credit | credit | cred / Ncr | table | table | tabl / Ntb |
| your | your | your / Nyr | conference | conferenc | e conf / Ncf |
| font | order | font / Nft | CsemiCol | ; | Cscl / NCs |
| X000 | 000 | 000 / N00 | Cpar | ( | Cpar / NCp |
| money | money | mone/ Nmn | Ccroch | [ | Ccro / NCc |
| hp | hp | hp / Nhp | Cexclam | ! | Cexc / NCe |
| hpl | hpl | hpl / Nhl | Cdollar | \$ | Cdol / NCd |
| george | george | geor / Nge | Cdiese | # | Cdie / NCi |
Un deuxième tableau liste 4 variables dont celle dénombrant le nombre de lettres majuscules.
Tableau 2:Liste de 4 variables, de leur libellé et des modalités après recodage.
| Code variable | Libellé | Modalités |
|---|---|---|
| Spam | Type de message pourriel ou non | Spam / Nsp |
| CapLM | Nombre moyen de capitales par mot | Mm1 / Mm2 / Mm3 |
| CapLsup | Nombre max de capitales par mot | Ms1 / Ms2 / Ms3 |
| CapLtot | Nombre totale de lettres capitales | Mt1 / Mt2 / Mt3 |
Cette étude nécessite une série d'étapes afin de mettre en oeuvre les comparaisons recherchées.
XGBoost).
In [ ]:
spamq=read.table("spamq.dat")
summary(spamq)
In [ ]:
spamr=read.table("spam.dat",header=T)
summary(spamr)
Extractions des échantillons d'apprentissage (spamq.appt) et de test (spamq.test) dans le cas qualitatif, spamr.appt et spamr.test dans le cas quantitatif.
In [ ]:
# choisir une valeur entière pour "xx"
xx=...
set.seed(xx)
npop=nrow(spamq)
test=sample(1:npop,1000)
spamq.test=spamq[test,]
appt=setdiff(1:npop,test)
spamq.appt=spamq[appt,]
spamr.test=spamr[test,]
spamr.appt=spamr[appt,]
Les échantillons tests sont identiques afin de pouvoir mieux comparer les résultats.
La taille de l'échantillon initial relativement importante permettrait d'extraire un autre échantillon de validation surlequel serait minimiser les erreurs de prévision afin d'optimiser la complexité des modèles. Cette stratégie mise en oeuvre par défaut dans SAS Enterprise Miner s'avère dans R finalement plus compliquée à exécuter, que celle utilisant systématiquement la validation croisée, sans bénéfice notable sur la qualité des optimisations (cf. le tutoriel). C'est donc la validation croisée qui est privilégiée dans la recherche de modèle optimaux.
Plusieurs approches peuvent être comparées selon le type de données ou de transformation des données pris en compte:
Q Quel intérêt pour la régression logistique de considérer des variables qualitatives obtenues par découpage en classe de variables quantitatives?
Q Corrélativement, pourquoi les méthodes basées sur des arbres de décision ont pus intérêt à prendre en compte les variables quantitatives initiales.
La régression logistique est détaillée en utilisant deux algorithmes de sélections de variables basés sur une pénalisation AIC (backward et stepwise).
Ces algorihtmes de sélection seraient à comparer avec ceux utilisant une pénalisation Lasso qui sont systématiquement utilisés dans la librairie Scikit-learn de Python.
Les modèles trop complexes génèrent de très nombreux warnings dus à des estimations de probabilités 0 ou 1 pour certaines observations.
les warnings sont temporairement supprimés:
In [ ]:
oldw <- getOption("warn")
options(warn = -1)
In [ ]:
spamq.logit1=glm(spamf~.,data=spamq.appt,family=binomial)
spamq.log1=step(spamq.logit1,direction='backward', trace=0)
spamq.logit2=glm(spamf~1,data=spamq.appt,family=binomial)
spamq.log2=step(spamq.logit2,direction='both',
scope=list(lower=~1,upper=~make + address +
all + X3d + our + over + remove + internet +
order + mail + receive + will + people +
report + addresses + free + business +
email + you + credit + your + font + X000 +
money + hp + hpl + george + X650 + lab +
labs + telnet + X857 + data + X415 + X85 +
technology + X1999 + parts + pm + direct +
cs + meeting + original + project + re +
edu + table + conference + CsemiCol +
Cpar + Ccroch + Cexclam + Cdollar +
Cdiese + CapLMq + CapLsupq + CapLtotq), trace=0)
L'option trace=0 supprime les sorties très volumineuses lorsque le nombre de variables est important.
Comparer les performances des deux modèles par validation croisée.
In [ ]:
library(boot)
cv.glm(spamq.appt,spamq.log1,K=10)$delta[1]
cv.glm(spamq.appt,spamq.log2,K=10)$delta[1]
Retenir le meilleur modèle disons spamq.log2 mais cela dépend de l'échantillon d'apprentissage tiré.
La commande:
In [ ]:
anova(spamq.log2,test="Chisq")
permet d'identifier la variable la moins significative du modèle. Il serait alors possible de retirer cette variable du modèle, le ré-estimer avant de vérifier si l'erreur de prévision sur l'échantillon de validation est améliorée ou pas, ... itérer ou non cette démarche.
Même chose avec les variables quantitatives.
In [ ]:
spamr.logit1=glm(spam~.,data=spamr.appt,family=binomial)
spamr.log1=step(spamr.logit1,direction='backward', trace=0)
spamr.logit2=glm(spam~1,data=spamr.appt,family=binomial)
spamr.log2=step(spamr.logit2,direction='both',
scope=list(lower=~1,upper=~make + address +
all + X3d + our + over + remove + internet +
order + mail + receive + will + people +
report + addresses + free + business +
email + you + credit + your + font + X000 +
money + hp + hpl + george + X650 + lab +
labs + telnet + X857 + data + X415 + X85 +
technology + X1999 + parts + pm + direct +
cs + meeting + original + project + re +
edu + table + conference + CsemiCol +
Cpar + Ccroch + Cexclam + Cdollar +
Cdiese + CapLM + CapLsup + CapLtot),trace=0)
Comparer les modèles.
In [ ]:
library(boot)
cv.glm(spamr.appt,spamr.log1,K=10)$delta[1]
cv.glm(spamr.appt,spamr.log2,K=10)$delta[1]
In [ ]:
# retour des warnings
options(warn = oldw)
Optimiser éventuellement le choix opéré en tâchant de réduire l'erreur par validation croisée par le retrait des variables les moins significatives du modèle.
Appliquer le meilleur modèle obtenu à l'achantillon test.
In [ ]:
pred.log=predict(spamq.log2,newdata=spamq.test)
table(pred.log>0.5,spamq.test[,"spamf"]) # 0.065
In [ ]:
pred.log=predict(spamr.log2,newdata=spamr.test)
table(pred.log>0.5,spamr.test[,"spam"]) # 0.065
Q Quelle ensemble de variables retenir: qualitatives ou quantitatives, pour la régression logistique?
La suite du travail se propose d'utiliser la librairie caret(Kuhn, 2008) afin de faciliter et industrialiser les traitements. Attention, cette librairie fait appel à de nombreuses autres librairies qui doivent être installées à la demande. La version "artisanale" des commandes, sans l'utilisation de caret est proposée en annexe.
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 validation croisée pour l'optimisation des paramètres le sont avec un gain de temps très appréciable fonction 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 [ ]:
require(caret)
# extraction des données
# remplacer les niveaux "0" et "1"
# par "nospam" et "spam"
Y=rep("spam",nrow(spamr))
Y[spamr[,1]==0]="nospam"
Y=as.factor(Y)
summary(Y)
In [ ]:
X=spamr[,-1]
summary(X)
In [ ]:
# indices de l'échantillon d'apprentissage
set.seed(xx)
inTrain = createDataPartition(X[,1],p = 78/100, list = FALSE)
# Extraction des échantillons
trainDescr=X[inTrain,]
testDescr=X[-inTrain,]
trainY=Y[inTrain]
testY=Y[-inTrain]
In [ ]:
summary(trainY)
Il est recommandé de centrer et réduire les variables dans plusieurs méthodes. C'est fait systématiquement et simplement en utilisant évidemment les mêmes transformations sur l'échantillon test que celles mises en place sur l'échantillon d'apprentissage.
In [ ]:
# Normalisation
xTrans=preProcess(trainDescr)
trainDescr=predict(xTrans,trainDescr)
testDescr=predict(xTrans,testDescr)
# Choix de la validation croisée
cvControl=trainControl(method="cv",number=10)
Il peut arriver sur certains jeux de données qu'après partitionnement de l'échantillon par validation croisée ou bootstrap, certaines variables deviennent strictement constantes sur une classe. Ceci pose des problèmes à certaines méthodes comme l'analyse discriminante ou les classifieurs bayésien naïf car la variance est nulle pour ces variables au sein d'une classe.
La librairie caret inclut une fonction permettant d'itentifier ces variables en vue de leur suppression.
In [ ]:
# Consulter la documentation
?nearZeroVar
nearZeroVar(spamr, saveMetrics= TRUE)
Le problème de ces données textuelles très creuses : beaucoup de mots ne sont pas présents dans la plupart des messages, est que la grande majorité des variables est considérée comme near zero variance. Le seuil de décision peut être remonté mais cela n'élimine alors plus le risque de variance constante!
La librairie intègre beaucoup plus de méthodes et celles sélectionnées ci-dessous semblent les plus pertinentes. Certaines méthodes ne sont pas utilisables à cause du nombre de variables. Il pourrait être utile et intéressant d'en tester encore d'autres. Exécuter successivement les blocs de commandes pour tracer séparamment chacun des graphes afin de contrôler le bon comportement de l'optimisation du paramètre de complexité de chaque modèle.
Consulter la liste des méthodes disponibles dans l'aide de la fonction : ?train et tester les principales. L'utilisation de certaines comme la régression logistique est moins flexible qu'en utilisation "manuelle" en particulier dans le choix de l'algorithme de sélection de variables. Il faut se montrer (très) patient pour certaines optimisations alors que d'autres sont immédiates, voire inutiles.
L'expérimenter fait partie de l'expertise à acquérir.
Q Que signifie: trControl = cvControl dans lacommande ci-dessous?
In [ ]:
#1 Régression logistique
set.seed(2)
rlogFit = train(trainDescr, trainY,method = "glmStepAIC",trControl = cvControl)
rlogFit
In [ ]:
#2 analyse discriminante linéaire
set.seed(2)
ldaFit = train(trainDescr, trainY,method = "lda2",trControl = cvControl)
ldaFit
Q Que signifie le paramètre tuneLength pour les méthodes suivantes? pourquoi n'est-il pas présent avec la méthode glmStepAIC? Quel est son impact?
In [ ]:
#3 K plus proches voisins
set.seed(2)
knnFit = train(trainDescr, trainY,method = "knn", tuneLength = 10,trControl = cvControl)
knnFit
plot(knnFit)
Q Comment interpréter le graphe précédent?
Q Quel est le paramètre optimisé pour la méthode suivante?
In [ ]:
#4 Arbre de décision
set.seed(2)
rpartFit = train(trainDescr, trainY, method = "rpart", tuneLength = 10,trControl = cvControl)
rpartFit
plot(rpartFit)
Q Quel est le modèle ci-desosus ? Quels sont les paramètres optimisés?
In [ ]:
#5 Réseau de neurones
set.seed(2)
nnetFit = train(trainDescr, trainY,method = "nnet", tuneLength = 6,trControl = cvControl)
nnetFit
plot(nnetFit)
Q Quel est le paramètre optimisé de random forest?
In [ ]:
#6 Random forest
set.seed(2)
rfFit = train(trainDescr, trainY,method = "rf", tuneLength = 6,trControl = cvControl)
rfFit
plot(rfFit)
Q Quel algorihtme de boosting est utilisé? Quels sont ses paramètres? Quels sont ceux optimisés?
In [ ]:
#7 Boosting d'arbres
set.seed(2)
gbmFit = train(trainDescr, trainY,method = "gbm", tuneLength = 6,trControl = cvControl)
gbmFit
plot(gbmFit)
Q Quels autres paramètres viennent s'ajouter avec l'extrême boosting?
In [ ]:
#8 xgboost ou extrem gradient boosting
set.seed(2)
xgboostFit = train(trainDescr, trainY,method = "xgbTree", tuneLength = 8,trControl = cvControl)
plot(xgboostFit)
Q Quel noyau est utilisé? Quels sont les paramètres?
In [ ]:
#9 Support vector machine avec noyau gaussien
set.seed(2)
svmgFit = train(trainDescr, trainY,method = "svmRadial", tuneLength = 8,trControl = cvControl)
svmgFit
plot(svmgFit)
In [ ]:
models=list(knn=knnFit,cart=rpartFit,rf=rfFit,gbm=gbmFit,svmg=svmgFit,xgb=xgboostFit)
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. L'analyse discriminante qui ne fournit pas de probabilités d'occurence de la classe ``spam'' n'est pas adaptée au tracé d'une courbe ROC. Comme elle fait par ailleurs partie des moins performantes sur l'échantillon test, elles est retirée de la liste des modèles.
In [ ]:
# Courbes ROC
library(ROCR)
models=list(cart=rpartFit,rf=rfFit,gbm=gbmFit,xgb=xgboostFit)
testProb=predict(models, newdata = testDescr,type="prob")
predroc=lapply(testProb,function(x)prediction(x[,1],testY=="nospam"))
perfroc=lapply(predroc,function(x)performance(x, "tpr", "fpr"))
plot(perfroc$cart,col=1)
plot(perfroc$rf,col=2,add=TRUE)
plot(perfroc$gbm,col=3,add=TRUE)
plot(perfroc$xgb,col=4,add=TRUE)
legend("bottomright",legend=c("CART","RF","boost","XGBoost"),col=c(1:4),pch="_")
Q Interpréter les résultats obtenus. Quelles sont les méthodes les plus performantes à ce niveau de l'analyse?
Les meilleures méthodes (forêts aléatoires, boosting...) sont très peu explicites, car de véritables "boîtes noires", quant au rôle et impact des variables sur la prévision des observations. Néanmoins, des indicateurs d'importance sont proposés pour les forêts aléatoires ou XGBoost.
Q comment son définis les indicateurs d'importance?
In [ ]:
rfFit2=randomForest(trainDescr,trainY,importance=T)
imp.mdrr=sort(rfFit2$importance[,3],decreasing=T)[1:20]
par(xaxt="n")
plot(imp.mdrr,type="h",ylab="Importance",xlab="Variables")
points(imp.mdrr,pch=20)
par(xaxt="s")
axis(1,at=1:20,labels=names(imp.mdrr),cex.axis=0.8, las=3)
L'échantillon est de taille raisonnable, mais 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 comme sur la réalité des différences entre les méthodes. Il est donc important d'itérer le processus sur plusieurs échantillons tests.
Exécuter la fonction en annexe en choisissant les méthodes semblant les plus performantes. Attention au temps de calcul !
In [ ]:
# Choisir la liste des méthodes
# l'effort d'optimisation
# Initialiser le générateur et
# fixer le nombre d'itérations
models=c("gbm","rf","xgbTree")
noptim=c(6,6,6)
Niter=10 ; Init=11
pred.spam=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.spam$obs
prev.spam=pred.spam$pred
res.spam=lapply(prev.spam,function(x)apply((x>0.5)==(obs==1),2,mean))
# Moyennes des taux de bien classés par méthode
lapply(res.spam,mean)
In [ ]:
# distributions des taux de bien classés
boxplot(data.frame(res.spam))
Les commandes suivantes concernent les courbes ROC.
In [ ]:
# Apprécier la faible dispersion des
# courbes ROC associées aux forêts aléatoires
plot(perfroc.spam$rf,col=2)
plot(perfroc.spam$rf,col=1,add=TRUE,lwd=2,avg="vertical")
In [ ]:
# Comparaison des méthodes par le
# tracer des courbes ROC moyennes
# Problème pas identifié avec rlogit !
predroc.spam=lapply(prev.spam,function(x)prediction(x,obs==1))
perfroc.spam=lapply(predroc.spam,function(x)performance(x,"tpr","fpr"))
plot(perfroc.spam$gbm,col=1,lwd=2,avg="vertical")
plot(perfroc.spam$rf,col=2,add=TRUE,lwd=2,avg="vertical")
plot(perfroc.spam$xgbTree,col=3,add=TRUE,lwd=3,avg="vertical")
legend("bottomright",legend=c("boost","RF","xgbTree"),col=c(1:3),pch="_")
Q Quelle méthode semble finalement la plus pertinente?
Q Que suggérer à Georges pour améliorer son détecteur de pourriel? Comment éviter que vos messages ne soient "mangés" par les anti-spams ?
La plupart des détecteurs de spams sont des classifieurs bayésiens "améliorés". Consulter le site Wikipedia à ce sujet.
Q Quelle est la raison principale de ce choix de classifieur ? Quel est le problème rencontré ici, comment est-il résolu par ces détecteurs.
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
}