Avant propos

Après avoir manipulé en détail les méthodes d'analyse factorielle que sont l'ACP et l'AFC, nous allons maintenant nous intéresser à la classification automatique. Là où les analyses graphiques nous permettent parfois de détecter de probables structures de groupes, nous n'avons pour l'instant pas utilisé d'outils pour définir automatiquement ces groupes (appelés aussi clusters). Les méthodes destinées à accomplir cette tâche sont appelés des algorithmes de clustering et nous allons en étudier quelques uns, à commencer par la Classification Ascendante Hiérarchique (CAH) dans l'exercice 4 ci-dessous.

Exercice 4

Dans cet exercice, nous utilisons la base de données très ((beaucoup ?) trop ?) classique iris, présente par défaut dans la version de base de R. Ce dataset regroupe un total de 150 iris de 3 espèces différentes (50 chacune) observées sur 4 variables relatives à leur pétales et sépales. La vraie classification étant connue d'avance, nous allons pouvoir tester nos méthodes de clustering, pour ensuite comparer les groupes obtenus automatiquement avec les vraies espèces d'iris. Les algorithmes de clustering sont des méthodes dites d'apprentissage non-supervisé, un sous-ensemble des méthodes d'apprentissage automatique (machine learning en anglais). Comme habituellement, importons packages et données, et visualisons les données brutes.

#rm(list = ls())
library(tidyverse)
library(gridExtra)
library(FactoMineR)

raw_db = iris %>% as_tibble

raw_db
## # A tibble: 150 x 5
##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
##           <dbl>       <dbl>        <dbl>       <dbl> <fct>  
##  1          5.1         3.5          1.4         0.2 setosa 
##  2          4.9         3            1.4         0.2 setosa 
##  3          4.7         3.2          1.3         0.2 setosa 
##  4          4.6         3.1          1.5         0.2 setosa 
##  5          5           3.6          1.4         0.2 setosa 
##  6          5.4         3.9          1.7         0.4 setosa 
##  7          4.6         3.4          1.4         0.3 setosa 
##  8          5           3.4          1.5         0.2 setosa 
##  9          4.4         2.9          1.4         0.2 setosa 
## 10          4.9         3.1          1.5         0.1 setosa 
## # ... with 140 more rows

Pour pouvoir faire l'analyse dans des conditions 'réelles', nous devons appliquer les méthodes de clustering sur les 4 premières variables seulement, qui constitueront notre db d'étude. La variable Species va être stockée à part et nous servira seulement au moment d'évaluer la qualité des résultats. Jetons un premier regard également à un summary de la base de données ainsi obtenu.

Species = raw_db %>% pull(Species)
db = raw_db %>% select(- Species)
summary(db)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500

Pour une analyse plus détaillée et détecter d'éventuelles structures de groupes, nous allons commencer par appliquer une première ACP, avec laquelle nous sommes maintenant familliers.

res_pca = PCA(db, ncp = 4, graph = F)
 
vp = res_pca$eig %>% data.frame %>% rownames_to_column(var = 'Composante') %>% as_tibble

ggplot(vp) + geom_col(aes(x = Composante, y = percentage.of.variance)) +
             geom_point(aes(x = Composante, y = cumulative.percentage.of.variance)) +
             ggtitle("Pourcentage de variance expliquée cumulée") + 
             geom_hline(aes(yintercept = 90), color = 'red')

Nous voyons que deux axes suffisent à expliquer plus de 95% de la variance des données, nous allons donc retenir 2 composantes principales et, après avoir vérifié la contribution et qualité de représentation des individus et variables, les projeter sur ce premier plan factoriel.

summary(res_pca)
## 
## Call:
## PCA(X = db, ncp = 4, graph = F) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4
## Variance               2.918   0.914   0.147   0.021
## % of var.             72.962  22.851   3.669   0.518
## Cumulative % of var.  72.962  95.813  99.482 100.000
## 
## Individuals (the 10 first)
##                  Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## 1            |  2.319 | -2.265  1.172  0.954 |  0.480  0.168  0.043 | -0.128
## 2            |  2.202 | -2.081  0.989  0.893 | -0.674  0.331  0.094 | -0.235
## 3            |  2.389 | -2.364  1.277  0.979 | -0.342  0.085  0.020 |  0.044
## 4            |  2.378 | -2.299  1.208  0.935 | -0.597  0.260  0.063 |  0.091
## 5            |  2.476 | -2.390  1.305  0.932 |  0.647  0.305  0.068 |  0.016
## 6            |  2.555 | -2.076  0.984  0.660 |  1.489  1.617  0.340 |  0.027
## 7            |  2.468 | -2.444  1.364  0.981 |  0.048  0.002  0.000 |  0.335
## 8            |  2.246 | -2.233  1.139  0.988 |  0.223  0.036  0.010 | -0.089
## 9            |  2.592 | -2.335  1.245  0.812 | -1.115  0.907  0.185 |  0.145
## 10           |  2.249 | -2.184  1.090  0.943 | -0.469  0.160  0.043 | -0.254
##                 ctr   cos2  
## 1             0.074  0.003 |
## 2             0.250  0.011 |
## 3             0.009  0.000 |
## 4             0.038  0.001 |
## 5             0.001  0.000 |
## 6             0.003  0.000 |
## 7             0.511  0.018 |
## 8             0.036  0.002 |
## 9             0.096  0.003 |
## 10            0.293  0.013 |
## 
## Variables
##                 Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## Sepal.Length |  0.890 27.151  0.792 |  0.361 14.244  0.130 | -0.276 51.778
## Sepal.Width  | -0.460  7.255  0.212 |  0.883 85.247  0.779 |  0.094  5.972
## Petal.Length |  0.992 33.688  0.983 |  0.023  0.060  0.001 |  0.054  2.020
## Petal.Width  |  0.965 31.906  0.931 |  0.064  0.448  0.004 |  0.243 40.230
##                cos2  
## Sepal.Length  0.076 |
## Sepal.Width   0.009 |
## Petal.Length  0.003 |
## Petal.Width   0.059 |
plot(res_pca, axes = c(1,2), choix = 'var')

On trouve donc une structure où la variable Sepa.width contribue fortement au deuxième axe, et est globalement orthogonale aux trois autres variables, qui elles contribuent majoritairement au premier axe. Voyons maintenant si la projection des individus dans le plan factoriel fait émerger une structure de groupe par rapport aux caractéristiques que nous venons de détailler.

plot(res_pca, axes = c(1,2), choix = 'ind')

Si le nombre exact de groupes apparents n'est pas parfaitement établi (probablement 2 ou 3 à première vue), un séparation apparait tout de même évidente et il est légitime d'ensuite chercher à appliquer un algorithme de clustering. Gardons à l'esprit que, de même que pour un précédent exercice, nous avons écarté deux axes factoriels, qui auraient pu nous aider à discriminer les individus d'une façon différente et peut être informative.

Mettons de côté mainenant ces résultats pour nous concentrer sur la mise en place de la CAH, à l'aide de la fonction hclust() du package FactoMineR. Il a été vu en cours que la CAH permet de construire un dendrogramme, un graphique qui résume la structure des groupes qui sont successivement construits à chaque étape de regroupement. Sur ce grapique, la taille des segments horizontaux sont proportionels à la valeur de dissimilarité inter-groupe entre les deux groupes rassemblés. Puisque l'on cherche à regrouper ensemble des individus qui se ressemblent et séparer les autres, il faut éviter justement de rassembler des sous groupes avec un fort score de dissimilarité inter-groupe. C'est pourquoi nous cherchons à stopper à 'couper' le dendrogramme avant que de grandes 'branches' horizontales n'apparaissent. Voyons tout d'abord un exemple d'utilisation de la fonction hclust(), après avoir calculé les distances euclidiennes entre les individus (centrés réduits avec scale()) à l'aide de la fonction dist():

db %>% scale() %>% 
       dist(method = 'euclidean') %>% 
       hclust() %>% 
       plot()

Tout ceci est très intéressant, mais n'oublions pas qu'il existe divers critères de regroupement avec leur caractéristiques propres, qui sont proposés comme argument de la fonction hclust. Affichons donc ci-dessous les résultats en fonction du critère choisi:

dist = db %>% scale() %>% dist(method = 'euclidean')

par(mfrow=c(2,2))
hclust(dist, method="ward.D2") %>% plot(main = "Dendrogramme avec la methode de Ward")
hclust(dist, method = "single") %>% plot(main = "Dendrogramme avec le lien simple")
hclust(dist, method = "complete") %>% plot(main = "Dendrogramme avec le lien complet")
hclust(dist, method = "average") %>% plot(main = "Dendrogramme avec le lien moyen")