6.1 K-moyennes

Les k-moyennes (ou “k-means” en anglais) représentent une autre façon de regrouper les individus d’un tableau multivarié. Par rapport à la CAH, cette technique est généralement moins efficace, mais elle a l’avantage de permettre le regroupement d’un très grand nombre d’individus (gros jeu de données), là où la CAH nécessiterait trop de temps de calcul et de mémoire vive. Il est donc utile de connaître cette seconde technique à utiliser comme solution de secours lorsque le dendrogramme de la CAH devient illisible sur de très gros jeux de données.

À vous de jouer !

Le principe des k-moyennes est très simple¹⁴ :

L’utilisateur choisi le nombre de groupes k qu’il veut obtenir à l’avance.
La position des k centres est choisie au hasard au début.
Les individus sont attribués aux k groupes en fonction de leurs distances aux centres (attribution au groupe de centre le plus proche).
Les k centres sont replacés au centre de gravité des groupes ainsi obtenus.
Les individus sont réaffectés en fonction de leurs distances à ces nouveaux centres.
Si au moins un individu a changé de groupe, le calcul est réitéré. Sinon, nous considérons avoir atteint la configuration finale.

La technique est illustrée dans la vidéo suivante :

Essayez par vous même via l’application ci-dessous qui utilise le célèbre jeu de données iris. Notez que vous devez utiliser des variables numériques. Par exemple, Species étant une variable qualitative, vous verrez que cela ne fonctionne pas dans ce cas.

6.1.1 Exemple simple

R propose la fonction kmeans() et différents packages la supplémente pour enrichir votre boite à outils. Nous rassemblons tout cela sous une interface cohérente dans SciViews autour d’une fonction légèrement différente : k_means(). Cependant, ce code n’est pas encore inclus dans un package. Vous aurez donc à intégrer le contenu du chunk suivant dans vos scripts et documents R Markdown avant de pouvoir utiliser nos fonctions (cliquez sur “voir le code” pour le dérouler).

# k_means for SciViews, version 1.0.0
# Copyright (c) 2021, Philippe Grosjean (phgrosjean@sciviews.org)

SciViews::R
library(broom)

# scale() is a generic function, but it does not provide a method for data
# frames. As such, data frames and tibbles are converted into matrices by the
# default method, which is not what we want
scale.data.frame <- function(x, center = TRUE, scale = TRUE)
  as.data.frame(scale(as.matrix(x)))
# This is for tibbles
scale.tbl_df <- function(x, center = TRUE, scale = TRUE)
  tibble::as_tibble(scale(as.matrix(x)))
# This is for data.tables
scale.data.table <- function(x, center = TRUE, scale = TRUE)
  data.table::as.data.table(scale(as.matrix(x)))

# This is a reworked version of factoextra::fviz_nbclust() to help chosing the
# number of clusters for kmeans()
profile_k <- function(x, FUNcluster = kmeans, method = "wss", k.max = NULL, ...) {
  if (NROW(x) < 2)
    stop("You must provide an data frame or matrix with at least two rows")
  if (is.null(k.max))
    k.max <- min(nrow(x) - 1, 10) # Avoid error with very small datasets in fviz_nbclust()
  factoextra::fviz_nbclust(x, FUNcluster = FUNcluster, method = method,
    k.max = k.max, ...)
}

# Traditional kmeans does not store the data... and this is a problem for plot
# later on. So, we define k_means() which store the original data by default
k_means <- function(x, k, centers = k, iter.max = 10L, nstart = 1L,
algorithm = c("Hartigan-Wong", "Lloyd", "Forgy", "MacQueen"), trace = FALSE,
keep.data = TRUE) {
  # k and centers are synonyms
  res <- kmeans(x, centers = centers, iter.max = iter.max, nstart = nstart,
    algorithm = algorithm, trace = trace)
  if (isTRUE(keep.data))
    res$data <- as.data.frame(x)
  class(res) <- unique(c("k_means", class(res)))
  res
}

# broom::augment.kmeans() seems buggy when data is called 'x'
augment.kmeans <- function(object, data) {
  res <- broom::fix_data_frame(data, newcol = ".rownames")
  res$.cluster <- factor(object$cluster)
  res
}

# No predict() method for kmeans, but we add one for k_means
predict.k_means <- function(object, ...)
  factor(object$cluster)

# There is no plot, autoplot and chart methods for kmeans objects => make them
# for k_means objects, because we have both the k-means results and the data
# Since data is no contained in the kmeans object, one has to provide it also
plot.k_means <- function(x, y, data = x$data, choices = 1L:2L,
  col = NULL, c.shape = 8, c.size = 3, ...) {
  nclust <- nrow(x$centers)
  if (is.null(col))
    col <- 1:nclust
  plot(as.data.frame(data)[, choices], col = col[x$cluster], ...)
  points(as.data.frame(x$centers)[, choices], col = col[1:nclust],
    pch = c.shape, cex = c.size)
}

autoplot.k_means <- function(object, data = object$data, choices = 1L:2L,
alpha = 1, c.shape = 8, c.size = 3, theme = NULL, use.chart = FALSE, ...) {
  data <- as.data.frame(data)
  vars <- choices
  if (is.numeric(choices))
    vars <- colnames(data)[choices]
  var_x <- as.name(vars[1])
  var_y <- as.name(vars[2])
  centers <- broom::tidy(object, col.names = colnames(data))
  cluster <- factor(object$cluster)
  if (isTRUE(use.chart)) {
    fun <- chart::chart
  } else {
    fun <- ggplot2::ggplot
  }
  res <- fun(data = data, mapping = aes(x = {{var_x}}, y = {{var_y}},
    col  = cluster)) +
    geom_point(alpha = alpha) +
    geom_point(data = centers, size = c.size, shape = c.shape)
  if (!is.null(theme))
    res <- res + theme
  res
}

chart.k_means <- function(data, ..., type = NULL, env = parent.frame())
  autoplot(data, type = type, theme = theme_sciviews(), use.chart = TRUE, ...)

Vous avez les méthodes print(), plot(), autoplot(), chart(), augment(), tidy(), glance(), predict() et fitted(). Nous verrons leur utilisation au fur et à mesure des explications dans cette partie. Enfin, la fonction profile_k() permet de recherche k, le nombre optimal de clusters à réaliser.

Afin de comparer la classification par k-moyennes à celle par CAH, nous reprendrons ici le même jeu de données zooplankton.

zoo <- read("zooplankton", package = "data.io")
zoo

# # A tibble: 1,262 x 20
#      ecd  area perimeter feret major minor  mean  mode   min   max std_dev range
#    <dbl> <dbl>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl>
#  1 0.770 0.465      4.45 1.32  1.16  0.509 0.363 0.036 0.004 0.908   0.231 0.904
#  2 0.700 0.385      2.32 0.728 0.713 0.688 0.361 0.492 0.024 0.676   0.183 0.652
#  3 0.815 0.521      4.15 1.33  1.11  0.598 0.308 0.032 0.008 0.696   0.204 0.688
#  4 0.785 0.484      4.44 1.78  1.56  0.394 0.332 0.036 0.004 0.728   0.218 0.724
#  5 0.361 0.103      1.71 0.739 0.694 0.188 0.153 0.016 0.008 0.452   0.110 0.444
#  6 0.832 0.544      5.27 1.66  1.36  0.511 0.371 0.02  0.004 0.844   0.268 0.84 
#  7 1.23  1.20      15.7  3.92  1.37  1.11  0.217 0.012 0.004 0.784   0.214 0.78 
#  8 0.620 0.302      3.98 1.19  1.04  0.370 0.316 0.012 0.004 0.756   0.246 0.752
#  9 1.19  1.12      15.3  3.85  1.34  1.06  0.176 0.012 0.004 0.728   0.172 0.724
# 10 1.04  0.856      7.60 1.89  1.66  0.656 0.404 0.044 0.004 0.88    0.264 0.876
# # … with 1,252 more rows, and 8 more variables: size <dbl>, aspect <dbl>,
# #   elongation <dbl>, compactness <dbl>, transparency <dbl>, circularity <dbl>,
# #   density <dbl>, class <fct>

Commençons par l’exemple simplissime de la réalisation de deux groupes à partir de six individus issus de ce jeu de données, comme nous l’avons fait avec la CAH :

zoo %>.%
  select(., -class) %>.% # Élimination de la colonne class
  slice(., 13:18) -> zoo6      # Récupération des lignes 13 à 18

zoo6_kmn <- k_means(zoo6, k = 2)
zoo6_kmn

# K-means clustering with 2 clusters of sizes 3, 3
# 
# Cluster means:
#         ecd      area perimeter    feret    major     minor      mean      mode
# 1 0.6292647 0.3188667  3.224133 1.159200 1.096433 0.4023333 0.1871667 0.1026667
# 2 1.1926500 1.1279667 10.346667 2.201133 1.677067 0.8596333 0.3217333 0.3533333
#          min       max   std_dev     range      size    aspect elongation
# 1 0.01066667 0.5400000 0.1166667 0.5293333 0.7493833 0.4753843   6.333315
# 2 0.00400000 0.8986667 0.2620000 0.8946667 1.2683500 0.5149422  23.046713
#   compactness transparency circularity    density
# 1    2.727708   0.14732060   0.4900333 0.06943333
# 2    7.987806   0.06173831   0.1357000 0.37630000
# 
# Clustering vector:
# [1] 1 2 1 2 2 1
# 
# Within cluster sum of squares by cluster:
# [1]  54.18647 200.03837
#  (between_SS / total_SS =  68.1 %)
# 
# Available components:
# 
#  [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
#  [6] "betweenss"    "size"         "iter"         "ifault"       "data"

Nous voyons que la fonction k_means() effectue notre classification¹⁵. Nous lui fournissons le tableau de départ et spécifions le nombre k de groupes souhaités via l’argument k =. Ne pas oublier d’assigner le résultat du calcul à une nouvelle variable, ici zoo6_kmn, pour pouvoir l’inspecter et l’utiliser par la suite. L’impression du contenu de l’objet nous donne plein d’information dont :

le nombre d’individus dans chaque groupe (ici 3 et 3),
la position des centres pour les k groupes dans Cluster means,
l’appartenance aux groupes dans Clustering vector (dans le même ordre que les lignes du tableau de départ),
la sommes des carrés des distances entre les individus et la moyenne au sein de chaque groupe dans Within cluster sum of squares ; le calcul between_SS / total_SS est à mettre en parallèle avec le $R^2$ de la régression linéaire : c’est une mesure de la qualité de regroupement des données (plus la valeur est proche de 100% mieux c’est, mais attention que cette valeur augmente d’office en même temps que k),
et enfin, la liste des composants accessibles via l’opérateur $ ; par exemple, pour obtenir la taille de chaque groupe (en nombre d’individus), nous ferons :

zoo6_kmn$size

# [1] 3 3

Le package {broom} contient trois fonctions complémentaires qui nous seront utiles : tidy(), augment() et glance(). broom::glance() retourne un data.frame avec les statistiques permettant d’évaluer la qualité de la classification obtenue :

glance(zoo6_kmn)

# # A tibble: 1 x 4
#   totss tot.withinss betweenss  iter
#   <dbl>        <dbl>     <dbl> <int>
# 1  796.         254.      542.     1

Si nous voulons déterminer la valeur optimale de k, pous pouvons utiliser profile_k() appliqué au jeu de données initial, ou à la composante data de notre objet k_means (spécifier éventuellemznt une valeur différente de celle par défaut pour l’argument k.max = de la fonction, voir l’aide en ligne de ?factoextra::fviz_nbclust) :

profile_k(zoo6) # ou zoo6_kmn$data

Le graphique obtenu montre la décroissance de la somme des carrés des distances intra-groupes en fonction de k. Avec k = 1, nous considérons toutes les données dans leur ensemble et nous avons simplement la somme des carrés des distances euclidiennes entre tous les individus et le centre de gravité du nuage de points dont les coordonnées sont les moyennes de chaque variable. C’est le point de départ qui nous indique de combien les données sont dispersées (la valeur absolue de ce nombre n’est pas importante).

Ensuite, avec k croissant, notre objectif est de faire des regroupement qui diminuent la variance intra-groupe autant que possible, ce que nous notons par la diminution de la somme des carrés intra-groupes (la variance du groupe est, en effet, la somme des carrés des distances euclidiennes entre les points et le centre du groupe, divisée par les degrés de liberté).

Nous recherchons ici des sauts importants dans la décroissance de la somme des carrés, tout comme dans le dendrogramme obtenu par la CAH nous recherchions des sauts importants dans les regroupements (hauteur des barres verticales du dendrogramme). Nous observons ici un saut important pour k = 2, puis une diminution un peu moins forte pour k = 3 et ensuite une stagnation. Ceci suggère que nous pourrions considérer deux voire trois groupes.

Le nombre de groupes proposé par profile_k() n’est qu’indicatif ! Si vous avez par ailleurs d’autres informations qui vous suggèrent un regroupement différent, ou si vous voulez essayer un regroupement plus ou moins détaillé par rapport à ce qui est proposé, c’est tout aussi correct.

La fonction profile_k() propose d’ailleurs deux autres méthodes pour déterminer le nombre optimal de groupes k, avec method = "silhouette" ou method = "gap_stat". Voyez l’aide en ligne de cette fonction ?factoextra::fviz_nbclust. Ces différentes méthodes peuvent d’ailleurs suggérer des regroupements différents pour les mêmes données… preuve qu’il n’y a pas une et une seule solution optimale !

A ce stade, nous pouvons collecter les groupes et les ajouter à notre tableau de données avec augment(). Notez que si vous voulez juste récupérer les groupes, vous pouvez utiliser alors predict(zoo6_kmn). augment() crée une nouvelle colonne nommée .cluster que nous renommons ici immédiatement en cluster, et ensuite, nous enregistrons le tout dans un nouveau jeu de données nommé zoo6b¹⁶.

augment(zoo6_kmn, zoo6) %>.%
  rename(., cluster = .cluster) -> zoo6b
names(zoo6b)

#  [1] "ecd"          "area"         "perimeter"    "feret"        "major"       
#  [6] "minor"        "mean"         "mode"         "min"          "max"         
# [11] "std_dev"      "range"        "size"         "aspect"       "elongation"  
# [16] "compactness"  "transparency" "circularity"  "density"      "cluster"

La nouvelle variable cluster contient ceci :

zoo6b$cluster

# [1] 1 2 1 2 2 1
# Levels: 1 2

C’est le contenu de zoo6_kmn$cluster, mais transformé en variable factor.

class(zoo6b$cluster)

# [1] "factor"

Nous pouvons enfin utiliser tidy() pour obtenir un tableau avec les coordonnées des k centres. Nous l’enregistrerons dans la variable zoo6_centers, en ayant bien pris soin de nommer les variables du même nom que dans le tableau original zoo6 (argument col.names = names(zoo6), afin de conserver le noms de nos variables initiales dans ce nouveau tableau :

zoo6_centers <- tidy(zoo6_kmn, col.names = names(zoo6))
zoo6_centers

# # A tibble: 2 x 21
#     ecd  area perimeter feret major minor  mean  mode    min   max std_dev range
#   <dbl> <dbl>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>   <dbl> <dbl>
# 1 0.629 0.319      3.22  1.16  1.10 0.402 0.187 0.103 0.0107 0.540   0.117 0.529
# 2 1.19  1.13      10.3   2.20  1.68 0.860 0.322 0.353 0.004  0.899   0.262 0.895
# # … with 9 more variables: size <int>, aspect <dbl>, elongation <dbl>,
# #   compactness <dbl>, transparency <dbl>, circularity <dbl>, density <dbl>,
# #   withinss <dbl>, cluster <fct>

La dernière colonne de ce tableau est également nommée cluster. C’est le lien entre le tableau zoo6b augmenté et zoo6_centers. Nous avons maintenant tout ce qu’il faut pour représenter graphiquement les regroupements effectués par les k-moyennes en colorant les points en fonction de la nouvelle variable cluster.

chart(data = zoo6b, area ~ circularity %col=% cluster) +
  geom_point() + # Affiche les points représentant les individus
  geom_point(data = zoo6_centers, size = 5, shape = 17) # Ajoute les centres

Comparez avec le graphique équivalent au module précédent consacré à la CAH. Outre que l’ordre des groupes est inversé et que les données n’ont pas été standardisées ici, un point est classé dans un groupe différent par les deux méthodes. Il s’agit du point ayant environ 0.25 de circularité et 0.5 de surface. Comme nous connaissons par ailleurs la classe à laquelle appartient chaque individu, nous pouvons la récupérer comme colonne supplémentaire du tableau zoo6b et ajouter cette information sur notre graphique.

zoo6b$class <- zoo$class[13:18]
zoo6_centers$class <- "" # Ceci est nécessaire pour éviter le label des centres
chart(data = zoo6b, area ~ circularity %col=% cluster %label=% class) +
  geom_point() +
  ggrepel::geom_text_repel() + # Ajoute les labels intelligemment
  geom_point(data = zoo6_centers, size = 5, shape = 17)

Nous constatons que le point classé différemment est un “Poecilostomatoïd”. Or, l’autre groupe des k-moyennes contient aussi un individu de la même classe. Donc, CAH a mieux classé notre plancton que les k-moyennes dans le cas présent. Ce n’est pas forcément toujours le cas, mais souvent.

Maintenant que nous comprenons bien la logique de création de ce graphique en l’ayant réalisé par l’extraction des groupes avec augment() et des centres de ces groupes à l’aide de tidy(), nous pouvons nous simplifier la vie en utilisant chart() simplement sur notre objet k_means pour obtenir directement ce graphique. Dans ce cas, nous n’utilisons pas une formule pour indiquer les variables à utiliser, mais un vecteur de deux nombre qui indique l’index de ces variables, ou leur noms directement dans choices = :

chart(zoo6_kmn, choices = c("circularity", "area"),
  alpha = 0.6, c.size = 5, c.shape = 17)

La taille et la forme des centres sont indiqués respectivement par c.size = et c.shape = , tandis que alpha = modifie la transparence des points. Ce dernier argument est surtout utile pour des gros jeux de données avec beaucoup de points à afficher sur le graphique.

Comme les k-moyennes partent d’une position aléatoire des k centres, le résultat final peut varier et n’est pas forcément optimal. Pour éviter cela, nous pouvons indiquer à k_means() d’essayer différentes situations de départ via l’argument nstart =. Par défaut, nous prenons une seule situation aléatoire de départ nstart = 1, mais en indiquant une valeur plus élevée pour cet argument, il est possible d’essayer plusieurs situations de départ et ne garder que le meilleur résultat final. Cela donne une analyse plus robuste et plus reproductible… mais le calcul est naturellement plus long.

k_means(zoo6, k = 2, nstart = 50) # 50 positions de départ différentes

# K-means clustering with 2 clusters of sizes 3, 3
# 
# Cluster means:
#         ecd      area perimeter    feret    major     minor      mean      mode
# 1 0.6292647 0.3188667  3.224133 1.159200 1.096433 0.4023333 0.1871667 0.1026667
# 2 1.1926500 1.1279667 10.346667 2.201133 1.677067 0.8596333 0.3217333 0.3533333
#          min       max   std_dev     range      size    aspect elongation
# 1 0.01066667 0.5400000 0.1166667 0.5293333 0.7493833 0.4753843   6.333315
# 2 0.00400000 0.8986667 0.2620000 0.8946667 1.2683500 0.5149422  23.046713
#   compactness transparency circularity    density
# 1    2.727708   0.14732060   0.4900333 0.06943333
# 2    7.987806   0.06173831   0.1357000 0.37630000
# 
# Clustering vector:
# [1] 1 2 1 2 2 1
# 
# Within cluster sum of squares by cluster:
# [1]  54.18647 200.03837
#  (between_SS / total_SS =  68.1 %)
# 
# Available components:
# 
#  [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
#  [6] "betweenss"    "size"         "iter"         "ifault"       "data"

Dans ce cas simple, cela ne change pas grand chose. Mais avec un plus gros jeu de données plus complexe, cela peut être important.

À vous de jouer !

6.1.2 Classification du zooplancton

Maintenant que nous savons utiliser k_means() et les fonctions annexes, nous pouvons classer le jeu de données zoo tout entier.

zoo %>.%
  select(., -class) %>.%
  profile_k(., k.max = 15)

Nous observons un saut maximal pour k = 2, mais le saut pour k = 3 est encore conséquent. Afin de comparer avec ce que nous avons fait par CAH, nous utiliserons donc k = 3. Enfin, comme un facteur aléatoire intervient, qui définira au final le numéro des groupes, nous utilisons set.seed() pour rendre l’analyse reproductible. Pensez à donner une valeur différente à cette fonction pour chaque utilisation ! Et pensez aussi à éliminer les colonnes non numériques à l’aide de select().

set.seed(562)
zoo_kmn <- k_means(select(zoo, -class), k = 3, nstart = 50)
zoo_kmn

# K-means clustering with 3 clusters of sizes 786, 385, 91
# 
# Cluster means:
#         ecd     area perimeter    feret     major     minor      mean
# 1 0.6664955 0.431915  3.575374 1.134705 0.9744768 0.4780780 0.2388065
# 2 0.9715857 1.009902  9.197299 2.668022 1.8468984 0.6194652 0.1723774
# 3 1.3774670 1.998097 19.653860 4.063837 2.1465758 0.9602846 0.1488495
#         mode         min       max   std_dev     range      size    aspect
# 1 0.09256997 0.007094148 0.7269109 0.1842660 0.7198168 0.7262774 0.5372808
# 2 0.04455065 0.004207792 0.6315844 0.1512922 0.6273766 1.2331818 0.5349924
# 3 0.02470330 0.004000000 0.7013187 0.1472286 0.6973187 1.5534302 0.5362249
#   elongation compactness transparency circularity    density
# 1   7.184451    3.002093   0.07385014  0.42917214 0.09349338
# 2  27.898079    9.529156   0.11719954  0.11197351 0.16938468
# 3  61.837019   20.325398   0.09737903  0.05186813 0.31140879
# 
# Clustering vector:
#    [1] 1 1 1 1 1 1 3 1 3 1 1 2 1 2 1 2 2 1 3 1 1 1 1 2 2 1 1 2 1 1 1 1 1 1 2 1 2
#   [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 2 1 2 2 1 2 1 1 1 1 3
#   [75] 1 1 1 1 1 2 1 1 1 1 1 3 1 1 2 1 1 2 1 3 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1
#  [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 3 2 3 2
#  [149] 1 1 1 1 1 2 2 2 2 2 1 3 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
#  [186] 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1
#  [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1
#  [260] 1 1 1 1 1 1 1 1 1 3 3 3 2 2 3 1 1 2 1 2 2 1 2 1 1 3 1 1 3 2 3 2 2 1 1 1 1
#  [297] 2 2 3 1 2 1 2 2 1 2 2 2 3 1 2 2 3 1 2 3 3 1 3 2 2 2 1 1 2 1 1 2 1 1 1 1 1
#  [334] 2 1 3 1 1 1 2 1 1 1 1 1 1 2 2 2 1 2 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 2 1 2
#  [371] 1 1 1 1 2 1 1 1 1 1 1 2 2 1 1 1 1 1 2 2 1 1 2 1 2 2 1 2 2 1 1 1 1 1 1 1 1
#  [408] 2 1 1 1 2 2 1 1 1 1 1 1 1 2 2 2 3 1 2 2 1 1 1 2 2 2 1 1 1 1 3 1 1 2 2 3 1
#  [445] 2 1 2 1 2 1 2 1 1 2 1 2 1 2 2 1 1 1 1 1 2 1 1 2 1 1 1 2 1 1 2 1 1 1 2 1 1
#  [482] 1 1 1 1 1 1 2 1 2 2 2 1 2 1 2 1 2 2 2 1 1 1 1 1 1 1 1 2 2 1 2 3 2 1 3 1 1
#  [519] 2 2 2 1 2 1 1 1 3 1 3 3 2 1 2 1 1 1 2 1 1 1 1 2 2 1 1 1 2 2 1 1 2 1 2 1 1
#  [556] 2 3 2 1 1 1 1 3 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 2 2 2 3 1 3 2 2 3 1 2 2
#  [593] 2 1 1 2 2 3 1 1 1 1 1 1 1 1 1 2 1 2 2 1 1 1 1 1 1 1 3 1 1 1 1 1 2 2 3 2 1
#  [630] 2 3 2 1 2 1 2 2 1 1 1 1 2 2 2 2 2 2 2 1 1 2 1 1 1 1 1 1 1 2 2 2 2 2 2 1 3
#  [667] 2 2 2 2 1 1 2 2 1 1 1 1 2 2 2 1 1 1 1 1 1 2 1 1 1 2 1 1 2 2 2 2 2 1 1 3 1
#  [704] 1 3 2 1 2 1 1 2 3 3 2 1 1 3 2 1 3 2 2 2 3 2 1 2 3 1 2 1 1 3 2 1 3 1 3 2 2
#  [741] 2 1 3 1 2 2 3 2 2 1 1 2 1 2 2 1 1 1 1 1 2 2 2 1 1 1 1 2 2 2 1 1 1 1 2 1 1
#  [778] 1 2 1 1 2 1 2 1 1 1 2 1 1 2 1 1 2 1 1 2 1 1 1 2 1 2 3 2 1 1 2 2 3 1 2 1 1
#  [815] 2 2 2 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 2 2 1 1 1 2 1 2 2 1 2 1 1 1 1 1 1 1 2
#  [852] 2 2 1 1 2 2 2 1 1 1 1 1 1 2 1 2 3 2 1 2 1 2 1 1 3 1 3 1 2 2 2 2 2 2 2 1 1
#  [889] 2 1 1 1 2 3 1 2 1 1 1 1 2 1 2 1 1 1 2 2 1 2 1 1 2 1 1 2 1 3 1 2 2 1 1 2 2
#  [926] 1 1 1 3 2 3 2 3 2 2 2 3 1 2 3 2 2 2 1 2 2 2 2 2 2 3 2 2 1 2 1 1 1 1 1 1 2
#  [963] 1 1 1 1 1 1 2 1 1 3 1 3 1 1 1 1 2 2 2 2 1 1 1 3 3 1 2 2 1 3 1 2 2 3 2 2 2
# [1000] 2 1 2 3 2 3 2 3 1 2 1 1 1 1 2 2 1 2 1 3 2 1 1 2 2 2 2 2 1 2 2 2 1 1 2 2 1
# [1037] 3 2 2 2 3 1 2 3 1 2 2 3 1 1 2 1 1 1 2 1 1 1 1 1 1 1 1 2 2 2 2 2 3 2 2 2 2
# [1074] 1 2 2 2 1 2 3 1 2 2 2 3 2 1 1 1 1 1 3 1 2 2 1 3 3 2 1 1 2 2 1 1 3 2 2 2 2
# [1111] 2 1 1 1 3 1 1 1 1 1 1 1 2 1 1 1 2 1 2 2 1 1 2 1 1 1 1 1 1 2 2 1 2 2 2 2 3
# [1148] 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 1 1 2 2 1 1 1 2 1 1 1 1
# [1185] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1
# [1222] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
# [1259] 2 2 2 2
# 
# Within cluster sum of squares by cluster:
# [1] 24356.96 39813.31 43792.25
#  (between_SS / total_SS =  77.0 %)
# 
# Available components:
# 
#  [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
#  [6] "betweenss"    "size"         "iter"         "ifault"       "data"

Récupérons les groupes dans zoob

augment(zoo_kmn, zoo) %>.%
  rename(., cluster = .cluster) -> zoob

Et enfin, effectuons un graphique similaire à celui réalisé pour la CAH au module précédent. À noter que nous pouvons ici choisir n’importe quelle paire de variables quantitatives pour représenter le nuage de points. Nous ajoutons des ellipses pour matérialiser les groupes à l’aide de stat_ellipse(). Elles contiennent 95% des points du groupe à l’exclusion des extrêmes. Enfin, comme il y a beaucoup de points, nous choisissons de les rendre semi-transparents avec l’argument alpha = 0.2 pour plus de lisibilité du graphique.

chart(zoo_kmn, choices = c("ecd", "compactness"), alpha = 0.2) +
  stat_ellipse()

Nous observons ici un regroupement beaucoup plus simple qu’avec la CAH, essentiellement stratifié de bas en haut en fonction de la compacité des points (Compactness). La tabulation des groupes en fonction des classes connues par ailleurs montre aussi que les k-moyennes les séparent moins bien que ce qu’à pu faire la CAH :

table(zoob$class, zoob$cluster)

#                   
#                      1   2   3
#   Annelid           38   6   6
#   Appendicularian   21  15   0
#   Calanoid          82 165  41
#   Chaetognath        6  45   0
#   Cirriped          14   8   0
#   Cladoceran        50   0   0
#   Cnidarian         13   6   3
#   Cyclopoid          5  40   5
#   Decapod          117   9   0
#   Egg_elongated     50   0   0
#   Egg_round         49   0   0
#   Fish              50   0   0
#   Gastropod         50   0   0
#   Harpacticoid       1  29   9
#   Malacostracan     54  41  26
#   Poecilostomatoid 143  14   1
#   Protist           43   7   0

Le groupe numéro 3 n’est pas vraiment défini en terme des classes de plancton car aucune classe ne s’y trouve de manière majoritaire. Le groupe numéro 1 contient la majorité des items de diverses classes, alors que le groupe 2 a une majorité de calanoïdes, de cyclopoïdes et d’harpacticoïdes (différents copépodes). Globalement, le classement a un sens, mais est moins bien corrélé avec les classes de plancton que ce que la CAH nous a fourni. Notez que, si nous avions standardisé les données avant d’effectuer les k-moyennes comme nous l’avons fait pour la CAH, nous aurions obtenu d’autres résultats. La transformation des variables préalablement à l’analyse reste une approche intéressante pour moduler l’importance des différentes variables entre elles dans leur impact sur le calcul des distances, et donc, des regroupements réalisés. Nous vous laissons réaliser les k-moyennes sur les données zoo standardisées à l’aide de la fonction scale() comme pour la CAH comme exercice.

À vous de jouer !

Pour en savoir plus

Il existe une approche mixte qui combine la CAH et les k-moyennes. Cette approche est intéressante pour les gros jeux de données. Elle est expliquée ici, et l’implémentation dans la fonction factoextra::hkmeans() est détaillée ici (en anglais).

Cet article explique dans le détail kmeans() et hclust() dans R, et montre aussi comment on peut calculer les k-moyennes à la main pour bien en comprendre la logique (en anglais).

À vous de jouer !

Effectuez maintenant les exercices du tutoriel B06La_kmeans (Partitionnement par k-moyennes).

BioDataScience2::run("B06La_kmeans")

Réalisez l’assignation B06Ia_fish_market, partie I.

Assignation GitHub Classroom pour les étudiants inscrits au cours de Science des Données Biologiques II à l’UMONS

Assignation GitHub Classroom pour les étudiants inscrits au cours de Bioinformatique et Science des Données II à Charleroi

Si vous êtes un utilisateur non enregistré ou que vous travaillez en dehors d’un cours, faites un “fork” de ce dépôt.

Voyez les explications dans le fichier README.md, partie I.

En pratique, différents algorithmes avec diverses optimisations existent. Le plus récent et le plus sophistiqué est celui de Hartigan-Wong. Il est utilisé par défaut par la fonction kmeans(). En pratique, il y a peu de raison d’en changer.↩︎
Utilisez l’aide en ligne de ?kmeans pour connaître les arguments. Seul centers = est changé en k = dans k_means(), mais avec centers = aussi accepté.↩︎
De manière générale, éviter de rajouter des données calculées dans le jeu de données initial. Cela peut amener à des erreurs particulièrement délicates si vous relancer ensuite l’analyse sur ce tableau.↩︎