k-means

Importation des librairies et données

library(tidyverse)
library(factoextra)
library(igraph)
library(cccd)
library(mclust)
library(umap)
load(url("https://lbmc.gitbiopages.ens-lyon.fr/hub/formations/ens_m1_ml/practical_b.Rdata"), verbose = T)

Commentaires

• Pour l’initialisation des centroïdes, nous avons choisi de tirer au hasard des individus pour définir nos centroïdes. Une autre solution aurait été de prendre aléatoirement des points dans l’espace de dimension n (nombre de composantes principales). Nous avons choisi la première pour deux raisons:
  • aucun cluster ne peut être vide.
  • dans un espace de dimension n > 3, les résultats étaient meilleurs.
• Pour la fonction metric_example, nous avons choisi d’autres arguments que ceux demandés car cela nous semblait plus direct pour calculer I_w.
• Pour la fonction kmeans_example, nous avons choisi d’arrêter la fonction lorsque le rapport de I_w à l’étape \([h - 1]\) avec I_w à l’étape \([h]\) est inférieur à une valeur empiriquement choisie (i.e. 1.0001). Nous avons rajouté une sécurité qui arrête la fonction après 100 itérations dans le cas où I_w ne se stabiliserait pas.

Initialisation de la position des centroïdes

kmeans_initiation <- function(x, k) {
  centroids <- list()
  points <- sample(seq_len(nrow(x)), k)
  for (i in 1:k) {
    centroids[[i]] <- x[points[i], ]
  }
  return(centroids)
}

Calcul de la distance euclidienne aux centroïdes

compute_distance <- function(x, centroid) {
  size <- length(centroid)
  distance <- rep(NA, nrow(x))
  for (i in 1:size) {
    distance <- cbind(distance, apply(x,
                                      1,
                                      function(x1, x2) {
                                        sum((x1 - x2)^2)
                                      },
                                      x2 = centroid[[i]]))
  }
  return(distance[, 2:(size + 1)])
}

Assigner les individus aux différents centroïdes

cluster_assignment <- function(distance) {
  return(apply(distance, 1, which.min))
}

Mise à jour de la position des centroïdes

centroid_update <- function(x, cluster, k) {
  updated_centroids <- list()
  for (i in 1:k) {
    if (sum(cluster == i) > 1) {
      updated_centroids[[i]] <- colMeans(x[cluster == i, ])
    } else {
      updated_centroids[[i]] <- x[cluster == i, ]
    }
  }
  return(updated_centroids)
}

Calcul de I_w

metric_example <- function(x, cluster, centroid) {
  iw <- 0
  for (i in seq_along(centroid)) {
    iw <- iw + sum(apply(x[cluster == i, ],
                         1,
                         function(x1, x2) {
                           sum((x1 - x2)^2)
                         },
                         x2 = centroid[[i]]))
  }
  return(iw)
}

Association des différentes fonctions pour faire le clustering

kmeans_example <- function(x, k) {
  centroid <- kmeans_initiation(x, k)
  iw1 <- 0
  iw2 <- 0
  i <- 0
  while (i < 3 || iw1 / iw2 > 1.0001) {
    iw1 <- iw2
    distance <- compute_distance(x, centroid)
    cluster <- cluster_assignment(distance)
    centroid <- centroid_update(x, cluster, k)
    iw2 <- metric_example(x, cluster, centroid)
    i <- i + 1
    if (i > 100) {
      break
    }
  }
  return(cluster)
}

Test de notre fonction

data_pca <- data[var_gene_2000[1:600], ] %>%
  t() %>%
  prcomp()

data_pca %>%
  fviz_pca_ind(
    geom = "point",
    col.ind = as.factor(kmeans_example(data_pca$x[, 1:2], k = 9))
  )

Comparaison avec la fonction originelle

data_pca <- data[var_gene_2000[1:600], ] %>%
  t() %>%
  prcomp()

data_dist <- dist(data_pca$x[, 1:2])

data_kmeans <- kmeans(data_dist, centers = 9)

data_pca %>%
  fviz_pca_ind(
    geom = "point",
    col.ind = as.factor(data_kmeans$cluster)
  )