:encoding: ISO-8859-1 Visualisation et classification multidimensionnelles ==================================================== Christophe Lalanne Mars 2008 .Avertissement ****************************************************************************** Ce document est supposé couvrir un ensemble cohérent de techniques de visualisation de données multidimensionnelles. En l'état actuel, il s'agit surtout d'un *document de travail*, peu structuré et avec de nombreuses erreurs à corriger. Il est conçu de manière itérative, au fil de mes lectures et de mes expériences dans ce domaine. Tant que ce message sera visible, il convient donc de rester prudent quant aux éventuelles anomalies ou erreurs, et surtout de pardonner à l'auteur le désordre intime de cette page html. L'historique des modifications de ce document se trouve dans le fichier link:History[History]. À terme, ce document devrait être disponible en version html, pdf et docbook. *La version consultée a été mise à jour en septembre 2008.* ****************************************************************************** [[introduction,Introduction]] Introduction ------------ La visualisation de données multidimensionnelles constitue une étape préalable importante à des techniques de modélisation plus classiques (régression, analyse discriminante, etc.). Selon la nature et le type de variables, certains graphiques sont naturellement plus appropriés que d'autres. Ils sont généralement présentés dans la plupart des cours de statistiques. On pourra consulter, par exemple, http://www.lsp.ups-tlse.fr/Besse/pub/sdm2.pdf[Statistique descriptive multidimensionnelle] (A. Baccini et P. Besse, Laboratoire de Statistique et Probabilités). Toutefois, la visualisation ne se limite pas au cas bi- ou tri-varié, mais s'étend également à la représentation de données plus complexes. Pour s'en convaincre, le lecteur peut consulter les différentes http://en.wikipedia.org/wiki/Visualization_(graphic)[rubriques de Wikipedia]. Nous nous intéresserons principalement aux techniques d'analyses factorielles et de classification, ainsi qu'au moyen de représenter des matrices de données de dimensions élevées (données quantitatives, qualitatives, distances, etc.) dans des espaces à 2 ou 3 dimensions. Avant toute chose, il est intéressant de remarquer que le domaine de la classification (_clustering_ dans la littérature anglo-saxonne) partage des liens étroits avec celui de l'apprentissage automatique (_machine learning_). De ce point de vue, on peut alors distinguer deux types de classes de méthodes, selon qu'elles relèvent plutôt de l'apprentissage supervisé ou non supervisé. Évidemment, une autre typologie s'impose d'elle-même lorsque l'on songe simplement au type et au statut des variables considérées. La première distinction rejoint le problème de la métrique éventuellement associée aux variables d'intérêt : sont-elles de nature qualitative (nominale, ordinale, binaire) ou quantitative (échelle d'intervalle ou de rapport) ? D'autre part, le rôle des variables est-il symétrique ou, au contraire, certaines variables sont-elles explicatives (ou prédictrices) et d'autres dépendantes (variables réponses) ? Le tableau ci-dessous adresse ces différentes questions en regroupant certaines techniques d'analyse multivariée classiques. .Aperçu des techniques d'analyses multivariées [grid="all"] .------------------------.-----------------------------.----------------------- Méthode Variables prédictrices Variable réponse ------------------------------------------------------------------------------- Régression linéaire p var. quant. ou qual. var. quant. ANOVA p var. qual. var. quant. ANCOVA 1 var. qual. / 1 var. qual. var. quant. Régression logistique p var. quant. ou qual. var. bin. ou ord. Analyse discriminante p var. quant. / 1 var. qual. CART p var. quant. ou qual. var. bin. ------------------------------------------------------------------------------- D'autre part, lorsque l'on considère une matrice n x p, on peut chercher à expliquer les corrélations entre les variables (colonnes) ou chercher les ressemblances, ou plutôt les dissimilarités, entre les individus (lignes). Cela amène bien évidemment à une nouvelle distinction possible entre les techniques évoquées ci-dessus. L'ensemble des exemples présentés dans cet article sont analysés avec R version 2.6.2 (2008-02-08). Les principaux ``packages'' utilisés sont les suivants : - +MASS+ (à présent, intégré sous forme de bundle avec +class+, +nnet+ et +spatial+) : Functions and datasets to support Venables and Ripley, 'Modern Applied Statistics with S' (4th edition). - +tree+ : Classification and Regression Trees. - +Vegan+ : Ordination methods, diversity analysis and other functions for community and vegetation ecologists. - +clv+ : Package contains most of the popular internal and external cluster validation methods ready to use for the most of the outputs produced by functions coming from package "cluster". Package contains also functions and examples of usage for cluster stability approach applied to algorithms implemented in "cluster" package. - +flexclust+ : The main function kcca implements a general framework for k-centroids cluster analysis supporting arbitrary distance/similarity measures and centroid computation. Further cluster methods include hard competitive learning, neural gas and QT clustering. - +cluster+ : Cluster Analysis, extended original from Peter Rousseeuw, Anja Struyf and Mia Hubert. - +mclust+ : Model-based clustering and normal mixture modeling including Bayesian regularization - +LLAhclust+ : The likelihood linkage analysis is a general agglomerative hierarchical clustering method developed in France by Lerman in a long series of research articles and books. Initially proposed in the framework of variable clustering, it has been progressively extended to allow the clustering of very general object descriptions. The approach mainly consists in replacing the value of the estimated similarity coefficient by the probability of finding a lower value under the hypothesis of 'absence of link'. The package LLAhclust contains routines for computing various types of probablistic similarity coefficients between variables or object descriptions. Once the similarity values between variables/objects are computed, a hierarchical clustering can be performed using several probabilistic and non-probabilistic aggregation criteria, and indices measuring the quality of the partitions compatible with the resulting hierarchy can be computed. - +ade4+ : Multivariate data analysis and graphical display. Il peut être intéressant de jeter un oeil aux ``task views'' sur le site CRAN, notamment *Multivariate* et *MachineLeraning*, pour avoir une vue cohérente des différents ``packages''. Nous ne traitons pas de la modélisation de données, de type de régression ou modèles mixtes, et n'utilisons pas le package ``lattice''. Pour le lecteur intéressé par les graphiques en treillis, il existe une abondante documentation sur le site http://netlib.bell-labs.com/cm/ms/departments/sia/project/trellis/[Treillis displays], ainsi qu'un ouvrage à paraître chez Springer : http://www.springer.com/statistics/computational/book/978-0-387-75968-5[Lattice. Multivariate Data Visualization with R] (D. Sarkar, approx. 290 p.). En revanche, le logiciel http://www.ggobi.org/[GGobi] est intéressant pour visualiser des jeux de données multivariées. Un ouvrage a d'ailleurs été publié récemment aux éditions Springer (toujours !) : http://www.springer.com/statistics/statistical+theory+and+methods/book/978-0-387-71761-6[Interactive and Dynamic Graphics for Data Analysis. With R and GGobi] (D. Cook and D. F. Swayne, 2007). Des http://www.ggobi.org/book/index.html[ressources] sont disponibles en ligne. L'intérêt de ce logiciel est, d'une part, d'être multiplateforme (Mac, Linux, Windows), et, d'autre part, d'offrir la possibilité de visualiser dynamiquement les données (exemple de http://www.ggobi.org/book/chap-toolbox/toolbox-PP2D.mov[poursuite projective en 2D], au format +mov+). Entrée en matière ----------------- Quelques rappels d'algèbre ~~~~~~~~~~~~~~~~~~~~~~~~~~ Le cas de trois variables ~~~~~~~~~~~~~~~~~~~~~~~~~ Cette section est (très) largement inspirée du tutoriel http://pbil.univ-lyon1.fr/R/fichestd/ter1.pdf[Représentations triangulaires] (D. Chessel, A. B. Dufour et J. R. Lobry, 2005). Les méthodes d'analyse factorielle ---------------------------------- .Remarque ****************************************************************************** Nous n'abordons les techniques d'analyses factorielles que dans le contexte d'une analyse exploratoire des données. Le but de ces dernières est souvent de proposer une interprétation ou un modèle explicatif en relation avec les données observées. Le lecteur intéressé par ces aspects trouvera toutes les ressources nécessaires sur le site du http://pbil.univ-lyon1.fr/R/enseignement.html[Laboratoire de Biométrie] de Lyon. ****************************************************************************** L'Analyse en Composantes Principales ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ L'objet de l'Analyse en Composantes Principales (ACP) est de réduire, si possible, l'ensemble initial des p variables à un ensemble de k variables (k < p) expliquant une part significative de la variance originale et non-corrélées entre elles. On retrouvera également dans la littérature le terme d'_ordination en espace réduit_ puisque l'objectif principal de l'ACP est de ``projeter'' les observations dans un espace de faible dimension tout en conservant le maximum de l'information originale, et cette information n'est autre que la variabilité des mesures observées sur chaque variable. Dans le contexte de l'ACP, mais également des techniques d'analyse de données en général, on désignera par _individu_ une unité statistique sur laquelle on été recueillie diverses mesures, et _variable_ les différents supports de mesure. On se donne donc pour commencer une matrice n x p (n individus, p variables), et l'on cherche à réduire l'espace à p dimensions en un espace à k dimensions permettant de reconstruire les relations entre les p variables de manière aussi précises que possible. // TODO // détailler aspects mathématiqes + schéma Lorsqu'aucune inférence n'est effectuée sur les paramètres de ce modèle, il n'y a aucune condition de validité. Le cas échéant (e.g. estimation des intervalles de confiance associés aux coefficients), on pose l'hypothèse, à vérifier, d'une distribution multinormale des données. On a dit en <> que cette technique s'applique préférentiellement sur un ensemble de variables continues. À cela, deux exceptions peuventêtre mentionnées : - il est tout à fait possible d'appliquer une ACP à des variables ordinales ou binaires, - des variables catégorielles supplémentaires peuvent être utilisées en même temps qu'un ensemble de variables quantitatives lors de la représentation graphique des relations entre variables ou de l'espace des individus. Une bonne référence concernant ce sujet est <>, mais on trouvera dans <> une présentation de l'ACP en lien avec l'<> ou l'<>. L'ACP est si fondamentale en analyse de données et dans la plupart des étapes préalables d'exploration et de synthèse numérique des données qu'il est difficile de ne pas y consacrer une large partie de ce document. On peut mentionner deux raisons à l'importance de l'ACP dans le domaine de l'analyse de données. D'une part, on y retrouve l'essence même de la plupart des méthodes statistiques explicatives, à savoir expliquer la variabilité observée dans une ou plusieurs variables réponses, en fonction d'autres variables. D'autre part, le schéma de l'ACP et ses propriétés mathématiques permettent de dériver la plupart des autres modèles d'analyse factorielle. On notera toutefois que cet aspect est plutôt lié à la tradition française de l'analyse des données car les chercheurs anglo-saxons limitent généraement l'usage de l'ACP à la construction de l'espace des variables, négligeant souvent la représentation des individus dans l'espace réduit. Les fonctions R permettant de réaliser ce type d'analyse sont +prcomp+ et +princomp+. La différence essentielle entre ces deux fonctions tient au fait que +princomp+ ne permet pas d'analyser des jeux de données où p > n, alors que +prcomp+ repose sur une décomposition en valeurs singulières. footnote:[La décomposition en valeurs singulières (voir +svd+ sous R) est en fait la généralisation mathématique de l'ACP.] On préférera donc cette dernière solution par souci de précision numérique ou dans le cas de matrices déséquilibrées (plus de variables que d'individus). Dernière remarque, on peut utiliser cette technique pour réduire l'espace des variables explicatives et utiliser ensuite une régression linéaire (ou non linéaire) sur les composantes principales : on parle ainsi de régression sur composantes principales. Intéressante par la simplicité de sa mise en oeuvre, cette technique nécessite de la part de l'utilisateur de prendre en considération le fait que les conditions d'application (hypothèse de multi-normalité, hypothèses classiques associées au modèle linéaire) ainsi que l'interpréation des coefficients de régression s'en trouve sensiblement compliquées. [[analyse factorielle des correspondances]] L'Analyse Factorielle des Correspondances ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ [[analyse des correspondances multiples]] L'Analyse des Correspondances Multiples ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ L'Échelonnement Multi-Dimensionnel ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Encore appelé Analyse de Co-Inertie, ou Multdimensional Scaling dans la littérature anglo-saxonne, l'échelonnement multi-dimensionnel vise à positionner des individus dans un espace dont la métrique est liée aux variables mesurées. Souvent, celles-ci sont des mesures de proximité, de distance, de préférence. L'Analyse Factorielle Exploratoire ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Contrairement à l'analyse en composantes principales dans laquelle on s'intéresse aux relations existant entre différentes variables quantitatives (combinaisons linéaires), l'analyse factorielle permet d'analyser des variables non directement observables. On parle alors de _variables latentes_ (e.g. intelligence, compétence linguistique). On suppose généralement que celles-ci sont ``liées'' à des variables qui sont, elles, observables ou mesurables. Ces dernières sont appelées _variables manifestes_. Il existe de très bonnes introductions à ce type de modèle dans <> et <>. Nous ne nous intéresserons dans cette section qu'à l'analyse factorielle exploratoire, en laissant de côté l'approche confirmatoire souvent associée aux modèles d'équations structurelles. Sur ce sujet, on pourra également consulter <> ou <>. Le principe général des modèles factoriels est similaire à celui appliqué en régression linéaire multiple, sauf que dans ce cadre on ne peut pas estimer directement les coefficients de régression. Dans l'analyse factorielle, les variables manifestes sont régressées sur les variables latentes, appelées également _facteurs communs_. L'idée est que les corrélations (ou covariances) observées entre les variables manifestes résultent de leurs liens avec la ou les les variables latente(s). L'un des ouvrages de référence sur le sujet est <>. En fait, l'idée de base repose sur le modèle en facteur commun de Thurstone <> qui postule que chaque indicateur (variable manifeste) dans un ensemble de données observées est une fonction linéaire d'un ou plusieurs facteurs _communs_ et d'un facteur _unique_. L'analyse factorielle revient à partitionner la variance de cet indicateur, à partir de la matrice de variance/covariance, en deux parts : - la _variance commune_, ou variance ``expliquée'' par le facteur latent, qui est estimée à partir de la variance partagée avec les autres indicateurs ; - la _variance spécifique_ (ou _unique_), qui est une combinaison de la variance réelle spécifique de l'indicateur et d'une variance résiduelle assimilable à un terme d'erreur (par analogie avec les modèles de régression). Ce dernier terme est également appelé erreur de mesure ou imprécision de l'indicateur. Sous R, on dispose de la fonction +factanal()+ qui permet d'effectuer une analyse factorielle de base. L'estimation des paramètres se fait par le _maximum de vraisemblance_. footnote:[La fonction +factanal()+ de S-PLUS offre quant à elle le choix entre MV et ACP.] SAS, avec +PROC FACTOR+, propose deux types d'estimation : par MV et sur composantes principales. SPSS propose également un ensemble de possibilités pour l'extraction des facteurs. Celles-ci sont accessibles dans le menu *Data Reduction*. Par défaut, la méthode d'extraction des facteurs est effectuée sur la base des composantes principales, mais les autres méthodes incluent : moindres carrés ordinaires, moindres carrés pondérés, maximum de vraisemblance, factorisation de l'axe principal, factorisation alpha, factorisation de l'image. Les options associées à chacune des méthodes d'extraction, ainsi que la facilité d'utilisation, explique sans doute sa popularité auprès des psychologues et des chercheurs en sciences sociales. image::spss_factor.png[] .Remarque ****************************************************************************** Depuis la rédaction de cette section, il existe une fonction R réalisant l'extraction des facteurs sur la base des composantes principales. Il s'agit de la fonction +principal+ dans le package +psych+. ****************************************************************************** Voici un petit exemple de la mise en oeuvre de ces principes (tiré de ???). ------------------------------------------------------------------------------ > # A little demonstration, v2 is just v1 with noise, > # and same for v4 vs. v3 and v6 vs. v5 > # Last four cases are there to add noise > # and introduce a positive manifold (g factor) > v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6) > v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5) > v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6) > v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4) > v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5) > v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4) > m1 <- cbind(v1,v2,v3,v4,v5,v6) > cor(m1) v1 v2 v3 v4 v5 v6 v1 1.0000000 0.9393083 0.5128866 0.4320310 0.4664948 0.4086076 v2 0.9393083 1.0000000 0.4124441 0.4084281 0.4363925 0.4326113 v3 0.5128866 0.4124441 1.0000000 0.8770750 0.5128866 0.4320310 v4 0.4320310 0.4084281 0.8770750 1.0000000 0.4320310 0.4323259 v5 0.4664948 0.4363925 0.5128866 0.4320310 1.0000000 0.9473451 v6 0.4086076 0.4326113 0.4320310 0.4323259 0.9473451 1.0000000 > factanal(m1, factors=3) ------------------------------------------------------------------------------ La méthode de rotation par défaut est +varimax+, et donne les résultats suivants : ------------------------------------------------------------------------------ Call: factanal(x = m1, factors = 3) Uniquenesses: v1 v2 v3 v4 v5 v6 0.005 0.101 0.005 0.224 0.084 0.005 Loadings: Factor1 Factor2 Factor3 v1 0.944 0.182 0.267 v2 0.905 0.235 0.159 v3 0.236 0.210 0.946 v4 0.180 0.242 0.828 v5 0.242 0.881 0.286 v6 0.193 0.959 0.196 Factor1 Factor2 Factor3 SS loadings 1.893 1.886 1.797 Proportion Var 0.316 0.314 0.300 Cumulative Var 0.316 0.630 0.929 The degrees of freedom for the model is 0 and the fit was 0.4755 ------------------------------------------------------------------------------ On peut également utiliser la méthode +promax+, en modifiant l'option +rotation+: +factanal(m1, factors=3, rotation="promax")+. On obtient alors : ------------------------------------------------------------------------------ Call: factanal(x = m1, factors = 3, rotation = "promax") Uniquenesses: v1 v2 v3 v4 v5 v6 0.005 0.101 0.005 0.224 0.084 0.005 Loadings: Factor1 Factor2 Factor3 v1 0.985 v2 0.951 v3 1.003 v4 0.867 v5 0.910 v6 1.033 Factor1 Factor2 Factor3 SS loadings 1.903 1.876 1.772 Proportion Var 0.317 0.313 0.295 Cumulative Var 0.317 0.630 0.925 The degrees of freedom for the model is 0 and the fit was 0.4755 ------------------------------------------------------------------------------ Enfin, on utilise +prcomp()+ pour afficher les composantes : ------------------------------------------------------------------------------ > # The following shows the g factor as PC1 > prcomp(m1) Standard deviations: [1] 3.0368683 1.6313757 1.5818857 0.6344131 0.3190765 0.2649086 Rotation: PC1 PC2 PC3 PC4 PC5 PC6 v1 0.4168038 -0.52292304 0.2354298 -0.2686501 0.5157193 -0.39907358 v2 0.3885610 -0.50887673 0.2985906 0.3060519 -0.5061522 0.38865228 v3 0.4182779 0.01521834 -0.5555132 -0.5686880 -0.4308467 -0.08474731 v4 0.3943646 0.02184360 -0.5986150 0.5922259 0.3558110 0.09124977 v5 0.4254013 0.47017231 0.2923345 -0.2789775 0.3060409 0.58397162 v6 0.4047824 0.49580764 0.3209708 0.2866938 -0.2682391 -0.57719858 ------------------------------------------------------------------------------ On peut visualiser le diagramme des valeurs propres et projeter les variables dans le plan factoriel (1,2) image::afe_scree.png[] image::afe_biplot.png[] // est-ce interprétable comme pour l'ACP ? Si l'on poursuit l'analyse, on peut afficher les facteurs à l'aide de la commande : ------------------------------------------------------------------------------ > (m1.scores <- factanal(~v1+v2+v3+v4+v5+v6,factors=3,scores="Bartlett")\$scores) Factor1 Factor2 Factor3 1 -0.9039949 -0.9308984 0.9475392 2 -0.8685952 -0.9328721 0.9352330 3 -0.9082818 -0.9320093 0.9616422 4 -1.0021975 -0.2529689 0.8178552 5 -0.9039949 -0.9308984 0.9475392 6 -0.7452711 0.7273960 -0.7884733 7 -0.7098714 0.7254223 -0.8007795 8 -0.7495580 0.7262851 -0.7743704 9 -0.8080740 1.4033517 -0.9304636 10 -0.7452711 0.7273960 -0.7884733 11 0.9272282 -0.9307506 -0.8371538 12 0.9626279 -0.9327243 -0.8494600 13 0.9229413 -0.9318615 -0.8230509 14 0.8290256 -0.2528211 -0.9668378 15 0.9272282 -0.9307506 -0.8371538 16 0.4224366 2.0453079 1.2864761 17 1.4713902 1.2947716 0.5451562 18 1.8822320 0.3086244 1.9547752 ------------------------------------------------------------------------------ On peut visualiser les scores à l'aide de différentes représentations, comme le montrent les figures ci-après. // pas forcément le plus adapté... image::afe_contour.png[] image::afe_levelplot.png[] On peut voir l'AFE comme une certaine forme de classification des variables, en groupes homogènes du point de vue de la mesure utilisée. On réduit ainsi le nombre de variables initiales à un nombre restreint (généralement, 2 à 5 variables) de facteurs exprimés sous forme de combinaison linéaire des variables initiales. Attention toutefois ! Contrairement à l'ACP, les relations entre les variables au sein de chaque combinaison linéaire n'ont pas la même interprétation. Les méthodes de classification ------------------------------ Tout d'abord, il convient de remarquer que les techniques de classification renvoient généralement à deux domaines, qui ne se recouvrent pas entièrement : le domaine de l'intelligence artificielle et de l'apprentissage automatique et les statistiques exploratoires multidimensionnelles. À la première catégorie, on peut rattacher des ouvrages comme <> ou <> (je cite des ouvrages que j'ai sous les yeux mais il en existe _vraiment_ beaucoup d'autres), tandis que l'école de Lebart et coll. <> se concentre plutôt sur la deuxième approche. Attention, cet ouvrage est à certains égards assez technique. Le lecteur débutant préférera de loin <> dont l'approche est beaucoup plus pédagogique. Entre les deux approches, on trouve un bon compromis avec le livre de <>. Les techniques de classification, ou l'analyse en cluster, visent à identifier des groupes homogènes (dans un sens à définir plus loin) parmi les n observations décrites par p variables. Il existe plusiseurs façons d'aborder le problème de la classification. On distinguera essentiellement des _méthodes agglomératives_ (e.g. CAH) ou _divisives_, des _méthodes de partitionnement optimal_ (nuées dynamiques), ainsi que des méthodes reposant sur la maximisation de la vraisemblance <>, <>, <>. Ce dernier type de méthode évite les inconvénients liés aux autres techniques, à savoir une évaluation heuristique de la solution optimale, puisqu'elle utilise des critères statistiques (BIC) pour décider du nombre de clusters à retenir. En revanche, elle suppose un modèle multivarié sous-jacent. Pour les méthodes agglomératives, les fonctions +hclust+, +agnes \{cluster\}+, +hcluster \{amap\}+ permettent de réaliser les analyses de base. Pour les méthodes de partitionnement dynamique, on dispose de la fonction +kmeans+, mais là encore le ``package'' +cluster+ fournit des outils de base additionnels : +pam+, +clara+ et +fanny+. Quant aux méthodes divisives, on dispose, toujours dans le package +cluster+, des fonctions +diana+ et +mona+. L'article de <> décrit la plupart de ces fonctions. Concernant la procédure de classification par MLE, le ``package'' +mclust+ fournit les outils nécessaires. On notera que l'utilisation de ces fonctions nécessite un http://www.stat.washington.edu/fraley/mclust/license.txt[accord de licence]. Par ailleurs, si dans le cadre des techniques de classification présentées plus haut, on se situait résolument du côté d'une ``approche par les variables'', les techniques de classification présentées dans les paragaphes suivants s'orientent plutôt sur une classification par les individus. En fin de compte, on peut voir ces deux classes de techniques dans une optique de dualité. Parfois, on s'intéresse plutôt à regrouper des variables ensemble car elles sous-tendent une même dimension -- c'est le cas de l'analyse factorielle exploratoire --, alors que d'autre fois, ce sont plutôt des groupes homogènes d'individus ou d'observations que l'on souhaite constituer. Le choix est laissé à la libre appréciation de l'analyste, mais il faut bien garder à l'esprit que la portée interprétative dépend fortement du type d'approche retenu. Par exemple, dans le domaine biomédical où les classifications sont assez souvent utilisées, la question du choix de l'approche se pose légitimement : - en effectuant une classification par les variables, on décide de constituer des ensembles de variables (e.g. des items ou des scores obtenus sur plusieurs instruments de mesure), qui peuvent être des variables manifestes de _symptômes_ particuliers ; une collection de telles variables est donc assimilable à un _syndrôme_. - en effectuant une classification par les individus, on cherche à partitionner la population (à partir de l'échantillon) en groupes homogènes présentant une typologie de réponse spécifique. Classification hiérarchique ~~~~~~~~~~~~~~~~~~~~~~~~~~~ Les méthodes agglomératives consistent à regrouper, de manière itérative, les individus en paquets. On commence par exemple par n clusters (sachant qu'il y a n individus au total), puis on cherche à partitionner cet ensemble en n-1 clusters, etc. Différents critères d'optimisation sont utilisés (généralement la maximisation de la variance inter-cluster et la minimisation de la variance intra-cluster). // TODO // start here... image::classification_fond_75pc.png[] Le principe général des méthodes agglomératives et des nuées dynamiques est illustré dans la figure ci-dessus. Dans la classification hiérarchique ascendante, les principales étapes se résument à 1. sélectionner les deux unités dont la distance est minimale 2. calculer le nouveau point résultant de la fusion de ces deux unités 3. itérer avec les n-1 unités restantes Le dendrogramme traduit l'évolution de la procédure, l'axe vertical représentant la mesure de distance inter-unités retenue. Pour la procédure k-means, on part de l'ensemble des observations avec un nombre de clusters déterminé a priori, puis 1. on sélectionne aléatoirement 2 unités et on calcule les distances entre les autres unités et ces deux unités (ou ``graine'') 2. on assimile dans un même cluster les unités qui minimisent la variance intra-cluster (distance entre unités d'un même cluster) et maximisent la variance inter-cluster (distance entre unités appartenant à des clusters différents) 3. on calcule les centroïdes de chaque cluster 4. itérer les étapes 2 à 4 en considérant comme ``graines'' les points moyens calculés à l'étape 3 Avant d'aborder la classification hiérarchique ascendante (CAH), il est nécessaire de décrire les différents critères d'agglomération utilisés. Tout repose généralement sur une mesure de dissimilarité, analogue à la notion de distance au sens mathématique : on impose d(A,B) = d(B,A) >=0 et $$`d(A,A)!=0`$$. Comme premier exemple, nous utiliserons les données de B. Everitt <>, contenues dans le dataframe +life+. Ces données portent sur l'espérance de vie de différents pays, stratifiés par âge et par sexe. Les différentes méthodes -- "ward", "single", "complete", "average", "mcquitty", "median" et "centroid" -- sont utilisées succesivement afin de comparer les résultats auxquels elles permettent d'aboutir. Commençons par la méthode de Ward : ------------------------------------------------------------------------------ > source("chap4lifeexp.dat") > head(life) m0 m25 m50 m75 w0 w25 w50 w75 Algeria 63 51 30 13 67 54 34 15 Cameroon 34 29 13 5 38 32 17 6 Madagascar 38 30 17 7 38 34 20 7 Mauritius 59 42 20 6 64 46 25 8 Reunion 56 38 18 7 62 46 25 10 Seychelles 62 44 24 7 69 50 28 14 > life.hc <- hclust(dist(life), "ward") > plot(life.hc) ------------------------------------------------------------------------------ [[hc_plot1]] .Dendogramme obtenu avec la méthode de Ward. image::hc1.png[hc1.png] Pour identifier le nombre de groupes optimal, la première méthode est visuelle et consiste à examiner le dendrogramme afin de déterminer, sur la base d'arguments ou de critères sémantiques, un nombre de classes _a priori_ acceptable. La fonction +cutree+ permet de ``couper'' le dendogramme soit à une hauteur spécifiée (+h+), soit sur la base d'un nombre déterminé de groupes (+k+). Supposons que nous choisissons de retenir 5 groupes d'individus. ------------------------------------------------------------------------------ > life.hc.grp <- cutree(life.hc, k=5) > country <- rownames(life) > lapply(1:5, function(x) country[life.hc.grp==x]) [[1]] [1] "Algeria" "Tunisia" "Costa Rica" "Dominican Rep" [5] "Nicaragua" "Panama" "Trinidad(62)" [[2]] [1] "Cameroon" "Madagascar" [[3]] [1] "Mauritius" "Reunion" "Seychelles" [4] "El Salvador" "Greenland" "Grenada" [7] "Honduras" "Jamaica" "Mexico" [10] "Trinidad (67)" "United States (NW66)" "Chile" [13] "Columbia" "Ecuador" [[4]] [1] "South Africa(C)" "Guatemala" [[5]] [1] "South Africa(W)" "Canada" "United States (66)" [4] "United States (W66)" "United States (67)" "Argentina" ------------------------------------------------------------------------------ Si l'on considère les autres critères, on peut comparer les solutions obtenues. ------------------------------------------------------------------------------ > meth <- c("ward","single","complete","average","mcquitty","median","centroid") > res <- matrix(NA,nr=nrow(life),nc=length(meth)) > for (i in 1:length(meth)) + res[,i] <- cutree(hclust(dist(life), meth[i]), k=5) > res.tab <- apply(res,2,function(x) sort(table(x))) > dimnames(res.tab) <- list(1:5,meth) > res.tab ward single complete average mcquitty median centroid 1 2 1 2 1 1 1 1 2 2 1 2 2 2 1 1 3 6 2 7 2 2 2 2 4 7 2 8 2 13 6 2 5 14 25 12 24 13 21 25 > mean(res.tab[5,])/nrow(life) [1] 0.6175115 ------------------------------------------------------------------------------ Les méthodes *single* et *centroid* fournissent exactement les mêmes solutions. On retrouve toujours un cluster dominant, avec en moyenne 62 % des individus. Idéalement, on voudrait représenter les n(n-1)/2 matrices de confusion correspondantes. En voici deux exemples, sous forme graphique : .----------------------------------------.------------------------------------ image::hc2.png[hc2.png] image::hc3.png[hc3.png] ------------------------------------------------------------------------------ Il existe naturellement d'autres critères que l'examen du dendrograme. Mentionnons par exemple : - les indices t^2^, pseudo-F, ccc - l'indice gamma de Hubert [TIP] ============================================================================== Du côté logiciel, SAS et SPSS proposent chacun des procédures de classification incluant les techniques décrites ci-dessus (à l'exception des techniques de visualisation dynamiques...). Sous SAS, la proédure est +PROC CLUSTER+. Associée à +PROC DISTANCE+ permettant de créer une matrice de distance selon méthode voulue, elle permet également de représenter le dendrograme lorsqu'elle est couplée à +PROC TREE+. On peut voir un exemple d'enchaînement de ces trois procédures dans l'aide en ligne. Voici le script : link:consommation.sas[consommation.sas]. Sous SPSS (v. 16), les procédures de classification sont regroupées dans le menu *Classify*. En particulier, la classification ascendante hiérarchique s'effecteura en sélectionnant *Classify > Hierarchical Cluster*. Ce http://www2.chass.ncsu.edu/garson/PA765/cluster.htm[tutoriel] est assez bien documenté. Il existe également un article détaillant les règles appliquées par SPSS pour gérer les ``ex-aequos'' : http://www.leidenuniv.nl/fsw/mtlab/software/SPSS%20rules.pdf[The rules of SPSS' Hierarchical Cluster Analysis for processing ties]. image::spss_cluster.png[spss_cluster.png] ============================================================================== Nuées dynamiques ~~~~~~~~~~~~~~~~ Sur le plan informatique, il faut bien garder à l'esprit que ce type de méthode repose sur la détermination du nombre possible de partitions d'un ensemble à n éléments. On peut montrer (http://www.lsp.ups-tlse.fr/Carlier/.Hyper/polyclass/node7.html[Cours de classification automatique de A. Carlier, par exemple]) que le nombre le nombre total de partition, N~n~, augmente plus qu'exponentiellement en fonction de n ; il vaut environ 5.10^13^ pour n=20. Le nombre de partitions contenant exactement k parties d'un ensemble E~n~ à n éléments vaut en fait image::eq_partitions.png[eq_partitions.png] À titre d'illustration, le programme http://www.aliquote.org[+partitions.c+] énumère le nombre de partitions pour un ensemble donné d'éléments de l'ensemble des entiers naturels. Un exemple d'algorithme pour la procédure ``k-means'' est fournit ci-dessous. Le script est programmé avec Octave. .kmeans2.m ------------------------------------------------------------------------------ function out = kmeans2(m,k) %% kmeans.m %% a quick implementation of k-means algorithm %% %% input: %% - m, data matrix %% - k, number of clusters %% %% output: %% - [m,g], data matrix plus group membership (1..k) %% [nr,nc] = size(m) p = randperm(size(m,1)); % random initialization for i=1:k c(i,:)=m(p(i),:) end tmp = zeros(nr,1) while 1, d = dist(m,c); % object centroïd distances [z,g] = min(d,[],2); % group matrix if g == tmp, break; else tmp = g; end for i=1:k f = find(g == i) if f % compute centroïd if not empty set c(i,:) = mean(m(find(g == i),:),1); end end end out = [m,g] ------------------------------------------------------------------------------ Il existe naturellement des algorithmes plus puissants, et surtout un peu mieux optimisés, que celui-ci, e.g. http://www.cs.ucsd.edu/~elkan/fastkmeans.html[Fast k-means code for Matlab] ou http://www.inference.phy.cam.ac.uk/itprnn/code/kmeans/[K-mean clustering and Soft K-means clustering ] (site de D. MacKay, auteur du fameux livre _Information Theory, Inference, and Learning Algorithms_) Utilisation des centroïdes de classe ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ On peut comparer les résultats précédents avec ceux obtenus par une méthode ``k-means''. Dans la pratique, il est plus courant de débuter avec une CAH pour isoler un nombre de clusters qui servira à initialiser l'algorithme de sélection/fusion de la procédure ``k-means'' (section suivante). ------------------------------------------------------------------------------ > set.seed(134) > life.km <- kmeans(life, 5) > life.km K-means clustering with 5 clusters of sizes 2, 13, 7, 7, 2 Cluster means: m0 m25 m50 m75 w0 w25 w50 w75 1 36.00000 29.50000 15.00000 6.000000 38.00000 33.00000 18.50000 6.500000 2 59.38462 43.00000 22.53846 7.923077 63.92308 46.84615 25.46154 9.615385 3 66.14286 45.28571 23.28571 7.857143 72.85714 51.14286 28.14286 10.428571 4 63.14286 50.57143 26.14286 10.571429 66.71429 50.85714 29.14286 12.428571 5 49.50000 39.50000 21.00000 8.000000 53.00000 42.00000 23.00000 8.000000 Clustering vector: Algeria Cameroon Madagascar 4 1 1 Mauritius Reunion Seychelles 2 2 3 South Africa(C) South Africa(W) Tunisia 5 3 4 Canada Costa Rica Dominican Rep 3 4 4 El Salvador Greenland Grenada 2 2 2 Guatemala Honduras Jamaica 5 2 2 Mexico Nicaragua Panama 2 4 4 Trinidad(62) Trinidad (67) United States (66) 4 2 3 United States (NW66) United States (W66) United States (67) 2 3 3 Argentina Chile Columbia 3 2 2 Ecuador 2 Within cluster sum of squares by cluster: [1] 25.5000 410.1538 102.8571 546.0000 15.0000 Available components: [1] "cluster" "centers" "withinss" "size" > plot(life, col=life.km\$cluster) ------------------------------------------------------------------------------ [[km_plot1]] .Matrice de dispersion pour les données +life+ avec les solutions ``k-means''. image::km1.png[km1.png] On peut donc directement comparer les solutions issues des deux méthodes : ------------------------------------------------------------------------------ > table(as.numeric(life.hc.grp),life.km\$cluster) ------------------------------------------------------------------------------ ce qui donne la matrice de confusion suivante : .............................................................................. 1 2 3 4 5 1 0 0 0 7 0 2 2 0 0 0 0 3 0 13 1 0 0 4 0 0 0 0 2 5 0 0 6 0 0 .............................................................................. On voit clairement que les solutions diffèrent entre les deux méthodes : si l'on retrouve bien un groupe dominant (13 individus), dans les deux cas il existe une incertitude sur le classement d'un individu (ligne 3, colonne 3). Pour évaluer la qualité de la solution proposée, on peut vérifier comment évolue la SS intra-groupe pour différentes solutions (k=2..10). ------------------------------------------------------------------------------ > kmax <- 10 > ssw <- NULL > for (i in 2:kmax) + ssw <- append(ssw, sum(kmeans(life, i)\$withinss)) > ssw.1 <- (nrow(life)-1)*sum(apply(life, 2, var)) > ssw <- c(ssw.1,ssw) > plot(1:kmax, ssw, type="l", xlab="Group #") > abline(h=ssw[4], lty=2) ------------------------------------------------------------------------------ [[km_plot2]] .Évolution de la SS intra-groupe en fonction du nombre de clusters. image::km2.png[km2.png] À la lecture de graphique, il apparaît qu'une solution à 4 clusters semble suffisante. Il est tout à fait possible d'associer la procédure ``k-means'' à une analyse exploratoire avec `GGobi`, comme le propose <>. Par exemple, nous avons reproduit ci-dessous l'analyse des données +fleabeetles \{DPpackage\}+ que ces auteurs proposent dans leur article (p. 4). Ce code fonctionne sous Mac OS X, avec GGobi 2.1.7 (sept. 2007). ------------------------------------------------------------------------------ > data(fleabeetles, package="DPpackage") > head(fleabeetles) fjft sjft mwhbee mwafp faa awfs species 1 191 131 53 150 15 104 1 2 185 134 50 147 13 105 1 3 200 137 52 144 14 102 1 4 173 127 50 144 16 97 1 5 171 118 49 153 13 106 1 6 160 118 47 140 15 99 1 > library(clusterfly) > ref <- as.numeric(fleabeetles\$species) > fl <- scale(as.matrix(fleabeetles[,1:6])) > x <- ggobi(fleabeetles)\$fleabeetles > glyph_colour(x) <- ref > glyph_colour(x) <- clarify(kmeans(fl,3)\$cluster, ref) ------------------------------------------------------------------------------ Les trois dernières lignes permettent de transférer les données à GGobi, puis de refléter la solution obtenue par la procédure ``k-means'' dans les données. Puisque l'algorithme est pour partie stochastique, la solution peut différer d'un essai à l'autre, à moins de fixer le générateur de nombres aléatoires (e.g. +set.seed(134)+). La séquence vidéo ci-dessous illustre une projection des 74 individus dans l'espace engendré par les 6 variables, en faisant apparaître en couleur les différentes espèces (mode *2D Tour*). // border: 1px solid silver; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ En plus ``exotique'', on peut utiliser http://www.xlispstat.8rf.com/[xlispstat] pour effectuer des visualisations dynamiques de données. Bien que son utilisation soit largement inférieure de nos jours par rapport à l'essor connu par ce logiciel dans les années 90, il reste utile pour un certain nombre d'opérations. Son installation sous Mac OS X est assez simple, puisque si l'on a installé le gestionnaire de paquets +fink+, il suffit de taper en ligne de commande +$ fink install xlispstat+. Normalement, depuis Mac OS X dans sa version 10.5, le serveur X11 est pré-installé. Si ce n'est pas le cas, on veillera à l'installer car xlispstat fonctionne sous X11 (eh oui, on ne change pas une équipe qui gagne). Il existe de nombreux tutoriels disponibles sur internet, et des ressources diverses, par exemple sur le site maintenu par Jan de Leeuw (éditeur du _Journal of Statistical Software_, entre autres) : http://www.stat.ucla.edu/xlispstat/code/[Département de statistiques de l'UCLA]. Pour visualiser les données de l'exemple précédent, on procèderai de la manière suivante : ------------------------------------------------------------------------------ \$ xlispstat XLISP-PLUS version 3.04 Portions Copyright (c) 1988, by David Betz. Modified by Thomas Almy and others. XLISP-STAT Release 3.52.20 (Beta). Copyright (c) 1989-1999, by Luke Tierney. > (def fb (read-data-columns "fleabeetles.data")) FB > (def fjft (select fb 0)) FJFT > (def sjft (select fb 1)) SJFT > (def mwhbee (select fb 2)) MWHBEE > (def mwafp (select fb 3)) MWAFP > (def faa (select fb 4)) FAA > (def awfs (select fb 5)) AWFS > (scatterplot-matrix (list fjft sjft mwhbee mwafp faa awfs) :variable-labels (list "fjft" "sjft" "mwhbee" "mwafp" "faa" "awfs")) # > (spin-plot (list fjft sjft mwhbee mwafp faa awfs) :title "fleabeetles" ) # ------------------------------------------------------------------------------ La commande +scatterplot-matrix+ affiche une matrice de dispersion, alors que la commande +spin-plot+ permet de naviguer dans l'espace à 6 dimensions défini par les variables. La variable +species+ pourrait être utilisée pour colorier les points selon la modalité. On obtient un résultat qui ressemble aux deux captures d'écran ci-dessous. [[xlisp1]] .Matrice de dispersion avec xlispstat. image::xlisp1.png[xlisp1.png] [[xlisp2]] .Spin-plot avec xlispstat. image::xlisp2.png[xlisp2.png] Le package +clusterfly+ facilite également la mise en correspondance des solutions obtenues à partir de différentes techniques. En reprenant l'exemple de classification CAH (avec la distance de Ward) et ``k-means'' des données +life+, on peut reproduire très facilement une matrice de confusion sous forme graphique. ------------------------------------------------------------------------------ > life.cf <- clusterfly(life) > life.cf <- cfly_cluster(life.cf, kmeans, 5) > life.cf <- cfly_cluster(life.cf, hierarchical, "ward") > life.cf Data: m0, m25, m50, m75, w0, w25, w50, w75 [31x8] Clusters: kmeans, hierarchical > cfly_fluct(life.cf, "kmeans","hierarchical", clarify=TRUE) ------------------------------------------------------------------------------ [[km_plot3]] .Diagramme de fluctuation (CAH vs. ``k-means''). image::km3.png[km3.png] Analyse par les médoïdes ^^^^^^^^^^^^^^^^^^^^^^^^ Plutôt que d'utiliser les centroïdes pour évaluer le critère d'optimalité d'appartenance à un cluster, on peut utiliser à la place ce que l'on appelle des médoïdes. La fonction +pam+ implémente cette dernière approche. Analyse combinée ~~~~~~~~~~~~~~~~ Face à de gros volumes de données, les techniques précédentes peuvent se révéler assez lourdes, voire inefficaces. Le problème tient essentiellement à la détermination _a priori_ du nombre de classe pour la classification par nuées dynamiques. Une solution généralement adoptée consiste en une approche mixant à la fois la classification hiérarchique et l'analyse en cluster, en suivant les étapes suivantes : 1. Utilisation procédure kmeans pour restreindre le nombre d'individus à classer (on prend généralement un nombre de classe k = 1/10 de l'effectif total). 2. Réalisation d'une CAH sur les barycentres des classes obtenues à l'étape 1. 3. Détermination du nombre de classes à partir des résultats de la CAH, et représentation du dendogramme. 4. Réallaocation itérative des individus sur les barycentres obtenus à l'étape 3 (le poids des barycentres est fixé à la somme des poids des individus de la classe correspondante). Dans le cas où la taille de l'effectif reste raisonnable, l'étape 1 peut être omise. Analyse sur scores ~~~~~~~~~~~~~~~~~~ On peut également utiliser directement des coordonnées issues d'une analyse en composantes principales. Dans ce cas, par construction, les variables composites (composantes) sont orthogonales mais rien n'enmpêche de poursuivre par une classification. ------------------------------------------------------------------------------ > life.acp <- princomp(life) > life.acp.km <- kmeans(life.acp\$scores, 4) > life.acp.km K-means clustering with 4 clusters of sizes 13, 2, 2, 14 Cluster means: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 1 10.2023451 -0.2348864 0.46560329 -0.09725453 0.05446874 0.1147774 2 -42.5479375 0.1814746 -0.08688113 0.43175704 0.19716032 0.2812810 3 -17.1627365 -1.3694832 1.00584535 -0.68731409 -0.19598963 0.5302735 4 -0.9435098 0.3878243 -0.56362652 0.12681593 -0.05074536 -0.2225154 Comp.7 Comp.8 1 0.031219670 -0.1343382 2 0.383423624 -0.3286510 3 -0.650278255 0.1237019 4 0.009132397 0.1540210 Clustering vector: Algeria Cameroon Madagascar 1 2 2 Mauritius Reunion Seychelles 4 4 1 South Africa(C) South Africa(W) Tunisia 3 1 4 Canada Costa Rica Dominican Rep 1 1 1 El Salvador Greenland Grenada 4 4 4 Guatemala Honduras Jamaica 3 4 4 Mexico Nicaragua Panama 4 1 1 Trinidad(62) Trinidad (67) United States (66) 1 4 1 United States (NW66) United States (W66) United States (67) 4 1 1 Argentina Chile Columbia 1 4 4 Ecuador 4 Within cluster sum of squares by cluster: [1] 752.4615 25.5000 15.0000 622.7857 Available components: [1] "cluster" "centers" "withinss" "size" > plot(life.acp\$scores, xlab="Comp. 1", ylab="Comp. 2", type="n") > text(life.acp\$scores, abbreviate(rownames(life),6), col=life.acp.km\$cluster) > abline(h=0); abline(v=0) ------------------------------------------------------------------------------ [[km_plot4]] .Projection des individus (données +life+) dans le plan factoriel (1,2). image::km4.png[km4.png] On retrouve les 4 clusters distribués sur l'axe horizontal, de droite à gauche : en d'autres termes, le premier axe factoriel synthétise la variance inter-cluster. La variance intra-cluster s'exprime quant à elle sur le second axe factoriel. Les valeurs propres correspondates sont d'ailleurs très constrastées, comme on peut le voir dans le diagramme des VP ci-dessous. [[km_acp]] .Diagramme des valeurs propres pour l'ACP des données +life+. image::km_acp.png[km_acp.png] La variance intra-cluster vaut : ------------------------------------------------------------------------------ > life.acp.km\$withinss [1] 752.4615 25.5000 15.0000 622.7857 ------------------------------------------------------------------------------ alors que la variance totale (identique à +ssw.1+ calculé plus haut) vaut 7012.5. La variance entre les clusters peut être estimée à 151283.5 : la variance intra représente ainsi moins de 3 % de la variance inter. ------------------------------------------------------------------------------ > (nrow(life.acp\$scores)-1)*sum(apply(life.acp\$scores, 2, var)) [1] 7012.516 > var(life.acp.km\$withinss) [1] 151283.5 > mean(life.acp.km\$withinss)/var(life.acp.km\$withinss) [1] 0.002339559 ------------------------------------------------------------------------------ Si on travaille avec des données binaires, il est préférable de partir sur les coordonnées d'une AFC. ------------------------------------------------------------------------------ > life.afc <- corresp(life, nf=6) > life.afc.km <- kmeans(life.afc\$rscore, 4) > life.afc.km K-means clustering with 4 clusters of sizes 2, 3, 7, 19 Cluster means: [,1] [,2] [,3] [,4] [,5] [,6] 1 -1.6425851 0.60173703 0.35313145 -0.9393376 -0.99442954 -2.12659516 2 -0.5701006 -0.03258134 -2.40742059 0.4177319 0.98256649 0.05012366 3 -0.5688215 0.47073539 0.05067545 -0.5359058 -0.55028939 1.34216687 4 0.4674511 -0.19462197 0.31167916 0.1381490 0.02185289 -0.20649978 Clustering vector: Algeria Cameroon Madagascar 1 3 1 Mauritius Reunion Seychelles 4 4 2 South Africa(C) South Africa(W) Tunisia 3 4 2 Canada Costa Rica Dominican Rep 4 4 4 El Salvador Greenland Grenada 3 4 4 Guatemala Honduras Jamaica 3 4 4 Mexico Nicaragua Panama 4 4 4 Trinidad(62) Trinidad (67) United States (66) 2 4 4 United States (NW66) United States (W66) United States (67) 4 4 4 Argentina Chile Columbia 4 3 3 Ecuador 3 Within cluster sum of squares by cluster: [1] 5.188307 32.732206 30.393023 54.804409 Available components: [1] "cluster" "centers" "withinss" "size" > plot(life.afc\$rscore,col=life.afc.km\$cluster,pch=19, + xlab="Comp. 1",ylab="Comp. 2") > abline(h=0); abline(v=0) ------------------------------------------------------------------------------ [[km_plot5]] .Projection des individus dans le plan factoriel (1,2). image::km5.png[km5.png] Une autre option pour les données binaires consiste à adapter la métrique utilisée. C'est ce que propose la fonction +daisy+. Les arbres de classification ---------------------------- À la différence des méthodes de classificaiton, où l'objectif est clairement d'identifier des groupes homogènes d'individus du point de vue de leurs réponses sur un ensemble de variables, l'objet des méthodes de classification par arbres se concentre plutôt sur la sélection des variables. On a une représentation classique en arbre binaire, avec des noeuds terminaux. Ceux-ci forment une partition de l'espace des variables. Idéalement, cette partition permet de classer correctement les individus en fonction des réponses observées. La construction d'un arbre binaire peut être effectuée sous R avec les ``packages'' +tree+ et +rpart+. Notons que ce type d'approche s'utilise également dans un contexte de modélisation : il s'agit alors d'_arbres de régression_, où les noeuds terminaux fournissent une valeur prédite. Globalement, l'ensemble de ces techniques est également regroupé sous le terme méthodes CART. Un exemple ~~~~~~~~~~ On charge dans un premier un jeu de données contenu dans le ``package'' +mlbench+. Il s'agit d'un ``package'' contenant exclusivement des données réelles ou artificielles utilisées dans le domaine de l'apprentissage automatique. La description du jeu de données est reproduite ci-dessous : .............................................................................. BreastCancer -- Wisconsin Breast Cancer Database The objective is to identify each of a number of benign or malignant classes. Samples arrive peri- odically as Dr. Wolberg reports his clinical cases. The database therefore reflects this chronological grouping of the data. This grouping information appears immediately below, having been removed from the data itself. Each variable except for the first was converted into 11 primitive numerical attributes with values ranging from 0 through 10. There are 16 missing attribute values. See cited below for more details. .............................................................................. On commence par charger les données et produire un résumé numérique du dataframe. ------------------------------------------------------------------------------ > data(BreastCancer, package="mlbench") > summary(BreastCancer) Id Cl.thickness Cell.size Cell.shape Marg.adhesion Length:699 1 :145 1 :384 1 :353 1 :407 Class :character 5 :130 10 : 67 2 : 59 2 : 58 Mode :character 3 :108 3 : 52 10 : 58 3 : 58 4 : 80 2 : 45 3 : 56 10 : 55 10 : 69 4 : 40 4 : 44 4 : 33 2 : 50 5 : 30 5 : 34 8 : 25 (Other):117 (Other): 81 (Other): 95 (Other): 63 Epith.c.size Bare.nuclei Bl.cromatin Normal.nucleoli Mitoses 2 :386 1 :402 2 :166 1 :443 1 :579 3 : 72 10 :132 3 :165 10 : 61 2 : 35 4 : 48 2 : 30 1 :152 3 : 44 3 : 33 1 : 47 5 : 30 7 : 73 2 : 36 10 : 14 6 : 41 3 : 28 4 : 40 8 : 24 4 : 12 5 : 39 (Other): 61 5 : 34 6 : 22 7 : 9 (Other): 66 NA's : 16 (Other): 69 (Other): 69 (Other): 17 Class benign :458 malignant:241 ------------------------------------------------------------------------------ Pour mettre en oeuvre l'arbre binaire de décision (la variable de classe considérée est ici la nature de la tumeur -- bénigne ou maligne), nous sélectionnons un échantillon aléatoire qui nous servira de base. Le taux de sondage est fixé à 0.2. ------------------------------------------------------------------------------ > set.seed(124) > bc.ech <- BreastCancer[sample(1:nrow(BreastCancer),ceiling(nrow(BreastCancer)*0.2)),-1] > summary(bc.ech) Cl.thickness Cell.size Cell.shape Marg.adhesion Epith.c.size 1 :32 1 :86 1 :74 1 :82 2 :80 5 :28 10 :13 2 :16 2 :13 4 :13 3 :25 2 : 8 3 :12 10 : 9 1 :10 4 :19 3 : 8 10 : 9 3 : 8 6 :10 8 :11 6 : 6 7 : 8 5 : 8 3 : 9 2 :10 8 : 6 4 : 6 6 : 7 5 : 9 (Other):15 (Other):13 (Other):15 (Other):13 (Other): 9 Bare.nuclei Bl.cromatin Normal.nucleoli Mitoses Class 1 :90 2 :39 1 :95 1 :115 benign :100 10 :29 3 :34 8 : 9 2 : 10 malignant: 40 2 : 6 1 :27 10 : 9 3 : 6 5 : 5 7 :14 2 : 8 4 : 3 4 : 2 4 : 7 6 : 5 7 : 2 (Other): 6 5 : 6 9 : 5 10 : 2 NA's : 2 (Other):13 (Other): 9 (Other): 2 ------------------------------------------------------------------------------ Nous pouvons ensuite procéder à la détermination de l'arbre. Attention, si l'on réalise l'analyse sur les variables telles quelles (variables ordinales, cf. +class(bc.ech$Cell.size)+), la procédure n'utilisera pas toutes les variables pour effectuer le partitionnement. ------------------------------------------------------------------------------ > tmp <- data.frame(matrix(as.numeric(as.matrix(bc.ech[,-10])),nc=9),bc.ech[,10]) > colnames(tmp) <- colnames(bc.ech) > bc.ech2 <- tmp > pairs(bc.ech2[,-10]) > summary(bc.ech2) Cl.thickness Cell.size Cell.shape Marg.adhesion Min. : 1.000 Min. : 1.00 Min. : 1.000 Min. : 1.000 1st Qu.: 2.000 1st Qu.: 1.00 1st Qu.: 1.000 1st Qu.: 1.000 Median : 4.000 Median : 1.00 Median : 1.000 Median : 1.000 Mean : 3.929 Mean : 2.95 Mean : 2.907 Mean : 2.764 3rd Qu.: 5.000 3rd Qu.: 4.00 3rd Qu.: 4.000 3rd Qu.: 4.000 Max. :10.000 Max. :10.00 Max. :10.000 Max. :10.000 Epith.c.size Bare.nuclei Bl.cromatin Normal.nucleoli Min. : 1.000 Min. : 1.000 Min. : 1.000 Min. : 1.00 1st Qu.: 2.000 1st Qu.: 1.000 1st Qu.: 2.000 1st Qu.: 1.00 Median : 2.000 Median : 1.000 Median : 3.000 Median : 1.00 Mean : 3.086 Mean : 3.391 Mean : 3.386 Mean : 2.80 3rd Qu.: 4.000 3rd Qu.: 5.000 3rd Qu.: 4.000 3rd Qu.: 3.25 Max. :10.000 Max. :10.000 Max. :10.000 Max. :10.00 NA's : 2.000 Mitoses Class Min. : 1.0 benign :100 1st Qu.: 1.0 malignant: 40 Median : 1.0 Mean : 1.5 3rd Qu.: 1.0 Max. :10.0 > library(rpart) > bc.tree <- rpart(Class ~ ., data=bc.ech2, parms=list(split='information'), cp=0.001) > summary(bc.tree) ------------------------------------------------------------------------------ Pour visualiser l'arbre binaire, on utlise la commande +plot+, comme suit ------------------------------------------------------------------------------ > opar <- par(xpd=NA) > ------------------------------------------------------------------------------ Références ---------- + [[[Banfield1993]]] Banfield, J. D. and Raftery, A. E. (1993). Model-basd Gaussian and non-Gaussian clustering. _Biometrics_, _49_, 803-821. + [[[Bartholomew1999]]] Bartholomew, D. J. and Knott, M. (1999). _Latent Variable Models and Factor Analysis_. Arnold. + [[[Berthold2003]]] M. Berthold and D. J. Hand (2003). _Intelligent Data Analysis_. Springer. + [[[Brown2006]]] T. A. Brown (2006). _Confirmatory Factor Analysis for Applied Research_. Guilford. [http://people.bu.edu/tabrown/cfabook.html[www]] + [[[Carroll1997]]] J. D. Carroll, P. E. Green, and A. Chaturvedi (1997). _Mathematical Tools for Applied Multivariate Analysis_. Academic Press. + [[[Everitt2005]]] B. Everitt (2005). _An R and S-PLUS Companion to Multivariate Analysis_. Springer. + [[[Falissard2001]]] B. Falissard (2001). _Comprendre et utiliser les statistiques dans les sciences de la vie_. Masson + [[[Fraley2003]]] Fraley, C. and Raftery, A. E. (2003). Enhanced Software for Model-Based Clustering, Density Estimation, and Discriminant Analysis : MCLUST. _Journal of Classification_, _20_, 263-286. + [[[Hastie2001]]] T. Hastie, R. Tibshirani, and J. Friedman (2001). _The Elements of Statistical Learning. Data Mining, Inference, and Prediction_. Springer. + [[[Jolliffe2004]]] I. T. Jolliffe (2004). _Principal Component Analysis_. Springer. + [[[Kline2005]]] R. B. Kline (2005). _Principles and Practice of Structural Equation Modeling_. Guilford. + [[[Lebart2006]]] L. Lebart, M. Piron, and A. Morineau (2006). _Statistique exploratoire multidimensionnelle. Visualisation et inférence en fouille de données_. Dunod. + [[[Michie1994]]] D. Michie, D. J. Spiegelhalter, and C. C. Taylor (eds) (1994). _Machine Learning, Neural and Statistical Classification_. Ellis Horwood. [http://www.amsta.leeds.ac.uk/~charles/statlog/[pdf]] + [[[Ripley1994]]] Ripley, B. D. (1994) Flexible non-linear approaches to classification. In _From Statistics to Neural Networks. Theory and Pattern Recognition Applications_ eds V. Cherkassky, J. H. Friedman and H. Wechsler, Springer, pp. 105-126. + [[[Kaufman1990]]] Kaufman, L. and Rousseeuw, P.J. (1990). _Finding Groups in Data: An Introduction to Cluster Analysis_. Wiley, New York. + [[[Stevens2002]]] J. P. Stevens (2002). _Applied Multivariate Statistics for the Social Sciences_ (4e ed.). Lawrence Erlbaum Associates. + [[[Struyf1997]]] Struyf, A., Hubert, M., and Rousseuw, P. (1997). Clustering in an Object-Oriented Environment. Journal of Statistical Software, 1(4). [http://www.jstatsoft.org/v01/i04[pdf]] + [[[Thurstone1947]]] L. L. Thurstone (1947). _Multiple-factor Analysis_. Chicago: University of Chicago Press. + [[[Venables2002]]] W. N. Venables and B. D. Ripley (2002). _Modern Applied Statistics with S_ (4th Ed.). Springer. + [[[Webb2002]]] A. R. Webb (2002). _Statistical Pattern Recognition_. Wiley + [[[Wickham2006]]] Wickham, H., Hofmann, H., and Cook, D. (2007). Exploring cluster analysis. (Paper submitted to CMV) [http://had.co.nz/clusterfly/clusters.pdf[pdf]]