确定聚类分析中的类别个数的方法
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了确定聚类分析中的类别个数的方法相关的知识，希望对你有一定的参考价值。
参考技术A 以k-means聚类为例，对于一个聚类个数k，首先利用k-means聚类将样本聚成k类，然后计算k类中各类内各点与类中心的距离加和W(ki)，进而计算k类的距离加和W(k)=sum(W(k1),…,W(ki),…,W(kk))；根据原始数据的特点产生B个均匀分布的参考数据集，对于每个数据集都计算W(sk)，计算B个数据集的平均E.W(k)=mean(W(1k),…,W(sk),…,W(bk))；那么对于每个k就有：gap(k)=log(E.W(k))-log(W(k))；然后选取最小的k，使得gap(k)为局部最大值，并且超出了其邻居1个标准差，即gap(k)-gap(k+1)>0.25*sd(W(s(k+1)))
参考文献：Single-cell messenger RNA sequencing reveals rare intestinal celltypes. Nature 2015
这里需要注意的是，gap statistic适用于可以直接设定聚类个数的聚类方法，如k均值和层次聚类，而类似密度聚类和PhenoGraph聚类方法，这两者是用参数(半径，resolution)来划分样本，无法直接设定类别个数。[谨慎!!!可以先调整参数使得得出的聚类个数和该方法得出的结果相同!!!]
轮廓系数，是聚类效果好坏的一种评价方式。最早由 Peter J. Rousseeuw 在 1986 提出。它结合内聚度和分离度两种因素。可以用来在相同原始数据的基础上评价不同算法、或者算法不同运行方式对聚类结果所产生的影响。
1. 计算样本i到同簇其他样本的平均距离ai。ai 越小，说明样本i越应该被聚类到该簇。将ai 称为样本i的簇内不相似度。簇C中所有样本的a i 均值称为簇C的簇不相似度。
2. 计算样本i到其他某簇Cj 的所有样本的平均距离bij，称为样本i与簇Cj 的不相似度。定义为样本i的簇间不相似度：bi =minbi1, bi2, ..., bik。bi越大，说明样本i越不属于其他簇。
3. 根据样本i的簇内不相似度ai和簇间不相似度bi，定义样本i的轮廓系数：
4. 判断：
a) si接近1，则说明样本i聚类合理；
b) si接近-1，则说明样本i更应该分类到另外的簇；
c) 若si近似为0，则说明样本i在两个簇的边界上。
所有样本的si的均值称为聚类结果的轮廓系数，si越大，说明聚类效果越好。
参考自: CSDN
引用文献: Cluster analysis of novel isometric strength measures produces a valid and evidence-based classification structure for wheelchair track racing. (Br J Sports Med. 2018 11.645)
mclust使用高斯混合模型对数据进行聚类分析。具体算法过于复杂，这里不详细解释。mclust包方法有点“暴力”，聚类数目可以自定义，比如选取从1到20，然后一共有14种模型，每一种模型都计算聚类数目从1到20的BIC值，最终确定最佳聚类数目。该方法弊端在于时间消耗特别高。使用的数据为R自带数据集wine。
这些分别代表着相关性(完全正负相关——对角线、稍强正负相关——椭圆、无关——圆)等参数的改变对应的模型。简单的看，BIC值越大则说明所选取的变量集合拟合效果越好。上图中除了两个模型一直递增，其他的12模型数基本上都是在聚类数目为3的时候达到峰值，所以该算法由此得出最佳聚类数目为3的结论。
这里的BIC的定义为：
和贝叶斯信息准则不是一回事!!!
相关文献: mclust version 4 for R: normal mixture modeling for model-based clustering, classification, and density estimation. (2012,引用次数877)
类似mclust，也是自己定义了几十种评估指标，然后遍历每一个设定的聚类数目，然后通过这些指标看分别在聚类数为多少时达到最优，最后选择指标支持数最多的聚类数目就是最佳聚类数目。
大体过程是用某种已有的聚类算法或者是划分类别的方法Kmeans，Ward(最小化类内方差)，Single(最小距离)，Complete(最大距离)和Average(平均距离)等，对每一个设定的类别个数进行聚类，得出聚类结果后用评估指标评估。
评价指标有:
"kl", "ch", "hartigan", "ccc", "scott", "marriot", "trcovw", "tracew", "friedman", "rubin", "cindex", "db", "silhouette"(轮廓系数), "duda", "pseudot2", "beale", "ratkowsky", "ball", "ptbiserial", "gap", "frey", "mcclain", "gamma", "gplus", "tau", "dunn", "hubert", "sdindex", "dindex", "sdbw"
"all" (all indices except GAP, Gamma, Gplus and Tau), "alllong" (all indices with Gap, Gamma, Gplus and Tau included).
这里不详细解释每一种评估指标。
相关文献: NbClust Package: finding the relevant number of clusters in a dataset. (2012,引用次数 25); Package 'nbclust'. (2014,引用次数 676)
一个简单可行的指标，SSE(sum of squared error组内平方误差)
library(factoextra)
library(ggplot2)
set.seed(1234)
fviz_nbclust(scale(wine[,-1]), kmeans, method = "wss") +
geom_vline(xintercept = 3, linetype = 2)
选用簇中位置最中心的对象，试图对n个对象给出k个划分；代表对象也被称为是中心点，其他对象则被称为非代表对象；最初随机选择k个对象作为中心点，该算法反复地用非代表对象来代替代表对象，试图找出更好的中心点，以改进聚类的质量；在每次迭代中，所有可能的对象对被分析，每个对中的一个对象是中心点，而另一个是非代表对象。对可能的各种组合，估算聚类结果的质量；一个对象Oi可以被使最大平方-误差值减少的对象代替；在一次迭代中产生的最佳对象集合成为下次迭代的中心点。
对比kmeans：k-means是每次选簇的均值作为新的中心，迭代直到簇中对象分布不再变化。其缺点是对于离群点是敏感的，因为一个具有很大极端值的对象会扭曲数据分布。而PAM考虑新的簇中心不选择均值而是选择簇内的某个对象，只要使总的代价降低就可以。kmedoids算法比kmenas对于噪声和孤立点更鲁棒，因为它最小化相异点对的和（minimizes a sum of pairwise dissimilarities ）而不是欧式距离的平方和（sum of squared Euclidean distances.）。一个中心点（medoid）可以这么定义：簇中某点的平均差异性在这一簇中所有点中最小。
R包fpc中的pamk函数可以确定最佳聚类个数，原理是通过尝试每一个设定的聚类数，用评价标准来看哪一个聚类数最佳，评价标准有:
1."asw" : average silhouette(平均轮廓系数) width given out by pam/clara
2."multiasw": average silhouette width computed by distcritmulti
3."ch": Calinski-Harabasz
Calinski-Harabasz公式:
其中，n表示聚类的数目 ,k 表示当前的类, trB(k)表示类间离差矩阵的迹, trW(k) 表示类内离差矩阵的迹。
相关文献: A new and efficient k-medoid algorithm for spatial clustering. (2005 引用次数: 87)
引用文献: A simple and fast algorithm for K-medoids clustering. (2009 引用数: 993)
Calinski-Harabasz准则有时称为方差比准则 (VRC)，它可以用来确定聚类的最佳K值。
Calinski Harabasz 指数定义为：
其中，K是聚类数，N是样本数，SSB是组与组之间的平方和误差，SSw是组内平方和误差。因此，如果SSw越小、SSB越大，那么聚类效果就会越好，即Calinsky criterion值越大，聚类效果越好。
相关文献: A dendrite method for cluster analysis. (1974 引用数4447)
相关文献: Clustering using flower pollination algorithm and Calinski-Harabasz index. (2016 引用数19)
这个本质上是类似kmeans或者层次聚类一样，是一种聚类方法。
AP算法的基本思想是将全部样本看作网络的节点，然后通过网络中各条边的消息传递计算出各样本的聚类中心。聚类过程中，共有两种消息在各节点间传递，分别是吸引度(responsibility)和归属度(availability)。AP算法通过迭代过程不断更新每一个点的吸引度和归属度值，直到产生m个高质量的Exemplar（类似于质心），同时将其余的数据点分配到相应的聚类中。
相关文献: APCluster: an R package for affinity propagation clustering. (2011 引用数240)
以上内容撰写时间较早，如有参考其他文章而未标注引用的部分，实在抱歉，请及时告知，我立即修改。
R中的聚类分析：确定最佳聚类数
【中文标题】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。以上是关于确定聚类分析中的类别个数的方法的主要内容，如果未能解决你的问题，请参考以下文章