Análisis de conglomerados en R: determine la cantidad óptima de conglomerados

Siendo novato en R, no estoy muy seguro de cómo elegir el mejor número de clusters para hacer un análisis de k-medias. Después de trazar un subconjunto de los datos a continuación, ¿cuántos conglomerados serán apropiados? ¿Cómo puedo realizar el análisis de dendro en clúster?

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 su pregunta es how can I determine how many clusters are appropriate for a kmeans analysis of my data? , aquí hay algunas opciones. El artículo de wikipedia sobre determinación de números de clusters tiene una buena revisión de algunos de estos métodos.

Primero, algunos datos reproducibles (los datos en la Q son … confusos para mí):

 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) 

enter image description here

Uno . Busque un doblez o un codo en la sum del diagtwig de sedimento de error cuadrado (SSE). Consulte http://www.statmethods.net/advstats/cluster.html y http://www.mattpeeples.net/kmeans.html para obtener más información. La ubicación del codo en la gráfica resultante sugiere una cantidad adecuada de grupos para los 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") 

Podríamos concluir que 4 clústeres estarían indicados por este método: enter image description here

Dos . Puede hacer particiones alrededor de medoids para estimar el número de clusters usando la función pamk en el paquete 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)) 

enter image description hereenter image description here

 # 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 

Tres . Criterio Calinsky: otro enfoque para diagnosticar cuántos clústeres se ajustan a los datos. En este caso, intentamos de 1 a 10 grupos.

 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! 

enter image description here

Cuatro . Determine el modelo óptimo y el número de conglomerados de acuerdo con el Criterio de información bayesiano para la maximización de las expectativas, inicializado por agrupamiento jerárquico para modelos de mezclas gaussianas parametrizadas.

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

enter image description hereenter image description hereenter image description here

Cinco . Agrupación de afinidad por propagación (AP), ver 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) 

enter image description hereenter image description here

Seis . Estadística de brecha para estimar el número de conglomerados. Vea también algunos códigos para una buena salida gráfica . Probando 2-10 clusters aquí:

 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 

Aquí está el resultado de la implementación de Edwin Chen de la estadística de brecha: enter image description here

Siete . También puede resultarle útil explorar sus datos con clustergrams para visualizar la asignación de clústeres, consulte http://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r- código / para más detalles.

Ocho . El paquete NbClust proporciona 30 índices para determinar el número de clústeres en un conjunto de datos.

 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! 

enter image description here

Si su pregunta es how can I produce a dendrogram to visualize the results of my cluster analysis , entonces debe comenzar con estos: 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/ Y mira aquí para ver más exótico métodos: http://cran.r-project.org/web/views/Cluster.html

Aquí están algunos ejemplos:

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

enter image description here

 # 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 

enter image description here

También para datos de gran dimensión se encuentra la biblioteca pvclust que calcula los valores p para la agrupación jerárquica a través del remuestreo de bootstrap multiescala. Aquí está el ejemplo de la documentación (no funcionará con datos dimensionales tan bajos como en mi ejemplo):

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

enter image description here

¿Algo de eso ayuda?

Es difícil agregar algo también una respuesta tan elaborada. Aunque creo que deberíamos mencionar aquí, especialmente porque @Ben muestra muchos ejemplos de dendrogtwig.

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

identity le permite elegir interactivamente clusters de un dendrogtwig y almacena sus elecciones en una lista. Presione Esc para salir del modo interactivo y regresar a la consola R. Tenga en cuenta que la lista contiene los índices, no los nombres de fila (a diferencia de cutree ).

Para determinar el k-cluster óptimo en métodos de agrupamiento. Usualmente uso el método de Elbow para acompañar el procesamiento en paralelo para evitar el consumo de tiempo. Este código puede muestrear de esta manera:

Método de codo

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

Ejecución de codo paralelo

 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) 

Funciona bien.

Espléndida respuesta de Ben. Sin embargo, estoy sorprendido de que el método de Propagación de afinidad (AP) haya sido sugerido aquí solo para encontrar el número de clúster para el método k-means, donde en general AP hace un mejor trabajo al agrupar los datos. Por favor, consulte el artículo científico que apoya este método en Science aquí:

Frey, Brendan J. y Delbert Dueck. “Agrupación al pasar mensajes entre puntos de datos”. science 315.5814 (2007): 972-976.

Por lo tanto, si no está predispuesto a k-means, le sugiero que use AP directamente, lo que agrupará los datos sin necesidad de conocer la cantidad de clústeres:

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

Si las distancias euclidianas negativas no son apropiadas, entonces puede usar otras medidas de similitud proporcionadas en el mismo paquete. Por ejemplo, para similitudes basadas en correlaciones de Spearman, esto es lo que necesita:

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

Tenga en cuenta que esas funciones para las similitudes en el paquete AP se proporcionan simplemente por simplicidad. De hecho, la función apcluster () en R aceptará cualquier matriz de correlaciones. Lo mismo antes con corSimMat () se puede hacer con esto:

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

o

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

dependiendo de lo que quieras agrupar en tu matriz (filas o columnas).

Estos métodos son geniales, pero al tratar de encontrar k para conjuntos de datos mucho más grandes, estos pueden ser muy lentos en R.

Una buena solución que he encontrado es el paquete “RWeka”, que tiene una implementación eficiente del algoritmo X-Means, una versión extendida de K-Means que se escala mejor y determinará la cantidad óptima de clústeres para usted.

En primer lugar, querrá asegurarse de que Weka esté instalado en su sistema y tenga instalado XMeans a través de la herramienta de gestión de paquetes 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 

Las respuestas son geniales Si desea dar la oportunidad a otro método de agrupamiento, puede usar la agrupación jerárquica y ver cómo se dividen los datos.

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

enter image description here

Dependiendo de cuántas clases necesite, puede cortar su dendrogtwig como;

 > 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 escribe ?cutree , verá las definiciones. Si su conjunto de datos tiene tres clases, será simplemente cutree(hc.complete, k = 3) . El equivalente para cutree(hc.complete,k = 2) es cutree(hc.complete,h = 4.9) .