Analyse des clusters dans R: déterminer le nombre optimal de clusters

Étant un débutant en R, je ne sais pas trop 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 dendro de cluster?

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 )

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

Tout d'abord, certaines données reproductibles (les données du Q ne sont pas 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 un coude dans la somme des éboulis d'erreurs au carré (SSE). Voir http://www.statmethods.net/advstats/cluster.html et 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 grappes 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 de medoids pour estimer le nombre de clusters à l'aide de la pamkfonction dans le package 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 Calinsky: une autre approche pour diagnostiquer le nombre de clusters adaptés 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 des attentes, initialisé par le clustering hiérarchique pour les modèles de mélange gaussiens 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 . Mise en cluster de la 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 . Écart statistique pour estimer le nombre de grappes. Voir aussi du code pour une sortie graphique agréable . Essayer 2 à 10 clusters 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 la mise en œuvre par Edwin Chen de la statistique de l'écart:

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

Huit . Le package 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 celles-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 / vues / 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

La pvclustbibliothèque qui calcule les valeurs de p pour le clustering hiérarchique via un rééchantillonnage bootstrap à plusieurs échelles est également pour les données de grande dimension . Voici l'exemple de la documentation (ne fonctionnera pas sur des données de dimension aussi faible que dans mon exemple):

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

Est-ce que tout cela aide?

Il est difficile d'ajouter quelque chose d'une réponse aussi élaborée. Bien que je pense que nous devrions mentionner identifyici, en particulier parce que @Ben montre de nombreux exemples de dendrogrammes.

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

identifyvous permet de choisir de manière interactive des clusters dans un dendrogramme et stocke vos choix dans une liste. Appuyez sur Échap pour quitter le mode interactif et revenir à la console R. Notez que la liste contient les indices, pas les noms de domaine (par opposition à cutree).

Afin de déterminer le k-cluster optimal dans les méthodes de clustering. J'utilise généralement la Elbowméthode accompagnée d'un traitement parallèle pour éviter de prendre du temps. Ce code peut échantillonner 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)
}

Courir le coude en 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)

Ça marche bien.

Splendide 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, où en général 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 en passant des messages entre les points de données." science 315,5814 (2007): 972-976.

Donc, si vous n'êtes pas biaisé vers k-means, je suggère d'utiliser directement AP, qui regroupera les données sans avoir besoin 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 utiliser une autre mesure de similitude fournie dans le même package. Par exemple, pour les similitudes basées sur les corrélations de Spearman, voici ce dont vous avez besoin:

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

Veuillez noter que ces fonctions de similitudes dans le package AP sont fournies à des fins de simplicité. En fait, la fonction apcluster () dans R acceptera toute matrice de corrélations. La même chose avec corSimMat () peut être effectuée avec ceci:

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

ou

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

en fonction de ce que vous souhaitez regrouper sur votre matrice (lignes ou cols).

Ces méthodes sont excellentes, mais lorsque vous essayez de trouver k pour des ensembles de données beaucoup plus grands, celles-ci peuvent être folles et lentes en R.

Une bonne solution que j'ai trouvée est le package "RWeka", qui a une implémentation efficace de l'algorithme X-Means - une version étendue de K-Means qui évolue mieux et déterminera le nombre optimal de clusters pour vous.

Vous devez d'abord vous assurer que Weka est installé sur votre système et que XMeans est installé via l'outil de gestion de paquets de Weka.

library(RWeka)

# Print a list of available options for the X-Means algorithm
WOW("XMeans")

# Create a Weka_control object which will specify our parameters
weka_ctrl <- Weka_control(
    I = 1000,                          # max no. of overall iterations
    M = 1000,                          # max no. of iterations in the kMeans loop
    L = 20,                            # min no. of clusters
    H = 150,                           # max no. of clusters
    D = "weka.core.EuclideanDistance", # distance metric Euclidean
    C = 0.4,                           # cutoff factor ???
    S = 12                             # random number seed (for reproducibility)
)

# Run the algorithm on your data, d
x_means <- XMeans(d, control = weka_ctrl)

# Assign cluster IDs to original data set
d$xmeans.cluster <- x_means$class_ids

Une solution simple est la bibliothèque factoextra. Vous pouvez modifier la méthode de clustering et la méthode de calcul du meilleur nombre de groupes. Par exemple, si vous voulez connaître le meilleur nombre de clusters pour un k- signifie:

Données: mtcars

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

Enfin, nous obtenons un graphique comme:

Les réponses sont excellentes. Si vous souhaitez donner une chance à une autre méthode de clustering, vous pouvez utiliser le clustering hiérarchique et voir comment les données sont fractionnées.

> set.seed(2)
> x=matrix(rnorm(50*2), ncol=2)
> hc.complete = hclust(dist(x), method="complete")
> plot(hc.complete)

Selon le nombre de classes dont vous avez besoin, vous pouvez couper votre dendrogramme comme;

> cutree(hc.complete,k = 2)
 [1] 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 2 1 1 1
[26] 2 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 2 1 1 1 1 1 1 1 2

Si vous tapez, ?cutreevous verrez les définitions. Si votre ensemble de données comprend trois classes, ce sera simplement cutree(hc.complete, k = 3). L'équivalent de cutree(hc.complete,k = 2)est cutree(hc.complete,h = 4.9).