雷鋒網(wǎng) AI科技評(píng)論按，本文作者貝爾塔，原文載于知乎專欄數(shù)據(jù)分析與可視化，雷鋒網(wǎng) AI科技評(píng)論獲其授權(quán)發(fā)布。

在聚類分析的時(shí)候確定最佳聚類數(shù)目是一個(gè)很重要的問題，比如kmeans函數(shù)就要你提供聚類數(shù)目這個(gè)參數(shù)，總不能兩眼一抹黑亂填一個(gè)吧。之前也被這個(gè)問題困擾過，看了很多博客，大多泛泛帶過。今天把看到的這么多方法進(jìn)行匯總以及代碼實(shí)現(xiàn)并盡量弄清每個(gè)方法的原理。
數(shù)據(jù)集選用比較出名的wine數(shù)據(jù)集進(jìn)行分析

library(gclus)
data(wine)
head(wine)

Loading required package: cluster

因?yàn)槲覀円乙粋€(gè)數(shù)據(jù)集進(jìn)行聚類分析，所以不需要第一列的種類標(biāo)簽信息，因此去掉第一列。
同時(shí)注意到每一列的值差別很大，從1到100多都有，這樣會(huì)造成誤差，所以需要?dú)w一化，用scale函數(shù)

dataset <- wine[,-1] #去除分類標(biāo)簽
dataset <- scale(dataset)

去掉標(biāo)簽之后就可以開始對(duì)數(shù)據(jù)集進(jìn)行聚類分析了，下面就一一介紹各種確定最佳聚類數(shù)目的方法

判定方法

1.mclust包

mclust包是聚類分析非常強(qiáng)大的一個(gè)包，也是上課時(shí)老師給我們介紹的一個(gè)包，每次導(dǎo)入時(shí)有一種科技感 :) 幫助文檔非常詳盡，可以進(jìn)行聚類、分類、密度分析
Mclust包方法有點(diǎn)“暴力”，聚類數(shù)目自定義，比如我選取的從1到20，然后一共14種模型，每一種模型都計(jì)算聚類數(shù)目從1到20的BIC值，最終確定最佳聚類數(shù)目，這種方法的思想很直接了當(dāng)，但是弊端也就顯然易見了——時(shí)間復(fù)雜度太高，效率低。

library(mclust)
m_clust <- Mclust(as.matrix(dataset), G=1:20) #聚類數(shù)目從1一直試到20
summary(m_clust)

Gaussian finite mixture model fitted by EM algorithm


Mclust EVE (ellipsoidal, equal volume and orientation) model with 3 components:


log.likelihood   n  df       BIC       ICL
      -3032.45 178 156 -6873.257 -6873.549


Clustering table:
1  2  3

63 51 64
可見該函數(shù)已經(jīng)把數(shù)據(jù)集聚類為3種類型了。數(shù)目分別為63、51、64。再畫出14個(gè)指標(biāo)隨著聚類數(shù)目變化的走勢(shì)圖

plot(m_clust, "BIC")

下表是這些模型的意義

它們應(yīng)該分別代表著相關(guān)性(完全正負(fù)相關(guān)——對(duì)角線、稍強(qiáng)正負(fù)相關(guān)——橢圓、無關(guān)——圓)等參數(shù)的改變對(duì)應(yīng)的模型，研究清楚這些又是非常復(fù)雜的問題了，先按下表，知道BIC值越大則說明所選取的變量集合擬合效果越好。上圖中除了兩個(gè)模型一直遞增，其他的12模型數(shù)基本上都是在聚類數(shù)目為3的時(shí)候達(dá)到峰值，所以該算法由此得出最佳聚類數(shù)目為3的結(jié)論。
mclust包還可以用于分類、密度估計(jì)等，這個(gè)包值得好好把玩。

注意：此BIC并不是貝葉斯信息準(zhǔn)則！??！

最近上課老師講金融模型時(shí)提到了BIC值，說BIC值越小模型效果越好，頓時(shí)想起這里是在圖中BIC極大值為最佳聚類數(shù)目，然后和老師探討了這個(gè)問題，之前這里誤導(dǎo)大家了，Mclust包里面的BIC并不是貝葉斯信息準(zhǔn)則。

1.維基上的貝葉斯信息準(zhǔn)則定義

與log(likelihood)成反比，極大似然估計(jì)是值越大越好，那么BIC值確實(shí)是越小模型效果越好

2.Mclust包中的BIC定義[3]

這是Mclust包里面作者定義的“BIC值”，此BIC非彼BIC，這里是作者自己定義的BIC，可以看到，這里的BIC與極大似然估計(jì)是成正比的，所以這里是BIC值越大越好，與貝葉斯信息準(zhǔn)則值越小模型越好的結(jié)論并不沖突

2.Nbclust包

Nbclust包是我在《R語言實(shí)戰(zhàn)》上看到的一個(gè)包，思想和mclust包比較相近，也是定義了幾十個(gè)評(píng)估指標(biāo)，然后聚類數(shù)目從2遍歷到15（自己設(shè)定），然后通過這些指標(biāo)看分別在聚類數(shù)為多少時(shí)達(dá)到最優(yōu)，最后選擇指標(biāo)支持?jǐn)?shù)最多的聚類數(shù)目就是最佳聚類數(shù)目。

library(NbClust)
set.seed(1234) #因?yàn)閙ethod選擇的是kmeans，所以如果不設(shè)定種子，每次跑得結(jié)果可能不同
nb_clust <- NbClust(dataset,  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:                                                
* 5 proposed 2 as the best number of clusters
* 16 proposed 3 as the best number of clusters
* 1 proposed 10 as the best number of clusters
* 1 proposed 12 as the best number of clusters
* 1 proposed 14 as the best number of clusters
* 3 proposed 15 as the best number of clusters

                  ***** Conclusion *****                            

* According to the majority rule, the best number of clusters is  3

*******************************************************************

barplot(table(nb_clust$Best.nc[1,]),xlab = "聚類數(shù)",ylab = "支持指標(biāo)數(shù)")

可以看到有16個(gè)指標(biāo)支持最佳聚類數(shù)目為3，5個(gè)指標(biāo)支持聚類數(shù)為2，所以該方法推薦的最佳聚類數(shù)目為3.

3. 組內(nèi)平方誤差和——拐點(diǎn)圖

想必之前動(dòng)輒幾十個(gè)指標(biāo)，這里就用一個(gè)最簡單的指標(biāo)——sum of squared error (SSE)組內(nèi)平方誤差和來確定最佳聚類數(shù)目。這個(gè)方法也是出于《R語言實(shí)戰(zhàn)》，自定義的一個(gè)求組內(nèi)誤差平方和的函數(shù)。

wssplot <- function(data, nc=15, seed=1234){
   wss <- (nrow(data)-1)*sum(apply(data,2,var))
   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")}

wssplot(dataset)

隨著聚類數(shù)目增多，每一個(gè)類別中數(shù)量越來越少，距離越來越近，因此WSS值肯定是隨著聚類數(shù)目增多而減少的，所以關(guān)注的是斜率的變化，但WWS減少得很緩慢時(shí)，就認(rèn)為進(jìn)一步增大聚類數(shù)效果也并不能增強(qiáng)，存在得這個(gè)“肘點(diǎn)”就是最佳聚類數(shù)目，從一類到三類下降得很快，之后下降得很慢，所以最佳聚類個(gè)數(shù)選為三

另外也有現(xiàn)成的包(factoextra)可以調(diào)用

library(factoextra)
library(ggplot2)
set.seed(1234)
fviz_nbclust(dataset, kmeans, method = "wss") +
   geom_vline(xintercept = 3, linetype = 2)

Loading required package: ggplot2

選定為3類為最佳聚類數(shù)目
用該包下的fviz_cluster函數(shù)可視化一下聚類結(jié)果

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

4. PAM(Partitioning Around Medoids) 圍繞中心點(diǎn)的分割算法

k-means算法取得是均值，那么對(duì)于異常點(diǎn)其實(shí)對(duì)其的影響非常大，很可能這種孤立的點(diǎn)就聚為一類，一個(gè)改進(jìn)的方法就是PAM算法，也叫k-medoids clustering
首先通過fpc包中的pamk函數(shù)得到最佳聚類數(shù)目

library(fpc)
pamk.best <- pamk(dataset)
pamk.best$nc

3

pamk函數(shù)不需要提供聚類數(shù)目，也會(huì)直接自動(dòng)計(jì)算出最佳聚類數(shù)，這里也得到為3
得到聚類數(shù)提供給cluster包下的pam函數(shù)并進(jìn)行可視化

library(cluster)
clusplot(pam(dataset, pamk.best$nc))

5.Calinsky criterion

這個(gè)評(píng)估標(biāo)準(zhǔn)定義[5]如下：
其中，k是聚類數(shù)，N是樣本數(shù)，SSw是我們之前提到過的組內(nèi)平方和誤差， SSb是組與組之間的平方和誤差，SSw越小，SSb越大聚類效果越好，所以Calinsky criterion值一般來說是越大，聚類效果越好

library(vegan)
ca_clust <- cascadeKM(dataset, 1, 10, iter = 1000)
ca_clust$results

可以看到該函數(shù)把組內(nèi)平方和誤差和Calinsky都計(jì)算出來了，可以看到calinski在聚類數(shù)為3時(shí)達(dá)到最大值。

calinski.best <- as.numeric(which.max(ca_clust$results[2,]))
calinski.best

3

畫圖出來觀察一下

plot(fit, sortg = TRUE, grpmts.plot = TRUE)

注意到那個(gè)紅點(diǎn)就是對(duì)應(yīng)的最大值，自帶的繪圖橫軸縱軸取的可能不符合我們的直覺，把數(shù)據(jù)取出來自己單獨(dú)畫一下

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

Warning message:
"Removed 1 rows containing missing values (geom_path)."

這個(gè)看上去直觀多了。這就很清晰的可以看到在聚類數(shù)目為3時(shí)，calinski指標(biāo)達(dá)到了最大值，所以最佳數(shù)目為3

6.Affinity propagation (AP) clustering

這個(gè)本質(zhì)上是類似kmeans或者層次聚類一樣，是一種聚類方法，因?yàn)椴恍枰駅means一樣提供聚類數(shù)，會(huì)自動(dòng)算出最佳聚類數(shù)，因此也放到這里作為一種計(jì)算最佳聚類數(shù)目的方法。
AP算法的基本思想是將全部樣本看作網(wǎng)絡(luò)的節(jié)點(diǎn)，然后通過網(wǎng)絡(luò)中各條邊的消息傳遞計(jì)算出各樣本的聚類中心。聚類過程中，共有兩種消息在各節(jié)點(diǎn)間傳遞，分別是吸引度( responsibility)和歸屬度(availability) 。AP算法通過迭代過程不斷更新每一個(gè)點(diǎn)的吸引度和歸屬度值，直到產(chǎn)生m個(gè)高質(zhì)量的Exemplar（類似于質(zhì)心），同時(shí)將其余的數(shù)據(jù)點(diǎn)分配到相應(yīng)的聚類中[7]

library(apcluster)
ap_clust <- apcluster(negDistMat(r=2), dataset)
length(ap_clust@clusters)

15

該聚類方法推薦的最佳聚類數(shù)目為15，再用熱力圖可視化一下

heatmap(ap_clust)

選x或者y方向看(對(duì)稱)，可以數(shù)出來“葉子節(jié)點(diǎn)”一共15個(gè)

7. 輪廓系數(shù)Average silhouette method

輪廓系數(shù)是類的密集與分散程度的評(píng)價(jià)指標(biāo)。

a(i)是測(cè)量組內(nèi)的相似度,b(i)是測(cè)量組間的相似度，s(i)范圍從-1到1，值越大說明組內(nèi)吻合越高，組間距離越遠(yuǎn)——也就是說，輪廓系數(shù)值越大，聚類效果越好[9]

require(cluster)
library(factoextra)
fviz_nbclust(dataset, kmeans, method = "silhouette")

可以看到也是在聚類數(shù)為3時(shí)輪廓系數(shù)達(dá)到了峰值，所以最佳聚類數(shù)為3

8. Gap Statistic

之前我們提到了WSSE組內(nèi)平方和誤差，該種方法是通過找“肘點(diǎn)”來找到最佳聚類數(shù)，肘點(diǎn)的選擇并不是那么清晰，因此斯坦福大學(xué)的Robert等教授提出了Gap Statistic方法，定義的Gap值為[9]

取對(duì)數(shù)的原因是因?yàn)閃k的值可能很大
通過這個(gè)式子來找出Wk跌落最快的點(diǎn)，Gap最大值對(duì)應(yīng)的k值就是最佳聚類數(shù)

library(cluster)
set.seed(123)
gap_clust <- clusGap(dataset, kmeans, 10, B = 500, verbose = interactive())
gap_clust

Clustering Gap statistic ["clusGap"] from call:
clusGap(x = dataset, FUNcluster = kmeans, K.max = 10, B = 500,     verbose = interactive())
B=500 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
--> Number of clusters (method 'firstSEmax', SE.factor=1): 3
         logW   E.logW       gap     SE.sim
[1,] 5.377557 5.863690 0.4861333 0.01273873
[2,] 5.203502 5.758276 0.5547745 0.01420766
[3,] 5.066921 5.697322 0.6304006 0.01278909
[4,] 5.023936 5.651618 0.6276814 0.01243239
[5,] 4.993720 5.615174 0.6214536 0.01251765
[6,] 4.962933 5.584564 0.6216311 0.01165595
[7,] 4.943241 5.556310 0.6130690 0.01181831
[8,] 4.915582 5.531834 0.6162518 0.01139207
[9,] 4.881449 5.508514 0.6270646 0.01169532
[10,] 4.855837 5.487005 0.6311683 0.01198264

library(factoextra)
fviz_gap_stat(gap_clust)

可以看到也是在聚類數(shù)為3的時(shí)候gap值取到了最大值，所以最佳聚類數(shù)為3

9.層次聚類

層次聚類是通過可視化然后人為去判斷大致聚為幾類，很明顯在共同父節(jié)點(diǎn)的一顆子樹可以被聚類為一個(gè)類

h_dist <- dist(as.matrix(dataset))
h_clust<-hclust(h_dist)
plot(h_clust, hang = -1, labels = FALSE)
rect.hclust(h_clust,3)

10.clustergram

最后一種算法是Tal Galili[10]大牛自己定義的一種聚類可視化的展示，繪制隨著聚類數(shù)目的增加，所有成員是如何分配到各個(gè)類別的。該代碼沒有被制作成R包，可以去Galili介紹頁面)里面的github地址找到源代碼跑一遍然后就可以用這個(gè)函數(shù)了，因?yàn)樵创a有點(diǎn)長我就不放博客里面了，直接放出運(yùn)行代碼的截圖。

clustergram(dataset, k.range = 2:8, line.width = 0.004)

Loading required package: colorspace
Loading required package: plyr

隨著K的增加，從最開始的兩類到最后的八類，圖肯定是越到后面越密集。通過這個(gè)圖判斷最佳聚類數(shù)目的方法應(yīng)該是看隨著K每增加1，分出來的線越少說明在該k值下越穩(wěn)定。比如k=7到k=8，假設(shè)k=7是很好的聚類數(shù)，那分成8類時(shí)應(yīng)該可能只是某一類分成了兩類，其他6類都每怎么變。反應(yīng)到圖中應(yīng)該是有6簇平行線，有一簇分成了兩股，而現(xiàn)在可以看到從7到8，線完全亂了，說明k=7時(shí)效果并不好。按照這個(gè)分析，k=3到k=4時(shí)，第一股和第三股幾本沒變，就第二股拆成了2類，所以k=3是最佳聚類數(shù)目

方法匯總與比較

wine數(shù)據(jù)集我們知道其實(shí)是分為3類的，以上10種判定方法中：

1. 層次聚類和clustergram方法、肘點(diǎn)圖法，需要人工判定，雖然可以得出大致的最佳聚類數(shù)，但算法本身不會(huì)給出最佳聚類數(shù)

2. 除了Affinity propagation (AP) clustering 給出最佳聚類數(shù)為15，剩下6種全都是給出最佳聚類數(shù)為3

3. 選用上次文本挖掘的矩陣進(jìn)行分析(667*1623)

• mclust效果很差，14種模型只有6種有結(jié)果

• bclust報(bào)錯(cuò)

• SSE可以運(yùn)行

• fpc包中的pamk函數(shù)聚成2類，明顯不行

• Calinsky criterion聚成2類

• Affinity propagation (AP) clustering 聚成28類，相對(duì)靠譜

• 輪廓系數(shù)Average silhouette聚類2類

• gap-Statistic跑不出結(jié)果

可見上述方法中有的因?yàn)閿?shù)據(jù)太大不能運(yùn)行，有的結(jié)果很明顯不對(duì)，一個(gè)可能是數(shù)據(jù)集的本身的原因（缺失值太多等），但是也告訴了我們?cè)诖_定最佳聚類數(shù)目的時(shí)候需要多嘗試幾種方法，并沒有固定的套路，然后選擇一種可信度較高的聚類數(shù)目。
最后再把這10種方法總結(jié)一下：

參考文獻(xiàn)
[1]R語言實(shí)戰(zhàn)第二版
[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]輪廓系數(shù)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)

雷峰網(wǎng)版權(quán)文章，未經(jīng)授權(quán)禁止轉(zhuǎn)載。詳情見轉(zhuǎn)載須知。