13.4 Nuées dynamiques
Les méthodes des nuées dynamiques regroupent plusieurs algorithmes, tous plus ou moins liés avec l’algorithme le plus connu : k-means, originalement proposé par James MacQueen (1967). Nous présentons également ici plusieurs variantes du k-means, soit le k-medians, le k-medioids, le c-means et le c-medians.
13.4.1 K-means
13.4.1.1 Fonctionnement de l’algorithme
Nous commençons ici par détailler le fonctionnement de cet algorithme afin de mieux le cerner. D’emblée, cet algorithme nécessite que certains éléments soient définis d’avance :
- Une matrice de données X comportant n lignes (nombre d’observations) et p colonnes (nombre de variables). Chaque variable de cette matrice doit être quantitative et continue et de préférence dans une échelle standardisée (par exemple des variables centrées réduites).
- Le nombre de groupes à identifier k doit être choisi par l’utilisateur ou l’utilisatrice.
- La distance d à utiliser entre les observations.
Le fonctionnement classique du k-means est le suivant :
- Définir k centres de groupes de façon aléatoire.
- Déterminer pour chaque observation le centre de son groupe le plus proche en utilisant la fonction de distance.
- Pour chacun des groupes ainsi formés, recalculer le centre du groupe en calculant le centroïde (moyennes le plus souvent) des observations appartenant à ce groupe.
- Répéter l’opération 2 avec les nouveaux centres.
- Calculer l’inertie expliquée par la nouvelle classification.
- Comparer cette inertie expliquée avec celle obtenue lors de l’itération précédente.
- Si la variation entre les deux valeurs est supérieure à une certaine limite, reprendre à l’étape 2, sinon, l’algorithme prend fin.
Ainsi, l’algorithme k-means part d’une première classification obtenue aléatoirement et la raffine jusqu’au point où l’amélioration de la classification devient négligeable. Du fait de ce point de départ aléatoire, cet algorithme est dit heuristique, car deux exécutions risquent de ne pas donner exactement le même résultat. Par conséquent, en relaçant l’algorithme, vous pourriez obtenir des résultats légèrement différents, avec par exemple des groupes similaires, mais obtenus dans un autre ordre, le groupe 1 étant devenu le groupe 3 et vice-versa. Il est aussi possible d’obtenir des résultats radicalement différents d’une tentative à l’autre, ce qui signifie que les groupes formés sont très instables et ne sont pas représentatifs de la population étudiée.
Réplicabilité des résultats dans R
Lorsqu’une méthode heuristique ou faisant appel au hasard est utilisée dans R, il est nécessaire de s’assurer que les résultats sont reproductibles. Cela permet notamment de relancer le même code et de réobtenir exactement les mêmes résultats : l’idée étant de figer le hasard.
Ultimement, un programme informatique est incapable de générer un résultat véritablement aléatoire, car il ne fait que suivre une suite d’opérations prédéterminées. Pour générer des résultats qui ressemblent au hasard, des algorithmes ont été proposés, partant d’une configuration initiale et appliquant une série d’opérations complexes permettant de générer des nombres semblant se distribuer aléatoirement. Si nous connaissons le point de départ de la suite d’opérations et que nous réappliquons ces dernières, alors nous sommes certains d’obtenir le même résultat. Il est possible, dans R, de définir un état initial de hasard à l’aide de la fonction set.seed
. Avec ce point de départ défini, nous sommes certains d’obtenir les mêmes résultats en relançant les mêmes opérations.
Prenons un exemple concret en sélectionnant aléatoire 3 chiffres dans un vecteur allant de 1 à 10.
<- 1:10
vec
# prenons un premier échantillon
sample(vec, size = 3)
## [1] 4 9 3
# et un second échantillon
sample(vec, size = 3)
## [1] 6 10 8
Nous obtenons bien deux échantillons différents. Recommençons en utilisant la fonction set.seed
pour obtenir cette fois-ci des résultats identiques.
<- 1:10
vec
# prenons un premier échantillon
set.seed(123)
sample(vec, size = 3)
## [1] 3 10 2
# et un second échantillon
set.seed(123)
sample(vec, size = 3)
## [1] 3 10 2
# prenons un troisème échantillon
set.seed(4568997)
sample(vec, size = 3)
## [1] 5 6 7
# et un quatrième échantillon
set.seed(4568997)
sample(vec, size = 3)
## [1] 5 6 7
Vous constatez que nous utilisons cette fonction plusieurs fois au cours de cette section. Elle nous permet de nous assurer que les résultats obtenus ne changent pas entre le moment où nous écrivons le livre et le moment où nous le formatons. Sinon, le texte pourrait ne plus être en phase avec les images ou les tableaux.
Pour mieux comprendre le fonctionnement du k-means, nous proposons ici une visualisation de ses différentes itérations.
Nous pouvons constater que, pour ce jeu de données relativement simple, l’algorithme converge très rapidement et que sa solution varie peu au-delà de la troisième itération. Si vous utilisez la version HTLM du livre, la figure 13.18 devrait être animée et illustrer pourquoi le k-means est aussi appelé algorithme de nuées dynamiques.
Centre de groupe et k-means
À nouveau, puisque chaque itération du k-means nécessite de recalculer les centres des groupes formés, des problèmes peuvent être rencontrés avec certains types de distance. C’est pourquoi il est recommandé d’utiliser la distance euclidienne avec le k-means original. Si des distances plus complexes doivent être utilisées, il est préférable d’utiliser la classification ascendante hiérarchique.
13.4.1.2 Choix du nombre optimal de groupes
Comme pour la CAH, le principal enjeu avec le k-means est de déterminer le nombre idéal de groupes pour effectuer la classification. Si ce nombre n’est pas connu à l’avance et qu’aucune forte justification théorique n’existe, il est possible d’utiliser les mêmes techniques que pour la CAH, soit la méthode du coude, l’indicateur de silhouette ou la méthode GAP.
13.4.2 K-médianes
Le k-medians est une variante du k-means. Contrairement au k-means privilégiant la distance euclidienne, le k-medians est à utiliser en priorité avec une distance de Manhattan. En effet, le centre d’un groupe n’est pas déterminé comme la moyenne des variables des observations appartenant à ce groupe (k-means), mais comme la médiane pour chaque variable (k-medians). En dehors de ces deux spécificités, il reprend le fonctionnement décrit plus haut pour le k-means. Il est particulièrement pertinent de l’utiliser quand un jeu de données comprend un très grand nombre de colonnes, car dans ce contexte, la distance euclidienne peine à représenter les différences entre les observations. De plus, l’utilisation de la médiane le rend moins sensible aux valeurs extrêmes.
13.4.3 K-médoïds
Le k-médoïds est également une variante du k-means. Le k-means crée des groupes en cherchant les centres de ces groupes dans l’espace multidimensionnel des données. Ces centres de groupes peuvent très bien ne pas correspondre à un point du jeu de données, au même titre que la moyenne d’une variable ne coïncide que rarement avec une observation réelle de cette variable. Pour le k-médoïds, les groupes sont formés en cherchant les centres de ces groupes parmi les observations du jeu de données. Ainsi, chaque groupe est centré sur une observation réelle, la plus similaire à l’ensemble des observations du groupe.
L’algorithme effectue les opérations suivantes :
- Sélectionner aléatoirement k observations du jeu de données, elles constituent les centres des groupes initiaux.
- Attribuer chaque observation au centre du groupe le plus proche.
- Tant que la nouvelle solution est plus efficace, effectuer les opérations suivantes :
- pour chaque centre m et pour chaque observation o,
- considérer l’inversion de m et o
- si cette permutation est meilleure que les précédentes, la conserver en mémoire
- effectuer la meilleure permutation retenue si elle améliore la classification, sinon l’algorithme prend fin.
Le k-médoïds est moins utilisé que le k-means, mais il est plus performant quand des distances autres que la distance euclidienne sont utilisées ou encore que des valeurs aberrantes/extrêmes sont présentes dans les données.
13.4.4 Mise en oeuvre dans R
Pour cet exemple, nous proposons d’utiliser le jeu de données spatiales LyonIris
du package geocmeans
. Ce jeu de données spatiales pour l’agglomération lyonnaise (France) comprend dix variables, dont quatre environnementales (EN) et six socioéconomiques (SE), pour les îlots regroupés pour l’information statistique (IRIS) de l’agglomération lyonnaise (tableau 13.7 et figure 12.4). Nous proposons de réaliser une analyse similaire à celle de l’article de Gelb et Apparicio (2021b), soit de classer les IRIS de Lyon selon ces caractéristiques pour déterminer si certains groupes d’IRIS combinent des situations désavantageuses sur les plans sociaux et environnementaux, dans une perspective d’équité environnementale.
Notez ici que la fonction st_drop_geometry
provenant du package sf
permet de retirer l’information géographique du jeu de données LyonIris
pour obtenir un simple dataframe
.
Nom | Intitulé | Type | Moy. | E.-T. | Min. | Max. |
---|---|---|---|---|---|---|
Lden | Bruit routier (Lden dB(A)) | EN | 55,6 | 4,9 | 33,9 | 71,7 |
NO2 | Dioxyde d’azote (ug/m3) | EN | 28,7 | 7,9 | 12,0 | 60,2 |
PM25 | Particules fines (PM\(_{2,5}\)) | EN | 16,8 | 2,1 | 11,3 | 21,9 |
VegHautPrt | Canopée (%) | EN | 18,7 | 10,1 | 1,7 | 53,8 |
Pct0_14 | Moins de 15 ans (%) | SE | 18,5 | 5,7 | 0,0 | 54,0 |
Pct_65 | 65 ans et plus (%) | SE | 16,2 | 5,9 | 0,0 | 45,1 |
Pct_Img | Immigrants (%) | SE | 14,5 | 9,1 | 0,0 | 59,8 |
TxChom1564 | Taux de chômage | SE | 14,8 | 8,1 | 0,0 | 98,8 |
Pct_brevet | Personnes à faible scolarité (%) | SE | 23,5 | 12,6 | 0,0 | 100,0 |
NivVieMed | Médiane du niveau de vie (Euros) | SE | 21 804,5 | 4 922,5 | 11 324,0 | 38 707,0 |
13.4.4.1 Préparation des données
La première étape consiste donc à charger les données et à les préparer pour l’analyse. Toutes les variables que nous utilisons sont des variables continues. Cependant, elles ne sont pas exprimées dans la même échelle, nous proposons donc de les standardiser ici en les centrant (moyenne = 0) et en les réduisant (écart-type = 1). Cette opération peut être effectuée simplement dans R en utilisant la fonction scale
.
# Chargement des données
library(geocmeans)
library(sf)
data(LyonIris)
# NB : LyonIris est un objet spatial, il faut donc extraire uniquement son DataFrame
<- st_drop_geometry(LyonIris[c("Lden","NO2","PM25","VegHautPrt","Pct0_14","Pct_65","Pct_Img",
X "TxChom1564","Pct_brevet","NivVieMed")])
# Centrage et réduction de chaque colonne du DataFrame
for (col in names(X)){
<- scale(X[[col]], center = TRUE, scale = TRUE)
X[[col]] }
13.4.4.2 Choix du nombre de groupes optimal
La seconde étape consiste à déterminer le nombre de groupes optimal. Pour cela, nous comparons les résultats des trois méthodes proposées : la méthode du coude, l’indice de silhouette et la méthode GAP. Pour chaque méthode, nous testons les nombres de groupes de 2 à 10.
13.4.4.2.1 Méthode du coude
Commençons par appliquer la méthode du coude. Nous calculons donc l’inertie expliquée par la classification pour différentes valeurs de k (nombre de groupes) avant de construire la figure 13.19.
<- 2:10
ks
## ---- Méthode du coude ---- ##
<- sapply(ks, function(k){
inertie_exps # calcul du kmeans avec k
<- kmeans(X, centers = k)
resultat # calcul de l'inertie expliquée (1 - inertie intragroupe / inertie totale)
<- 1-(sum(resultat$withinss) / resultat$totss)
inertie_exp return(inertie_exp)
})
<- data.frame(
df k = ks,
inertie_exp = inertie_exps
)
ggplot(df) +
geom_line(aes(x = k, y = inertie_exp)) +
geom_point(aes(x = k, y = inertie_exp), color = "red") +
labs(x = "nombre de groupes", y = "inertie expliquée (%)")
Dans l’article original, quatre groupes avaient été retenus. Nous pouvons constater ici qu’un coude fort se situe à k = 3 et qu’au-delà de cette limite, l’ajout d’un groupe supplémentaire contribue à expliquer une plus petite partie de l’inertie supplémentaire comparativement au précédent.
13.4.4.2.2 Indice de silhouette
Poursuivons avec l’indice de silhouette calculé de nouveau avec des valeurs de k allant de 2 à 10. Notez que nous devons au préalable créer une matrice de distances entre les observations du jeu de données pour construire notre indice de silhouette. Puisque nous utilisons l’algorithme k-means, nous utilisons la distance euclidienne.
<- 2:10
ks
# calcul d'une matrice de distance euclidienne entre les observations
<- dist(X, method = "euclidean")
dist_mat
## ---- indice de silhouette ---- ##
<- sapply(ks, function(k){
values <- kmeans(X,centers = k)
resultat <- resultat$cluster
groupes # calcul des valeurs de silhouette
<- silhouette(groupes, dist = dist_mat)
sil # extraction de l'indice global (moyenne des moyennes)
<- mean(summary(sil)$clus.avg.widths)
idx return(idx)
})
<- data.frame(
df k = ks,
silhouette = values
)
ggplot(df) +
geom_line(aes(x = k, y = silhouette)) +
geom_point(aes(x = k, y = silhouette), color = "red") +
labs(x = "nombre de groupes", y = "Indice de silhouette")
À nouveau, la figure 13.20 indique que le nombre de groupes optimal est trois selon l’indice de silhouette.
13.4.4.2.3 Méthode GAP
Pour appliquer la méthode GAP, nous proposons d’utiliser la fonction clusGap
du package NbClust
. Pour l’utiliser, il est nécessaire de définir une fonction renvoyant pour le nombre de groupes k et le jeu de données x une liste comprenant un vecteur attribuant chaque observation à chaque groupe. Il est possible de considérer ce type de fonction comme un « adaptateur ».
library(NbClust)
# définition de la fonction adaptateur
<- function(x,k){
adaptor <- kmeans(x,k)
clust return(list(
"cluster" = clust$cluster
))
}
# calcul de la méthode GAP
<- clusGap(X, adaptor, K.max = 10, verbose = FALSE)
vals <- data.frame(vals$Tab)
tab $k <- 1:nrow(tab)
tab
# détermination des valeurs de k retenues par la méthode (1ere et 2e)
<- sapply(2:nrow(tab), function(i){
is_valid -1,"gap"] >= (tab[i,"gap"] - tab[i,"SE.sim"])
tab[i
})<- subset(tab,is_valid)[1,]
valids <- subset(tab,is_valid)[2,]
valids2
# réalisation du graphique
ggplot(tab) +
geom_line(aes(x = k, y = gap)) +
geom_segment(x = valids$k, xend = valids$k, y = min(tab$gap), yend = valids$gap,
linetype = "dashed") +
geom_segment(x = valids2$k, xend = valids2$k, y = min(tab$gap), yend = valids2$gap,
linetype = "dashed") +
geom_point(aes(x = k, y = gap), color = 'red') +
scale_x_continuous(breaks = 1:10) +
labs(x = "nombre de groupes", y = "GAP")
La figure 13.21 indique également que le nombre de groupes à retenir est trois. Nous retenons cependant quatre groupes pour pouvoir plus facilement comparer nos résultats avec ceux de l’article original.
13.4.4.3 Application l’algorithme du k-means
Maintenant que nous avons choisi le nombre de groupes à former, nous pouvons simplement appliquer la fonction kmeans
présente de base dans R.
set.seed(145)
<- kmeans(X, centers = 4) resultats
13.4.4.4 Interprétation des résultats
Une fois les groupes obtenus, l’étape la plus importante est de parvenir à interpréter ces groupes. Pour cela, il est nécessaire de les explorer en profondeur au travers des variables utilisées pour les constituer. Dans notre cas, le jeu de données LyonIris
est spatialisé, nous pouvons donc commencer par cartographier les groupes.
library(tmap)
$groupes <- paste("groupe", resultats$cluster, sep = " ")
LyonIris
tm_shape(LyonIris) +
tm_polygons(col = "groupes", palette =
c("#EFBE89", "#4A6A9F", "#7DB47C", "#FAF29C"),lty = 1, lwd = 0.1)
Il est ainsi possible de constater que le groupe 3 forme un ensemble assez compact d’IRIS au centre de Lyon. Le groupe 4 correspond quant à lui à des IRIS situés en périphérie plutôt éloignée, essentiellement à l’ouest. Le groupe 1 correspond à une périphérie proche du groupe 2 et apparaît comme un ensemble d’enclaves dispersées.
Pour distinguer rapidement les profils des différents groupes, il est possible d’utiliser un graphique en radar. La construction d’un tel graphique peut être un peu fastidieuse dans R, cependant le package geocmeans
propose une fonction assez pratique : spiderPlots
.
library(geocmeans)
# création d'une matrice d'appartenance binaire des groupes
<- fastDummies::dummy_cols(resultats$cluster, remove_selected_columns = TRUE)
matrice_gp
# réalisation du graphique
par(mfrow=c(3,2), mai = c(0.1,0.1,0.1,0.1))
<- spiderPlots(X, matrice_gp,
plots chartcolors = c("#EFBE89", "#4A6A9F", "#7DB47C", "#FAF29C"))
Il est ainsi possible de constater, à la figure 13.23, que le groupe 3 est caractérisé par un niveau de vie élevé, mais par des niveaux de concentration de pollution atmosphérique plus élevés également. Le groupe 4 en revanche est caractérisé par un important couvert végétal, un niveau de vie médian élevé et une plus forte proportion de personnes de plus de 65 ans. Le groupe 1 est quant à lui marqué par des niveaux sonores plus élevés. Enfin, le groupe 2 se caractérise par une plus grande proportion de population ayant obtenu comme diplôme le plus élevé le brevet des collèges, d’immigrants, de jeunes de moins de 15 ans et un taux de chômage plus élevé.
Notez que ces graphiques nous permettent rapidement de nous faire une idée des caractéristiques des groupes, mais uniquement sur une échelle relative. En effet, ils ne nous indiquent à aucun moment la taille des écarts entre les groupes. Pour cela, il est nécessaire de réaliser des graphiques en violon pour chaque variable. Pour ce type de graphique, il est préférable d’utiliser les données originales non transformées pour pouvoir mieux appréhender si les différences entre les groupes sont importantes ou négligeables.
<- st_drop_geometry(LyonIris[c("Lden","NO2","PM25","VegHautPrt","Pct0_14",
X2 "Pct_65","Pct_Img","TxChom1564","Pct_brevet",
"NivVieMed")])
<- violinPlots(X2, as.character(resultats$cluster))
plots ggarrange(plotlist = plots, ncol = 2, nrow = 5)
Il est également recommandé de calculer des statistiques descriptives par groupe et de les rapporter dans un tableau.
# obtention d'un tableau par groupe
<- summarizeClusters(X2, matrice_gp, dec = 1, silent = TRUE)
tableaux
# concaténation des tableaux
<- do.call(rbind, tableaux) tableau_tot
Lden | NO2 | PM25 | VegHautPrt | Pct014 | Pct65 | PctImg | TxChom1564 | Pctbrevet | NivVieMed | |
---|---|---|---|---|---|---|---|---|---|---|
groupe 1 | ||||||||||
Q5 | 53,8 | 25,2 | 15,6 | 6,6 | 12,4 | 10,0 | 8,1 | 7,6 | 17,4 | 15 822,7 |
Q10 | 54,5 | 26,4 | 15,9 | 8,0 | 13,9 | 11,4 | 9,2 | 9,8 | 18,1 | 16 976,0 |
Q25 | 56,3 | 29,3 | 17,0 | 11,0 | 16,3 | 13,6 | 11,8 | 11,5 | 20,8 | 18 454,0 |
Q50 | 58,9 | 32,3 | 18,2 | 15,3 | 18,3 | 16,6 | 15,9 | 13,7 | 24,1 | 19 559,0 |
Q75 | 62,5 | 36,5 | 18,7 | 22,3 | 20,5 | 19,5 | 18,6 | 16,9 | 30,0 | 21 575,2 |
Q90 | 64,7 | 39,3 | 19,0 | 30,3 | 22,7 | 22,8 | 21,1 | 19,8 | 33,0 | 23 544,7 |
Q95 | 66,7 | 40,7 | 19,3 | 35,6 | 24,9 | 25,1 | 23,0 | 21,9 | 37,7 | 24 765,7 |
Mean | 59,6 | 32,7 | 17,8 | 17,2 | 18,2 | 16,7 | 15,7 | 14,1 | 25,5 | 19 948,0 |
Std | 4,2 | 5,2 | 1,2 | 8,5 | 4,0 | 4,7 | 6,5 | 4,3 | 7,6 | 2 637,2 |
groupe 2 | ||||||||||
Q5 | 50,8 | 19,3 | 13,9 | 6,1 | 18,3 | 6,4 | 17,7 | 16,5 | 30,5 | 12 323,5 |
Q10 | 52,0 | 20,2 | 14,3 | 7,8 | 20,0 | 8,7 | 20,2 | 16,8 | 32,9 | 12 747,0 |
Q25 | 53,9 | 23,0 | 15,3 | 11,2 | 22,8 | 10,7 | 23,5 | 19,6 | 36,2 | 13 530,5 |
Q50 | 56,4 | 25,1 | 16,2 | 14,5 | 25,2 | 13,6 | 28,0 | 24,5 | 39,9 | 15 340,0 |
Q75 | 58,4 | 29,4 | 17,0 | 18,1 | 27,8 | 16,9 | 33,5 | 32,5 | 46,0 | 16 342,0 |
Q90 | 63,1 | 33,7 | 18,5 | 24,4 | 31,3 | 20,2 | 38,2 | 36,3 | 50,2 | 18 363,9 |
Q95 | 64,9 | 39,0 | 19,1 | 28,1 | 32,8 | 21,4 | 41,1 | 38,0 | 60,2 | 20 128,6 |
Mean | 56,8 | 26,3 | 16,2 | 15,4 | 25,4 | 13,8 | 28,5 | 26,6 | 42,0 | 15 401,4 |
Std | 4,1 | 6,3 | 1,5 | 6,8 | 6,2 | 4,8 | 8,0 | 10,5 | 11,3 | 2 340,4 |
groupe 3 | ||||||||||
Q5 | 50,2 | 28,3 | 17,0 | 5,0 | 6,8 | 5,2 | 5,8 | 7,7 | 6,7 | 19 029,9 |
Q10 | 51,2 | 29,5 | 17,6 | 6,6 | 9,5 | 7,2 | 6,7 | 8,5 | 7,7 | 19 546,4 |
Q25 | 53,2 | 31,3 | 18,6 | 9,2 | 11,4 | 9,5 | 7,9 | 11,0 | 9,5 | 21 652,0 |
Q50 | 55,2 | 35,4 | 19,3 | 12,6 | 14,1 | 12,4 | 11,0 | 12,9 | 12,0 | 23 342,0 |
Q75 | 58,0 | 39,6 | 19,8 | 16,0 | 16,2 | 16,0 | 13,1 | 15,1 | 14,8 | 25 954,1 |
Q90 | 60,0 | 43,0 | 20,1 | 21,0 | 18,2 | 18,9 | 16,3 | 18,0 | 16,6 | 28 951,0 |
Q95 | 61,2 | 44,4 | 20,4 | 29,3 | 19,5 | 20,9 | 17,7 | 19,3 | 19,1 | 31 922,5 |
Mean | 55,6 | 35,8 | 19,0 | 13,6 | 13,8 | 12,6 | 11,1 | 13,1 | 12,1 | 23 999,7 |
Std | 3,7 | 5,6 | 1,0 | 6,8 | 4,2 | 4,7 | 4,6 | 4,4 | 4,0 | 3 870,3 |
groupe 4 | ||||||||||
Q5 | 44,9 | 14,7 | 12,7 | 8,2 | 13,1 | 12,7 | 3,8 | 6,7 | 10,8 | 19 496,0 |
Q10 | 46,6 | 15,7 | 13,0 | 11,8 | 14,9 | 13,5 | 4,4 | 7,0 | 12,4 | 20 257,0 |
Q25 | 50,3 | 19,0 | 13,8 | 16,7 | 17,1 | 16,1 | 6,0 | 8,0 | 15,7 | 21 963,0 |
Q50 | 52,6 | 22,0 | 14,7 | 24,9 | 18,9 | 19,3 | 8,1 | 9,8 | 20,0 | 24 109,8 |
Q75 | 54,9 | 25,2 | 15,5 | 32,6 | 20,8 | 22,8 | 10,9 | 12,0 | 25,6 | 26 698,0 |
Q90 | 57,7 | 27,5 | 16,3 | 41,2 | 22,3 | 27,4 | 14,5 | 14,7 | 30,6 | 29 947,0 |
Q95 | 59,1 | 28,7 | 16,8 | 43,6 | 22,7 | 29,7 | 18,1 | 16,6 | 33,1 | 32 386,5 |
Mean | 52,3 | 21,8 | 14,7 | 25,5 | 18,6 | 20,1 | 8,8 | 10,4 | 21,1 | 24 761,6 |
Std | 4,3 | 4,4 | 1,2 | 11,0 | 3,1 | 5,6 | 4,2 | 3,4 | 7,6 | 4 008,7 |
Les constats que nous avons faits précédemment sont confirmés par la figure 13.24 et le tableau 13.8. Nous retrouvons ici les groupes originaux décrits dans l’article de Gelb et Apparicio (2021b) :
- Groupe 1 : les espaces interstitiels, formant une périphérie proche du centre et relativement hétérogène sur les variables étudiées, mais caractérisée par des niveaux de bruit importants.
- Groupe 2 : les banlieues jeunes et défavorisées, avec des niveaux d’exposition aux pollutions atmosphérique et sonore relativement élevés comparativement à l’ensemble de la région.
- Groupe 3 : les quartiers centraux aisés, mais marqués par les plus hauts niveaux de pollution atmosphérique.
- Groupe 4 : les communes rurales, aisées et vieillissantes.
Interprétation interactive
Si, comme dans notre exemple, vos données comportent une dimension spatiale, le package geocmeans
propose une fonction intéressante appelée sp_clust_explorer
démarrant une application permettant d’explorer les résultats de votre classification. Le seul enjeu est de créer un objet de la classe FCMres
. Voici un court exemple :
# création d'une matrice binaire d'appartenance
<- dummy_cols(resultats$cluster, remove_selected_columns = TRUE)
kmeans_mat
# extraction des centres de notre classification
<- resultats$centers
centres
# création de l'objet FCMres
<- FCMres(list(
kmeansres "Centers" = centres,
"Belongings" = kmeans_mat,
"Data" = X2,
"m" = 1,
"algo" = "kmeans"
))
# démarrage de l'application shiny
sp_clust_explorer(object = kmeansres, spatial = LyonIris)
13.4.4.5 K-médianes et K-médoides
Nous présentons simplement ici comment effectuer la même analyse en utilisant les variantes du k-means, soit le k-medians et le k-mediods.
Il existe relativement peu d’implémentation du k-medians dans R, nous optons donc ici pour la fonction kGmedian
du package Gmedian
. Pour le k-mediods, nous avons retenu la fonction pam
du package cluster
.
library(Gmedian)
<- kGmedian(X, 4)
k_median_res
library(cluster)
<- pam(X,4) k_mediods_res
Juste pour le plaisir des yeux, nous pouvons cartographier les trois classifications obtenues en nous assurant au préalable de faire coïncider les groupes les plus similaires de nos trois classifications.
<- dummy_cols(resultats$cluster,
matrice_gp_kmeans remove_selected_columns = TRUE)
<- dummy_cols(as.vector(k_median_res$cluster),
matrice_gp_kmedians remove_selected_columns = TRUE)
<- dummy_cols(k_mediods_res$cluster,
matrice_gp_kmedioids remove_selected_columns = TRUE)
# Appariement des groupes du k-medians avec ceux du kmeans
<- geocmeans::groups_matching(as.matrix(matrice_gp_kmeans),
matrice_gp_kmedians as.matrix(matrice_gp_kmedians))
# Appariement des groupes du k-medioids avec ceux du kmeans
<- geocmeans::groups_matching(as.matrix(matrice_gp_kmeans),
matrice_gp_kmedioids as.matrix(matrice_gp_kmedioids))
# ajouts des colonnes nécessaires à LyonIris
colnames(matrice_gp_kmeans) <- paste0("groupe_", 1:4)
colnames(matrice_gp_kmedians) <- paste0("groupe_", 1:4)
colnames(matrice_gp_kmedioids) <- paste0("groupe_", 1:4)
$kmeans <- colnames(matrice_gp_kmeans)[max.col(matrice_gp_kmeans)]
LyonIris$kmedians <- colnames(matrice_gp_kmedians)[max.col(matrice_gp_kmedians)]
LyonIris$kmedioids <- colnames(matrice_gp_kmedioids)[max.col(matrice_gp_kmedioids)]
LyonIris
# construction de la figure
<- c("#EFBE89", "#4A6A9F", "#7DB47C", "#FAF29C")
couleurs
<- tm_shape(LyonIris) +
map1 tm_polygons(col = "kmeans",palette = couleurs,lty = 1, lwd = 0.1)
<- tm_shape(LyonIris) +
map2 tm_polygons(col = "kmedians",palette = couleurs,lty = 1, lwd = 0.1)
<- tm_shape(LyonIris) +
map3 tm_polygons(col = "kmedioids",palette = couleurs,lty = 1, lwd = 0.1)
tmap_arrange(map1, map2, map3,
ncol = 2, nrow = 2)
Les trois cartes sont très similaires (figure 13.25), ce qui signifie que les trois algorithmes tendent à attribuer les observations aux mêmes groupes. Cependant, nous observons des différences, notamment au nord avec des observations alternant entre les groupes 2 et 3 selon la méthode employée. Cela peut notamment signifier que ces observations sont « indécises », qu’il est difficile de les attribuer définitivement à une catégorie en particulier. Pour prendre en compte cette forme d’incertitude, il est possible d’opter pour des méthodes de classification en logique floue.
13.4.5 Extensions en logique floue : c-means, c-medoids
Comme nous l’avons mentionné en introduction de cette section, les méthodes de classification floues ont pour objectif d’évaluer le degré d’appartenance de chaque observation à chaque groupe plutôt que d’attribuer chaque observation à un seul groupe. Il est ainsi possible de repérer des observations incertaines, à cheval entre plusieurs groupes. Nous présentons ici deux algorithmes appartenant à cette famille : le c-means et le c-medoids. Il s’agit dans les deux cas d’extensions des k-means et k-medoids vus précédemment.
Pour ces deux méthodes, comme pour le k-means, le nombre de groupes k doit être spécifié. Elles comprennent cependant un paramètre supplémentaire : m, appelé paramètre de floutage qui contrôle à quel point le résultat obtenu sera flou ou strict. Une valeur de 1 produit une classification stricte (chaque observation appartient à un seul groupe) et une valeur plus grande conduit à des classifications de plus en plus floues, jusqu’à ce que chaque observation appartienne à un degré identique à chacun des groupes. Il est recommandé de sélectionner m en même temps que k, car ces deux valeurs influencent simultanément la qualité de la classification. La meilleure approche consiste à tester un ensemble de combinaisons de m et de k et à comparer les valeurs obtenues pour différents indicateurs de qualité de classification floue. Parmi ces indicateurs, il est notamment recommandé d’utiliser le pourcentage de l’inertie expliquée, l’indice de silhouette pour classification floue, l’indice de Xie et Beni (1991), et de Fukuyama et Sugeno (Fukuyama 1989).
13.4.5.1 Mise en oeuvre du c-means dans R
Le package fclust
comprend un très grand nombre de méthodes pour effectuer des classifications floues, nous l’utilisons donc en priorité ici en combinaison avec des fonctions d’interprétation du package geocmeans
.
13.4.5.1.1 Préparation des données
Comme pour le k-means, cette méthode nécessite de disposer d’un jeu de données ne comprenant que des variables quantitatives dans la même échelle. Nous commençons donc à nouveau par standardiser nos données. Pour varier les plaisirs, nous optons cette fois-ci pour une transformation des variables dans une échelle allant de 0 à 100.
library(fclust)
data(LyonIris)
# NB : LyonIris est un objet spatial, il faut donc extraire uniquement son DataFrame
<- st_drop_geometry(LyonIris[c("Lden","NO2","PM25","VegHautPrt","Pct0_14",
X "Pct_65","Pct_Img",
"TxChom1564","Pct_brevet","NivVieMed")])
# changement d'échelle des données (0 à 100)
<- function(x){
to_0_100 return((x-min(x)) / (max(x) - min(x)) * 100)
}
for (col in names(X)){
<- to_0_100(X[[col]])
X[[col]] }
13.4.5.1.2 Sélection de k et de m
La seconde étape consiste à sélectionner les valeurs optimales pour k et m. Nous testons ici toutes les valeurs de k de 2 à 7, et les valeurs de m de 1,5 à 2,5 (avec des écarts de 0,1).
library(e1071)
set.seed(123)
<- seq(1.5,2.5,by = 0.1)
ms <- 2:7
ks
# calcul de toutes les combinaisons
<- expand.grid(ms,ks)
combinaisons
<- c("Explained.inertia", "Silhouette.index", "FukuyamaSugeno.index")
eval_indices
<- apply(combinaisons, MARGIN = 1, FUN = function(row){
values <- row[[1]]
m <- row[[2]]
k <- FKM(X, k, m)
resultats <- geocmeans::calcqualityIndexes(as.matrix(X),
idx as.matrix(resultats$U),
m = m,
indices = eval_indices)
return(c(k,m,unlist(idx)))
})
<- data.frame(t(values))
df_scores names(df_scores) <- c("k", "m", "inertie", "silhouette", "FukuyamaSugeno")
# changer l'échelle de l'indice pour un graphique plus joli
$FukuyamaSugeno <- round(df_scores$FukuyamaSugeno/10000,2)
df_scores
# création de trois figures pour représenter les trois indicateurs
library(viridis)
<- ggplot(df_scores) +
plot1 geom_raster(aes(x = k, y = m, fill = inertie)) +
scale_fill_viridis() +
scale_x_continuous(breaks = c(2,3,4,5,6,7)) +
coord_fixed(ratio=4) +
guides(fill = guide_colourbar(barwidth = 5, barheight = 0.5)) +
labs(fill = "Inertie expliquée") +
theme(legend.position = "bottom", legend.box = "horizontal",
legend.title = element_text( size=9), legend.text=element_text(size=8))
<- ggplot(df_scores) +
plot2 geom_raster(aes(x = k, y = m, fill = silhouette)) +
scale_fill_viridis() +
scale_x_continuous(breaks = c(2,3,4,5,6,7)) +
coord_fixed(ratio=4) +
guides(fill = guide_colourbar(barwidth = 5, barheight = 0.5)) +
labs(fill = "Indice de silhouette") +
theme(legend.position = "bottom", legend.box = "horizontal",
legend.title = element_text( size=9), legend.text=element_text(size=8))
<- ggplot(df_scores) +
plot3 geom_raster(aes(x = k, y = m, fill = FukuyamaSugeno)) +
scale_fill_viridis() +
scale_x_continuous(breaks = c(2,3,4,5,6,7)) +
coord_fixed(ratio=4) +
guides(fill = guide_colourbar(barwidth = 5, barheight = 0.5)) +
labs(fill = "Indice de Fukuyama et Sugeno") +
theme(legend.position = "bottom", legend.box = "horizontal",
legend.title = element_text( size=9), legend.text=element_text(size=8))
ggarrange(plot1, plot2, plot3, ncol = 2, nrow = 2)
Les trois graphiques à la figure 13.26 semblent indiquer des solutions différentes. Sans surprise, augmenter le niveau de flou (m) réduit l’inertie expliquée, alors qu’augmenter le nombre de groupes (k) augmente l’inertie expliquée. L’indice de silhouette indique assez clairement que trois groupes serait le meilleur choix, suivi par deux ou quatre groupes, si m est inférieur à 1,8. Cependant, ne retenir que trois groupes ne permet d’expliquer que 30% de l’inertie. Afin de nous rapprocher des résultats de l’article original (Gelb et Apparicio 2021b), nous retenons m = 1,5
et k = 4
.
13.4.5.1.3 Application l’algorithme c-means
Avec k et m définis, il ne reste plus qu’à appliquer l’algorithme à nos observations.
set.seed(123)
<- FKM(X, 4, 1.5) cmeans_resultats
L’objet obtenu cmeans_resultats
contient les résultats de la classification. Plus spécifiquement, cmeans_resultats$U
est la matrice d’appartenance, soit une matrice de taille n x k, dont chaque case \(U_{ij}\) indique la probabilité pour l’observation i d’appartenir au groupe j. cmeans_resultats$H
contient le centre des groupes, et cmeans_resultats$Clus
, le groupe auquel chaque observation à le plus de chances d’appartenir. Pour comparer plus facilement nos résultats avec ceux du k-means, nous pouvons changer l’ordre des groupes obtenus pour les faire coïncider avec les groupes les plus similaires obtenus avec la méthode k-means.
# changeons l'ordre des groupes
<- cmeans_resultats$U
U <- geocmeans::groups_matching(as.matrix(matrice_gp_kmeans), as.matrix(U))
U2
# mais aussi du centre des classes
<- as.integer(gsub("Clus ","",colnames(U2), fixed = TRUE))
idx <- cmeans_resultats$H[idx,]
H2
# et recalcul du groupe le plus probable
<- data.frame(
Clus2 "Cluster" = (1:4)[max.col(U2,ties.method="first")],
"Membership degree" = apply(U2, MARGIN = 1, max)
)
colnames(U2) <- paste("Clus",1:4, sep = " ")
rownames(H2) <- paste("Clus",1:4, sep = " ")
$U <- U2
cmeans_resultats$H <- H2
cmeans_resultats$Clus <- Clus2 cmeans_resultats
13.4.5.1.4 Interprétation des résultats
Globalement, les approches pour interpréter les résultats issus d’une classification obtenue par c-means sont les mêmes que pour une classification obtenue par k-means.
Commençons donc par créer plusieurs cartes des probabilités d’appartenir aux différents groupes.
<- mapClusters(LyonIris, cmeans_resultats$U)
maps
tmap_arrange(maps$ProbaMaps, ncol = 2, nrow = 2)
Sur les cartes de la figure 13.27, l’intensité de bleu correspond à la probabilité pour chaque IRIS d’appartenir aux différents groupes. Nous retrouvons les principales structures spatiales que nous avons identifiées avec le k-means; cependant, nous pouvons à présent constater que le groupe 1 est bien plus incertain que les autres. Nous pouvons une fois encore générer un graphique en radar pour comparer les profils des quatre groupes.
par(mfrow=c(3,2), mai = c(0.1,0.1,0.1,0.1))
spiderPlots(X, cmeans_resultats$U,
chartcolors = c("#EFBE89", "#4A6A9F", "#7DB47C", "#FAF29C"))
Sans surprise, nous retrouvons essentiellement les profils que nous avons obtenus avec le k-means dans la figure 13.28. Pour compléter la lecture des résultats, il est nécessaire de se pencher sur le tableau des statistiques descriptives des différents groupes. Une fois encore, nous proposons d’utiliser la fonction summarizeClusters
du package geocmeans
. Notez que cette fonction calcule les statistiques descriptives pondérées en fonction de l’appartenance des observations aux groupes. Ainsi, une observation ayant une faible chance d’appartenir à un groupe ne contribue que faiblement aux statistiques descriptives de ce groupe.
<- st_drop_geometry(LyonIris[c("Lden","NO2","PM25","VegHautPrt","Pct0_14",
df "Pct_65","Pct_Img", "TxChom1564",
"Pct_brevet","NivVieMed")])
<- summarizeClusters(data = df,
tableaux belongmatrix = cmeans_resultats$U,
weighted = TRUE, dec = 1)
<- do.call(rbind, tableaux) tableau_tot
Lden | NO2 | PM25 | VegHautPrt | Pct014 | Pct65 | PctImg | TxChom1564 | Pctbrevet | NivVieMed | |
---|---|---|---|---|---|---|---|---|---|---|
groupe 1 | ||||||||||
Q5 | 48,7 | 15,8 | 13,5 | 6,2 | 12,2 | 9,5 | 4,3 | 6,7 | 10,7 | 16 191,5 |
Q10 | 50,5 | 18,7 | 13,9 | 7,4 | 14,2 | 11,4 | 5,7 | 7,6 | 13,7 | 17 616,0 |
Q25 | 52,4 | 20,9 | 14,7 | 11,4 | 16,8 | 14,2 | 7,8 | 9,3 | 18,4 | 19 559,0 |
Q50 | 54,7 | 25,0 | 15,6 | 15,8 | 19,0 | 17,7 | 11,2 | 12,1 | 24,0 | 21 947,5 |
Q75 | 57,6 | 28,6 | 16,8 | 22,3 | 21,3 | 21,1 | 15,6 | 15,1 | 29,8 | 24 069,2 |
Q90 | 60,7 | 33,0 | 18,5 | 29,8 | 23,2 | 24,4 | 20,3 | 18,5 | 34,7 | 26 164,0 |
Q95 | 63,1 | 37,8 | 19,0 | 34,3 | 25,4 | 27,5 | 23,5 | 22,3 | 39,0 | 28 931,0 |
Mean | 55,1 | 25,4 | 15,8 | 17,5 | 18,9 | 17,9 | 12,4 | 13,0 | 24,7 | 21 951,6 |
Std | 4,4 | 6,3 | 1,7 | 8,6 | 4,5 | 5,7 | 6,8 | 6,8 | 10,3 | 3 744,1 |
groupe 2 | ||||||||||
Q5 | 50,8 | 19,6 | 14,0 | 6,4 | 12,9 | 7,5 | 8,4 | 9,4 | 15,5 | 12 392,0 |
Q10 | 52,1 | 21,4 | 14,6 | 7,7 | 16,6 | 8,8 | 12,9 | 13,0 | 22,0 | 12 920,0 |
Q25 | 54,5 | 23,2 | 15,8 | 10,8 | 20,1 | 11,1 | 19,6 | 16,8 | 30,5 | 14 078,0 |
Q50 | 57,1 | 26,6 | 16,6 | 14,5 | 23,6 | 13,9 | 25,8 | 22,0 | 37,5 | 15 985,0 |
Q75 | 59,5 | 31,9 | 17,8 | 19,0 | 27,0 | 17,4 | 32,1 | 30,1 | 44,2 | 18 566,2 |
Q90 | 63,3 | 37,4 | 18,8 | 25,4 | 29,7 | 20,7 | 37,8 | 33,9 | 48,6 | 21 146,0 |
Q95 | 64,7 | 39,6 | 19,3 | 30,7 | 32,3 | 23,3 | 39,4 | 37,7 | 52,9 | 23 712,5 |
Mean | 57,3 | 28,1 | 16,7 | 15,8 | 23,4 | 14,4 | 25,5 | 23,0 | 36,9 | 16 658,7 |
Std | 4,4 | 6,6 | 1,6 | 7,6 | 6,3 | 5,1 | 9,7 | 9,3 | 12,3 | 3 611,9 |
groupe 3 | ||||||||||
Q5 | 50,6 | 27,5 | 16,3 | 6,2 | 8,6 | 6,8 | 6,0 | 7,6 | 7,5 | 16 976,0 |
Q10 | 51,9 | 28,8 | 17,1 | 7,7 | 10,4 | 8,4 | 7,0 | 8,6 | 8,7 | 18 434,0 |
Q25 | 53,8 | 30,9 | 18,3 | 10,2 | 12,6 | 10,8 | 9,1 | 11,1 | 11,0 | 19 785,0 |
Q50 | 56,5 | 35,0 | 18,9 | 13,4 | 15,2 | 14,1 | 12,3 | 13,1 | 14,8 | 22 289,0 |
Q75 | 59,4 | 38,7 | 19,6 | 17,6 | 17,9 | 17,5 | 15,9 | 15,5 | 21,5 | 24 522,0 |
Q90 | 62,7 | 41,1 | 20,0 | 24,8 | 20,5 | 20,5 | 18,9 | 18,4 | 27,6 | 27 600,0 |
Q95 | 64,3 | 44,3 | 20,2 | 30,6 | 21,8 | 23,5 | 21,0 | 20,6 | 32,2 | 29 732,8 |
Mean | 56,8 | 34,9 | 18,7 | 14,8 | 15,3 | 14,2 | 12,8 | 13,7 | 16,9 | 22 595,1 |
Std | 4,4 | 5,9 | 1,3 | 7,2 | 4,6 | 5,2 | 5,8 | 5,5 | 9,1 | 3 969,2 |
groupe 4 | ||||||||||
Q5 | 44,2 | 14,7 | 12,4 | 11,3 | 12,5 | 11,1 | 4,0 | 6,7 | 9,3 | 18 686,5 |
Q10 | 45,7 | 15,6 | 12,9 | 16,4 | 14,1 | 13,0 | 4,5 | 7,2 | 11,2 | 20 220,0 |
Q25 | 49,3 | 18,7 | 13,7 | 24,7 | 16,5 | 15,8 | 5,8 | 7,9 | 14,2 | 22 622,0 |
Q50 | 52,3 | 22,0 | 14,6 | 30,4 | 18,6 | 19,0 | 7,5 | 9,7 | 17,9 | 24 843,0 |
Q75 | 55,1 | 26,3 | 15,9 | 37,9 | 20,8 | 22,3 | 10,1 | 12,1 | 22,8 | 28 590,0 |
Q90 | 58,9 | 30,9 | 17,1 | 42,3 | 22,4 | 27,2 | 14,4 | 15,3 | 30,1 | 31 294,0 |
Q95 | 60,8 | 34,5 | 18,3 | 45,7 | 23,2 | 28,9 | 19,3 | 18,6 | 33,2 | 34 574,0 |
Mean | 52,3 | 22,8 | 14,9 | 30,3 | 18,4 | 19,4 | 8,8 | 10,9 | 19,5 | 25 547,4 |
Std | 5,1 | 6,3 | 1,7 | 10,0 | 4,1 | 5,9 | 5,5 | 6,2 | 9,2 | 4 603,5 |
13.4.5.2 Mise en oeuvre du c-medoids dans R
La méthode du c-medoids dans R peut être mise en oeuvre avec la fonction FKM.med
du package fclust
. Le processus d’analyse et de validation est identique à celui présenté ci-dessus pour le c-means. Nous ne donnons donc pas un exemple complet de la méthode.
Stabilité des groupes obtenus par les méthodes de nuées dynamiques :
Puisque la méthode k-means et ses variantes reposent sur une initialisation aléatoire de leur algorithme, les résultats peuvent varier en fonction de cet état de départ. Dans certains contextes, il est possible que les résultats obtenus varient significativement, ce qui signifie que les groupes obtenus ne sont pas représentatifs de la population étudiée. Une solution pour vérifier si ce problème se pose est simplement de relancer l’algorithme un grand nombre de fois (généralement 1000) et de comparer les résultats obtenus au cours de ces réplications.
Cette méthode est rarement implémentée directement et requiert d’écrire sa propre fonction. geocmeans
dispose d’une fonction déjà existante, mais ne pouvant être appliquée qu’avec l’algorithme c-means. Nous proposons ici une implémentation pour la méthode k-means qui peut facilement être adaptée aux autres méthodes de classifications heuristiques.
La démarche à suivre est la suivante :
- Appliquer l’algorithme une première fois pour obtenir une classification de référence à laquelle nous comparerons toutes les réplications.
- Effectuer 1000 itérations au cours desquelles :
- Une nouvelle classification est calculée.
- Les groupes obtenus sont comparés à ceux de la classification de référence.
- L’indice de Jacard est calculé entre les groupes des deux classifications.
- Les valeurs de l’indice de Jacard sont enregistrées.
- Les centres des groupes sont enregistrés.
Ainsi, nous obtenons 1000 valeurs de l’indice de Jacard pour chaque groupe. Cet indice permet de mesurer le degré d’accord entre deux variables (ici les probabilités d’appartenance des observations au même groupe pour deux classifications différentes.) Une valeur moyenne en dessous de 0,5 indique qu’un groupe est très instable, car nous obtenons des résultats très différents lors des réplications. Une valeur entre 0,6 et 0,75 indique qu’un groupe semble bien exister dans les données, bien que marqué par une certaine incertitude. Une valeur au-dessus de 0,8 indique un groupe bien identifié et stable.
Nous obtenons également les centres des groupes des 1000 classifications. Il est ainsi possible de représenter leurs histogrammes et de déterminer si les centres des groupes sont stables (unimodalité et faible variance) ou incertains (plusieurs modes et/ou forte variance).
# X sera le jeu de données pour la classification
# clust_ref sera le vecteur indiquant le groupe de chaque observation obtenu par kmeans
# nsim sera le nombre de simulations à effectuer
<- function(X, clust_ref, nsim, verbose = TRUE){
kmeans_stability
# définition de la matrice d'appartenance originale
<- length(unique(clust_ref))
k <- dummy_cols(clust_ref, remove_selected_columns = TRUE)
ref_mat colnames(ref_mat) <- paste0("groupe_",1:k)
# lancement des itérations
<- lapply(1:nsim, function(i){
sim_resultats
# afficher la progression si requis
if(verbose){
print(paste0("iteration numeros : ",i,"/",nsim))
}# calculer le kmeans
<- kmeans(X,k)
km_res <- dummy_cols(km_res$cluster, remove_selected_columns = TRUE)
sim_mat
# ajustement de l'ordre des groupes avec geocmeans
<- groups_matching(as.matrix(ref_mat), as.matrix(sim_mat))
sim_mat
# calcul des indices de jacard
<- sapply(1:k, function(j){
jac_idx calc_jaccard_idx(sim_mat[,j], ref_mat[,j])
})
# recuperation des centres des groupes
<- as.integer(gsub(".data_","",colnames(sim_mat),fixed = TRUE))
idx <- data.frame(km_res$centers)
centers <- centers[idx,]
centers $groupe <- 1:k
centers
return(list(
"jac_idx" = jac_idx,
"centers" = centers
))
})
# les simulations sont finies, nous pouvons combiner les resultats
<- do.call(rbind, lapply(sim_resultats, function(x){x$jac_idx}))
all_jac_values <- do.call(rbind, lapply(sim_resultats, function(x){x$centers}))
all_centers return(list(
"jacard_values" = all_jac_values,
"centers" = all_centers
)) }
Il ne nous reste plus qu’à utiliser notre nouvelle fonction pour déterminer si les groupes obtenus avec notre k-means sont stables.
data(LyonIris)
set.seed(123)
# NB : LyonIris est un objet spatial, il faut donc extraire uniquement son DataFrame
<- st_drop_geometry(LyonIris[c("Lden","NO2","PM25","VegHautPrt","Pct0_14","Pct_65","Pct_Img",
X "TxChom1564","Pct_brevet","NivVieMed")])
# Centrage et réduction de chaque colonne du DataFrame
for (col in names(X)){
<- scale(X[[col]], center = TRUE, scale = TRUE)
X[[col]]
}
# calcul du kmeans de référence
<- kmeans(X, 4)
kmeans_ref
# application de la fonction de stabilité
<- kmeans_stability(X, kmeans_ref$cluster, nsim = 1000, verbose = FALSE) stab_results
Nous pouvons à présent vérifier la stabilité de nos quatre groupes.
<- data.frame(stab_results$jacard_values)
jacard_values names(jacard_values) <- paste("groupe", 1:4, sep = "_")
<- reshape2::melt(jacard_values)
df $groupes <- as.factor(df$variable)
df
ggplot(df) +
geom_histogram(aes(x = value), bins = 30) +
facet_wrap(vars(groupes), ncol=2) +
labs(x = "", y = "Indice de Jacard")
La figure 13.29 indique clairement que les groupes 1 et 2 sont très stables, car les valeurs de Jacard obtenues sont le plus souvent supérieures à 0,75. Le groupe 3 a le plus souvent des valeurs légèrement inférieures aux deux premiers groupes, mais tout de même bien supérieures à 0,5. En revanche, le groupe 4 a un grand nombre de valeurs inférieures à 0,5 indiquant une tendance du groupe à se dissoudre lors des réplications.
Considérant que le dernier groupe est le plus instable, nous décidons d’observer les valeurs des centres qu’il obtient pour les différentes réplications.
<- subset(stab_results$centers, stab_results$centers$groupe == 4)
centers_groupe4 $groupe <- NULL
centers_groupe4
<- reshape2::melt(centers_groupe4)
df $variable <- as.factor(df$variable)
df
ggplot(df) +
geom_histogram(aes(x = value), bins = 30) +
facet_wrap(vars(variable), ncol=3, scales="free") +
labs(x = "", y = "")
Les différents histogrammes de la figure 13.30 indiquent clairement que pour plusieurs variables (Lden
, NO2
, PM25
, Pct_Img
, et NivVieMed
) les caractéristiques du groupe 4 varient grandement sur l’ensemble des réplications. Nous pouvons comparer ce graphique à celui du groupe 2 qui est bien plus stable.
<- subset(stab_results$centers, stab_results$centers$groupe == 2)
centers_groupe2 $groupe <- NULL
centers_groupe2
<- reshape2::melt(centers_groupe2)
df $variable <- as.factor(df$variable)
df
ggplot(df) +
geom_histogram(aes(x = value), bins = 30) +
facet_wrap(vars(variable), ncol=3, scales="free") +
labs(x = "", y = "")
Nous pouvons constater une plus faible variance des résultats obtenus (en regardant notamment l’axe horizontal) pour les centres des groupes à la figure 13.31.