Exercice 1¶

Commençons par récupérer le jeu de données, puis regardons ce qu'il contient.

Note : toute l'analyse ci-dessous peut être réalisée dans un autre langage (Python, Perl, ...). À vous de choisir ! Il y aurait un travail de recherche à effectuer ceci dit pour obtenir les équivalents de chaque commande (ou si c'est impossible, des contournements / autres façons de procéder).

In [1]:
data <- read.csv("../data/Pokemon.csv")
dim(data)
head(data)
  1. 1072
  2. 13
A data.frame: 6 × 13
numbernametype1type2totalhpattackdefensesp_attacksp_defensespeedgenerationlegendary
<int><chr><chr><chr><int><int><int><int><int><int><int><int><lgl>
11Bulbasaur GrassPoison31845 49 49 65 65451FALSE
22Ivysaur GrassPoison40560 62 63 80 80601FALSE
33Venusaur GrassPoison52580 82 83100100801FALSE
43Mega Venusaur GrassPoison62580100123122120801FALSE
53Gigantamax VenusaurGrassPoison52580 82 83100100801FALSE
64Charmander Fire 30939 52 43 60 50651FALSE

Il paraît clair que les colonnes "name", "type1", "type2" et "legendary" ne peuvent pas être utilisées, n'étant pas numériques. Elles peuvent être intéressantes ceci dit pour la visualisation ultérieure (variables supplémentaires, comme avec l'ACP/AFC). Ensuite, la génération d'un pokemon semble peu pertinente car elle définit en quelque sorte déjà un clustering : génération 1 = un groupe, génération 2 = un groupe etc. Finalement, la colonne "number" correspond probablement à un identifiant, et doit être ignorée. Vérifions cela rapidement :

In [2]:
length(unique(data$number))
levels(as.factor(data$generation)) #on peut aussi écrire simplement unique(data$generation)
898
  1. '0'
  2. '1'
  3. '2'
  4. '3'
  5. '4'
  6. '5'
  7. '6'
  8. '7'
  9. '8'

Bon, "number" n'est pas identifiant sur ce jeu de données. On le voit en fait dans l'extrait ci-dessus : Venusaur. Cette exemple semble indiquer que les numéros en double correspondent à des variantes d'un même pokemon. Vérifions.

In [3]:
# Recherche des numéros en double:
dbl <- unlist(lapply(unique(data$number), function(n) { if (sum(data$number == n) >= 2) return(n); return(NULL) }))
# Note: astuce, NULL dans une liste est ignoré par unlist(). Il y a plein d'autres façons de faire hein.
# Note2: attention length(data$number == n) renverrait 1072 ; il faut utiliser sum() ici.
dbl
  1. 3
  2. 6
  3. 9
  4. 12
  5. 15
  6. 18
  7. 19
  8. 20
  9. 25
  10. 26
  11. 27
  12. 28
  13. 37
  14. 38
  15. 50
  16. 51
  17. 52
  18. 53
  19. 65
  20. 68
  21. 74
  22. 75
  23. 76
  24. 77
  25. 78
  26. 79
  27. 80
  28. 83
  29. 88
  30. 89
  31. 94
  32. 99
  33. 103
  34. 105
  35. 110
  36. 115
  37. 122
  38. 127
  39. 130
  40. 131
  41. 133
  42. 142
  43. 143
  44. 144
  45. 145
  46. 146
  47. 150
  48. 181
  49. 199
  50. 208
  51. 212
  52. 214
  53. 222
  54. 229
  55. 248
  56. 254
  57. 257
  58. 260
  59. 263
  60. 264
  61. 282
  62. 302
  63. 303
  64. 306
  65. 308
  66. 310
  67. 319
  68. 323
  69. 334
  70. 354
  71. 359
  72. 362
  73. 373
  74. 376
  75. 380
  76. 381
  77. 382
  78. 383
  79. 384
  80. 386
  81. 413
  82. 428
  83. 445
  84. 448
  85. 460
  86. 475
  87. 479
  88. 487
  89. 492
  90. 531
  91. 554
  92. 555
  93. 562
  94. 569
  95. 618
  96. 641
  97. 642
  98. 645
  99. 646
  100. 647
  101. 648
  102. 658
  103. 678
  104. 681
  105. 710
  106. 711
  107. 718
  108. 719
  109. 720
  110. 741
  111. 745
  112. 746
  113. 774
  114. 800
  115. 809
  116. 812
  117. 815
  118. 818
  119. 823
  120. 826
  121. 834
  122. 839
  123. 841
  124. 842
  125. 844
  126. 849
  127. 851
  128. 858
  129. 861
  130. 869
  131. 875
  132. 876
  133. 877
  134. 879
  135. 884
  136. 888
  137. 889
  138. 890
  139. 892
  140. 893
  141. 898
In [4]:
data[data$number == 27,]
data[data$number == 445,]
data[data$number == 844,]
A data.frame: 2 × 13
numbernametype1type2totalhpattackdefensesp_attacksp_defensespeedgenerationlegendary
<int><chr><chr><chr><int><int><int><int><int><int><int><int><lgl>
4127Sandshrew Ground 3005075852030401FALSE
4227Alolan SandshrewIce Steel3005075901035407FALSE
A data.frame: 2 × 13
numbernametype1type2totalhpattackdefensesp_attacksp_defensespeedgenerationlegendary
<int><chr><chr><chr><int><int><int><int><int><int><int><int><lgl>
539445Garchomp DragonGround600108130 95 80851024FALSE
540445Mega GarchompDragonGround70010817011512095 924FALSE
A data.frame: 2 × 13
numbernametype1type2totalhpattackdefensesp_attacksp_defensespeedgenerationlegendary
<int><chr><chr><chr><int><int><int><int><int><int><int><int><lgl>
997844Sandaconda Ground510721071256570718FALSE
998844Gigantamax SandacondaGround510721071256570718FALSE

On ne va pas tous les faire : la conclusion est claire. Notons que les deux formes du pokemon Sandshrew n'ont pas les mêmes types, et que les caractéristiques de Sandaconda sont exactement les mêmes que celles de sa forme "Gigantamax". Après vérification il n'y a pas d'erreur : https://www.pokebip.com/pokedex/pokemon/dunaconda-gigamax. Bizarre tout de même. Monde étrange que celui des Pokemons...

Je ne sais pas vous, mais je suis curieux de savoir parmi tous ceux-là lesquels sont légendaires. Regardons.

In [5]:
dim(data[data$legendary,]) #Awai, 118 quand-même...
head(data[data$legendary,])
  1. 118
  2. 13
A data.frame: 6 × 13
numbernametype1type2totalhpattackdefensesp_attacksp_defensespeedgenerationlegendary
<int><chr><chr><chr><int><int><int><int><int><int><int><int><lgl>
195144Articuno Ice Flying58090 85100 95125 851TRUE
196144Galarian ArticunoPsychic Flying58090 85 85125100 958TRUE
197145Zapdos ElectricFlying58090 90 85125 901001TRUE
198145Galarian Zapdos FightingFlying58090125 90 85 901008TRUE
199146Moltres Fire Flying58090100 90125 85 901TRUE
200146Galarian Moltres Dark Flying58090 85 90100125 908TRUE

Ma curiosité n'étant pas encore satisfaite, demandons à Poképédia : "Un Pokémon légendaire est un Pokémon particulier, associé généralement à un mythe concernant la création et l'organisation du monde, et dont la rareté n'a d'égale que sa puissance au combat. L'obtention d'un Pokémon légendaire se fait donc en effectuant une quête particulière ou par le biais d'une distribution lors d'événements organisés par Nintendo."

In [6]:
# Pokemon légendaire = plus puissant ?
options(repr.plot.width=15, repr.plot.height=10)
hist((1:nrow(data))[ data[order(data$total, decreasing=TRUE), "legendary"] ], xlab="Rang", ylab="Compte", main="")
# Répartition des pokemons légendaires dans la liste par puissances décroissantes:
# Conclusion = oui, plutôt.
No description has been provided for this image

Il est temps de revenir à nos m... pokemons. De ce qui précède on conclut que les colonnes à retenir seraient "total, hp, attack, defense, sp_attack, sp_defense, speed". "total" est cependant redondant = somme des autres colonnes. Il risque d'aplanir les différences dans une certaine mesure, deux pokemons différents mais ayant le même total se verraient rapprochés artificiellement. On conserve donc finalement 6 colonnes.

In [7]:
data_clust <- subset(data, select=c("hp", "attack", "defense", "sp_attack", "sp_defense", "speed"))
rownames(data_clust) <- data$name #utile pour les graphes

Clustering hiérarchique¶

In [8]:
distances <- dist(data_clust)
h <- hclust(distances, method="ward.D") #distance de Ward: en général bons résultats
plot(h, labels=substr(data$name, 1, 5)) #substr(x, 1, 5): 5 premiers caractères (très moche sinon)
No description has been provided for this image

Le bas du dendrogramme est illisible, normal avec tant d'individus. La hauteur avant la dernière fusion est importante, indiquant qu'il faut garder au moins deux groupes. On peut regarder aussi le résultat pour 3, 4 ou 5 groupes par exemple (choix arbitraire). Je ne trace ici que $K=2$ et $K=3$ pour ne pas surcharger : on comparera à PAM et kmeans.

In [9]:
library(factoextra)
cl2_h <- cutree(h, 2)
fviz_cluster(list(data=data_clust, cluster=cl2_h))
Loading required package: ggplot2

Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

No description has been provided for this image

Sans surprise, on retrouve un groupe de pokemons plutôt puissant à droite, et un autre plus modeste à gauche. Vérifions cette affirmation:

In [10]:
par(mfrow=c(1,2))
hist(data[cl2_h==1, "total"], main="Groupe 1 à gauche")
hist(data[cl2_h==2, "total"], main="Groupe 2 à droite")
No description has been provided for this image
In [11]:
cl3_h <- cutree(h, 3)
fviz_cluster(list(data=data_clust, cluster=cl3_h))
No description has been provided for this image

Le cluster de droite se subdivise en moitiés haute et basse : comparons les deux individus extrêmes pour tenter de comprendre ce second axe ACP.

In [12]:
data[which(data$name %in% c("Shuckle", "Deoxys Attack Forme")),]
A data.frame: 2 × 13
numbernametype1type2totalhpattackdefensesp_attacksp_defensespeedgenerationlegendary
<int><chr><chr><chr><int><int><int><int><int><int><int><int><lgl>
273213Shuckle Bug Rock50520 10230 10230 52FALSE
475386Deoxys Attack FormePsychic 60050180 20180 201503 TRUE

Donc plus on est haut (resp. bas) sur le second axe plus on a de points de défense (resp. attaque).

PAM, kmeans¶

Ces deux algorithmes procédant de manière similaire, je les regroupe dans une section.

In [13]:
library(cluster)
cl2_p <- pam(data_clust, 2)
fviz_cluster(cl2_p)
cl2_k <- kmeans(data_clust, 2)
fviz_cluster(cl2_k, data=data_clust)
No description has been provided for this image
No description has been provided for this image
In [14]:
cl3_p <- pam(data_clust, 3)
fviz_cluster(cl3_p)
cl3_k <- kmeans(data_clust, 3)
fviz_cluster(cl3_k, data=data_clust)
No description has been provided for this image
No description has been provided for this image
In [15]:
par(mfrow=c(1,2))
hist(cl3_p$clustering)
hist(cl3_k$cluster)
No description has been provided for this image

À quelques exceptions près (Shuckle, Kadabra, ...), et modulo les tailles des clusters (cf. ci-dessus), les résultats sont les mêmes qu'avec hclust : pokemons puissants à droite, offensifs en bas, défensifs en haut.

Stabilité¶

Afin d'évaluer la confiance que l'on peut accorder en un clustering, comme évoqué en cours il y a trois pistes principales :

  • calculer un critère géométrique basé sur la forme d'un groupe,
  • comparer la classification obtenue à une référence,
  • chercher à apprendre le label d'un individu (mode supervisé).

Nous parlerons de classification supervisée plus tard dans le cours, donc je ne développe pas ce point ici. En gros le principe serait :

  1. lancer l'algorithme de clustering à $K$ groupes (de manière stable),
  2. déterminer l'erreur d'apprentissage $x \mapsto k$,
  3. si cette erreur est en-dessous d'un certain seuil, $K \leftarrow K+1$, revenir en 1.

Je choisis assez arbitrairement ici le critère géométrique Dunn, et le critère par comparison de partition Rand. Leur expression est résumée ci-dessous ; voir la vignette https://cran.r-project.org/web/packages/clusterCrit/vignettes/clusterCrit.pdf pour les détails.

"Let us denote by $d_{min}$ the minimal distance between points of different clusters and $d_{max}$ the largest within-cluster distance."
Then the Dunn index is $\frac{d_{min}}{d_{max}}$.

Rand index = taux de paires de points classés de façon similaire dans les deux partitions (même cluster, ou clusters différents). Pourquoi pas.

In [16]:
# D'abord, augmenter nstart pour améliorer le résultat du kmeans:
cl2_k <- kmeans(data_clust, 2, nstart=10)
cl3_k <- kmeans(data_clust, 3, nstart=10)

# Concernant PAM, l'aide (?pam) indique
# "By default, when ‘medoids’ are not specified, the algorithm first looks for a good initial set of medoids".
# On fait donc confiance disons.
In [17]:
library(clusterCrit)
int_crit <- "Dunn"

# intCriteria demande une matrice de réels en entrée (peu flexible...):
mdata <- apply(data_clust, 2, as.numeric)
intcrit_data <- matrix(
  c(intCriteria(mdata, cl2_k$cluster, int_crit),
    intCriteria(mdata, cl2_p$clustering, int_crit),
    intCriteria(mdata, cl2_h, int_crit),
    intCriteria(mdata, cl3_k$cluster, int_crit),
    intCriteria(mdata, cl3_p$clustering, int_crit),
    intCriteria(mdata, cl3_h, int_crit)),
  nrow=3, ncol=2)
barplot(intcrit_data, main="Index values", xlab="Number of clusters",
        col=c("darkblue","red","darkgreen"), legend = c("kmeans", "pam"," hclust"))
No description has been provided for this image

Au premier coup d'oeil, on a envie de dire "victoire de hclust". En y regardant de plus près, les valeurs de l'indice sont ridiculement basses : moins de 0.05. Cela signifie que les clusters sont assez grands, et très proches (en terme de distance single-linkage). C'est en effet ce qui ressort des graphes précédents. Ce critère géométrique dit donc "tu t'es planté y'a pas de clusters dans ce jeu de données". Il a peut-être raison.

Continuons tout de même avec le critère externe Rand (assez populaire, Google "rand index clustering"). À défaut d'avoir une partition de référence, je compare ici les partitions deux à deux. On obtient donc une mesure de l'adéquation entre les trois algorithmes.

In [18]:
ext_crit <- "Rand"

extcrit_data <- matrix(
  c(extCriteria(cl2_k$cluster, cl2_p$clustering, ext_crit),
    extCriteria(cl2_p$clustering, cl2_h, ext_crit),
    extCriteria(cl2_h, cl2_k$cluster, ext_crit),
    extCriteria(cl3_k$cluster, cl3_p$clustering, ext_crit),
    extCriteria(cl3_p$clustering, cl3_h, ext_crit),
    extCriteria(cl3_h, cl3_k$cluster, ext_crit)),
  nrow=3, ncol=2)
barplot(extcrit_data, main="Index values", xlab="Number of clusters",
        col=c("darkblue","red","darkgreen"), legend = c("kmeans", "pam"," hclust"))
No description has been provided for this image

Conclusion : les partitions sont plutôt semblables pour $K=2$. Un peu moins marqué pour $K=3$, mais on reste à + de 0.7 (le maximum étant 1). C'est pas mal, et donc les algorithmes retournent à peu près la même chose, conformément à l'impression visuelle précédente.

Exercice 2¶

J'ai choisi assez arbitrairement quatre jeux de données dont 3 sont quasi sûrs de poser problème.

In [19]:
data1 <- read.table("https://raw.githubusercontent.com/deric/clustering-benchmark/master/src/main/resources/datasets/artificial/cluto-t7-10k.arff", skip=13, sep=",")
target1 <- data1[,3] #0, 1, ..., 8, noise
target1[target1 == "noise"] <- "9" #je préfère un vecteur d'entiers
target1 <- as.integer(target1)
data1 <- data1[,-3]

data2 <- read.table("https://raw.githubusercontent.com/deric/clustering-benchmark/master/src/main/resources/datasets/artificial/3-spiral.arff", skip=12, sep=",")
target2 <- data2[,3] #1, 2, 3
data2 <- data2[,-3]

data3 <- read.table("https://raw.githubusercontent.com/deric/clustering-benchmark/master/src/main/resources/datasets/artificial/diamond9.arff", skip=9, sep=",")
target3 <- data3[,3] #0, 1, ..., 8
data3 <- data3[,-3]

data4 <- read.table("https://raw.githubusercontent.com/deric/clustering-benchmark/master/src/main/resources/datasets/artificial/target.arff", skip=18, sep=",")
target4 <- data4[,3] #1, 2 (3, 4, 5, 6: points isolés)
data4 <- data4[,-3]
In [20]:
# Regroupements à retrouver :
par(mfrow=c(2,2))
plot(data1, col=rainbow(10)[target1+1])
plot(data2, col=rainbow(3)[target2])
plot(data3, col=rainbow(9)[target3+1])
plot(data4, col=rev(rainbow(6))[target4])
No description has been provided for this image
In [21]:
# C'est parti !
cl <- function(data, k) kmeans(data, k, nstart=10)$cluster
plotAll <- function() {
  par(mfrow=c(2,2))
  plot(data1, col=cl(data1, 10))
  plot(data2, col=cl(data2, 3))
  plot(data3, col=cl(data3, 9))
  plot(data4, col=cl(data4, 6))
}
plotAll()
No description has been provided for this image

Comme prévu, seuls les losanges sont bien classés : les clusters sont alors de forme suffisamment "patatoïdale" pour que kmeans s'en sorte. Les trois autres situations mènent à des catastrophes, l'algorithme étant incapable de retrouver des formes arbitraires - comme dit en cours.

In [22]:
cl <- function(data, k) pam(data, k)$clustering
plotAll() #temps de calcul déjà beaucoup (beaucoup) plus long ! ...comme évoqué en cours
No description has been provided for this image

Résultats similaires à ceux du kmeans. Les clusters en haut à gauche changent légèrement mais restent mauvais.

In [23]:
cl <- function(data, k) cutree(hclust(dist(data), method="ward.D"), k)
plotAll()
No description has been provided for this image

C'est aussi mauvais qu'avant, car la distance de Ward correspond au critère que le kmeans minimise.

Essayons avec le single-linkage, qui intuitivement pourrait fonctionner dans certains cas ici :

In [24]:
cl <- function(data, k) cutree(hclust(dist(data), method="single"), k)
plotAll()
No description has been provided for this image

En présence de bruit cette méthode est vouée à l'échec, comme illustré en haut à gauche. De même, quand les groupes sont collés les uns aux autres l'effet de chaînage risque de les connecter dans un même cluster. En revanche, dans des cas idéaux en l'absence de bruit à droite, l'algorithme retrouve (logiquement) les classes attendues.