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

data <- read.csv("../data/Pokemon.csv")
dim(data)
head(data)

1. 1072
2. 13

A data.frame: 6 × 13

	number	name	type1	type2	total	hp	attack	defense	sp_attack	sp_defense	speed	generation	legendary
	<int>	<chr>	<chr>	<chr>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<lgl>
1	1	Bulbasaur	Grass	Poison	318	45	49	49	65	65	45	1	FALSE
2	2	Ivysaur	Grass	Poison	405	60	62	63	80	80	60	1	FALSE
3	3	Venusaur	Grass	Poison	525	80	82	83	100	100	80	1	FALSE
4	3	Mega Venusaur	Grass	Poison	625	80	100	123	122	120	80	1	FALSE
5	3	Gigantamax Venusaur	Grass	Poison	525	80	82	83	100	100	80	1	FALSE
6	4	Charmander	Fire		309	39	52	43	60	50	65	1	FALSE

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 :

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.

# 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

data[data$number == 27,]
data[data$number == 445,]
data[data$number == 844,]

A data.frame: 2 × 13

	number	name	type1	type2	total	hp	attack	defense	sp_attack	sp_defense	speed	generation	legendary
	<int>	<chr>	<chr>	<chr>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<lgl>
41	27	Sandshrew	Ground		300	50	75	85	20	30	40	1	FALSE
42	27	Alolan Sandshrew	Ice	Steel	300	50	75	90	10	35	40	7	FALSE

A data.frame: 2 × 13

	number	name	type1	type2	total	hp	attack	defense	sp_attack	sp_defense	speed	generation	legendary
	<int>	<chr>	<chr>	<chr>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<lgl>
539	445	Garchomp	Dragon	Ground	600	108	130	95	80	85	102	4	FALSE
540	445	Mega Garchomp	Dragon	Ground	700	108	170	115	120	95	92	4	FALSE

A data.frame: 2 × 13

	number	name	type1	type2	total	hp	attack	defense	sp_attack	sp_defense	speed	generation	legendary
	<int>	<chr>	<chr>	<chr>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<lgl>
997	844	Sandaconda	Ground		510	72	107	125	65	70	71	8	FALSE
998	844	Gigantamax Sandaconda	Ground		510	72	107	125	65	70	71	8	FALSE

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.

dim(data[data$legendary,]) #Awai, 118 quand-même...
head(data[data$legendary,])

1. 118
2. 13

A data.frame: 6 × 13

	number	name	type1	type2	total	hp	attack	defense	sp_attack	sp_defense	speed	generation	legendary
	<int>	<chr>	<chr>	<chr>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<lgl>
195	144	Articuno	Ice	Flying	580	90	85	100	95	125	85	1	TRUE
196	144	Galarian Articuno	Psychic	Flying	580	90	85	85	125	100	95	8	TRUE
197	145	Zapdos	Electric	Flying	580	90	90	85	125	90	100	1	TRUE
198	145	Galarian Zapdos	Fighting	Flying	580	90	125	90	85	90	100	8	TRUE
199	146	Moltres	Fire	Flying	580	90	100	90	125	85	90	1	TRUE
200	146	Galarian Moltres	Dark	Flying	580	90	85	90	100	125	90	8	TRUE

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

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

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.

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

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)

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.

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

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

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")

cl3_h <- cutree(h, 3)
fviz_cluster(list(data=data_clust, cluster=cl3_h))

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.

data[which(data$name %in% c("Shuckle", "Deoxys Attack Forme")),]

A data.frame: 2 × 13

	number	name	type1	type2	total	hp	attack	defense	sp_attack	sp_defense	speed	generation	legendary
	<int>	<chr>	<chr>	<chr>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<lgl>
273	213	Shuckle	Bug	Rock	505	20	10	230	10	230	5	2	FALSE
475	386	Deoxys Attack Forme	Psychic		600	50	180	20	180	20	150	3	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.

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)

cl3_p <- pam(data_clust, 3)
fviz_cluster(cl3_p)
cl3_k <- kmeans(data_clust, 3)
fviz_cluster(cl3_k, data=data_clust)

par(mfrow=c(1,2))
hist(cl3_p$clustering)
hist(cl3_k$cluster)

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

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

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"))

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.

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"))

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.

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]

# 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])

# 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()

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.

cl <- function(data, k) pam(data, k)$clustering
plotAll() #temps de calcul déjà beaucoup (beaucoup) plus long ! ...comme évoqué en cours

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

cl <- function(data, k) cutree(hclust(dist(data), method="ward.D"), k)
plotAll()

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 :

cl <- function(data, k) cutree(hclust(dist(data), method="single"), k)
plotAll()

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.