Análisis de conglomerados en R: determine el número óptimo de conglomerados

428

Al ser un novato en R, no estoy muy seguro de cómo elegir el mejor número de grupos para hacer un análisis de k-medias. Después de trazar un subconjunto de datos a continuación, ¿cuántos grupos serán apropiados? ¿Cómo puedo realizar el análisis dendro de 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 )
usuario2153893
fuente
44
Si no está completamente casado con kmeans, puede probar el algoritmo de agrupación DBSCAN, disponible en el fpcpaquete. Es cierto, luego debe establecer dos parámetros ... pero descubrí que fpc::dbscanhace un trabajo bastante bueno para determinar automáticamente un buen número de clústeres. Además, en realidad puede generar un solo clúster si eso es lo que le dicen los datos: algunos de los métodos en las excelentes respuestas de @ Ben no lo ayudarán a determinar si k = 1 es realmente el mejor.
Stephan Kolassa
Ver también stats.stackexchange.com/q/11691/478
Richie Cotton

Respuestas:

1020

Si su pregunta es how can I determine how many clusters are appropriate for a kmeans analysis of my data?, entonces aquí hay algunas opciones. El artículo de Wikipedia sobre la determinación de la cantidad de grupos tiene una buena revisión de algunos de estos métodos.

Primero, algunos datos reproducibles (los datos en la Q son ... poco claros 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)

ingrese la descripción de la imagen aquí

Uno . Busque una curva o codo en la gráfica de suma de error al 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 un número adecuado de grupos para los medios:

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 grupos serían indicados por este método: ingrese la descripción de la imagen aquí

Dos . Puede realizar particiones alrededor de medoides para estimar el número de clústeres utilizando la pamkfunción 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))

ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí

# 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 de Calinsky: otro enfoque para diagnosticar cuántos grupos se adaptan a los datos. En este caso tratamos 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!

ingrese la descripción de la imagen aquí

Cuatro . Determine el modelo óptimo y el número de grupos según el Criterio de información bayesiano para la maximización de expectativas, inicializado por agrupamiento jerárquico para modelos de mezcla gaussiana parametrizados

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

ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí

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

ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí

Seis . Estadística de brecha para estimar el número de grupos. Vea también algo de código para una buena salida gráfica . Tratando de 2 a 10 grupos 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 la estadística de brecha de Edwin Chen: ingrese la descripción de la imagen aquí

Siete . También puede resultarle útil explorar sus datos con clustergramas 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 la cantidad de clústeres en un conjunto de datos.

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!

ingrese la descripción de la imagen aquí

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 vea aquí para obtener métodos más exóticos: 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

ingrese la descripción de la imagen aquí

# 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

ingrese la descripción de la imagen aquí

También para datos de alta dimensión está la pvclustbiblioteca que calcula los valores p para la agrupación jerárquica a través del remuestreo bootstrap multiescala. Aquí está el ejemplo de la documentación (no funcionará en datos de dimensiones tan bajas como en mi ejemplo):

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

ingrese la descripción de la imagen aquí

¿Algo de eso ayuda?

Ben
fuente
Para el último dendograma (clúster dendograma con AU / BP) a veces es conveniente dibujar rectángulos alrededor de los grupos con valores p relativamente altos: pvrect (fit, alpha = 0.95)
Igor Elbert
Esto es exactamente lo que estaba buscando. Soy nuevo en R y me habría llevado mucho tiempo encontrar esto. Gracias @Ben por responder con tanto detalle. ¿Puede guiarme sobre dónde puedo encontrar la lógica detrás de cada uno de estos métodos, como qué métrica o criterio están utilizando para determinar el número óptimo de grupos, o cómo es cada uno de ellos diferente el uno del otro? Mi jefe quiere que le diga eso, para que podamos decidir cuál de los métodos usar. Gracias por adelantado.
nasia jaffri
1
@ Aleksandr Blekh También podría intentar convertir cualquier método gráfico en analítico. Por ejemplo, uso el método "codo" (mencionado por primera vez en la respuesta), pero trato de encontrarlo analíticamente. El punto del codo podría ser el punto con curvatura máxima. Para datos discretos, es el punto con la diferencia central máxima de segundo orden (análogo a la derivada máxima de segundo orden para datos continuos). Consulte stackoverflow.com/a/4473065/1075993 y stackoverflow.com/q/2018178/1075993 . Supongo que otros métodos gráficos también podrían convertirse a analíticos.
Andrey Sapegin
1
@AndreySapegin: podría, pero: 1) francamente, no lo considero una solución elegante (en mi humilde opinión, en la mayoría de los casos, los métodos visuales deben permanecer visuales, mientras que los analíticos deben seguir siendo analíticos); 2) He descubierto una solución analítica para esto, usando uno o varios Rpaquetes (está en mi GitHub, puedes echar un vistazo); 3) mi solución parece funcionar lo suficientemente bien, además, ha pasado un tiempo y ya he finalizado mi software de tesis, informe de tesis (tesis) y actualmente me estoy preparando para la defensa :-). De todos modos, agradezco mucho su comentario y enlaces. ¡Todo lo mejor!
Aleksandr Blekh
1
2,2 millones de filas están en mi conjunto de datos de agrupamiento actual. Ninguno de estos paquetes R funciona, espero. Simplemente abren mi computadora y luego se cae por mi experiencia. Sin embargo, parece que el autor conoce sus cosas para datos pequeños y para el caso general sin tener en cuenta la capacidad del software. No se deducen puntos debido al obvio buen trabajo del autor. Ustedes, por favor, sepan que el viejo R es horrible con 2.2 millones de filas, pruébelo usted mismo si no confía en mí. H2O ayuda pero se limita a un pequeño jardín amurallado de felicidad.
Geoffrey Anderson
21

Es difícil agregar algo también una respuesta tan elaborada. Aunque creo que deberíamos mencionarlo identifyaquí, particularmente porque @Ben muestra muchos ejemplos de dendrogramas.

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

identifyle permite elegir interactivamente grupos de un dendrograma y almacena sus elecciones en una lista. Pulsa Esc para salir del modo interactivo y volver a la consola R. Tenga en cuenta que la lista contiene los índices, no los nombres de fila (en oposición a cutree).

Matt Bannert
fuente
10

Para determinar el k-cluster óptimo en los métodos de agrupamiento. Usualmente uso el Elbowmétodo acompañado por el procesamiento en paralelo para evitar el consumo de tiempo. Este código puede muestrear así:

Método del 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.

VanThaoNguyen
fuente
2
Las funciones de codo y css provienen del paquete GMD: cran.r-project.org/web/packages/GMD/GMD.pdf
Rohan
6

Espléndida respuesta de Ben. Sin embargo, me sorprende que el método de Propagación de afinidad (AP) se haya 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 agrupando los datos. Consulte el documento científico que respalda 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.

Entonces, si no está predispuesto hacia k-means, sugiero usar AP directamente, que agrupará los datos sin requerir conocer el número de grupos:

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

Si las distancias euclidianas negativas no son apropiadas, 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 similitudes en el paquete AP solo se proporcionan 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 desee agrupar en su matriz (filas o columnas).

zsram
fuente
6

Estos métodos son geniales, pero cuando se trata 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á el número óptimo de grupos para usted.

Primero querrá asegurarse de que Weka esté instalado en su sistema y que XMeans esté instalado a través de la herramienta de administració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
RDRR
fuente
6

Una solución simple es la biblioteca factoextra. Puede cambiar el método de agrupación y el método para calcular la mejor cantidad de grupos. Por ejemplo, si desea conocer el mejor número de clústeres para un k- significa:

Datos: mtcars

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

Finalmente, obtenemos un gráfico como:

ingrese la descripción de la imagen aquí

Cro-Magnon
fuente
2

Las respuestas son geniales. Si desea dar una oportunidad a otro método de agrupación, puede utilizar 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)

ingrese la descripción de la imagen aquí

Dependiendo de cuántas clases necesite, puede cortar su dendrograma 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).

boyaronur
fuente
Prefiero Wards sobre completo.
Chris