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

Posted 2023-03-12

【中文标题】R中的聚类分析：确定最佳聚类数

【英文标题】：Cluster analysis in R: determine the optimal number of clusters

【发布时间】：2013-02-28 20:13:10

【问题描述】：

作为 R 的新手，我不太确定如何选择最佳数量的集群来进行 k-means 分析。在绘制以下数据的子集后，有多少个集群是合适的？如何进行聚类树状分析？

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 )

【问题讨论】：

如果您还没有完全依赖 kmeans，您可以尝试 DBSCAN 聚类算法，该算法在 fpc 包中提供。是的，你必须设置两个参数......但我发现fpc::dbscan 在自动确定大量集群方面做得很好。此外，如果数据告诉您，它实际上可以输出单个集群 - @Ben 的优秀答案中的一些方法无法帮助您确定 k=1 是否实际上是最好的。

另见stats.stackexchange.com/q/11691/478

【参考方案1】：

如果您的问题是how can I determine how many clusters are appropriate for a kmeans analysis of my data?，那么这里有一些选项。 wikipedia article 关于确定集群数量对其中一些方法进行了很好的回顾。

首先，一些可重现的数据（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 个集群：

两个。您可以使用 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 标准：另一种诊断有多少集群适合数据的方法。在这种情况下，我们尝试 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.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)

六个。用于估计聚类数量的差距统计。另见some code for a nice graphical output。在这里尝试 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 实施差距统计的输出：

七。您可能还会发现使用 clustergrams 探索数据以可视化集群分配很有用，请参阅http://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r-code/ 了解更多详细信息。

八。 NbClust package 提供 30 个索引来确定数据集中的集群数量。

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!

如果您的问题是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 值。这是文档中的示例（不会像我的示例那样处理低维数据）：

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

这些有帮助吗？

【讨论】：

对于最后一个树状图（带有 AU/BP 的集群树状图），有时在 p 值相对较高的组周围绘制矩形很方便：pvrect(fit, alpha=0.95)

这正是我想要的。我是 R 新手，我需要很长时间才能找到它。感谢@Ben 如此详细地回答。您能否指导我在哪里可以找到这些方法背后的逻辑，例如它们使用什么指标或标准来确定最佳集群数量，或者它们中的每一个之间有何不同。我的老板想让我说出来，这样我们就可以决定使用哪一种方法。提前致谢。

@Aleksandr Blekh 您也可以尝试将任何图形方法转换为分析方法。例如，我使用“肘部”方法（在答案中首先提到），但尝试分析地找到它。肘点可以是具有最大曲率的点。对于离散数据，它是具有最大二阶中心差的点（类似于连续数据的最大二阶导数）。请参阅 ***.com/a/4473065/1075993 和 ***.com/q/2018178/1075993。我想其他图形方法也可以转换为分析方法。

@AndreySapegin：我可以，但是：1）坦率地说，我不认为这是一个优雅的解决方案（恕我直言，在大多数情况下，视觉方法应该保持视觉，而分析方法应该保持分析）； 2）我已经找到了分析解决方案，使用一个或几个R包（它在我的GitHub上-欢迎你看看）； 3）我的解决方案似乎运行得很好，另外，已经有一段时间了，我已经完成了我的论文软件、论文报告（论文），目前我正在为答辩做准备:-)。无论如何，我非常感谢您的评论和链接。万事如意！

220 万行在我当前的聚类数据集中。我希望这些 R 包都不能用于它。他们只是弹出我的电脑，然后它就从我的经验中消失了。然而，看起来作者知道他的小数据和一般情况下的知识，而不考虑软件容量。不因作者明显的好作品而扣分。你们都知道普通的旧 R 在 220 万行时很糟糕——如果你不相信我，你自己试试。 H2O 有帮助，但仅限于快乐的小围墙花园。

【参考方案2】：

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

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

identify 让您以交互方式从树状图中选择集群并将您的选择存储到列表中。按 Esc 退出交互模式并返回 R 控制台。请注意，该列表包含索引，而不是行名（与 cutree 相反）。

【参考方案3】：

为了确定聚类方法中的最优 k-cluster。我通常使用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)

效果很好。

【讨论】：

肘部和 css 函数来自 GMD 包：cran.r-project.org/web/packages/GMD/GMD.pdf

GMD 不再适用于最新版本的 R，是否有替代品？

【参考方案4】：

一个简单的解决方案是库factoextra。您可以更改聚类方法和计算最佳组数的方法。例如，如果您想知道 k 均值的最佳聚类数：

数据：mtcars

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

最后，我们得到如下图：

【参考方案5】：

本的精彩回答。然而，令我感到惊讶的是，这里建议使用亲和传播 (AP) 方法只是为了找到 k-means 方法的聚类数，通常 AP 在聚类数据方面做得更好。请在此处查看支持此方法的科学论文：

弗雷、布伦丹 J. 和德尔伯特杜克。 “通过在数据点之间传递消息进行聚类。”科学 315.5814 (2007): 972-976。

因此，如果您不偏向于 k-means，我建议直接使用 AP，它可以对数据进行聚类，而无需知道聚类的数量：

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

如果负欧几里得距离不合适，那么您可以使用同一包中提供的其他相似性度量。例如，对于基于 Spearman 相关性的相似性，这就是您所需要的：

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

请注意，AP 包中的那些相似功能只是为了简单起见而提供的。事实上，R 中的 apcluster() 函数将接受任何相关矩阵。与之前的 corSimMat() 相同，可以这样做：

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

或

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

取决于您要在矩阵上聚类的内容（行或列）。

【参考方案6】：

这些方法很棒，但是当尝试为更大的数据集找到 k 时，这些方法在 R 中可能会非常缓慢。

我发现的一个很好的解决方案是“RWeka”包，它有效地实现了 X-Means 算法 - K-Means 的扩展版本，可以更好地扩展，并为您确定最佳的集群数量。

首先，您需要确保在您的系统上安装了 Weka，并通过 Weka 的包管理器工具安装了 XMeans。

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

【参考方案7】：

答案很棒。如果您想给其他聚类方法一个机会，您可以使用层次聚类并查看数据是如何拆分的。

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

根据您需要多少类，您可以将您的树状图切割为；

> 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,h = 4.9)。

【讨论】：

比起完整版，我更喜欢 Wards。