1. 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, 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 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.

Table: Aperçu des techniques d'analyses multivariées
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 :

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 Treillis displays, ainsi qu'un ouvrage à paraître chez Springer : Lattice. Multivariate Data Visualization with R (D. Sarkar, approx. 290 p.).

En revanche, le logiciel 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 !) : Interactive and Dynamic Graphics for Data Analysis. With R and GGobi (D. Cook and D. F. Swayne, 2007). Des 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 poursuite projective en 2D, au format mov).

2. Entrée en matière

2.1. Quelques rappels d'algèbre

2.2. Le cas de trois variables

Cette section est (très) largement inspirée du tutoriel Représentations triangulaires (D. Chessel, A. B. Dufour et J. R. Lobry, 2005).

3. Les méthodes d'analyse factorielle

3.1. 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.

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 [introduction] que cette technique s'applique préférentiellement sur un ensemble de variables continues. À cela, deux exceptions peuventêtre mentionnées :

Une bonne référence concernant ce sujet est [Jollife2004], mais on trouvera dans [Lebart2006] une présentation de l'ACP en lien avec l'[analyse factorielle des correspondances] ou l'[analyse des correspondances multiples]. 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.
[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).

3.2. L'Analyse Factorielle des Correspondances

3.3. L'Analyse des Correspondances Multiples

3.4. 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.

3.5. 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 [Brown2006] et [Stevens2002]. 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 [Brown2006] ou [Kline2005].

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 [Bartholomew1999]. En fait, l'idée de base repose sur le modèle en facteur commun de Thurstone [Thurstone1947] 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 :

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.
[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.

spss_factor.png

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)

afe_scree.png
afe_biplot.png

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.

afe_contour.png
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.

4. 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 [Michie1994] ou [Webb2002] (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. [Lebart2006] 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 [Falissard2001] dont l'approche est beaucoup plus pédagogique. Entre les deux approches, on trouve un bon compromis avec le livre de [Hastie2001].

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 [Kaufman1990], [Banfield1993], [Fraley2003]. 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 [Struyf1997] 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 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 :

4.1. 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).

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 [Everitt2005], 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)
hc1.png
Figure: Dendogramme obtenu avec la méthode de Ward.

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 :

hc2.png hc3.png

Il existe naturellement d'autres critères que l'examen du dendrograme. Mentionnons par exemple :

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 : 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 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” : The rules of SPSS' Hierarchical Cluster Analysis for processing ties.

spss_cluster.png

4.2. 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 (Cours de classification automatique de A. Carlier, par exemple) que le nombre le nombre total de partition, Nn, augmente plus qu'exponentiellement en fonction de n ; il vaut environ 5.1013 pour n=20. Le nombre de partitions contenant exactement k parties d'un ensemble En à n éléments vaut en fait

eq_partitions.png

À titre d'illustration, le programme 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.

Example: 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. Fast k-means code for Matlab ou K-mean clustering and Soft K-means clustering (site de D. MacKay, auteur du fameux livre Information Theory, Inference, and Learning Algorithms)

4.2.1. 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)
km1.png
Figure: Matrice de dispersion pour les données life avec les solutions “k-means”.

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)
km2.png
Figure: Évolution de la SS intra-groupe en fonction du nombre de clusters.

À 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 [Wickham2006]. 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).

En plus “exotique”, on peut utiliser 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) : 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"))
#<Object: 2748d8, prototype = SCATMAT-PROTO>
> (spin-plot
(list fjft sjft mwhbee mwafp faa awfs)
:title "fleabeetles"
)
#<Object: 277b78, prototype = SPIN-PROTO>

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.png
Figure: Matrice de dispersion avec xlispstat.
xlisp2.png
Figure: Spin-plot avec xlispstat.

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)
km3.png
Figure: Diagramme de fluctuation (CAH vs. “k-means”).

4.2.2. 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.

4.3. 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.

4.4. 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)
km4.png
Figure: Projection des individus (données life) dans le plan factoriel (1,2).

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.png
Figure: Diagramme des valeurs propres pour l'ACP des données life.

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)
km5.png
Figure: Projection des individus dans le plan factoriel (1,2).

Une autre option pour les données binaires consiste à adapter la métrique utilisée. C'est ce que propose la fonction daisy.

5. 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.

5.1. 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)
>

6. Références

  1. [Banfield1993] Banfield, J. D. and Raftery, A. E. (1993). Model-basd Gaussian and non-Gaussian clustering. Biometrics, 49, 803-821.

  2. [Bartholomew1999] Bartholomew, D. J. and Knott, M. (1999). Latent Variable Models and Factor Analysis. Arnold.

  3. [Berthold2003] M. Berthold and D. J. Hand (2003). Intelligent Data Analysis. Springer.

  4. [Brown2006] T. A. Brown (2006). Confirmatory Factor Analysis for Applied Research. Guilford. [www]

  5. [Carroll1997] J. D. Carroll, P. E. Green, and A. Chaturvedi (1997). Mathematical Tools for Applied Multivariate Analysis. Academic Press.

  6. [Everitt2005] B. Everitt (2005). An R and S-PLUS Companion to Multivariate Analysis. Springer.

  7. [Falissard2001] B. Falissard (2001). Comprendre et utiliser les statistiques dans les sciences de la vie. Masson

  8. [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.

  9. [Hastie2001] T. Hastie, R. Tibshirani, and J. Friedman (2001). The Elements of Statistical Learning. Data Mining, Inference, and Prediction. Springer.

  10. [Jolliffe2004] I. T. Jolliffe (2004). Principal Component Analysis. Springer.

  11. [Kline2005] R. B. Kline (2005). Principles and Practice of Structural Equation Modeling. Guilford.

  12. [Lebart2006] L. Lebart, M. Piron, and A. Morineau (2006). Statistique exploratoire multidimensionnelle. Visualisation et inférence en fouille de données. Dunod.

  13. [Michie1994] D. Michie, D. J. Spiegelhalter, and C. C. Taylor (eds) (1994). Machine Learning, Neural and Statistical Classification. Ellis Horwood. [pdf]

  14. [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.

  15. [Kaufman1990] Kaufman, L. and Rousseeuw, P.J. (1990). Finding Groups in Data: An Introduction to Cluster Analysis. Wiley, New York.

  16. [Stevens2002] J. P. Stevens (2002). Applied Multivariate Statistics for the Social Sciences (4e ed.). Lawrence Erlbaum Associates.

  17. [Struyf1997] Struyf, A., Hubert, M., and Rousseuw, P. (1997). Clustering in an Object-Oriented Environment. Journal of Statistical Software, 1(4). [pdf]

  18. [Thurstone1947] L. L. Thurstone (1947). Multiple-factor Analysis. Chicago: University of Chicago Press.

  19. [Venables2002] W. N. Venables and B. D. Ripley (2002). Modern Applied Statistics with S (4th Ed.). Springer.

  20. [Webb2002] A. R. Webb (2002). Statistical Pattern Recognition. Wiley

  21. [Wickham2006] Wickham, H., Hofmann, H., and Cook, D. (2007). Exploring cluster analysis. (Paper submitted to CMV) [pdf]