Clusteranalyse in R: Bestimmen Sie die optimale Anzahl von Clustern

428

Als Neuling in R bin ich mir nicht sicher, wie ich die beste Anzahl von Clustern für eine k-means-Analyse auswählen soll. Wie viele Cluster sind nach dem Zeichnen einer Teilmenge der folgenden Daten geeignet? Wie kann ich eine Cluster-Dendro-Analyse durchführen?

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 )
user2153893
quelle
4
Wenn Sie nicht vollständig mit kmeans verbunden sind, können Sie den im fpcPaket verfügbaren DBSCAN-Clustering-Algorithmus ausprobieren . Es ist wahr, Sie müssen dann zwei Parameter einstellen ... aber ich habe festgestellt, dass es fpc::dbscandann ziemlich gute Arbeit leistet, automatisch eine gute Anzahl von Clustern zu bestimmen. Außerdem kann tatsächlich ein einzelner Cluster ausgegeben werden, wenn die Daten dies aussagen. Einige der Methoden in @ Bens hervorragenden Antworten helfen Ihnen nicht dabei, festzustellen, ob k = 1 tatsächlich am besten ist.
Stephan Kolassa
Siehe auch stats.stackexchange.com/q/11691/478
Richie Cotton

Antworten:

1020

Wenn Ihre Frage ist how can I determine how many clusters are appropriate for a kmeans analysis of my data?, dann sind hier einige Optionen. Der Wikipedia-Artikel zur Bestimmung der Anzahl von Clustern bietet einen guten Überblick über einige dieser Methoden.

Erstens einige reproduzierbare Daten (die Daten im Q sind ... für mich unklar):

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)

Geben Sie hier die Bildbeschreibung ein

Eins . Suchen Sie nach einer Biegung oder einem Ellbogen in der Summe der SSE-Geröllkurven (Squared Error). Weitere Informationen finden Sie unter http://www.statmethods.net/advstats/cluster.html und http://www.mattpeeples.net/kmeans.html . Die Position des Ellbogens in der resultierenden Darstellung legt eine geeignete Anzahl von Clustern für die Kilometer nahe:

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")

Wir könnten daraus schließen, dass 4 Cluster durch diese Methode angezeigt würden: Geben Sie hier die Bildbeschreibung ein

Zwei . Sie können mithilfe der pamkFunktion im fpc-Paket eine Partitionierung um Medoide durchführen, um die Anzahl der Cluster zu schätzen .

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))

Geben Sie hier die Bildbeschreibung ein Geben Sie hier die Bildbeschreibung ein

# 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

Drei . Calinsky-Kriterium: Ein weiterer Ansatz zur Diagnose, wie viele Cluster zu den Daten passen. In diesem Fall versuchen wir 1 bis 10 Gruppen.

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!

Geben Sie hier die Bildbeschreibung ein

Vier . Bestimmen Sie das optimale Modell und die Anzahl der Cluster gemäß dem Bayes'schen Informationskriterium für die Erwartungsmaximierung, das durch hierarchisches Clustering für parametrisierte Gaußsche Mischungsmodelle initialisiert wird

# 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)

Geben Sie hier die Bildbeschreibung ein Geben Sie hier die Bildbeschreibung ein Geben Sie hier die Bildbeschreibung ein

Fünf . Affinity Propagation (AP) -Clustering, siehe 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)

Geben Sie hier die Bildbeschreibung ein Geben Sie hier die Bildbeschreibung ein

Sechs . Lückenstatistik zur Schätzung der Anzahl von Clustern. Siehe auch Code für eine schöne grafische Ausgabe . Versuchen Sie hier 2-10 Cluster:

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

Hier ist die Ausgabe von Edwin Chens Implementierung der Lückenstatistik: Geben Sie hier die Bildbeschreibung ein

Sieben . Es kann auch nützlich sein, Ihre Daten mit Clustergrammen zu untersuchen, um die Clusterzuweisung zu visualisieren. Siehe http://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r- Code / für weitere Details.

Acht . Das NbClust-Paket enthält 30 Indizes, um die Anzahl der Cluster in einem Dataset zu bestimmen.

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!

Geben Sie hier die Bildbeschreibung ein

Wenn Ihre Frage lautet how can I produce a dendrogram to visualize the results of my cluster analysis, sollten Sie mit folgenden Punkten beginnen: 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/ Weitere exotische Methoden finden Sie hier: http://cran.r-project.org/ web / views / Cluster.html

Hier einige Beispiele:

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

Geben Sie hier die Bildbeschreibung ein

# 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

Geben Sie hier die Bildbeschreibung ein

Ebenfalls für hochdimensionale Daten ist die pvclustBibliothek vorgesehen, die p-Werte für hierarchisches Clustering über Multiskalen-Bootstrap-Resampling berechnet. Hier ist das Beispiel aus der Dokumentation (funktioniert nicht mit so niedrigdimensionalen Daten wie in meinem Beispiel):

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

Geben Sie hier die Bildbeschreibung ein

Hilft irgendetwas davon?

Ben
quelle
Für das letzte Dendogramm (Cluster-Dendogramm mit AU / BP) ist es manchmal bequem, Rechtecke um die Gruppen mit relativ hohen p-Werten zu zeichnen: pvrect (fit, alpha = 0,95)
Igor Elbert
Genau das habe ich gesucht. Ich bin neu bei R und es hätte sehr lange gedauert, bis ich das gefunden hätte. Vielen Dank an @Ben für die ausführliche Antwort. Können Sie mich bitte anleiten, wo ich die Logik hinter jeder dieser Methoden finden kann, z. B. welche Metrik oder welches Kriterium sie verwenden, um die optimale Anzahl von Clustern zu bestimmen, oder wie sich die einzelnen Methoden voneinander unterscheiden. Mein Chef möchte, dass ich das sage, damit wir entscheiden können, welche der Methoden wir verwenden möchten. Danke im Voraus.
Nasia Jaffri
1
@Aleksandr Blekh Sie können auch versuchen, eine beliebige grafische Methode in die Analyse umzuwandeln. ZB verwende ich die "Ellbogen" -Methode (zuerst in der Antwort erwähnt), versuche aber, sie analytisch zu finden. Der Ellbogenpunkt könnte der Punkt mit maximaler Krümmung sein. Für diskrete Daten ist dies der Punkt mit der maximalen zentralen Differenz zweiter Ordnung (analog zur maximalen Ableitung zweiter Ordnung für kontinuierliche Daten). Siehe stackoverflow.com/a/4473065/1075993 und stackoverflow.com/q/2018178/1075993 . Ich denke, dass andere grafische Methoden auch in analytische umgewandelt werden könnten.
Andrey Sapegin
1
@AndreySapegin: Ich könnte, aber: 1) Ehrlich gesagt halte ich es nicht für eine elegante Lösung (IMHO sollten visuelle Methoden in den meisten Fällen visuell bleiben, während analytische analytisch bleiben sollten); 2) Ich habe eine analytische Lösung für dieses Problem mit einem oder mehreren RPaketen gefunden (es befindet sich auf meinem GitHub - Sie können es sich gerne ansehen). 3) Meine Lösung scheint gut genug zu funktionieren. Außerdem ist es eine Weile her und ich habe meine Dissertationssoftware und meinen Dissertationsbericht (Abschlussarbeit) bereits fertiggestellt und bereite mich derzeit auf die Verteidigung vor :-). Unabhängig davon freue ich mich sehr über Ihren Kommentar und Ihre Links. Alles Gute!
Aleksandr Blekh
1
In meinem aktuellen Clustering-Datensatz befinden sich 2,2 Millionen Zeilen. Ich gehe davon aus, dass keines dieser R-Pakete darauf funktioniert. Sie knallen einfach meinen Computer und dann fällt es aus meiner Erfahrung um. Es sieht jedoch so aus, als ob der Autor seine Kenntnisse für kleine Daten und für den allgemeinen Fall ohne Rücksicht auf die Softwarekapazität kennt. Keine Punkte wegen offensichtlich guter Arbeit des Autors abgezogen. Ihr wisst bitte nur, dass das alte R bei 2,2 Millionen Zeilen schrecklich ist - probiert es selbst aus, wenn ihr mir nicht vertraut. H2O hilft, beschränkt sich aber auf einen kleinen ummauerten Garten von Happy.
Geoffrey Anderson
21

Es ist schwer, etwas zu einer so ausführlichen Antwort hinzuzufügen. Obwohl ich identifydenke, wir sollten dies hier erwähnen , insbesondere weil @Ben viele Dendrogrammbeispiele zeigt.

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

identifyMit dieser Option können Sie Cluster aus einem Dendrogramm interaktiv auswählen und Ihre Auswahl in einer Liste speichern. Drücken Sie die Esc-Taste, um den interaktiven Modus zu verlassen und zur R-Konsole zurückzukehren. Beachten Sie, dass die Liste die Indizes enthält, nicht die Rownamen (im Gegensatz zu cutree).

Matt Bannert
quelle
10

Um den optimalen k-Cluster in Clustering-Methoden zu bestimmen. Normalerweise verwende ich eine ElbowMethode, die von der Parallelverarbeitung begleitet wird, um Zeitaufwand zu vermeiden. Dieser Code kann wie folgt aussehen:

Ellbogenmethode

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)
}

Ellbogen parallel laufen lassen

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)

Es funktioniert gut.

VanThaoNguyen
quelle
2
Die Ellbogen- und CSS-Funktionen stammen aus dem GMD-Paket: cran.r-project.org/web/packages/GMD/GMD.pdf
Rohan
6

Herrliche Antwort von Ben. Ich bin jedoch überrascht, dass hier die Affinity Propagation (AP) -Methode vorgeschlagen wurde, nur um die Anzahl der Cluster für die k-means-Methode zu ermitteln, bei der AP im Allgemeinen die Daten besser gruppiert. Bitte lesen Sie das wissenschaftliche Papier, das diese Methode in Science unterstützt, hier:

Frey, Brendan J. und Delbert Dueck. "Clustering durch Weiterleiten von Nachrichten zwischen Datenpunkten." science 315.5814 (2007): 972 & ndash; 976.

Wenn Sie also nicht auf k-means eingestellt sind, empfehle ich, AP direkt zu verwenden, wodurch die Daten geclustert werden, ohne dass die Anzahl der Cluster bekannt sein muss:

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

Wenn negative euklidische Abstände nicht angemessen sind, können Sie andere Ähnlichkeitsmaße verwenden, die im selben Paket enthalten sind. Für Ähnlichkeiten, die auf Spearman-Korrelationen basieren, benötigen Sie beispielsweise Folgendes:

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

Bitte beachten Sie, dass diese Funktionen für Ähnlichkeiten im AP-Paket der Einfachheit halber nur bereitgestellt werden. Tatsächlich akzeptiert die Funktion apcluster () in R jede Korrelationsmatrix. Das gleiche vorher mit corSimMat () kann damit gemacht werden:

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

oder

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

abhängig davon, was Sie in Ihrer Matrix gruppieren möchten (Zeilen oder Spalten).

zsram
quelle
6

Diese Methoden sind großartig, aber wenn Sie versuchen, k für viel größere Datensätze zu finden, können diese in R verrückt langsam sein.

Eine gute Lösung, die ich gefunden habe, ist das "RWeka" -Paket, das eine effiziente Implementierung des X-Means-Algorithmus enthält - eine erweiterte Version von K-Means, die besser skaliert und die optimale Anzahl von Clustern für Sie bestimmt.

Zuerst möchten Sie sicherstellen, dass Weka auf Ihrem System installiert ist und XMeans über das Paketmanager-Tool von Weka installiert ist.

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
RDRR
quelle
6

Eine einfache Lösung ist die Bibliothek factoextra. Sie können die Clustering-Methode und die Methode zur Berechnung der besten Anzahl von Gruppen ändern. Wenn Sie beispielsweise die beste Anzahl von Clustern für ein k-Mittel wissen möchten:

Daten: mtcars

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

Schließlich erhalten wir eine Grafik wie:

Geben Sie hier die Bildbeschreibung ein

Cro-Magnon
quelle
2

Die Antworten sind großartig. Wenn Sie einer anderen Clustering-Methode eine Chance geben möchten, können Sie hierarchisches Clustering verwenden und sehen, wie Daten aufgeteilt werden.

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

Geben Sie hier die Bildbeschreibung ein

Je nachdem, wie viele Klassen Sie benötigen, können Sie Ihr Dendrogramm wie folgt schneiden:

> 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

Wenn Sie eingeben ?cutree, werden die Definitionen angezeigt. Wenn Ihr Datensatz drei Klassen hat, ist es einfach cutree(hc.complete, k = 3). Das Äquivalent für cutree(hc.complete,k = 2)ist cutree(hc.complete,h = 4.9).

Boyaronur
quelle
Ich bevorzuge Wards gegenüber Complete.
Chris