Introduction au Cluster

1 Le Clustering

1.1 Chargement des packages

library(tidyverse)#manipulation de données + ggplot
library(MASS)#fonction statistiques (ex:kde2d)
library(mvtnorm)#simulation de lois normales multivariées

Introduction:

Nous allons initier à ces méthodes de clustering grace à l’article…

. Suivre scrupuleusement l’article

. Eventuellement appliquer à un jeu de données Non-supervisé…

Les Datas . Les données sont disponible dans ce lien : Télécharger le code R

.Ce sont des données simulés… Les autrices générent des données numériques (2 variables) et (4 variables) qualitatives.

Remarque: Le k-means se fait sur des données numériques … C’est la méthode que l’on verra en profondeur. Elles générent trois groupes d’individus, on les nomme black(177),red(237),green(186),soit en tout 600.

1:Proportions des groupes

set.seed(288) # La graine pour reproductibilité
pi = c(0.3, 0.4, 0.3) # Les proportions attendues de chaque

mu = matrix(c(0, 0, 3, 5, 0, 6), 2, 3, byrow=FALSE) #les moyennes
mu

     [,1] [,2] [,3]
[1,]    0    3    0
[2,]    0    5    6

Sigma.sph = diag(c(1, 1))#cercle(variance identique)
Sigma.sph

     [,1] [,2]
[1,]    1    0
[2,]    0    1

Sigma.ell = matrix(c(2, 1.3, 1.3, 1),2, 2, byrow=TRUE)#ellipse(corrélation entre variables)
Sigma.ell

     [,1] [,2]
[1,]  2.0  1.3
[2,]  1.3  1.0

cl = sample(c(1, 2, 3), 600, T, pi)#chaque individu reçoit un groupe selon les proportions
cl

  [1] 2 1 2 1 2 2 3 2 2 3 1 3 1 2 3 3 2 2 2 2 3 3 1 2 1 2 2 2 2 1 1 2 1 1 3 2 1
 [38] 1 2 1 2 3 2 1 2 3 1 2 1 3 3 1 2 3 1 1 2 3 2 2 1 1 3 3 3 2 3 2 1 1 3 2 2 3
 [75] 1 3 2 1 2 2 2 2 1 1 2 3 3 1 2 1 3 3 1 3 1 1 2 2 2 3 2 1 2 3 2 2 3 2 1 3 2
[112] 2 2 3 1 1 3 2 2 3 1 2 2 1 2 3 3 2 1 2 3 2 3 1 2 2 2 2 1 3 2 3 1 2 2 3 1 1
[149] 3 3 3 2 3 1 3 1 3 2 3 3 3 3 2 3 2 3 2 1 2 3 3 1 1 1 1 2 1 2 1 2 3 2 1 2 2
[186] 2 1 3 2 2 3 1 2 2 2 2 3 1 3 3 2 3 2 1 3 1 2 2 2 2 2 1 1 3 2 3 3 3 3 2 3 3
[223] 2 3 3 3 3 2 2 2 2 1 1 2 2 3 3 3 3 2 2 2 2 3 3 2 1 3 1 3 2 3 3 3 3 3 2 2 1
[260] 2 3 1 1 2 3 2 1 1 1 2 1 2 2 1 1 2 3 1 3 1 2 1 3 3 2 2 2 2 3 2 2 3 3 2 1 1
[297] 3 2 3 1 1 2 3 1 1 1 1 3 2 1 2 2 3 2 1 1 3 2 2 2 2 3 1 3 1 2 2 3 2 2 2 1 1
[334] 3 1 1 2 2 2 2 2 2 3 1 3 2 3 2 3 3 2 3 3 2 2 1 2 1 3 3 1 1 1 1 1 2 1 2 3 3
[371] 2 3 1 2 2 3 2 2 2 1 2 1 1 3 2 1 1 3 3 3 3 2 2 2 1 2 1 2 2 3 2 3 1 1 3 2 3
[408] 2 2 3 1 1 1 1 2 3 1 3 3 2 1 3 1 2 1 3 2 3 1 1 2 1 2 3 3 3 1 2 2 2 2 3 2 2
[445] 3 2 3 2 3 2 2 3 2 2 1 1 1 3 2 2 1 1 1 1 1 1 3 1 2 2 1 3 2 1 3 2 2 3 1 1 3
[482] 3 2 3 3 3 2 2 3 3 1 2 2 3 2 1 1 3 1 3 1 2 3 1 2 3 2 2 2 3 1 2 3 1 3 2 2 1
[519] 2 2 3 2 1 1 1 3 2 1 1 1 2 2 2 2 1 1 3 2 2 1 3 2 1 3 3 1 3 2 1 3 2 3 3 2 3
[556] 1 2 3 2 2 2 2 2 1 1 1 3 2 1 3 3 3 1 2 1 3 3 3 3 1 1 1 1 2 1 2 3 2 1 3 2 3
[593] 2 3 1 1 1 1 2 2

X = matrix(0, 600, 2) # X va recrevoir les données simulées et sont intialisées toutes à 0
#X

library(mvtnorm)

for(i in 1:600)
{
ifelse(cl[i] == 3, X[i,]<-rmvnorm(1, mu[,3], Sigma.ell), X[i,]<-rmvnorm(1, mu[,cl[i]], Sigma.sph))
}

2:Générer les données binaires et teraires

c1 = rbinom(600, 1, .5) + 1 # Génère une variable binaire qui peut valoir 1 ou 2


c2 = NULL

for(i in 1:600){

  ifelse(cl[i] == 1, c2[i]<-rbinom(1, 1, .1)+1, c2[i]<-rbinom(1, 1, .9)+1)

}

c3 = NULL

for(i in 1:600){

  if(cl[i] == 1) c3[i] = sample(1:3, 1, prob = c(.6, .2, .2))

  if(cl[i] == 2) c3[i] = sample(1:3, 1, prob = c(.1, .8, .1))

  if(cl[i] == 3) c3[i] = sample(1:3, 1, prob = c(.1, .1, .8))

}


c4 = NULL

for(i in 1:600){

  ifelse(cl[i] == 2, c4[i]<-rbinom(1, 1, .8)+1, c4[i]<-rbinom(1, 1, .2)+1)

}

vars = cbind(X, c1, c2, c3, c4)
#vars

library(tibble)
library(ggplot2)
vars|>class() # Pour connaîte le type d'objet

[1] "matrix" "array" 

vars|>as.tibble()|>ggplot(aes(x=V1,y=V2))+geom_point(col=cl)+
 ggtitle("Données continues: la couleur correspondant au groupe d'attribution")

Warning: `as.tibble()` was deprecated in tibble 2.0.0.
ℹ Please use `as_tibble()` instead.
ℹ The signature and semantics have changed, see `?as_tibble`.

Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.
ℹ The deprecated feature was likely used in the tibble package.
  Please report the issue at <https://github.com/tidyverse/tibble/issues>.

3:Distance euclidienne

#vars[1,1:2]
#vars[2,1:2]

(1.613755-0.3135524)^2+(4.678454 -(-0.3026124))^2

[1] 26.50155

dist(vars[1:2,1:2],method ="euclidean")^2

         1
2 26.50155

Représentation en 3D

result<-vars|>as.tibble()|>dplyr::select(1:2)
 result

# A tibble: 600 × 2
       V1     V2
    <dbl>  <dbl>
 1  1.61   4.68 
 2  0.314 -0.303
 3  3.14   3.77 
 4 -1.54  -0.153
 5  2.52   4.85 
 6  1.98   4.45 
 7  0.760  6.47 
 8  4.29   5.29 
 9  5.58   5.58 
10 -3.00   4.49 
# ℹ 590 more rows

 library(plotly)

Attachement du package : 'plotly'

L'objet suivant est masqué depuis 'package:MASS':

    select

L'objet suivant est masqué depuis 'package:ggplot2':

    last_plot

L'objet suivant est masqué depuis 'package:stats':

    filter

L'objet suivant est masqué depuis 'package:graphics':

    layout

library(plotly)

add_3Dbar <- function(p, x, y, z, width = 0.45){

  w <- width

  plotly::add_trace(

    p,

    type = "mesh3d",

    x = c(x-w, x-w, x+w, x+w, x-w, x-w, x+w, x+w),

    y = c(y-w, y+w, y+w, y-w, y-w, y+w, y+w, y-w),

    z = c(0, 0, 0, 0, z, z, z, z),

    i = c(7,0,0,0,4,4,2,6,4,0,3,7),

    j = c(3,4,1,2,5,6,5,5,0,1,2,2),

    k = c(0,7,2,3,6,7,1,2,5,5,7,6),

    opacity = 0.8,

    showscale = FALSE

  )

}

# exemple : histogramme 2D -> matrice de fréquences

x <- rnorm(500)

y <- rnorm(500)

bx <- seq(min(x), max(x), length.out = 12)

by <- seq(min(y), max(y), length.out = 12)

ix <- cut(x, breaks = bx, include.lowest = TRUE, labels = FALSE)

iy <- cut(y, breaks = by, include.lowest = TRUE, labels = FALSE)

tab <- table(ix, iy)

cx <- head(bx, -1) + diff(bx)/2

cy <- head(by, -1) + diff(by)/2

fig <- plot_ly()

for(i in seq_along(cx)){

  for(j in seq_along(cy)){

    z <- tab[i, j]

    if(z > 0){

      fig <- add_3Dbar(fig, cx[i], cy[j], z)

    }

  }

}

fig <- fig |>

  layout(title = "Histogramme 3D (Plotly) : effectifs par classe",

    scene = list(

      xaxis = list(title = "x"),

      yaxis = list(title = "y"),

      zaxis = list(title = "effectif")

    )

  )

fig

2 exercice deffinir les couleurs vert, noir et rouge

install.packages("MASS")

Warning: le package 'MASS' est en cours d'utilisation et ne sera pas installé

x_data <- result[[1]]

y_data <- result[[2]]
 
# 2. Calculer la densité de distribution (le "binning")

# n=20 définit la résolution (20x20 barres/points)
library(MASS)
densite <- kde2d(x_data, y_data, n = 30)
 
# 3. Créer le graphique 3D


plot_ly(x = densite$x, y = densite$y, z = densite$z) %>%

  add_surface(

    colorscale = 
list(
      list(0,   "green"),
      list(0.5, "black"),
      list(1,   "red")
    ) ,

    contours = list(

      z = list(show = TRUE, usecolormap = TRUE, highlightcolor = "#ff0000", project = list(z = TRUE))

    )

  ) %>%

  layout(

    title = "Histogramme de densité 3D",

    scene = list(

      xaxis = list(title = names(result)[1]),

      yaxis = list(title = names(result)[2]),

      zaxis = list(title = "Densité/Fréquence")

    )

  )

4:k-meams

result <- result |> mutate(couleur=cl)
result

# A tibble: 600 × 3
       V1     V2 couleur
    <dbl>  <dbl>   <dbl>
 1  1.61   4.68        2
 2  0.314 -0.303       1
 3  3.14   3.77        2
 4 -1.54  -0.153       1
 5  2.52   4.85        2
 6  1.98   4.45        2
 7  0.760  6.47        3
 8  4.29   5.29        2
 9  5.58   5.58        2
10 -3.00   4.49        3
# ℹ 590 more rows

kml = kmeans (result [,1:2], centers = 1, nstart = 50)

x <- result[,1]
y <- result[,2]
((x- mean(x$V1))^2+(y- mean(y$V2))^2)|>sum()

[1] 6481.454

km2 <- kmeans(result [,1:2], centers = 2, nstart = 50)

km3 <- kmeans(result [,1:2], centers = 3, nstart = 50)

km4 <- kmeans(result [,1:2], centers = 4, nstart = 50)

km5 <- kmeans(result [,1:2], centers = 5, nstart = 50)

km6 <- kmeans(result [,1:2], centers = 6, nstart = 50)

km7 <- kmeans(result [,1:2], centers = 7, nstart = 50)

km8 <- kmeans(result [,1:2], centers = 8, nstart = 50)

km9 <- kmeans(result [,1:2], centers = 9, nstart = 50)

5:Méthode du Coude

plot(
  1:9,
  c(kml$tot.withinss, km2$tot.withinss, km3$tot.withinss,
    km4$tot.withinss, km5$tot.withinss, km6$tot.withinss,
    km7$tot.withinss, km8$tot.withinss, km9$tot.withinss),
  type = 'b',
  pch = 19,
  xlab = 'Nombre de clusters (K)',
  ylab = 'Total Within-Cluster Sum of Squares',
  main = "Méthode du coude (Elbow Method)",
  xaxt = 'n',
  cex.axis = .75
)#tot.withinssc'est la somme des distances aux centres . Plus K augmante plus ça diminue

Une maniére plus moderne

library(factoextra)

Welcome to factoextra!

Want to learn more? See two factoextra-related books at https://www.datanovia.com/en/product/practical-guide-to-principal-component-methods-in-r/

#| label: fig-elbow
#| fig.height: 3.8
#| fig-cap: "Méthode du coude via factoextra. La ligne verticale pointillée marque l'inflexion (K = 3) suggérée."



library(factoextra)



## fviz_nbclust gère automatiquement la boucle de K=1 à 10 par défaut.
## On précise k.max = 9 
fviz_nbclust(result, kmeans, method = "wss", k.max = 9, nstart = 50) +
 # On ajoute manuellement la ligne verticale pour K=3
geom_vline(xintercept = 3, linetype = 2, color = "#B40000") +
 # Personnalisation des labels pour correspondre à votre style
 labs(
 title = "Inertie intra-classe selon le nombre de clusters K",
 x = "Nombre de clusters K",
 y = "WSS (Within-Cluster Sum of Squares)"
 ) +
theme_minimal()

3 les 3 clusters font réference aux 3 populations générées

question ? chercher une méthode moins rule of thumb

cl # les couleurs définis par la similation

  [1] 2 1 2 1 2 2 3 2 2 3 1 3 1 2 3 3 2 2 2 2 3 3 1 2 1 2 2 2 2 1 1 2 1 1 3 2 1
 [38] 1 2 1 2 3 2 1 2 3 1 2 1 3 3 1 2 3 1 1 2 3 2 2 1 1 3 3 3 2 3 2 1 1 3 2 2 3
 [75] 1 3 2 1 2 2 2 2 1 1 2 3 3 1 2 1 3 3 1 3 1 1 2 2 2 3 2 1 2 3 2 2 3 2 1 3 2
[112] 2 2 3 1 1 3 2 2 3 1 2 2 1 2 3 3 2 1 2 3 2 3 1 2 2 2 2 1 3 2 3 1 2 2 3 1 1
[149] 3 3 3 2 3 1 3 1 3 2 3 3 3 3 2 3 2 3 2 1 2 3 3 1 1 1 1 2 1 2 1 2 3 2 1 2 2
[186] 2 1 3 2 2 3 1 2 2 2 2 3 1 3 3 2 3 2 1 3 1 2 2 2 2 2 1 1 3 2 3 3 3 3 2 3 3
[223] 2 3 3 3 3 2 2 2 2 1 1 2 2 3 3 3 3 2 2 2 2 3 3 2 1 3 1 3 2 3 3 3 3 3 2 2 1
[260] 2 3 1 1 2 3 2 1 1 1 2 1 2 2 1 1 2 3 1 3 1 2 1 3 3 2 2 2 2 3 2 2 3 3 2 1 1
[297] 3 2 3 1 1 2 3 1 1 1 1 3 2 1 2 2 3 2 1 1 3 2 2 2 2 3 1 3 1 2 2 3 2 2 2 1 1
[334] 3 1 1 2 2 2 2 2 2 3 1 3 2 3 2 3 3 2 3 3 2 2 1 2 1 3 3 1 1 1 1 1 2 1 2 3 3
[371] 2 3 1 2 2 3 2 2 2 1 2 1 1 3 2 1 1 3 3 3 3 2 2 2 1 2 1 2 2 3 2 3 1 1 3 2 3
[408] 2 2 3 1 1 1 1 2 3 1 3 3 2 1 3 1 2 1 3 2 3 1 1 2 1 2 3 3 3 1 2 2 2 2 3 2 2
[445] 3 2 3 2 3 2 2 3 2 2 1 1 1 3 2 2 1 1 1 1 1 1 3 1 2 2 1 3 2 1 3 2 2 3 1 1 3
[482] 3 2 3 3 3 2 2 3 3 1 2 2 3 2 1 1 3 1 3 1 2 3 1 2 3 2 2 2 3 1 2 3 1 3 2 2 1
[519] 2 2 3 2 1 1 1 3 2 1 1 1 2 2 2 2 1 1 3 2 2 1 3 2 1 3 3 1 3 2 1 3 2 3 3 2 3
[556] 1 2 3 2 2 2 2 2 1 1 1 3 2 1 3 3 3 1 2 1 3 3 3 3 1 1 1 1 2 1 2 3 2 1 3 2 3
[593] 2 3 1 1 1 1 2 2

table(cl,km3$cluster)

   
cl    1   2   3
  1   0   0 177
  2 227  10   0
  3  33 153   0

Red<-result[1:5,1:2]

dist(Red)

         1        2        3        4
2 5.147966                           
3 1.776207 4.953916                  
4 5.769867 1.859623 6.103711         
5 0.919719 5.607630 1.252374 6.444347

# On a 5 individus 
# On procede de proche en proche 
# A chaque étape un seul individu est rattraché à un autre cluster 
# Il faut que ce rapprochement soit de distance minimale

#{1}{2}{3}{4}{5} commentaire
#{1,5}{2}{3}{4} commentaire
#{1,5,3}{2}{4} commentaire
#{1,5,3}{2,4} commentaire
#{1,5}{2}{3}{4} commentaire

hRed<-hclust(dist(Red),method="average")
hRed

Call:
hclust(d = dist(Red), method = "average")

Cluster method   : average 
Distance         : euclidean 
Number of objects: 5 

plot(hRed)

hc <- hclust(dist(result[1:2]),method="average")
plot(hc,main = "Dendrogramme de la classification hiérarchique ")

fviz_dend(hc,k=3,rect=TRUE,type="circular")+
  ggtitle("Dendrogramme circulaire des clusters")

cutree(hc,k=3) |> table() #coupe l'arbre en 3 groupes

  1   2   3 
419 177   4 

5:Modéles de mélanges gaussiens/finite mixture models

Ces méthodes consistent a ajuster les données observées avec des mélanges de distrubition normales (dites aussi gaussiennes).Nous allons ullustrer le cas bivarié ci-dessous .

\[ f(x,y)=\sum_{k=1}^{k}\tau_k f_k(x,y) \] Préciser tous ces élements de la formule ci-dessus

Algorithme

L’algorithme va estimer les parametres de cette distribution qui minimisent la vraisemblance de l’échantillon .il va le faire pour différente valeur de \(k\), par défaut de 1 à 9 . Le modéle retenu sera celui qui maximise le critére BIC

library(mclust)

Package 'mclust' version 6.1.2
Type 'citation("mclust")' for citing this R package in publications.

Attachement du package : 'mclust'

L'objet suivant est masqué depuis 'package:mvtnorm':

    dmvnorm

L'objet suivant est masqué depuis 'package:dplyr':

    count

L'objet suivant est masqué depuis 'package:purrr':

    map

ml <- Mclust (vars[, 1:2])
sml<-summary(ml)
ml$parameters

$pro
[1] 0.3878546 0.2945688 0.3175766

$mean
     [,1]        [,2]        [,3]
 2.992254  0.11328552 -0.03203617
 4.932755 -0.04727724  5.95321545

$variance
$variance$modelName
[1] "VVE"

$variance$d
[1] 2

$variance$G
[1] 3

$variance$sigma
, , 1

                    
 1.0722210 0.1869151
 0.1869151 0.9297125

, , 2

                          
 0.9665568283 0.0003461345
 0.0003461345 0.9662929272

, , 3

                  
 2.203563 1.439495
 1.439495 1.106058


$variance$scale
[1] 0.9807752 0.9664248 0.6042526

$variance$shape
          [,1]      [,2]      [,3]
[1,] 1.2245444 1.0003834 5.2881099
[2,] 0.8166302 0.9996168 0.1891035

$variance$orientation
                     
 0.8234704 -0.5673592
 0.5673592  0.8234704


$Vinv
NULL

ml$BIC

Bayesian Information Criterion (BIC): 
        EII       VII       EEI       VEI       EVI       VVI       EEE
1 -5448.591 -5448.591 -5377.939 -5377.939 -5377.939 -5377.939 -5313.801
2 -5034.370 -4943.649 -4912.760 -4855.237 -4896.074 -4832.905 -4918.559
3 -4809.654 -4790.897 -4813.893 -4785.884 -4808.834 -4795.366 -4769.950
4 -4750.342 -4735.761 -4739.474 -4741.631 -4744.989 -4752.285 -4708.504
5 -4751.013 -4735.634 -4746.494 -4732.024 -4761.346 -4727.257 -4725.861
6 -4744.222 -4718.074 -4742.745 -4716.427 -4763.885 -4743.413 -4730.691
7 -4749.425 -4742.531 -4739.313 -4738.624 -4766.467 -4773.469 -4667.456
8 -4765.684 -4760.166 -4750.299 -4764.463 -4789.868 -4798.474 -4686.120
9 -4759.392 -4721.664 -4760.715 -4727.735 -4785.921 -4769.058 -4697.642
        VEE       EVE       VVE       EEV       VEV       EVV       VVV
1 -5313.801 -5313.801 -5313.801 -5313.801 -5313.801 -5313.801 -5313.801
2 -4861.617 -4902.159 -4839.254 -4914.138 -4867.589 -4908.115 -4845.586
3 -4773.776 -4565.033 -4557.338 -4770.102 -4777.459 -4571.733 -4566.265
4 -4695.461 -4590.692 -4581.890 -4718.304 -4700.031 -4603.676 -4596.470
5 -4706.533 -4600.838 -4603.698 -4739.564 -4706.328 -4627.608 -4626.708
6 -4690.751 -4626.434 -4635.972 -4759.385 -4709.768 -4654.401 -4661.417
7 -4676.530 -4635.693 -4666.222 -4702.730 -4710.972 -4666.821 -4698.366
8 -4697.895 -4654.414 -4684.189 -4726.759 -4725.387 -4693.733 -4710.493
9 -4686.611 -4679.675 -4701.888 -4743.831 -4719.149 -4725.633 -4735.622

Top 3 models based on the BIC criterion: 
    VVE,3     EVE,3     VVV,3 
-4557.338 -4565.033 -4566.265 

ml$loglik

[1] -2230.692

mlv<-Mclust(vars[,1:2],verbose=FALSE)

Présenter les résultats de maniére plus esthétique en utulisant soit un kable soit stargazer

library(stargazer)

Please cite as: 

 Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.

 R package version 5.2.3. https://CRAN.R-project.org/package=stargazer 

stargazer(sml)

% Error: Unrecognized object type.

-Remarque: Verifier que les 600 observations ont été lobelées en fonction de leurs distances par rapport au vecteur moyen estimé.