Modélisation prédictive et apprentissage statistique avec R, S. Tufféry - Editions Ophrys by TO Groupe

Stéphane TUFFÉRY

Responsable statistique dans un grand groupe bancaire français Université de Rennes 1 stephane.tuffery@univ-rennes1.fr

Modélisation prédictive et Apprentissage statistique avec R

2017

Éditions TECHNIP

5 avenue de la République, 75011 PARIS

CHEZ LE MÊME ÉDITEUR • Data Mining et statistique décisionnelle. L’intelligence des données S. TUFFÉRY

• Étude de cas en statistique décisionnelle S. TUFFÉRY

• Probabilités, analyse des données et statistique G. SAPORTA

• Les techniques de sondage P. ARDILLY

• Économie générale O. HUEBER

• Approche pragmatique de la classification J.P. NAKACHE, J. CONFAIS

• Statistique explicative appliquée J.P. NAKACHE, J. CONFAIS

• Modèles statistiques pour données qualitatives J.-J. DROESBEKE, M. LEJEUNE, G. SAPORTA, Eds.

• Plans d’expériences. Applications à l’entreprise J.-J. DROESBEKE, J. FINE, G. SAPORTA, Eds.

• Méthodes bayésiennes en statistique J.-J. DROESBEKE, J. FINE, G. SAPORTA, Eds.

• Approches non paramétriques en régression J.-J. DROESBEKE, G. SAPORTA, Eds.

• Analyse statistique des données spatiales J.-J. DROESBEKE, M. LEJEUNE, G. SAPORTA, Eds.

Tous droits de traduction, de reproduction et d’adaptation réservés pour tous pays.

Toute représentation, reproduction intégrale ou partielle faite par quelque procédé que ce soit, sans le consentement de l’auteur ou de ses ayants cause, est illicite et constitue une contrefaçon sanctionnée par les articles 425 et suivants du Code pénal. Par ailleurs, la loi du 11 mars 1957 interdit formellement les copies ou les reproductions destinées à une utilisation collective.

Avant-propos

Cet ouvrage explique et met en œuvre les principales méthodes de modélisation d’une situation dans laquelle on recherche la prédiction d’une variable binaire à l’aide de variables qualitatives ou quantitatives. Cette situation est très fréquente en data mining et en statistique, et elle est à la base entre autres des algorithmes de scoring. Nous le ferons à l’aide du logiciel libre R sur la base d’une étude de cas. Nous avons choisi ce logiciel, parce qu’il est aujourd’hui le plus utilisé par les statisticiens et les data scientists, et aussi en raison de la richesse et de la variété sans égales de ses fonctionnalités, couvertes par des modules appelés « packages » (ou « paquets »), au nombre de plus de dix mille huit cents au moment où nous écrivons ces lignes. De plus, les ressources consacrées à R sur Internet et dans la littérature statistique sont innombrables et permettent de trouver la réponse à presque n’importe quelle question. Pour aider le lecteur, nous avons indiqué quelques références dans la bibliographie en fin d’ouvrage. Les méthodes présentées ici couvrent le champ de la statistique, mais aussi de ce qui dans la terminologie anglo-saxonne s’appelle « statistical learning » et « machine learning ». Chaque chapitre de l’ouvrage est consacré à une méthode, et comporte des rappels théoriques, suivis de larges développements mettant en œuvre R sur un jeu de données que nous allons décrire. Les paramètres et les sorties des fonctions de R sont commentés, et les résultats obtenus pour chaque méthode sont confrontés d’une part avec la théorie et d’autre part avec les résultats obtenus avec les autres méthodes. Le lecteur sera ainsi à même de choisir les méthodes les plus adaptées à ses problématiques, et de les programmer efficacement. Les textes de référence sont indiqués s’il veut approfondir la théorie. Les bases de la statistique sont supposées connues du lecteur, ainsi que les concepts d’estimateur, d’échantillonnage bootstrap, de courbe ROC et d’aire sous la courbe ROC, qui sont juste rappelés en temps utile. Si nécessaire, le lecteur pourra se référer aux indications bibliographiques fournies ou à notre ouvrage Data Mining et statistique décisionnelle, paru aux Éditions Technip et réédité en 2017. Les étapes de notre étude de cas sont les suivantes. Nous commençons par importer, échantillonner et décrire les données. Nous mesurons la liaison entre la variable à expliquer et chacune des variables explicatives. Une phase de mise en forme des données aboutit à des regroupements de modalités pour les variables qualitatives, et une discrétisation des variables explicatives quantitatives. Une deuxième méthode de discrétisation supervisée est proposée, avec une procédure basée sur la maximisation de l’aire sous la courbe ROC de la prédiction à l’aide de la variable discrétisée. Nous mettons en œuvre le modèle classique de régression logistique logit, avec recours à la sélection pas à pas, puis à la sélection globale (algorithme de Furnival et Wilson et méthode directe). Nous en déduisons une grille de score et nous établissons les règles de décision basées sur cette grille de score. Nous testons aussi les modèles probit et log-log.

Avant-propos

Nous voyons ensuite comment les méthodes de pénalisation s’appliquent à la régression logistique : régression ridge, lasso, group lasso, et l’elastic net qui est un compromis entre ridge et lasso. Nous évoquons la non-consistance de l’estimateur lasso et les solutions que sont le relaxed lasso et l’adaptive lasso. Nous testons la régression logistique PLS et faisons quelques rappels sur les intervalles de confiance qu’il est possible d’obtenir par bootstrap dans ce genre de situation. Nous rappelons ensuite les principes de l’arbre CART, que nous mettons en œuvre, en insistant particulièrement sur son mécanisme d’élagage. Après cela, nous présentons une méthode de Friedman et Fisher dérivée des arbres de décision mais moins connue : le « bump hunting » (ou algorithme PRIM). Suivent les méthodes d’agrégation de modèles, dites aussi « méthodes d’ensemble » : le bagging, les forêts aléatoires, les Extra-Trees et le boosting (Discrete AdaBoost et Real AdaBoost). Une variante des forêts aléatoires est présentée, consistant à appliquer le mécanisme de double randomisation (individus et variables) des forêts aléatoires au modèle logistique et non plus à des arbres. Nous appliquons aussi l’algorithme du boosting au modèle logistique. Ce panorama des méthodes d’apprentissage supervisé se poursuit avec les Support Vector Machines (séparateurs à vaste marge), et plusieurs noyaux possibles : linéaire, radial, sigmoïde et polynomial. Nous montrons le bénéfice que nous pouvons attendre de l’utilisation des SVM sur les coordonnées factorielles issues d’une analyse des correspondances multiples. Le panorama s’achève avec les réseaux de neurones, et le boosting et les forêts aléatoires de réseaux de neurones, plus spécifiquement de perceptrons à une couche cachée. Le dernier chapitre présente une synthèse des méthodes précédentes, avec la comparaison de leur pouvoir prédictif et autres propriétés. Une première annexe aborde ensuite la détection des règles d’association par l’algorithme Apriori. Nous verrons que cette technique présente de l’intérêt même hors de son cadre habituel du web mining et de l’analyse des tickets de caisse, et qu’elle nous fournit quelques règles très intéressantes. Elle est en quelque sorte hybride, midescriptive, mi-prédictive, surtout lorsque l’on fixe le conséquent des règles recherchées. Une deuxième annexe est consacrée à la question des performances de R en grande dimension, lorsque le nombre d’observations ou de variables est très grand. Ce n’est pas le cas de notre jeu de données, mais d’une part, certaines solutions peuvent déjà, même avec un volume de données limité, apporter une diminution appréciable des temps de calcul, telle la parallélisation des traitements que nous appliquons à l’algorithme des forêts aléatoires, et d’autre part, nous voulons présenter des solutions qui seront pleinement appréciables sur des bases beaucoup plus grandes. Pour illustrer ces solutions, nous générons à partir de notre jeu de données une base de près de 65 millions d’observations. Outre la diminution des temps de calcul, le principal axe d’amélioration des performances concerne en effet la gestion des données trop massives pour être chargées en mémoire vive, comme l’exige le mode de fonctionnement normal de R. Plusieurs packages ont été développés ces dernières années pour traiter les données sur le disque dur de la machine, au lieu de les charger en mémoire vive. Nous

Avant-propos

VII

présentons dans cette annexe les packages ff et ffbase. Il faut aussi signaler le package bigmemory, qui présente l’intérêt d’être bien interfacé avec un ensemble de fonctions matricielles classiques ainsi qu’avec les fonctions d’analyse des données rassemblées dans le package biganalytics. Le package bigmemory présente aussi la particularité de pouvoir travailler soit en mémoire vive (mais de façon optimisée), ce qui est généralement plus rapide, soit sur le disque dur lorsque les volumes de données excèdent les capacités de la mémoire. Ce package présente toutefois la limite de ne manipuler que des big matrix numériques, et de devoir par conséquent recoder numériquement les facteurs. Nous présentons aussi deux packages spécifiques aux méthodes de régression en grande dimension : biglm et speedglm. Le premier peut être utilisé conjointement ou non avec les packages précédents, donc soit sur des data frames traditionnels, soit sur des objets big matrix, soit sur des data frames ffdf. Comme l’indique son nom, le package speedglm permet d’obtenir un gain significatif de temps de calcul, même sur des données en volume modéré. Dans une troisième annexe, nous montrons comment le logiciel SAS, avec son module IML, est capable de lancer l’exécution de code R, et nous l’illustrons en utilisant le package rpart pour obtenir des arbres CART dans SAS. À la fin de l’ouvrage, nous aurons construit un outil de credit scoring qui peut être utilisé pour l’octroi de crédits à la consommation dans une banque ou un établissement spécialisé de crédit. Le jeu de données utilisé présente l’intérêt d’être public, complet et de bonne qualité, mais il est évident que les méthodes vues dans ce document, et leur mise en œuvre avec R, s’appliqueraient de façon entièrement similaire à toute autre situation et tout autre jeu de données comportant des variables explicatives (quantitatives et/ou qualitatives) et une variable à expliquer binaire.

Préliminaire R est un logiciel libre et son code est public, mais en même temps centralisé sur des plate-formes de développement (entre autres GitHub et Bioconductor) et une plate-forme de référence : le CRAN1 (Comprehensive R Archive Network). Le CRAN est géré par l’organisation responsable de R, la « R Foundation for Statistical Computing », qui est une organisation à but non lucratif enregistrée en Autriche avec son siège à Vienne. Elle possède le copyright sur le logiciel R et sa documentation. Elle organise et met à disposition des utilisateurs de R des infrastructures, des moyens d'informations et d'échanges, du support, de la documentation. Le noyau de R ne peut être modifié que par les membres de la « R Development Core Team », mais le CRAN est ouvert aux développeurs qui veulent y déposer des packages. Ceux-ci doivent toutefois être vérifiés et documentés. Le CRAN est un repository historisé : toutes les versions successives sont archivées et accessibles.

http://cran.r-project.org/

VIII

Avant-propos

Pour effectuer les calculs décrits dans cet ouvrage, vous devez avoir installé R sur votre machine. Vous pouvez le télécharger sur un site miroir du CRAN et son installation est simple et rapide : il suffit de sélectionner son système d’exploitation (Unix, Windows ou macOS) et de lancer l’installation en suivant les instructions propres à chaque système, en sélectionnant éventuellement la version 32 ou 64 bits. Tous les calculs décrits dans l’ouvrage ont été réalisés avec des versions de R au moins égales à 2.15. Il faudra en outre installer et charger certains packages supplémentaires, qui ne sont pas inclus dans l’installation de base de R et que nous indiquerons au fur et à mesure de leur utilisation. L’installation se fait par la commande install.packages, et n’est à faire qu’une seule fois. Elle nécessite d’être connecté à Internet ou d’avoir préalablement téléchargé le package sur le site du CRAN. On peut régulièrement vérifier si une mise à jour du package est disponible sur ce site. Pour ces actions, pour que le système puisse installer les packages dans des répertoires autorisés, il est possible qu’il pose la question « Would you like to use a personal library instead? » à laquelle il faudra répondre « oui ». Le chargement d’un package se fait par la commande library et est à refaire après chaque ouverture d’une nouvelle session de R. Pour une utilisation plus agréable de R, plutôt que la « console » de base, on peut aussi utiliser un éditeur de texte amélioré comme par exemple Notepad++ (et son plugin NpptoR2) ou l’interface de développement RStudio3. Notepad++ est très simple à installer et permet de disposer d’une numérotation des lignes de code, de leur coloriage et de la possibilité de continuer à saisir du code pendant que R effectue les calculs. De son côté, RStudio Server permet de travailler en client-serveur, avec R s’exécutant sur le serveur, et sur le client une interface Web utilisable avec un navigateur Internet. Le lecteur débutant en R trouvera dans la bibliographie à la fin de notre ouvrage des références utiles pour s’initier à R. Parmi les nombreuses ressources disponibles sur l’Internet, nous apprécions la « Short R Reference Card »4 ainsi que le cours d’Anne Philippe à l’Université de Nantes.

Remerciements Je remercie vivement Arthur Charpentier pour ses suggestions utiles, Brigitte Gelein pour sa relecture attentive et ses conseils avisés, Olivier Decourt pour son esprit critique et ses remarques stimulantes, Franck Berthuit pour ses observations pertinentes, Magalie Fromont pour nos échanges, et Jérémy Magnanensi pour ses renseignements sur le package plsRglm.

http://npptor.sourceforge.net/

http://www.rstudio.com/

http://cran.r-project.org/doc/contrib/Short-refcard.pdf

Table des matières

Avant-propos .................................................................................................................... V Présentation du jeu de données ....................................................................................... 1 Préparation des données .................................................................................................. 7 2.1. 2.2. 2.3. 2.4. 2.5.

Lecture des données dans un fichier texte........................................................................... 7 Lecture des données dans une table SAS ou SPSS............................................................ 8 Premières transformations ................................................................................................... 8 Générateurs de nombres aléatoires..................................................................................... 9 Échantillonnage ................................................................................................................. 10

Exploration des données ................................................................................................ 13 3.1. 3.2. 3.3. 3.4. 3.5.

Vue d’ensemble des données............................................................................................ 13 Analyse des variables continues ........................................................................................ 15 Discrétisation des variables continues ............................................................................... 19 Liaison des variables explicatives avec la variable à expliquer .......................................... 24 Liaisons entre les variables explicatives ............................................................................ 37

Discrétisation automatique supervisée des variables continues ................................ 43 La régression logistique ................................................................................................. 57 5.1. 5.2. 5.3. 5.4. 5.5. 5.6. 5.7. 5.8. 5.9. 5.10. 5.11.

Sélection pas à pas des variables ..................................................................................... 58 Sélection globale des variables ......................................................................................... 69 Recherche des prédicteurs maximisant l’aire sous la courbe ROC .................................... 83 Modélisation avec tous les croisements possibles de variables ......................................... 85 Regroupements définitifs des modalités ............................................................................ 88 Modèle logit retenu ............................................................................................................ 91 Grille de score ................................................................................................................... 96 Détermination des seuils de score ................................................................................... 104 Courbe ROC et courbe de lift........................................................................................... 107 Modèle probit ................................................................................................................... 110 Modèle log-log ................................................................................................................. 114

La régression logistique pénalisée ridge .................................................................... 117 6.1. 6.2. 6.3.

Principe de la pénalisation ridge ...................................................................................... 117 La régression ridge dans R .............................................................................................. 119 Grille de score du modèle ridge ....................................................................................... 132

X 6.4.

Table des matières Affinement du modèle ridge ............................................................................................. 135

La régression logistique pénalisée lasso .................................................................... 139 7.1. 7.2. 7.3. 7.4. 7.5. 7.6.

Principe de la pénalisation lasso ...................................................................................... 139 Le lasso dans R ............................................................................................................... 141 Non-consistance de l’estimateur lasso............................................................................. 149 Relaxed lasso, SCAD et adaptive lasso ........................................................................... 151 Le « group lasso » ........................................................................................................... 154 La pénalisation elastic net ............................................................................................... 159

La régression logistique PLS ....................................................................................... 161 L’arbre de décision CART ............................................................................................. 175 9.1. 9.2. 9.3. 9.4.

Introduction aux arbres de décision ................................................................................. 175 Les arbres de décision dans R ........................................................................................ 179 Complexité et élagage d’un arbre de décision ................................................................. 184 Prédiction et mesure du pouvoir discriminant d’un arbre de décision .............................. 192

L’algorithme PRIM ......................................................................................................... 199 10.1. 10.2. 10.3. 10.4. 10.5.

Principe de l’algorithme PRIM.......................................................................................... 199 PRIM en deux dimensions ............................................................................................... 201 PRIM sur l’ensemble des variables quantitatives ............................................................. 205 PRIM sur les variables quantitatives et qualitatives ......................................................... 207 PRIM sur les coordonnées factorielles des variables explicatives .................................... 211

Les forêts aléatoires ...................................................................................................... 219 11.1. 11.2. 11.3. 11.4. 11.5. 11.6. 11.7.

Introduction aux forêts aléatoires ..................................................................................... 219 Les forêts aléatoires dans R ............................................................................................ 222 Estimateurs de l’erreur de prédiction ............................................................................... 228 Le graphique de proximité ............................................................................................... 232 Mesure de l’importance des variables .............................................................................. 234 Choix du nombre de variables ......................................................................................... 240 Les Extra-Trees ............................................................................................................... 246

Le bagging...................................................................................................................... 249 Les forêts aléatoires de modèles logistiques ............................................................ 255 Le boosting .................................................................................................................... 271 14.1. 14.2. 14.3. 14.4. 14.5.

Introduction au boosting .................................................................................................. 271 Discrete AdaBoost ........................................................................................................... 275 Real AdaBoost ................................................................................................................ 280 Compléments sur le boosting .......................................................................................... 285 Boosting de modèles logistiques ..................................................................................... 288

Les Support Vector Machines ...................................................................................... 293 15.1. Introduction aux machines à vecteurs de support ............................................................ 293

Table des matières 15.2. 15.3. 15.4. 15.5. 15.6. 15.7. 15.8. 15.9. 15.10.

SVM à noyau linéaire....................................................................................................... 299 SVM à noyau linéaire : le calcul des coefficients des prédicteurs .................................... 305 SVM à noyau radial gaussien .......................................................................................... 308 SVM à noyau radial laplacien .......................................................................................... 310 SVM à noyau sigmoïde.................................................................................................... 314 SVM à noyau polynomial ................................................................................................. 315 SVM après réduction de dimension ................................................................................. 319 Expression explicite du modèle SVM sur composantes factorielles ................................. 326 Conclusion....................................................................................................................... 333

Les réseaux de neurones .............................................................................................. 335 16.1. 16.2. 16.3. 16.4. 16.5.

Introduction aux réseaux de neurones ............................................................................. 335 Modélisation avec le perceptron à une couche cachée .................................................... 339 Les poids d’un réseau de neurones ................................................................................. 345 Le boosting de perceptrons à une couche cachée ........................................................... 348 Les forêts aléatoires de perceptrons à une couche cachée ............................................. 352

Synthèse des méthodes prédictives ............................................................................ 357 Annexe 1 : La détection des règles d’association ...................................................... 361 Annexe 2 : La performance des calculs en grande dimension.................................. 377 19.1. 19.2. 19.3. 19.4.

La parallélisation dans R ................................................................................................. 377 La parallélisation des forêts aléatoires de modèles logistiques ........................................ 380 La parallélisation du package randomForest ................................................................ 383 Les packages R pour le traitement des données massives ............................................. 385

Annexe 3 : L’utilisation de R dans le logiciel SAS ...................................................... 395 Bibliographie .................................................................................................................. 403 Index des packages R utilisés ...................................................................................... 407

Chapitre 1

Présentation du jeu de données

Le jeu de données utilisé dans l’ouvrage est connu sous le nom de « German Credit Data » et il est disponible sur Internet5. Il contient 1 000 dossiers de crédit à la consommation, dont 300 avec des impayés. En réalité, il s’agit d’un échantillon stratifié et redressé sur la variable à expliquer, le taux d’impayés réel étant sûrement plus proche de 5 %. Ce jeu de données a d’abord été fourni par le professeur Hans Hofmann, de l’Université de Hambourg. Il a été utilisé par L. Fahrmeir et G. Tutz dans leur ouvrage Multivariate Statistical Modelling Based on Generalized Linear Models (1994). Il se compose de la variable à expliquer « cible », avec sa modalité 1 (pas d’impayé : 700 observations) et sa modalité 2 (présence d’impayés : 300 observations), et de 19 variables explicatives (nous en avons omise une sur l’origine du demandeur) : trois variables numériques continues : la durée du crédit demandé en mois, le montant du crédit en euros et l’âge du demandeur en années ; six variables numériques discrétisées : le solde moyen sur compte courant, l’encours d'épargne, le nombre de crédits déjà détenus dans la banque, le taux d’effort (ou taux d’endettement), l’ancienneté au domicile, le nombre de personnes à charge ; dix variables qualitatives : l’objet du crédit, l’historique de remboursement du demandeur, les autres crédits détenus (hors de la banque), les biens de valeur détenus par le demandeur, ses garanties, sa situation familiale, son ancienneté dans l’emploi, son statut au domicile, son type d’emploi occupé et la fourniture d’un numéro de téléphone. Certaines variables portent sur le crédit demandé (durée, montant, objet), d’autres sur le profil financier du demandeur (encours, autres créditsY) et les dernières sur le profil personnel du demandeur (âge, situation de familleY). Toutes ces données sont bien sûr observées au moment de la demande de crédit et non ultérieurement. Ce jeu de données présente l’intérêt d’être réel ou au moins réaliste, simple à décrire et à appréhender, d’avoir des variables diversifiées et d’être suffisamment riche pour permettre des analyses variées, sans être trop complexe. Il présente l’inconvénient d’être un peu petit. En particulier, avec seulement 300 cas positifs6, il n’est visiblement pas possible de découper ce jeu de données en trois échantillons indépendants de taille

Notamment à l’adresse https://archive.ics.uci.edu/ml/datasets/Statlog+(German+Credit+Data).

Nous nous conformons à l’usage en parlant de cas positifs, non pour indiquer un jugement, mais pour

signaler les cas que nous cherchons à détecter. Ici, les cas positifs sont les dossiers avec des impayés, de même que le résultat positif d’un test médical peut signifier la présence d’une pathologie.

1. Présentation du jeu de données

suffisante, pour l’apprentissage, la validation et le test. Or, la méthode rigoureuse d’élaboration d’un modèle consiste à l’ajuster sur un échantillon d’apprentissage, à utiliser l’échantillon de validation pour choisir entre plusieurs modèles possibles celui qui semble le mieux se généraliser, puis à utiliser l’échantillon de test7 pour mesurer avec le moins de biais possible la qualité de généralisation du modèle, par exemple en mesurant son aire sous la courbe ROC8. L’échantillon de validation qui permet d’optimiser les paramètres des modèles ne devrait pas être confondu avec l’échantillon de test qui permet de comparer les modèles entre eux9. Le même échantillon ne devrait pas intervenir dans la construction et la comparaison des modèles. C’est pourtant ce à quoi nous devrons nous résoudre ici, ne pouvant scinder en deux un échantillon de test qui ne contient déjà qu’une centaine de cas positifs. Notre échantillon de validation sera donc aussi celui de test. Que le lecteur veuille bien nous pardonner cette entorse à la méthodologie statistique. Les formats des variables autres que continues sont disponibles sur Internet à l’adresse indiquée. Nous les indiquons ci-après. Le solde sur compte courant : COMPTES 'A11' = 'CC < 0 euros ' 'A12' = 'CC ∈ [0-200 euros[' 'A13' = 'CC >= 200 euros' 'A14' = 'Pas de compte ' Nous notons que l’établissement financier n’est pas un pur établissement de crédit ne gérant pas les comptes courants, puisque dans certains cas il peut détenir le compte courant du client et aussi ses comptes d’épargne (voir plus bas). 7

Cet échantillon ne provient pas toujours de la même population, mais parfois d’un tirage effectué à une

autre date. On parle d’échantillon « hors-temps » (« out of time »). 8

Rappelons que l’aire sous la courbe ROC d’un modèle prédictif binaire est la probabilité que la prédiction du modèle vérifie p(x) > p(y) si x est un cas positif tiré au hasard et y est un cas négatif tiré au hasard. Si le modèle parvient à séparer parfaitement les positifs des négatifs, cette aire sera égale à 1. Elle sera en revanche égale à 0,5 si la prédiction du modèle n’est pas supérieure à une prédiction aléatoire, et égale à 0 si la prédiction du modèle est toujours fausse. Le nom de cette probabilité vient de ce qu’elle est égale à l’aire comprise entre la courbe ROC et l’axe des abscisses. La courbe ROC (« Receiver Operating Characteristic ») est l’ensemble des points (a(s),b(s)), où, pour s décrivant toutes les valeurs prédites possibles de la plus grande à la plus petite, a(s) est la proportion de négatifs x pour lesquels p(x) ≥ s, et b(s) est la proportion de positifs x pour lesquels p(x) ≥ s. On appelle les derniers les « vrais positifs » et les premiers les « faux positifs », puisque ce sont des négatifs considérés comme positifs en raison de la valeur p(x) ≥ s. La courbe ROC représente donc la proportion de vrais positifs en fonction de celle de faux positifs, pour les différents seuils s possibles. La courbe ROC répond à une problématique habituelle dans de nombreux domaines qui vont de l’analyse des signaux radars à la recherche médicale : détecter le plus possible de positifs tout en ayant le moins possible de fausses alertes. 9

C’est tout particulièrement vrai pour les réseaux de neurones, dont le « early stopping » largement

pratiqué consiste à interrompre l’apprentissage lorsque l’erreur commence à remonter sur l’échantillon de validation.

Chapitre 2

Préparation des données

La description des données étant connue, nous allons y accéder et les manipuler dans une structure de données appelée « data frame » dans R. Nous constituons aussi un échantillon d’apprentissage et un échantillon de test. L’utilisation d’une « graine » dans le générateur de nombres aléatoires permettra d’obtenir des échantillons reproductibles par le lecteur.

2.1. Lecture des données dans un fichier texte Le plus simple est de lire le jeu de données dans un fichier texte qui aura été récupéré et stocké sur votre disque dur, votre serveur, ou en accédant directement à Internet. > credit <- read.table("http://archive.ics.uci.edu/ml/machinelearning-databases/statlog/german/german.data", h=FALSE, col.names=myVariableNames) > credit <- read.table("http://blogperso.univrennes1.fr/stephane.tuffery/public/german.txt", h=FALSE, col.names=myVariableNames) Ici, il faut avoir préalablement défini le vecteur du nom des variables. > myVariableNames <-c ("Comptes", "Duree_credit", "Historique_credit", "Objet_credit", "Montant_credit", "Epargne", "Anciennete_emploi", "Taux_effort", "Situation_familiale", "Garanties", "Anciennete_domicile", "Biens", "Age", "Autres_credits", "Statut_domicile", "Nb_credits", "Type_emploi", "Nb_pers_charge", "Telephone", "Etranger", "Cible") Ou en anglais : > myVariableNamesE <-c ("checking_status", "duration", "credit_history", "purpose", "credit_amount", "savings", "employment", "installment_rate", "personal_status", "other_parties", "residence_since", "property_magnitude", "age", "other_payment_plans", "housing", "existing_credits", "job", "num_dependents", "telephone", "foreign_worker", "class")

2. Préparation des données

On ajoute un identifiant qui n’est autre que le numéro de la ligne. Il n’est pas utilisé par la suite, mais pourrait être utile en cas de tri de la base qui pourrait modifier l’ordre des observations et altérer la correspondance établie en section 2.3 avec les numéros des observations sélectionnées pour les échantillons d’apprentissage et de test qui seront utilisés dans la suite de cette étude. > credit$Cle <- seq(1, nrow(credit)) Nous supprimons une variable. > credit$Etranger <- NULL

2.2. Lecture des données dans une table SAS ou SPSS Le package sas7bdat permet à R de lire une table SAS à condition qu’elle ne soit pas compressée (option « compress » de SAS). > install.packages("sas7bdat") > library(sas7bdat) > credit <- read.sas7bdat("C:\\Users\\...\\Etude de cas\\credit.sas7bdat") Ce package ne permet pas de sélectionner les variables à lire dans la base. Quant à la lecture d’une table SPSS, elle peut se faire à l’aide du package foreign10. > library(foreign) > credit <- read.spss("C:\\Users\\...\\Etude de cas\\credit.sav", use.value.labels=TRUE, max.value.labels=Inf, to.data.frame=TRUE)

2.3. Premières transformations Nous commençons par recoder les valeurs 1 / 2 (« OK / KO ») de la variable à expliquer par 0 / 1, parce que la codification 0 / 1 est acceptée par tous les packages. Nous déclarons aussi cette variable comme qualitative par la commande factor. > credit$Cible[credit$Cible == 1] <- 0 # crédits OK > credit$Cible[credit$Cible == 2] <- 1 # crédits KO 10

Le package foreign ne permet toutefois pas d’importer un format SPSS compressé ZSAV. Ce

format apparu dans la version 21 de SPSS fait appel à une compression dont le taux est sensiblement plus élevé que la compression antérieure.

Chapitre 4

Discrétisation automatique supervisée des variables continues

Dans le chapitre 3, nous avons proposé une discrétisation des variables continues par cette méthode : découpage en déciles (ou d’autres quantiles, selon l’effectif), examen des taux de positifs (ici, d’impayés) par déciles, et regroupement des déciles dont les taux de positifs sont proches, en tenant compte également des effectifs par déciles et en évitant les tranches contenant trop peu de positifs, dont le taux serait assorti d’un large intervalle de confiance. Cette méthode rappelle bien sûr l’algorithme CHAID, qui fusionne les modalités adjacentes formant avec la variable à expliquer des tableaux 2 x 2 dont la probabilité du χ² est supérieure au seuil de significativité choisi. Nous avons ainsi obtenu les discrétisations suivantes : pour l’âge : moins de 25 ans et plus de 25 ans ; pour la durée de crédit : moins de 15 mois, entre 15 et 36 mois, et plus de 36 mois ; pour le montant de crédit : moins de 4000 euros et plus de 4000 euros. Nous allons maintenant exposer une autre méthode permettant d’obtenir automatiquement une discrétisation supervisée des variables. Sur notre jeu de données, il se trouve qu’elle fournira de nouvelles discrétisations des trois variables précédentes très proches de celles déjà obtenues. Cette méthode oblige d’abord à fixer le nombre k de classes souhaité puis consiste à rechercher le découpage de la variable en k classes qui maximise l’aire sous la courbe ROC d’un modèle logistique ajusté sur la variable ainsi discrétisée. Plus précisément, voici les étapes de la méthode pour une discrétisation supervisée de la variable X par rapport à une variable à expliquer Y : 1. Fixer le nombre de classes k (paramètre nbmod plus bas). 2. Calculer les k-tiles (par exemple, les quartiles si k = 4), en ayant préalablement isolé les éventuelles valeurs manquantes dans une k+1e classe. 3. Ajuster un modèle logit Y ~ Xd, Xd étant la variable catégorielle dont les facteurs sont les classes déterminées à l’étape précédente : ]– ∞ , q1], ]q1 , q2] , Y, ]qk-1, + ∞[. 4. Mesurer l’aire sous la courbe ROC du modèle logit, ou l’aire sous une portion de la courbe ROC comprise entre les spécificités 1 et un paramètre pAUC compris entre 0 et 1. 5. Modifier les seuils qi des classes et revenir à l’étape 3) de façon à obtenir une variable Xd modifiée et un modèle logit correspondant dont l’aire sous la courbe ROC soit plus élevée que l’aire précédente. 6. S’interrompre lorsque l’aire sous la courbe ROC n’augmente plus significativement.

4. Discrétisation automatique supervisée des variables continues

Les étapes 3 à 6 sont programmées dans R comme suit. Nous définissons une fonction fitauc réalisant les étapes 3 et 4 en acceptant un vecteur de seuils qi en entrée et en renvoyant en sortie l’aire sous la courbe ROC du modèle correspondant. Les étapes 5 et 6 sont réalisées à l’aide de la fonction optim qui recherche le minimum d’une fonction par un algorithme qui est par défaut celui de Nelder-Mead (méthode du simplexe) mais peut aussi être celui de quasi-Newton, de recuit simulé, de descente du gradient ou de Brent. Ici, nous recherchons le minimum de l’opposé de l’aire sous la courbe ROC, et l’algorithme de Nelder-Mead donne généralement de bons résultats. Les étapes 3 à 6 sont effectuées sur les vecteurs X et Y obtenus après exclusion des éventuelles valeurs manquantes de X (on suppose Y sans valeur manquante), et les découpages obtenus sont ensuite appliqués sur les vecteurs Xt et Yt incorporant les valeurs manquantes. Le nombre de classes final peut être strictement inférieur à k si deux seuils de discrétisation coïncident. Nous pouvons choisir de maximiser, non pas l’aire sous l’ensemble de la courbe ROC, mais l’aire sous une portion comprise entre les spécificités 1 et un paramètre pAUC compris entre 0 et 1, soit une abscisse comprise entre 0 et 1 – pAUC. Ce choix peut être légitime si l’on veut optimiser le pouvoir prédictif, non pas sur l’ensemble de la population, mais sur les individus aux scores les plus élevés. C’est parfois cet objectif qui est recherché, par exemple dans un score de détection de la fraude, ou dans un score d’appétence où il importe de savoir détecter les clients ayant la plus forte probabilité d’acheter, plus que de distinguer ceux qui ont une faible ou moyenne probabilité d’acheter. Le package pROC16 sait calculer cette aire sous la courbe ROC partielle. L’ensemble des étapes 1 à 6 est rassemblé dans la fonction decoup, fonction de discrétisation automatique dont les paramètres sont : le nom du fichier (base) ; le nom de la variable à découper (x) ; le nom de la variable à expliquer (y) ; le nombre de modalités à produire (nbmod), par défaut égal à 3 (en l’absence d’autre précision) ; un paramètre h utilisé pour renforcer l’homogénéité des effectifs des classes produites ; un paramètre k utilisé pour renforcer l’hétérogénéité des taux de positifs des classes produites ; le paramètre pAUC décrit plus haut ; le paramètre calcul qui vaut 1 par défaut pour déclencher les traitements ici décrits de recherche d’une discrétisation, et vaut 0 si l’on veut directement appliquer des seuils fournis dans un vecteur bornes et obtenir l’affichage des caractéristiques de ce découpage ; l’algorithme d’optimisation, par défaut celui de Nelder-Mead ; un paramètre graphe qui vaut par défaut 0 pour ne pas afficher de graphique, et dont la valeur 1 déclenche l’affichage du graphique décrit plus loin.

Xavier Robin, Natacha Turck, Alexandre Hainard, Natalia Tiberti, Frédérique Lisacek, Jean-Charles

Sanchez and Markus Müller (2011). pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics, 12, p. 77.

Chapitre 6

La régression logistique pénalisée ridge

Nous savons que la variance d’un modèle augmente et que son biais diminue avec l’accroissement de sa complexité. La tâche classique du modélisateur consiste à trouver le bon niveau de complexité, assurant un juste équilibre entre la variance et le biais, et évitant dans les modèles de régression les inconvénients liés à la multicolinéarité des prédicteurs. Mais il n’est pas toujours possible de simplifier le modèle aussi qu’on le souhaiterait : parfois ses utilisateurs, les experts du domaine, requièrent l’intervention simultanée de tous les critères disponibles, même s’ils sont très corrélés. La méthode de base pour utiliser un grand nombre de prédicteurs sans pâtir de l’accroissement de la complexité et de la colinéarité est la « pénalisation ». Une régression pénalisée est une régression dont les coefficients sont recherchés en optimisant un critère habituel, mais avec une contrainte supplémentaire qui limite la norme de ces coefficients. Du choix de la norme dépend la solution trouvée. La plus courante est la norme L², qui donne lieu à la régression ridge dont nous parlons maintenant.

6.1. Principe de la pénalisation ridge De même que dans la régression linéaire ridge l’erreur quadratique est pénalisée, dans la régression logistique ridge c’est la log-vraisemblance qui est pénalisée. On recherche donc les coefficients N? qui minimisent l’expression : $

W XY<∑Z XZ[Z

−2 E T' log V \ + F1 − ' G log V1 − \] + ^ E N? XY <∑Z XZ [Z9 XY <∑Z XZ [Z9 1+W 1+W _; ?_;

λ étant un paramètre ≥ 0. Il est équivalent de chercher à minimiser : $

W XY<∑Z XZ[Z 9

W XY<∑Z XZ [Z 9

−2 E T' log V \ + (1 − ' )log V1 − \] XY <∑Z XZ [Z9 XY <∑Z XZ [Z9 1+W 1+W _;

% sous la condition ∑?_; N? ≤ a (C ≥ 0 est lié à λ). Pour cette raison, on parle aussi en français de régression « bornée » et en anglais de « shrinkage method » (rétrécissement). Plus λ est grand (plus C est petit) et plus les coefficients sont « rétrécis ». L’idée est d’éviter ce qui se produit parfois en présence de forte colinéarité entre les prédicteurs :

118

6. La régression logistique pénalisée ridge

l’apparition de coefficients instables et de coefficients très grands se « compensant », au prix éventuellement d’incohérences de signes. La régression logistique ridge fait partie des méthodes dites de pénalisation (ou régularisation). Elle se généralise dans deux directions. D’une part, le principe de pénalisation de la log-vraisemblance permet de l’étendre au modèle linéaire généralisé. La variable à expliquer peut être quantitative (normale, Gamma ou de Poisson, par exemple), binaire, qualitative (nominale) ou même de durée (régression de Cox)22. D’autre part, la pénalisation peut être de n’importe quelle forme c

^ ∑ _;bN? b , %

avec δ ≥ 0, les valeurs de δ comprises entre 0 et 1 étant plus adaptées aux cas où de nombreux coefficients à estimer sont faibles. On reconnaît dans cette expression une norme Lδ. La difficulté calculatoire est de trouver suffisamment rapidement les valeurs optimales simultanées de λ et δ, comme l’a montré Jerome Friedman23. À chaque valeur de δ correspond un domaine de contrainte, une région définie par la pénalisation : en dimension 2, c’est un disque N; + N ≤ D lorsque δ = 2, un losange |N; | + |N | ≤ D lorsque δ = 1, une étoile lorsque δ < 1, étoile se contractant vers les deux axes lorsque δ se rapproche de 0. Les branches de l’étoile correspondent aux cas où l’un c

des coefficients à estimer est faible. Dans le cas limite où δ = 0, la somme ∑ _;bN? b n’est autre que le nombre de coefficients non nuls, et l’effet de l’augmentation de la pénalisation est la diminution du nombre de variables. Lorsque δ < 1, le domaine de contrainte n’est pas convexe : cela pose des problèmes d’optimisation et ce cas est rare hormis celui où δ = 0. Les valeurs δ > 2 sont également peu fréquentes, et les valeurs les plus usitées sont δ = 0, 1 et 2. Le cas δ = 0 est celui qui consiste à appliquer une pénalisation basée sur le nombre de coefficients, telle que l’AIC ou le BIC, et nous avons vu en section 5.1 qu’elle peut déjà donner des résultats satisfaisants, notamment avec le BIC plus sévère. L’augmentation progressive de la pénalisation conduit à éliminer une à une des variables du modèle dans une sélection pas à pas. Cette méthode rapide ne conduit néanmoins pas à des résultats optimaux, et ce d’autant moins que le nombre de variables est grand. Par ailleurs, une sélection globale consistant à rechercher le sous-ensemble optimal de variables, telle que celle basée sur l’algorithme leaps and bounds de Furnival et Wilson (voir la section 5.2), n’est en général pas envisageable lorsque le nombre de variables dépasse 40 ou 50. Pour une méthode calculable avec un grand nombre de variables, et néanmoins performante, on est amené à se placer dans la première situation convexe rencontrée en faisant croître δ, celle où δ = 1. Nous reviendrons dans le chapitre suivant sur cette pénalisation appelée lasso, et nous verrons qu’elle a aussi pour effet d’annuler des coefficients, comme dans le cas précédent. %

Voir notamment : Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika 80: 27-

38. Et : Heinze, G. (2006). A comparative investigation of methods for logistic regression with separated or nearly separated data. Statistics in Medicine 25: 4216-4226. 23

Voir

ses

présentations

http://www-stat.stanford.edu/~jhf/talks/GPS_R.pdf

stat.stanford.edu/~jhf/ftp/GPSpaper.pdf.

http://www-

Chapitre 16

Les réseaux de neurones

Les réseaux de neurones connurent un grand succès dans les années 1980, avec un prestige venant de leur structure réputée imiter le cerveau humain, organisé en neurones et synapses. Leur réglage est délicat et leur pouvoir prédictif parfois élevé n’est pas toujours facile à obtenir, mais ils offrent dans ce chapitre l’occasion de revenir sur les notions de early stopping et de weight decay, qui rappellent ce que l’on a vu en boosting.

16.1. Introduction aux réseaux de neurones79 Les réseaux de neurones se présentent comme des ensembles d’unités (encore appelées neurones formels ou nœuds) connectées entre elles, chaque unité étant activée en recevant une donnée et en en renvoyant une autre. Ces ensembles d’unités sont organisés en niveaux appelés couches, et chaque couche envoie ses données à la couche suivante, comme des neurones transmettant un signal par leurs synapses.

n1 X1

p1 q1

X2 X3 données

s(Σn ipi)

Σ niqi

...

n4 ...

Σ nipi

s [Π.s(Σnipi) + Θ.s(Σ niq i)]

s(Σniqi)

p5 n5

couche cachée

couche de sortie

couche d'entrée

Figure 16.1 – Réseau de neurones avec une couche cachée

Bishop, C. (1995). Neural Networks for Pattern Recognition, Clarendon Press, Oxford.

Ripley, B. D. (1996). Pattern Recognition and Neural Networks, Cambridge University Press.

336

16. Les réseaux de neurones

La Figure 16.1 représente l’un des modèles génériques de réseaux de neurones : le perceptron multicouches, ici un perceptron à une couche cachée. Dans le perceptron, la connexion entre deux unités i et j, c'est-à-dire l’équivalent d’une synapse, se fait en multipliant la sortie ni de l’unité i par un certain poids pij. L’unité j est activée par le signal ∑ ( ? ƒ , la somme étant prise sur l’ensemble des unités i de la couche précédente. On peut ajouter une « constante » au modèle avec une unité particulière (le « biais ») qui transmet toujours la valeur 1 (voir la Figure 16.4). On peut éventuellement appliquer une fonction d’activation s(∑ ( ? ƒ ) à ce signal. Le premier niveau, la couche d’entrée, est composé d’unités qui contiennent les données de l’échantillon d’apprentissage. Chaque variable quantitative ou binaire correspond à une unité, et chaque variable qualitative correspond à autant d’unités qu’elle a de modalités. Les réseaux de neurones peuvent être utilisés pour le classement ou la prédiction, et dans les réseaux utilisés pour la prédiction (régression ou classement), la variable à expliquer correspond à un niveau appelé couche de sortie, avec une seule unité dans le cas de la régression, et une unité pour chaque modalité de la variable à expliquer lorsque celle-ci en a au moins trois. Si la variable à expliquer est binaire, la couche de sortie peut avoir un ou deux unités, comme précisé plus bas. Entre la couche d’entrée et la couche de sortie sont parfois connectées des unités appartenant à un niveau intermédiaire : la couche cachée, ainsi dénommée car elle correspond à des données non apparentes. Il peut exister plusieurs couches cachées. Quand un signal arrive à une couche cachée, on lui applique généralement une fonction d’activation qui est la sigmoïde 1/(1+e-x) ou la tangente hyperbolique (ex–e-x)/(ex+e-x). Cela permet de transmettre au niveau suivant une donnée normalisée (dans [-1 ; +1] pour la tangente hyperbolique, ou dans [0 ; +1] pour la sigmoïde), dont on verra plus loin l’intérêt pour l’activation d’une unité. Les premiers réseaux de neurones80 utilisaient aussi comme fonction d’activation la fonction de seuil, valant 1 si x dépassait un certain seuil, et valant 0 sinon. Mais une fonction non dérivable pose des difficultés pour l’optimisation, et de plus les unités « saturent » immédiatement à la réception d’un signal dépassant le seuil, au lieu de réagir progressivement. Lorsque le signal arrive à la couche de sortie, on peut appliquer ou non une fonction d’activation s(x) différente de l’identité à ∑ ( ? ƒ . Pour les problèmes de régression, il est plutôt conseillé de prendre s(x) = x ou une fonction linéaire s(x) = ax+b, tandis que pour les problèmes de classement, il est plutôt conseillé, selon les cas : pour le classement binaire Y = 0/1 avec une seule unité de sortie, de prendre pour fonction d’activation la sigmoïde (fonction logistique), qui donne une sortie dans [0 ; 1], estimateur de P(Y=1|x) ; pour le classement multi-classes (et éventuellement aussi binaire), avec une unité par classe, d’appliquer la fonction « softmax »

W [¢ -@ (O) = £ [ ∑J_; W 6

à la kème composante xk de la sortie ; ceci assure d’avoir une somme ∑£ @_; -@ (O) = 1 sur l’ensemble des K classes, et chaque sortie est un estimateur P(Y=k|x). 80

McCulloch, S. W. and Pitts, W. (1943). A logical calculus of ideas immanent in neural activity. Bulletin

of Mathematical Biophysics 5, 115-133.

16. Les réseaux de neurones

337

Comme dans tout apprentissage supervisé, il faut définir une fonction d’erreur, et choisir l’algorithme utilisé pour l’optimiser. Pour les problèmes de régression, la fonction d’erreur est généralement l’erreur quadratique (L2) : •F , (O, ')G = ‖ (O) − '‖² = E( (O ) − ' ) .

On peut aussi recourir à l’erreur absolue (L1) :

•F , (O, ')G = E| (O ) − ' |.

L’erreur quadratique ou absolue est moins adaptée aux problèmes de classement, et on utilise généralement l’entropie comme dans certains arbres de décision. Dans le cas binaire Y = 0/1, elle vaut : •F , (O, ')G = − E¥' logF (O )G + (1 − ' ) logF1 − (O )G¦.

Dans le cas multi-classes, elle vaut :

•F , (O, ')G = − E E ' @ logF @ (O )G , @_;

yik étant la variable à expliquer 0/1 associée à la kème classe pour la ième observation. En posant f(x) = P(Y=1|X=x), on reconnaît dans cette expression la moitié de la déviance du modèle linéaire généralisé. La minimisation de l’erreur revient donc à la recherche du maximum de vraisemblance. Avec la fonction d’activation sigmoïde, on reconnaît d’ailleurs dans le perceptron une régression logistique sur les unités de la couche cachée. Le maximum de vraisemblance donnera les poids cherchés pour ces unités. Pour la recherche du minimum de l’erreur, les algorithmes d’optimisation sont le plus souvent la descente du gradient (conjugué jusqu’à quelques dizaines de milliers d’observations, et stochastique au-delà) ou la méthode de quasi-Newton. La méthode de Newton elle-même, avec son calcul de l’inverse de la matrice hessienne, est instable et c’est pourquoi on lui préfère la méthode de quasi-Newton, qui approxime l’inverse de la matrice hessienne de façon itérative sans le calculer explicitement. Pour les réseaux de neurones, on parle d’algorithmes de rétropropagation. Par nature, ils se prêtent à la parallélisation des calculs, puisque les échanges d’information ne se font pas globalement mais seulement entre unités connectées. L’apprentissage peut être en batch (la mise à jour des poids n’est effectuée qu’après traitement de toutes les observations) ou on-line (à chaque nouvelle observation traitée). Ce mécanisme d’apprentissage impose quelques contraintes. D’abord sur les poids initiaux : plus ils sont proches de 0, plus le modèle est proche d’un modèle linéaire. Des poids plus éloignés de 0 prennent en compte les effets nonlinéaires, ce qui est intéressant, mais s’ils sont trop grands, ils saturent le réseau qui renvoie alors toujours des valeurs proches. On considère que les poids initiaux doivent être initialisés selon la loi uniforme dans un intervalle autour de 0, et l’on fixe parfois §−o6¨(( + 1) ; o6¨(( + 1)ª, où p est le nombre d’unités de la couche précédente.