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

作者: 只怕再见是故人
发布时间: 2024-12-28 06:24:49 (1月前)
转自：

7 条回复

0#
回复此人
夏花 | 2019-08-31 10-32

<div class =“post-text”itemprop =“text”> 这些方法很棒但是当试图为更大的数据集找到k时，这些在R中可能会很慢。 我发现的一个很好的解决方案是“RWeka”软件包，它具有X-Means算法的有效实现 - K-Means的扩展版本可以更好地扩展，并将为您确定最佳的簇数。 首先，您需要确保在您的系统上安装了Weka，并通过Weka的软件包管理器工具安装了XMeans。 <pre> <code> 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 </code> </pre> </DIV>

编辑
1#
回复此人
生如夏花 | 2019-08-31 10-32

<div class =“post-text”itemprop =“text”> 为了确定聚类方法中的最优k-聚类。我经常使用 <code> Elbow </code> 方法伴随并行处理以避免时间消耗。此代码可以像这样样本： 的 肘法 </强> <pre> <code> 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) } </code> </pre> 的 跑肘平行 </强> <pre> <code> 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) </code> </pre> 它运作良好。 </DIV>

编辑
2#
回复此人
部落用户 | 2019-08-31 10-32

<div class =“post-text”itemprop =“text”> 如果你的问题是 <code> how can I determine how many clusters are appropriate for a kmeans analysis of my data? </code> ，那么这里有一些选择。该 <a href="http://en.wikipedia.org/wiki/Determining_the_number_of_clusters_in_a_data_set"> 维基百科文章 </A> 在确定群集数量时，对这些方法中的一些方法进行了很好的回顾。 首先，一些可重现的数据（Q中的数据......我不清楚）： <pre> <code> 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) </code> </pre> <img src =“https://i.stack.imgur.com/vlioB.png”alt =“在此处输入图片说明”/> 的 一 </强> 。在平方误差（SSE）碎石图中寻找弯曲或弯头。看到 <a href="http://www.statmethods.net/advstats/cluster.html"> http://www.statmethods.net/advstats/cluster.html </A> ＆安培; <a href="http://www.mattpeeples.net/kmeans.html"> http://www.mattpeeples.net/kmeans.html </A> 更多。在结果图中肘部的位置表明适合kmeans的簇的数量： <pre> <code> 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") </code> </pre> 我们可以得出结论，这种方法将表明4个集群： <img src =“https://i.stack.imgur.com/VVzNR.png”alt =“在此处输入图片说明”/> 的 二 </强> 。您可以围绕medoids进行分区，以估计使用的群集数量 <code> pamk </code> 功能在fpc包中。 <pre> <code> 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)) </code> </pre> <img src =“https://i.stack.imgur.com/lQrhi.png”alt =“在此处输入图片说明”/> <img src =“https://i.stack.imgur.com/qF2jI.png”alt =“在此处输入图片说明”/> <pre> <code> # 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 </code> </pre> 的 三 </强> 。 Calinsky准则：另一种诊断适合数据的簇数的方法。在这种情况下，我们尝试1到10组。 <pre> <code> 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! </code> </pre> <img src =“https://i.stack.imgur.com/yEbi3.png”alt =“在此处输入图片说明”/> 的 四 </强> 。根据期望最大化的贝叶斯信息准则确定最优模型和聚类数，通过参数化高斯混合模型的层次聚类初始化 <pre> <code> # 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) </code> </pre> <img src =“https://i.stack.imgur.com/ZdR5C.png”alt =“在此处输入图片说明”/> <img src =“https://i.stack.imgur.com/hDdKg.png”alt =“在此处输入图片说明”/> <img src =“https://i.stack.imgur.com/3zqaN.png”alt =“在此处输入图片说明”/> 的 五 </强> 。亲和传播（AP）聚类，参见 <a href="http://dx.doi.org/10.1126/science.1136800"> http://dx.doi.org/10.1126/science.1136800 </A> <pre> <code> 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) </code> </pre> <img src =“https://i.stack.imgur.com/XIZFU.png”alt =“在此处输入图片说明”/> <img src =“https://i.stack.imgur.com/ajx39.png”alt =“在此处输入图片说明”/> 的 六 </强> 。用于估计群集数量的差距统计。也可以看看 <a href="https://github.com/echen/gap-statistic"> 一些漂亮的图形输出的代码 </A> 。在这里尝试2-10个集群： <pre> <code> 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 </code> </pre> 这是Edwin Chen实施差距统计的输出： <img src =“https://i.stack.imgur.com/2kQHH.png”alt =“在此处输入图片说明”/> 的 七 </强> 。您还可以发现使用集群图浏览数据以查看集群分配是否有用，请参阅 <a href="http://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r-code/"> http://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r-code/ </A> 更多细节。 的 八 </强> 。该 <a href="http://cran.r-project.org/web/packages/NbClust/index.html"> NbClust包 </A> 提供30个索引来确定数据集中的簇数。 <pre> <code> 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! </code> </pre> <img src =“https://i.stack.imgur.com/ptSXf.png”alt =“在此处输入图片说明”/> 如果你的问题是 <code> how can I produce a dendrogram to visualize the results of my cluster analysis </code> 那么你应该从这些开始： <a href="http://www.statmethods.net/advstats/cluster.html"> http://www.statmethods.net/advstats/cluster.html </A> <a href="http://www.r-tutor.com/gpu-computing/clustering/hierarchical-cluster-analysis"> http://www.r-tutor.com/gpu-computing/clustering/hierarchical-cluster-analysis </A> <a href="http://gastonsanchez.wordpress.com/2012/10/03/7-ways-to-plot-dendrograms-in-r/"> http://gastonsanchez.wordpress.com/2012/10/03/7-ways-to-plot-dendrograms-in-r/ </A> 并在此处了解更多信息奇异方法： <a href="http://cran.r-project.org/web/views/Cluster.html"> http://cran.r-project.org/web/views/Cluster.html </A> 这里有一些例子： <pre> <code> d_dist <- dist(as.matrix(d)) # find distance matrix plot(hclust(d_dist)) # apply hirarchical clustering and plot </code> </pre> <img src =“https://i.stack.imgur.com/olsxB.png”alt =“在此处输入图片说明”/> <pre> <code> # 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 </code> </pre> <img src =“https://i.stack.imgur.com/4Iyh4.png”alt =“在此处输入图片说明”/> 对于高维数据也是如此 <code> pvclust </code> 通过多尺度引导程序重采样计算层次聚类的p值的库。这是文档中的示例（不会像我的示例那样处理如此低维数据）： <pre> <code> library(pvclust) library(MASS) data(Boston) boston.pv <- pvclust(Boston) plot(boston.pv) </code> </pre> <img src =“https://i.stack.imgur.com/OcVa3.png”alt =“在此处输入图片说明”/> 这有什么帮助吗？ </DIV>

编辑
3#
回复此人
妖邪 | 2019-08-31 10-32

<div class =“post-text”itemprop =“text”> 本的精彩回答。然而，我很惊讶这里建议的亲和传播（AP）方法只是为了找到k-means方法的簇数，其中一般来说AP可以更好地聚类数据。请参阅Science中支持此方法的科学论文： 的 Frey，Brendan J.和Delbert Dueck。 “通过在数据点之间传递消息进行聚类。” science 315.5814（2007）：972-976。 </强> 因此，如果您不偏向k-means，我建议直接使用AP，这将聚集数据，而无需知道群集的数量： <pre> <code> library(apcluster) apclus = apcluster(negDistMat(r=2), data) show(apclus) </code> </pre> 如果负欧氏距离不合适，则可以使用同一包中提供的其他相似性度量。例如，对于基于Spearman相关性的相似性，这就是您所需要的： <pre> <code> sim = corSimMat(data, method="spearman") apclus = apcluster(s=sim) </code> </pre> 请注意，AP包中的相似功能仅为简单起见而提供。实际上，R中的apcluster（）函数将接受任何相关矩阵。使用corSimMat（）之前也可以这样做： <pre> <code> sim = cor(data, method="spearman") </code> </pre> 要么 <pre> <code> sim = cor(t(data), method="spearman") </code> </pre> 取决于您想要在矩阵（行或列）上聚类的内容。 </DIV>

编辑
4#
回复此人
v-star*위위 | 2019-08-31 10-32

<div class =“post-text”itemprop =“text”> 答案很棒。如果您想有机会使用其他聚类方法，可以使用层次聚类并查看数据如何拆分。 <pre> <code> > set.seed(2) > x=matrix(rnorm(50*2), ncol=2) > hc.complete = hclust(dist(x), method="complete") > plot(hc.complete) </code> </pre> <a href="https://i.stack.imgur.com/9tsAY.png" rel="nofollow noreferrer"> <img src =“https://i.stack.imgur.com/9tsAY.png”alt =“在此处输入图片说明”/> </A> 根据您需要的课程数量，您可以将树形图剪切为; <pre> <code> > 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 </code> </pre> 如果你输入 <code> ?cutree </code> 你会看到定义。如果您的数据集有三个类，那就简单了 <code> cutree(hc.complete, k = 3) </code> 。相当于 <code> cutree(hc.complete,k = 2) </code> 是 <code> cutree(hc.complete,h = 4.9) </code> 。 </DIV>

编辑
5#
回复此人
荧惑 | 2019-08-31 10-32

<div class =“post-text”itemprop =“text”> 一个简单的解决方案是图书馆 <code> factoextra </code> 。您可以更改聚类方法和计算最佳组数的方法。例如，如果您想知道k-的最佳簇数： <H3> 数据：mtcars </H3> <pre> <code> library(factoextra) fviz_nbclust(mtcars, kmeans, method = "wss") + geom_vline(xintercept = 3, linetype = 2)+ labs(subtitle = "Elbow method") </code> </pre> 最后，我们得到一个图形： <a href="https://i.stack.imgur.com/gm0Xz.jpg" rel="nofollow noreferrer"> <img src =“https://i.stack.imgur.com/gm0Xz.jpg”alt =“在此处输入图片说明”/> </A> </DIV>

编辑

登录后才能参与评论