R中的聚类分析:确定聚类的最佳数量

作为R的新手,我不太确定如何select最佳数量的聚类来进行k均值分析。 在绘制以下数据的子集之后,将有多less个集群适合? 我怎样才能进行聚类dendro分析?

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 ) 

如果你的问题是how can I determine how many clusters are appropriate for a kmeans analysis of my data? ,那么这里有一些选项。 关于确定群集数量的维基百科文章对这些方法中的一些进行了很好的回顾。

首先,一些可重复的数据(Q中的数据是不清楚的):

 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) 

在这里输入图像描述

一个 。 在平方误差(SSE)碎石图的总和中寻找弯曲或肘部。 有关更多信息 ,请参阅http://www.statmethods.net/advstats/cluster.html&http://www.mattpeeples.net/kmeans.html 。 肘部在结果图中的位置为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") 

我们可以得出结论:这个方法将表示4个簇: 在这里输入图像描述

。 您可以对medoids进行分区,以使用fpc软件包中的pamk函数来估计簇的数量。

 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 

。 Calinsky准则:诊断有多less个群集适合数据的另一种方法。 在这种情况下,我们尝试1到10组。

 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! 

在这里输入图像描述

。 根据用于期望最大化的贝叶斯信息准则确定最优模型和数量的聚类,通过用于参数化高斯混合模型的分层聚类来初始化

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

在这里输入图像描述在这里输入图像描述在这里输入图像描述

。 亲和传播(AP)聚类,请参阅http://dx.doi.org/10.1126/science.11​​36800

 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) 

在这里输入图像描述在这里输入图像描述

。 用于估计群集数量的差距统计。 另请参阅一些很好的graphics输出的代码 。 在这里尝试2-10簇:

 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 

下面是Edwin Chen实施差距统计的结果: 在这里输入图像描述

。 您也可能发现使用群集探索数据以可视化群集分配很有用,请参阅http://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r-代码/更多细节。;

。 NbClust软件包提供30个索引来确定数据集中的聚类数量。

 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! 

在这里输入图像描述

如果你的问题是how can I produce a dendrogram to visualize the results of my cluster analysis ,那么你应该从这些开始: 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/在这里看到更多的异国情调方法:; http : //cran.r-project.org/web/views/Cluster.html

这里有几个例子:

 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 

在这里输入图像描述

对于高维数据也是通过多尺度自pvclust来计算用于分层聚类的p值的pvclust库。 以下是文档中的示例(不会像在我的示例中那样处理这种低维数据):

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

在这里输入图像描述

有帮助吗?

很难添加一些如此详尽的答案。 虽然我觉得我们应该在这里identify ,特别是因为@Ben显示了很多树状图例子。

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

identify让您可以交互式地从树形图中select聚类,并将您的select存储到列表中。 按Esc离开交互模式并返回到R控制台。 请注意,该列表包含的索引,而不是cutree (而不是cutree )。

为了确定聚类方法中的最优k-聚类。 我通常使用Elbow方法伴随着并行处理来避免时间混淆。 这段代码可以像这样采样:

弯头法

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

平行运行肘部

 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) 

它运作良好。

来自本的精彩答案。 然而,我很惊讶的是这里提出的亲和传播(Affinity Propagation,AP)方法只是为了findk-均值方法的聚类数量,一般而言,AP对数据进行更好的聚类。 这里请看Science中支持这种方法的科学论文:

Frey,Brendan J.和Delbert Dueck。 “通过在数据点之间传递消息进行分组”。 科学315.5814(2007):972-976。

所以,如果你不偏向k-means,我build议直接使用AP,这将聚集数据,而不需要知道簇的数量:

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

如果负欧氏距离不合适,则可以使用同一包中提供的其他相似性度量。 例如,基于Spearman相关性的相似性,这就是你所需要的:

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

请注意,AP包中的相似function仅仅是为了简单起见。 实际上,R中的apcluster()函数将接受任何相关matrix。 以前用corSimMat()可以做到这一点:

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

要么

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

取决于你想在你的matrix(行或列)上聚类。

答案很好。 如果你想给另一个聚类方法的机会,你可以使用层次聚类,看看数据是如何分裂的。

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

在这里输入图像描述

根据您需要的class级数量,您可以将树状图减less为

 > 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 

如果您键入?cutree您将看到定义。 如果你的数据集有三个类,它将是简单的cutree(hc.complete, k = 3)cutree(hc.complete,k = 2)等价于cutree(hc.complete,k = 2) cutree(hc.complete,h = 4.9)