究竟應該聚多少類？聚類分析
分析(Classification analysis) ，前者是無監督學習，而后者是有監督學習。
聚類分析(Cluster analysis)是一組將研究對象分為相對同質的群組(Clusters)的統計分析技術。 聚類分析區別于分類分析(Classification analysis) ，前者是無監督學習，而后者是有監督學習。
無監督學習也就是預先沒有類別的標準，K-means方法是最常見的聚類方法之一，在R語言的kmeans()函數中，聚類數目的參數centers是必須人為輸入的。那么究竟如何確定聚類的數目呢?這里引用一句話(參考文獻1)"While there are no best solutions for the problem of determining the number of clusters to extract..."，這里的潛臺詞是目前有很多方法，只是不知道選哪種。Anyway，大家知道有哪些方法嗎?
下面通過一個實例對幾種常見確定聚類數目的方法進行介紹。
數據
收集了96例患者的兩個生化指標(a、b)，研究者擬根據a、b兩個指標將96例患者分類若干組，進而對不同的組深入分析，制定不同的治療方案。a、b兩個指標的關系如下圖(下面是模擬數據，實際情況通常沒有那么明顯的距離)。那么，我們該聚多少類呢?
- set.seed(2018)
- n=100;g=6
- mydata=data.frame(a=unlist(lapply(1:g,function(i) rnorm(n/g, runif(1)*i^2))), b=unlist(lapply(1:g,function(i) rnorm(n/g, runif(1)*i^2))))
- mydata=scale(mydata)
1. 組內平方和
該方法是畫出不同聚類數目(通常是1到10)對應的組內平方和，以組內平方和基本不再明顯變化為標準(較為主觀)，從而確定聚類的數目。如下圖，我們可以看到組內平方和在聚類數目為5之后，基本不再變化，所以可以確定聚類數目為5。
- ss=(nrow(mydata)-1)*sum(apply(mydata,2,var))
- for(i in 2:10){ss[i]=kmeans(mydata,centers=i)$tot.withinss}
- plot(1:10,ss,type="b",xlab="聚類數目",ylab="組內平方和")
2. Medoids周圍分類法
這個方法的操作比較簡單，顧名思義就是通過partitioning around medoids來估計最優的聚類數目。從下圖的橢圓數目可以看出聚類數目應該為5。
- library(fpc)
- library(cluster)
- fit=pamk(mydata)
- plot(pam(mydata,fit$nc))
3. Calinsky準則
Calinsky是一個定量的評價指標，該值越大越好。下圖的右邊可以看出Calinsky較大值對應的聚類數目為5。
- library(vegan)
- fit=cascadeKM(mydata,inf.gr=1,sup.gr=10,iter=1000)
- plot(fit,sortg=TRUE,grpmts.plot=TRUE)
4. BIC
沒錯，此處的BIC就是貝葉斯信息準則(Bayesian Information Criterion)，但與其用于模型的擬合判斷(越小越好)不同的是，這里用于判斷聚類數目的判斷標準為越大越好。如下圖，右下角的圖例表示各種不同的模型，從折線可以看出，所有的模型都在聚類數目為5時，BIC取得較大。
- library(mclust)
- fit=Mclust(mydata)
- plot(fit)
5. AP法
AP表示Affinity propagation的意思，旨在通過Affinity propagation的方法估計最優的聚類數目。如下圖，左側和上側的彩色帶的顏色種類表示聚類的數目，為5。
- library(apcluster)
- fit=apcluster(negDistMat(r=2),mydata)
- heatmap(fit)
6. 多目標決策
NbClust函數提供30種的判斷指標，我們可以用類似投票的原理，看看通過30種指標選出的最優聚類數目。如下圖，最后一句表示，聚類數目應該取3，這個結果與前述所有結果不一致，根據原始數據的情況，這個方法不太可靠，建議選擇前述的方法。
- library(NbClust)
- fit=NbClust(mydata,min.nc=2,max.nc=10,method="kmeans",index="alllong")
參考文獻
https://www.statmethods.net/advstats/cluster.html
https://stackoverflow.com/questions/15376075/cluster-analysis-in-r-determine-the-optimal-number-of-clusters
