在聚类分析的时候确定最佳聚类数目是一个很重要的问题，比如kmeans函数就要你提供聚类数目这个参数，总不能两眼一抹黑乱填一个吧。之前也被这个问题困扰过，看了很多博客，大多泛泛带过。今天把看到的这么多方法进行汇总以及代码实现并尽量弄清每个方法的原理。

数据集选用比较出名的wine数据集进行分析

1. get data

(1) load data

#install.packages("gclus")
library(gclus) #聚类画图
data(wine)

dim(wine)
## [1] 178  14
head(wine)
##   Class Alcohol Malic  Ash Alcalinity Magnesium Phenols Flavanoids Nonflavanoid
## 1     1   14.23  1.71 2.43       15.6       127    2.80       3.06         0.28
## 2     1   13.20  1.78 2.14       11.2       100    2.65       2.76         0.26
## 3     1   13.16  2.36 2.67       18.6       101    2.80       3.24         0.30
## 4     1   14.37  1.95 2.50       16.8       113    3.85       3.49         0.24
## 5     1   13.24  2.59 2.87       21.0       118    2.80       2.69         0.39
## 6     1   14.20  1.76 2.45       15.2       112    3.27       3.39         0.34
##   Proanthocyanins Intensity  Hue OD280 Proline
## 1            2.29      5.64 1.04  3.92    1065
## 2            1.28      4.38 1.05  3.40    1050
## 3            2.81      5.68 1.03  3.17    1185
## 4            2.18      7.80 0.86  3.45    1480
## 5            1.82      4.32 1.04  2.93     735
## 6            1.97      6.75 1.05  2.85    1450
#View(wine)
# 共 14 列
colnames(wine) |> jsonlite::toJSON()
## ["Class","Alcohol","Malic","Ash","Alcalinity","Magnesium","Phenols","Flavanoids","Nonflavanoid","Proanthocyanins","Intensity","Hue","OD280","Proline"]
#"Class", 分类：1，2，3
#"Alcohol", 酒精含量
#"Malic", 苹果酸
#"Ash", 灰烬
#"Alcalinity", 灰的碱度
#"Magnesium", 镁
#"Phenols", 总酚
#"Flavanoids", 黄酮
#"Nonflavanoid", 非类黄酮
#"Proanthocyanins", 锦葵原花青素
#"Intensity", 色彩强度
#"Hue", 色相
#"OD280", 吸光度 OD280/OD315稀释葡萄酒
#"Proline" 脯氨酸

(2) pre-processing

# remove label
dat=wine[, -1]
dim(dat) #13 cols
## [1] 178  13
# scale to 100 of each column
dat.scale = apply(dat, 2, function(x){
  (x-mean(x)) / sd(x)
})
head(dat.scale)
##     Alcohol       Malic        Ash Alcalinity  Magnesium   Phenols Flavanoids
## 1 1.5143408 -0.56066822  0.2313998 -1.1663032 1.90852151 0.8067217  1.0319081
## 2 0.2455968 -0.49800856 -0.8256672 -2.4838405 0.01809398 0.5670481  0.7315653
## 3 0.1963252  0.02117152  1.1062139 -0.2679823 0.08810981 0.8067217  1.2121137
## 4 1.6867914 -0.34583508  0.4865539 -0.8069748 0.92829983 2.4844372  1.4623994
## 5 0.2948684  0.22705328  1.8352256  0.4506745 1.27837900 0.8067217  0.6614853
## 6 1.4773871 -0.51591132  0.3043010 -1.2860793 0.85828399 1.5576991  1.3622851
##   Nonflavanoid Proanthocyanins  Intensity        Hue     OD280     Proline
## 1   -0.6577078       1.2214385  0.2510088  0.3610679 1.8427215  1.01015939
## 2   -0.8184106      -0.5431887 -0.2924962  0.4048188 1.1103172  0.96252635
## 3   -0.4970050       2.1299594  0.2682629  0.3173170 0.7863692  1.39122370
## 4   -0.9791134       1.0292513  1.1827317 -0.4264485 1.1807407  2.32800680
## 5    0.2261576       0.4002753 -0.3183774  0.3610679 0.4483365 -0.03776747
## 6   -0.1755994       0.6623487  0.7298108  0.4048188 0.3356589  2.23274072
#head(scale(dat))
table((scale(dat) - dat.scale)<1e-5)
## 
## TRUE 
## 2314
dat=dat.scale

2. mclust

mclust包是聚类分析非常强大的一个包，也是上课时老师给我们介绍的一个包，每次导入时有一种科技感 :) 帮助文档非常详尽，可以进行聚类、分类、密度分析

Mclust包方法有点“暴力”，聚类数目自定义，比如我选取的从1到20，然后一共14种模型，每一种模型都计算聚类数目从1到20的BIC值，最终确定最佳聚类数目，这种方法的思想很直接了当，但是弊端也就显然易见了——时间复杂度太高，效率低。

# install.packages("mclust")
library(mclust) #聚类Clustering、分类Classification、密度分析Density estimation
m_clust <- Mclust(as.matrix(dat), G=1:20) #聚类数目从1一直试到20
summary(m_clust)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVE (ellipsoidal, equal orientation) model with 3 components: 
## 
##  log-likelihood   n  df       BIC       ICL
##       -2292.525 178 158 -5403.772 -5404.735
## 
## Clustering table:
##  1  2  3 
## 56 73 49
par(mar=c(4,4,1,1))
plot(m_clust, "BIC")

# 作者自己定义的BIC，值越大越好。并不是熟知的 贝叶斯信息准则

结论：除了2个指标外，其余都是n=3达到峰值。

3. Nbclust

Nbclust包是我在《R语言实战》上看到的一个包，思想和mclust包比较相近，也是定义了几十个评估指标，然后聚类数目从2遍历到15（自己设定），然后通过这些指标看分别在聚类数为多少时达到最优，最后选择指标支持数最多的聚类数目就是最佳聚类数目。

# install.packages("NbClust")
library(NbClust)
set.seed(2023) #因为method选择的是kmeans，所以如果不设定种子，每次跑得结果可能不同

par(mar=c(4,4,1,1))
nb_clust <- NbClust(dat, distance = "euclidean",
                    min.nc=2, max.nc=15, method = "kmeans",
                    index = "alllong", alphaBeale = 0.1)

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 3 proposed 2 as the best number of clusters 
## * 20 proposed 3 as the best number of clusters 
## * 1 proposed 13 as the best number of clusters 
## * 1 proposed 14 as the best number of clusters 
## * 2 proposed 15 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  3 
##  
##  
## *******************************************************************
#xlab="聚类数", ylab = "支持指标数"
par(mar=c(4,4,1,1))
barplot(table(nb_clust$Best.nc[1,]),xlab = "Cluster Number", ylab = "Supporting Index Number")

可以看到有20个指标支持最佳聚类数目为3，3个指标支持聚类数为2，所以该方法推荐的最佳聚类数目为3.

结论：分3类最优，支持的指标最多

4. [elbow]组内平方误差和——拐点图

想必之前动辄几十个指标，这里就用一个最简单的指标——sum of squared error (SSE)组内平方误差和来确定最佳聚类数目。这个方法也是出于《R语言实战》，自定义的一个求组内误差平方和的函数。

wssplot <- function(data, nc=15, seed=1234){
  # n=1时: 聚成一类的组内平方误差
  wss <- (nrow(data)-1)*sum(apply(data,2,var))
  # n=2:nc时
  for (i in 2:nc){
    set.seed(seed)
    wss[i] <- sum( kmeans(data, centers=i)$withinss )
  }
  plot(1:nc, wss, type="b", xlab="Number of Clusters",
       ylab="Within groups sum of squares")
}

par(mar=c(4,4,1,1))
wssplot(dat)

从一类到三类下降得很快，之后下降得很慢，所以最佳聚类个数选为3。

随着聚类数目增多，每一个类别中数量越来越少，距离越来越近，因此WSS值肯定是随着聚类数目增多而减少的，所以关注的是斜率的变化，但WWS减少得很缓慢时，就认为进一步增大聚类数效果也并不能增强，存在得这个“肘点”就是最佳聚类数目，从一类到三类下降得很快，之后下降得很慢，所以最佳聚类个数选为三

另外也有现成的包(factoextra)可以调用

5. factoextra

# install.packages("factoextra")
library(factoextra)
library(ggplot2)
set.seed(1234)
fviz_nbclust(dat, kmeans, method = "wss", k.max = 20) +
  geom_vline(xintercept = 3, linetype = 2)

选定为3类为最佳聚类数目

用该包下的fviz_cluster函数可视化一下聚类结果

km.res <- kmeans(dat,3)
fviz_cluster(km.res, data = dat)

6. PAM(Partitioning Around Medoids) 围绕中心点的分割算法

k-means算法取得是均值，那么对于异常点其实对其的影响非常大，很可能这种孤立的点就聚为一类，一个改进的方法就是PAM算法，也叫k-medoids clustering

首先通过fpc包中的pamk函数得到最佳聚类数目

#install.packages("fpc")
library(fpc)
pamk.best <- pamk(dat)
pamk.best$nc
## [1] 3

pamk函数不需要提供聚类数目，也会直接自动计算出最佳聚类数，这里也得到为3

得到聚类数提供给cluster包下的pam函数并进行可视化

clust_2 = pam(dat, pamk.best$nc)
plot(clust_2)

#library(cluster)
#par(mar=c(4,4,1,1))
# clusplot( clust_2, main="" ) # 同上图 左
# we could also do:
library(fpc)
asw <- numeric(20)
for (k in 2:20){
  asw[[k]] <- pam(dat, k)$silinfo$avg.width
}
k.best <- which.max(asw)
cat("silhouette-optimal number of clusters:", k.best, "\n")
## silhouette-optimal number of clusters: 3

7. Calinsky criterion

这个评估标准定义[5]如下：

其中，k是聚类数，N是样本数，SSw是我们之前提到过的组内平方和误差，SSb是组与组之间的平方和误差，SSw越小，SSb越大聚类效果越好，所以Calinsky criterion值一般来说是越大，聚类效果越好

#install.packages("vegan")
library(vegan)
ca_clust <- cascadeKM(dat, 1, 10, iter = 1000)
ca_clust$results
##          1 groups   2 groups   3 groups  4 groups   5 groups   6 groups
## SSE          2301 1649.45652 1270.72887 1168.5920 1095.13335 1032.77270
## calinski       NA   69.52087   70.94253   56.2041   47.62318   42.24261
##           7 groups  8 groups  9 groups 10 groups
## SSE      970.80549 918.97337 874.91280  834.6779
## calinski  39.05061  36.52283  34.43325   32.7927
par(mar=c(4,4,1,1))
plot(ca_clust$results[1,], type="o")

可以看到该函数把组内平方和误差和Calinsky都计算出来了，可以看到calinski在聚类数为3时达到最大值。

calinski.best <- as.numeric(which.max(ca_clust$results[2,]))
calinski.best
## [1] 3
#画图出来观察一下
plot(ca_clust, sortg = TRUE, grpmts.plot = TRUE)

# 注意到那个红点就是对应的最大值，自带的绘图横轴纵轴取的可能不符合我们的直觉，把数据取出来自己单独画一下
calinski<-as.data.frame(ca_clust$results[2,])
calinski$cluster <- c(1:10)
library(ggplot2)
ggplot(calinski,aes(x = calinski[,2], y = calinski[,1]))+geom_line()+theme_bw()
## Warning: Removed 1 row containing missing values (`geom_line()`).

这个看上去直观多了。这就很清晰的可以看到在聚类数目为3时，calinski指标达到了最大值，所以最佳数目为3

8.Affinity propagation (AP) clustering

  • 这个本质上是类似kmeans或者层次聚类一样，是一种聚类方法，因为不需要像kmeans一样提供聚类数，会自动算出最佳聚类数，因此也放到这里作为一种计算最佳聚类数目的方法。
  • AP算法的基本思想是将全部样本看作网络的节点，然后通过网络中各条边的消息传递计算出各样本的聚类中心。聚类过程中，共有两种消息在各节点间传递，分别是吸引度( responsibility)和归属度(availability) 。
  • AP算法通过迭代过程不断更新每一个点的吸引度和归属度值，直到产生m个高质量的Exemplar（类似于质心），同时将其余的数据点分配到相应的聚类中[7]
#install.packages("apcluster")
library(apcluster)
ap_clust <- apcluster(negDistMat(r=2), dat)
length(ap_clust@clusters) # 15
## [1] 15
#该聚类方法推荐的最佳聚类数目为15，再用热力图可视化一下
heatmap(ap_clust)

选x或者y方向看(对称)，可以数出来“叶子节点”一共15个

9. 轮廓系数 Average silhouette method

轮廓系数是类的密集与分散程度的评价指标。

s(i) = [b(i) - a(i)] / max(a(i), b(i))

a(i)是测量组内的相似度,b(i)是测量组间的相似度，s(i)范围从-1到1， 值越大说明组内吻合越高，组间距离越远——也就是说，轮廓系数值越大，聚类效果越好[9]

require(cluster)
library(factoextra)
fviz_nbclust(dat, kmeans, method = "silhouette", k.max = 20)

可以看到也是在聚类数为3时轮廓系数达到了峰值，所以最佳聚类数为3

10. Gap Statistic

之前我们提到了WSSE组内平方和误差，该种方法是通过找“肘点”来找到最佳聚类数，肘点的选择并不是那么清晰，因此斯坦福大学的Robert等教授提出了Gap Statistic方法，定义的Gap值为[9]

Gapn(k) = En( log(Wk)) - logWk

取对数的原因是因为Wk的值可能很

大通过这个式子来找出Wk跌落最快的点，Gap最大值对应的k值就是最佳聚类数

library(cluster)
set.seed(123)
gap_clust <- clusGap(dat, kmeans, K.max=15, B = 500, verbose = interactive())
gap_clust
## Clustering Gap statistic ["clusGap"] from call:
## clusGap(x = dat, FUNcluster = kmeans, K.max = 15, B = 500, verbose = interactive())
## B=500 simulated reference sets, k = 1..15; spaceH0="scaledPCA"
##  --> Number of clusters (method 'firstSEmax', SE.factor=1): 3
##           logW   E.logW       gap     SE.sim
##  [1,] 5.377557 5.863177 0.4856200 0.01198738
##  [2,] 5.203225 5.758357 0.5551325 0.01409569
##  [3,] 5.066921 5.696602 0.6296806 0.01219938
##  [4,] 5.026389 5.651454 0.6250652 0.01290231
##  [5,] 4.993046 5.615238 0.6221925 0.01197308
##  [6,] 4.961533 5.584118 0.6225850 0.01241987
##  [7,] 4.938117 5.556271 0.6181540 0.01185669
##  [8,] 4.921254 5.531437 0.6101830 0.01224035
##  [9,] 4.886470 5.508302 0.6218315 0.01218142
## [10,] 4.851747 5.487194 0.6354468 0.01162440
## [11,] 4.839963 5.466707 0.6267441 0.01145243
## [12,] 4.811363 5.447657 0.6362949 0.01150135
## [13,] 4.801165 5.429322 0.6281566 0.01125737
## [14,] 4.764846 5.412164 0.6473183 0.01165631
## [15,] 4.751258 5.394906 0.6436476 0.01122514
library(factoextra)
fviz_gap_stat(gap_clust)

可以看到也是在聚类数为3的时候gap值取到了最大值，所以最佳聚类数为3

11. 层次聚类 hclust

层次聚类是通过可视化然后人为去判断大致聚为几类，很明显在共同父节点的一颗子树可以被聚类为一个类

h_dist <- dist(as.matrix(dat), method = "euclidean")
h_clust<-hclust(h_dist, method="ward.D2")

par(mar=c(5,4,2,1))
plot(h_clust, hang = -1, labels = FALSE)
rect.hclust(h_clust,3)

12. clustergram

最后一种算法是Tal Galili[10]大牛自己定义的一种聚类可视化的展示，绘制随着聚类数目的增加，所有成员是如何分配到各个类别的。该代码没有被制作成R包，可以去Galili介绍页面)里面的 github地址 找到源代码跑一遍然后就可以用这个函数了，因为源代码有点长我就不放博客里面了，

# https://clustergram.readthedocs.io/en/stable/
# https://github.com/martinfleis/clustergram

# https://www.r-statistics.com/2010/06/clustergram-visualization-and-diagnostics-for-cluster-analysis-r-code/
# https://github.com/talgalili/R-code-snippets/blob/master/clustergram.r

#source("https://www.r-statistics.com/wp-content/uploads/2012/01/source_https.r.txt")
#source_https("https://raw.github.com/talgalili/R-code-snippets/master/clustergram.r")
par(cex.lab = 1.5, cex.main = 1.2)
set.seed(2023)
clustergram(as.matrix(dat), k.range = 2:8, line.width = 0.004)

随着K的增加，从最开始的两类到最后的八类，图肯定是越到后面越密集。

通过这个图判断最佳聚类数目的方法应该是看随着K每增加1，分出来的线越少说明在该k值下越稳定。比如k=7到k=8，假设k=7是很好的聚类数，那分成8类时应该可能只是某一类分成了两类，其他6类都没怎么变。反应到图中应该是有6簇平行线，有一簇分成了两股，而现在可以看到从7到8，线完全乱了，说明k=7时效果并不好。

按照这个分析，k=3到k=4时，第一股和第三股几本没变，就第二股拆成了2类，所以k=3是最佳聚类数目

方法汇总与比

wine数据集我们知道其实是分为3类的，以上10种判定方法中：

  • 层次聚类和clustergram方法、肘点图法，需要人工判定，虽然可以得出大致的最佳聚类数，但算法本身不会给出最佳聚类数
  • 除了Affinity propagation (AP) clustering 给出最佳聚类数为15，剩下6种全都是给出最佳聚类数为3

选用上次文本挖掘的矩阵进行分析(667*1623)

  • mclust效果很差，14种模型只有6种有结果
  • bclust报错
  • SSE可以运行
  • fpc包中的pamk函数聚成2类，明显不行
  • Calinsky criterion聚成2类
  • Affinity propagation (AP) clustering 聚成28类，相对靠谱
  • 轮廓系数Average silhouette聚类2类
  • gap-Statistic跑不出结果

可见上述方法中有的因为数据太大不能运行，有的结果很明显不对，一个可能是数据集的本身的原因（缺失值太多等），但是也告诉了我们在确定最佳聚类数目的时候需要多尝试几种方法，并没有固定的套路，然后选择一种可信度较高的聚类数目。

最后再把这10种方法总结一下：

方法 优点 缺点
mclust包 傻瓜式，强大 复杂度大
Nbclust包 傻瓜式，强大 复杂度大
WSSE 组内平方和误差肘点图 算法简单，复杂度小 不是很准
Calinsky criterion 比 WSSE 考虑更全面 复杂度高
轮廓系数Average silhouette method 比 WSSE考虑更全面
Gap Statistic 判定效果更好 复杂度高
PAM(Partitioning Around Medoids)围绕中心点的分割算法 对 kmeans 的改进
Affinity propagation (AP) clustering 热力图
clustergram 可视化 无量化指标判断
层次聚类 可视化 无量化指标判断

参考资料

  1. R语言实战第二版
  2. Partitioning cluster analysis: Quick start guide - Unsupervised Machine Learning
  3. BIC：http://www.stat.washington.edu/raftery/Research/PDF/fraley1998.pdf
  4. Cluster analysis in R: determine the optimal number of clusters
  5. Calinski-Harabasz Criterion：Calinski-Harabasz criterion clustering evaluation object
  6. Determining the optimal number of clusters: 3 must known methods - Unsupervised Machine Learning
  7. affinity-propagation：聚类算法Affinity Propagation(AP)
  8. 轮廓系数https://en.wikipedia.org/wiki/Silhouette(clustering))
  9. gap statistic-Tibshirani R, Walther G, Hastie T. Estimating the number of clusters in a data set via the gap statistic[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2001, 63(2): 411-423.
  10. ClustergramsClustergram: visualization and diagnostics for cluster analysis (R code)

I found the function pamk in fpc package to be most useful for my requirements.