Analyse par grappes en R : déterminer le nombre optimal de grappes

Question

Étant novice en R, je ne suis pas très sûr de savoir comment choisir le meilleur nombre de clusters pour faire une analyse k-means. Après avoir tracé un sous-ensemble de données ci-dessous, combien de clusters seront appropriés ? Comment puis-je effectuer une analyse dendrologique en grappes ?

n = 1000
kk = 10    
x1 = runif(kk)
y1 = runif(kk)
z1 = runif(kk)    
x4 = sample(x1,length(x1))
y4 = sample(y1,length(y1)) 
randObs <- function()
{
  ix = sample( 1:length(x4), 1 )
  iy = sample( 1:length(y4), 1 )
  rx = rnorm( 1, x4[ix], runif(1)/8 )
  ry = rnorm( 1, y4[ix], runif(1)/8 )
  return( c(rx,ry) )
}  
x = c()
y = c()
for ( k in 1:n )
{
  rPair  =  randObs()
  x  =  c( x, rPair[1] )
  y  =  c( y, rPair[2] )
}
z <- rnorm(n)
d <- data.frame( x, y, z )

Answer 1

Si votre question est how can I determine how many clusters are appropriate for a kmeans analysis of my data? alors voici quelques options. Le site article de wikipédia sur la détermination du nombre de clusters présente une bonne revue de certaines de ces méthodes.

Tout d'abord, quelques données reproductibles (les données dans le Q sont... peu claires pour moi) :

n = 100
g = 6 
set.seed(g)
d <- data.frame(x = unlist(lapply(1:g, function(i) rnorm(n/g, runif(1)*i^2))), 
                y = unlist(lapply(1:g, function(i) rnorm(n/g, runif(1)*i^2))))
plot(d)

Un . Recherchez un coude ou une courbe dans le diagramme de dispersion de la somme des erreurs au carré (SSE). Voir http://www.statmethods.net/advstats/cluster.html & http://www.mattpeeples.net/kmeans.html pour en savoir plus. L'emplacement du coude dans le graphique résultant suggère un nombre approprié de clusters pour les kmeans :

mydata <- d
wss <- (nrow(mydata)-1)*sum(apply(mydata,2,var))
  for (i in 2:15) wss[i] <- sum(kmeans(mydata,
                                       centers=i)$withinss)
plot(1:15, wss, type="b", xlab="Number of Clusters",
     ylab="Within groups sum of squares")

Nous pourrions conclure que 4 clusters seraient indiqués par cette méthode :

Deux . Vous pouvez effectuer un partitionnement autour des médoïdes pour estimer le nombre de clusters en utilisant la fonction pamk dans le paquet fpc.

library(fpc)
pamk.best <- pamk(d)
cat("number of clusters estimated by optimum average silhouette width:", pamk.best$nc, "\n")
plot(pam(d, pamk.best$nc))

# we could also do:
library(fpc)
asw <- numeric(20)
for (k in 2:20)
  asw[[k]] <- pam(d, k) $ silinfo $ avg.width
k.best <- which.max(asw)
cat("silhouette-optimal number of clusters:", k.best, "\n")
# still 4

Trois . Critère de Calinsky : Une autre approche pour diagnostiquer combien de clusters conviennent aux données. Dans ce cas, nous essayons de 1 à 10 groupes.

require(vegan)
fit <- cascadeKM(scale(d, center = TRUE,  scale = TRUE), 1, 10, iter = 1000)
plot(fit, sortg = TRUE, grpmts.plot = TRUE)
calinski.best <- as.numeric(which.max(fit$results[2,]))
cat("Calinski criterion optimal number of clusters:", calinski.best, "\n")
# 5 clusters!

Quatre . Déterminer le modèle optimal et le nombre de clusters selon le critère d'information bayésien pour la maximisation de l'espérance, initialisé par le clustering hiérarchique pour les modèles de mélange gaussien paramétrés.

# See http://www.jstatsoft.org/v18/i06/paper
# http://www.stat.washington.edu/research/reports/2006/tr504.pdf
#
library(mclust)
# Run the function to see how many clusters
# it finds to be optimal, set it to search for
# at least 1 model and up 20.
d_clust <- Mclust(as.matrix(d), G=1:20)
m.best <- dim(d_clust$z)[2]
cat("model-based optimal number of clusters:", m.best, "\n")
# 4 clusters
plot(d_clust)

Cinq . Regroupement par propagation d'affinité (AP), voir http://dx.doi.org/10.1126/science.1136800

library(apcluster)
d.apclus <- apcluster(negDistMat(r=2), d)
cat("affinity propogation optimal number of clusters:", length(d.apclus@clusters), "\n")
# 4
heatmap(d.apclus)
plot(d.apclus, d)

Six . Statistique d'écart pour l'estimation du nombre de clusters. Voir aussi un peu de code pour une sortie graphique agréable . J'essaie les groupes de 2 à 10 ici :

library(cluster)
clusGap(d, kmeans, 10, B = 100, verbose = interactive())

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 100)  [one "." per sample]:
.................................................. 50 
.................................................. 100 
Clustering Gap statistic ["clusGap"].
B=100 simulated reference sets, k = 1..10
 --> Number of clusters (method 'firstSEmax', SE.factor=1): 4
          logW   E.logW        gap     SE.sim
 [1,] 5.991701 5.970454 -0.0212471 0.04388506
 [2,] 5.152666 5.367256  0.2145907 0.04057451
 [3,] 4.557779 5.069601  0.5118225 0.03215540
 [4,] 3.928959 4.880453  0.9514943 0.04630399
 [5,] 3.789319 4.766903  0.9775842 0.04826191
 [6,] 3.747539 4.670100  0.9225607 0.03898850
 [7,] 3.582373 4.590136  1.0077628 0.04892236
 [8,] 3.528791 4.509247  0.9804556 0.04701930
 [9,] 3.442481 4.433200  0.9907197 0.04935647
[10,] 3.445291 4.369232  0.9239414 0.05055486

Voici le résultat de l'implémentation de la statistique de l'écart par Edwin Chen :

Sept . Vous pouvez également trouver utile d'explorer vos données avec des clustergrammes pour visualiser l'affectation des clusters, voir http://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r-code/ pour plus de détails.

Huit . Le site Paquet NbClust fournit 30 indices pour déterminer le nombre de clusters dans un ensemble de données.

library(NbClust)
nb <- NbClust(d, diss=NULL, distance = "euclidean",
        method = "kmeans", min.nc=2, max.nc=15, 
        index = "alllong", alphaBeale = 0.1)
hist(nb$Best.nc[1,], breaks = max(na.omit(nb$Best.nc[1,])))
# Looks like 3 is the most frequently determined number of clusters
# and curiously, four clusters is not in the output at all!

Si votre question est how can I produce a dendrogram to visualize the results of my cluster analysis alors vous devriez commencer par ceux-ci : http://www.statmethods.net/advstats/cluster.html http://www.r-tutor.com/gpu-computing/clustering/hierarchical-cluster-analysis http://gastonsanchez.wordpress.com/2012/10/03/7-ways-to-plot-dendrograms-in-r/ Et voyez ici pour des méthodes plus exotiques : http://cran.r-project.org/web/views/Cluster.html

Voici quelques exemples :

d_dist <- dist(as.matrix(d))   # find distance matrix 
plot(hclust(d_dist))           # apply hirarchical clustering and plot

# a Bayesian clustering method, good for high-dimension data, more details:
# http://vahid.probstat.ca/paper/2012-bclust.pdf
install.packages("bclust")
library(bclust)
x <- as.matrix(d)
d.bclus <- bclust(x, transformed.par = c(0, -50, log(16), 0, 0, 0))
viplot(imp(d.bclus)$var); plot(d.bclus); ditplot(d.bclus)
dptplot(d.bclus, scale = 20, horizbar.plot = TRUE,varimp = imp(d.bclus)$var, horizbar.distance = 0, dendrogram.lwd = 2)
# I just include the dendrogram here

De même, pour les données à haute dimension, le pvclust bibliothèque qui calcule les valeurs p pour le regroupement hiérarchique via un rééchantillonnage bootstrap multi-échelles. Voici l'exemple tiré de la documentation (ne fonctionnera pas sur des données de faible dimension comme dans mon exemple) :

library(pvclust)
library(MASS)
data(Boston)
boston.pv <- pvclust(Boston)
plot(boston.pv)

Est-ce que tout cela vous aide ?

Answer 2

Il est difficile d'ajouter quelque chose à une réponse aussi élaborée. Cependant, je pense que nous devrions mentionner identify ici, notamment parce que @Ben montre beaucoup d'exemples de dendrogrammes.

d_dist <- dist(as.matrix(d))   # find distance matrix 
plot(hclust(d_dist)) 
clusters <- identify(hclust(d_dist))

identify vous permet de choisir interactivement des clusters dans un dendrogramme et enregistre vos choix dans une liste. Appuyez sur Esc pour quitter le mode interactif et revenir à la console R. Notez que la liste contient les indices, pas les noms de rown (par opposition à cutree ).

Answer 3

Afin de déterminer le k-cluster optimal dans les méthodes de clustering. J'utilise généralement Elbow accompagné d'un traitement parallèle pour éviter les pertes de temps. Ce code peut être échantillonné comme ceci :

Méthode du coude

elbow.k <- function(mydata){
dist.obj <- dist(mydata)
hclust.obj <- hclust(dist.obj)
css.obj <- css.hclust(dist.obj,hclust.obj)
elbow.obj <- elbow.batch(css.obj)
k <- elbow.obj$k
return(k)
}

Coude parallèle

no_cores <- detectCores()
    cl<-makeCluster(no_cores)
    clusterEvalQ(cl, library(GMD))
    clusterExport(cl, list("data.clustering", "data.convert", "elbow.k", "clustering.kmeans"))
 start.time <- Sys.time()
 elbow.k.handle(data.clustering))
 k.clusters <- parSapply(cl, 1, function(x) elbow.k(data.clustering))
    end.time <- Sys.time()
    cat('Time to find k using Elbow method is',(end.time - start.time),'seconds with k value:', k.clusters)

Cela fonctionne bien.

Answer 4

Une solution simple est la bibliothèque factoextra . Vous pouvez changer la méthode de clustering et la méthode pour calculer le meilleur nombre de groupes. Par exemple, si vous voulez connaître le meilleur nombre de groupes pour une moyenne k :

Données : mtcars

library(factoextra)   
fviz_nbclust(mtcars, kmeans, method = "wss") +
      geom_vline(xintercept = 3, linetype = 2)+
      labs(subtitle = "Elbow method")

Finalement, nous obtenons un graphique comme :

Answer 5

Excellente réponse de Ben. Cependant je suis surpris que la méthode de Propagation d'Affinité (AP) ait été suggérée ici juste pour trouver le nombre de cluster pour la méthode k-means, alors qu'en général la méthode AP fait un meilleur travail de clustering des données. Veuillez consulter l'article scientifique soutenant cette méthode dans Science ici :

Frey, Brendan J., et Delbert Dueck. "Clustering by passing messages between data points". science 315.5814 (2007) : 972-976.

Donc, si vous n'avez pas de préjugés à l'égard de k-means, je vous suggère d'utiliser directement AP, qui regroupera les données sans qu'il soit nécessaire de connaître le nombre de clusters :

library(apcluster)
apclus = apcluster(negDistMat(r=2), data)
show(apclus)

Si les distances euclidiennes négatives ne sont pas appropriées, vous pouvez alors utiliser d'autres mesures de similarité fournies dans le même paquet. Par exemple, pour les similarités basées sur les corrélations de Spearman, c'est ce dont vous avez besoin :

sim = corSimMat(data, method="spearman")
apclus = apcluster(s=sim)

Veuillez noter que ces fonctions pour les similarités dans le paquet AP sont juste fournies pour la simplicité. En fait, la fonction apcluster() dans R acceptera n'importe quelle matrice de corrélations. La même chose que précédemment avec corSimMat() peut être faite avec ceci :

sim = cor(data, method="spearman")

ou

sim = cor(t(data), method="spearman")

en fonction de ce que vous voulez regrouper sur votre matrice (lignes ou colonnes).