Analisi del cluster in R: determina il numero ottimale di cluster

Essendo un novizio in R, non sono molto sicuro di come scegliere il miglior numero di cluster per fare un’analisi k-means. Dopo aver tracciato un sottoinsieme di dati sottostanti, quanti cluster saranno appropriati? Come posso eseguire l’analisi del dendro del 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 ) 

    Se la tua domanda è: how can I determine how many clusters are appropriate for a kmeans analysis of my data? , quindi ecco alcune opzioni. L’ articolo di wikipedia sulla determinazione del numero di cluster ha una buona recensione di alcuni di questi metodi.

    Innanzitutto, alcuni dati riproducibili (i dati nella Q sono … poco chiari per me):

     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) 

    inserisci la descrizione dell'immagine qui

    Uno . Cercare una piega o un gomito nella trama del complotto di sum di errore quadrato (SSE). Vedi http://www.statmethods.net/advstats/cluster.html e http://www.mattpeeples.net/kmeans.html per ulteriori informazioni. La posizione del gomito nella trama risultante suggerisce un numero adeguato di cluster per i kmea:

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

    Potremmo concludere che 4 cluster sarebbero indicati con questo metodo: inserisci la descrizione dell'immagine qui

    Due È ansible eseguire il partizionamento su medoids per stimare il numero di cluster utilizzando la funzione pamk nel pacchetto 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)) 

    inserisci la descrizione dell'immagine quiinserisci la descrizione dell'immagine qui

     # 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 

    Tre . Criterio di Calinsky: un altro approccio per diagnosticare quanti cluster soddisfano i dati. In questo caso proviamo da 1 a 10 gruppi.

     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! 

    inserisci la descrizione dell'immagine qui

    Quattro Determinare il modello e il numero di cluster ottimali in base al criterio di informazione bayesiano per massimizzare le aspettative, inizializzato dal clustering gerarchico per i modelli di miscela gaussiana parametrizzati

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

    inserisci la descrizione dell'immagine quiinserisci la descrizione dell'immagine quiinserisci la descrizione dell'immagine qui

    Cinque . Cluster di propagazione affinità (AP), vedi 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) 

    inserisci la descrizione dell'immagine quiinserisci la descrizione dell'immagine qui

    Sei . Gap Statistic for Estimating the Number of Clusters. Vedi anche qualche codice per un buon output grafico . Cercando 2-10 cluster qui:

     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 

    Ecco l'output dall'implementazione della statistica del divario da parte di Edwin Chen: inserisci la descrizione dell'immagine qui

    Sette . Potresti anche trovare utile esplorare i tuoi dati con clustergrams per visualizzare l'assegnazione dei cluster, vedi http://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r- codice / per maggiori dettagli.

    Otto . Il pacchetto NbClust fornisce 30 indici per determinare il numero di cluster in un set di dati.

     library(NbClust) nb < - NbClust(d, diss="NULL", distance = "euclidean", min.nc=2, max.nc=15, method = "kmeans", 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! 

    inserisci la descrizione dell'immagine qui

    Se la tua domanda è how can I produce a dendrogram to visualize the results of my cluster analysis , allora dovresti iniziare con questi: 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/ E vedi qui per altri esotici metodi: http://cran.r-project.org/web/views/Cluster.html

    Ecco alcuni esempi:

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

    inserisci la descrizione dell'immagine qui

     # 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 

    inserisci la descrizione dell'immagine qui

    Anche per i dati di dimensione elevata è la libreria pvclust che calcola i p-value per il clustering gerarchico tramite il ricampionamento bootstrap multiscala. Ecco l'esempio della documentazione (non lavorerai su dati di bassa dimensione come nel mio esempio):

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

    inserisci la descrizione dell'immagine qui

    Qualcuno di quello aiuta?

    È difficile aggiungere qualcosa di una risposta così elaborata. Anche se sento che dovremmo menzionare qui, in particolare perché @Ben mostra molti esempi di dendrogramma.

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

    Identify ti consente di selezionare in modo interattivo cluster da un dendrogramma e memorizza le tue scelte in un elenco. Premi Esc per uscire dalla modalità intertriggers e tornare alla console R. Si noti che la lista contiene gli indici, non i nomi dei giocatori (al contrario di un cutree ).

    Al fine di determinare il k-cluster ottimale nei metodi di clustering. Di solito uso il metodo Elbow accompagnato dall’elaborazione parallela per evitare il consumo di tempo. Questo codice può campionare in questo modo:

    Metodo del gomito

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

    Gomito in esecuzione parallelo

     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) 

    Funziona bene.

    Splendida risposta di Ben. Tuttavia sono sorpreso che il metodo Affinity Propagation (AP) sia stato qui suggerito solo per trovare il numero di cluster per il metodo k-means, dove in generale AP fa un lavoro migliore nel clustering dei dati. Si prega di consultare il documento scientifico che supporta questo metodo in Scienza qui:

    Frey, Brendan J. e Delbert Dueck. “Raggruppamento passando messaggi tra punti dati.” science 315.5814 (2007): 972-976.

    Quindi, se non sei prevenuto verso k-means, ti suggerisco di utilizzare direttamente l’AP, che raggrupperà i dati senza richiedere la conoscenza del numero di cluster:

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

    Se le distanze euclidee negative non sono appropriate, è ansible utilizzare altre misure di similarità fornite nello stesso pacchetto. Ad esempio, per le similitudini basate sulle correlazioni di Spearman, questo è ciò di cui hai bisogno:

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

    Si noti che queste funzioni per le somiglianze nel pacchetto AP sono solo fornite per semplicità. In effetti, la funzione apcluster () in R accetta qualsiasi matrice di correlazioni. Lo stesso prima con corSimMat () può essere fatto con questo:

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

    o

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

    a seconda di ciò che si desidera raggruppare sulla matrice (righe o colonne).

    Questi metodi sono grandiosi ma quando si cerca di trovare k per set di dati molto più grandi, questi possono essere follemente lenti in R.

    Una buona soluzione che ho trovato è il pacchetto “RWeka”, che ha un’efficiente implementazione dell’algoritmo X-Means – una versione estesa di K-Means che si adatta meglio e determinerà il numero ottimale di cluster per te.

    Prima di tutto, assicurati che Weka sia installato sul tuo sistema e che XMeans sia installato tramite lo strumento di gestione dei pacchetti di 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 

    Le risposte sono grandiose. Se vuoi dare una possibilità a un altro metodo di clustering puoi usare il clustering gerarchico e vedere come i dati si dividono.

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

    inserisci la descrizione dell'immagine qui

    A seconda del numero di classi di cui hai bisogno, puoi tagliare il tuo dendrogramma come;

     > 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 

    Se scrivi ?cutree vedrai le definizioni. Se il tuo set di dati ha tre classi sarà semplicemente cutree(hc.complete, k = 3) . L’equivalente per cutree(hc.complete,k = 2) è cutree(hc.complete,h = 4.9) .