如何聚类并挑选类

用来展示聚类三张图吧 这篇文章 提到了聚类展示的三张图，这里定义展示一下，也是对聚类概念的加强和理解； 原文地址 如何找到合适的聚类？

不过我们可以基于不同聚类过程中使用的相似性算法和模块划分参数，选择一个最合适的数目。在K-means，PAM和层次聚类中选择合适的聚类数目，这些方法包括直接方法和统计检验方法。 1.直接方法 设置一些适合的划分标准，比如elbow和average silhouette法 2.统计检验方法 就是常用的假设检验方法，比如gap statistic

# 如果因为缺少某些包而安装失败，请先安装其他依赖包再重新安装，安装时间比较久
# 可能需要依赖lme4，cowplot，ggpubr，FactoMineR
# 其中cowplot需要R3.3版本，其他版本可以试试用下面命令安装
devtools::install_url("https://github.com/wilkelab/cowplot/archive/1.0.0.zip")
if(!require(devtools)) install.packages('devtools')
if(!require(factoextra)) devtools::install_github('kassambara/factoextra')

#如果已经安装，跳过即可
pkgs <- c("cluster",  "NbClust")
install.packages(pkgs)

载入R包

library(factoextra)
library(cluster)
library(NbClust)

准备数据

data(iris)
head(iris)

# remove species column
iris.scaled = scale(iris[,-5]) #scale 注意行 列，scale是对列进行操作
head(iris.scaled)

1. 三种聚类方法

1.1 k-means() 和 pama()

这里演示了stat包中的k-means()，cluster包中的pama()的使用，把上面的归一化后的数据分成3个cluster。

# K-means clustering
set.seed(123)
km.res = kmeans(iris.scaled, 3, nstart=25)
# k-means  group number of each observation
str(km.res)
km.res$cluster

# visualize k-means result
fviz_cluster(km.res, data=iris.scaled, geom='point', stand=FALSE, ellipse.type='convex')
# convex / norm

# PAM clustering
pam.res = pam(iris.scaled, 3)
# str(pam.res)
pam.res$cluster

# visualize pam clusters
fviz_cluster(pam.res, stand=FALSE, geom='point', ellipse.type = 'norm')

1.2 hclust()

另一个是R中内建的方法hclust():

# 计算两两间的距离，计算方法比较多，这里选择欧几里德距离
dist.res = dist(iris.scaled, method='euclidean')
# 进行层次聚类，不同的算法说明可以查看函数帮助信息
hc = hclust(dist.res, method='complete')
# 展示聚类结果
plot(hc, labels=FALSE, hang=-1, xlab='hclust')
# 为3个group添加方框，原来还有这个函数，真神奇
rect.hclust(hc, k=3, border=2:4)

# 把层次聚类的结果分成3组 
hc.cut = cutree(hc,k=3)
hc.cut

2. 确定最佳分组数目的3种方法

这里介绍3种常用的方法：Elbow method（肘部法则），silhouette method 和 gap statistic 。

2.1 Elbow method

在聚类分割算法中，比如K-means聚类，为了确定不同的分类，需要保证每个类分组总变异量之和最小； 具体过程如下

1.对不同的k值，分别进行聚类。如K-means中k可以取从1到10 2.对每个k值，计算每个组的组内平方各（within-cluster sum of square）的和 3.绘制k值和组内平方和的总和的趋势图 4.从图上的转折点确定最佳分组数目

#用kmeans看看
set.seed(123)  
# k值从2到15
k.max = 15
data = iris.scaled
# 这里不必手动计算平方总和，kmeans中已经完成计算，直接调用
wss = sapply(1:k.max,
            function(k){kmeans(data,k,nstart=10)$tot.withinss})
plot(1:k.max, wss,
    type='b', pch=19, frame=FALSE,
    xlab = 'Number of cluster K',
    ylab = "Total within-clusters sum of squares")
abline(v=3, lty=2)

从上面的图，可以看出在k=3这个点上，曲线的变化率比较大，建议选择k=3作为最终的结果。当然你还可以看到k越大，组内平方和总和是越来越小，不过随着k变大，分类结果也更加分散，可能不能很好的表现数据聚类想要表达的信息。

也可以直接利用一个叫factoextra的R包来实现，使用它的提供的fviz_nbclust()函数

# fviz_nbclust(x,FUNcluster,method=c('silhourtte','wss'))
fviz_nbclust(iris.scaled, kmeans, method='wss') + geom_vline(xintercept = 3, linetype=2)

# 这里貌似不能直接给出一个推荐值，需要我们自己从图中寻找一个，也就是k=3

下面试试用PAM聚类的结果进行测试

fviz_nbclust(iris.scaled, pam, method='wss') + geom_vline(xintercept = 3, linetype=2)

最终结果也和k-means的聚类结果类似。最后再试试用层次聚类的结果来试试看。使用factoextra提供的 hcut对数据进行聚类并划分分组

fviz_nbclust(iris.scaled, hcut, method='wss') + geom_vline(xintercept = 3, linetype=2)

如果数据比较复杂的话，Elbow method可能会陷入局部最优的结果，随着k的增大，wss反而又增大的情况。

2.2 Average silhouette method

简单来说，该主方法用于评估聚类结果的质量。如果一个聚类的结果比较好，那么它的average silhouette就会比较高。计算过程类似Elbow method：

1.对不同的k值分别进行聚类 1.对每个k的聚类结果分别计算average silhouette（avg.sil）值 1.以k为x轴，avg.sil为y轴绘制连线图 1.avg.sil值最大处就是最优k值

具体的R执行代码可以利用cluster包中的silhouette()函数计算average silhouette值。先看看在K-means聚类中的使用

library(cluster)
k.max = 15
data = iris.scaled
sil = rep(0, k.max)
# k从2到15分别进行kmeans 
for (i in 2:k.max){
    km.res = kmeans(data, centers=i, nstart = 25)
    ss = silhouette(km.res$cluster, dist(data))
    sil[i] = mean(ss[,3])
}
# 画图
plot(1:k.max, sil, type='b', pch=19, frame=FALSE, xlab='Number of clusters k')
abline(v=which.max(sil), lty=2)

使用前面用到的fviz_nbclust()完成

fviz_nbclust(iris.scaled, kmeans, method='silhouette') #所以说这函数 标准化的数据、聚类方法和聚类评价方式

下面再介绍在PAM和层次聚类中使用

# for PAM clustering  
fviz_nbclust(iris.scaled, pam, method='silhouette')

# for hierarchical clustering
fviz_nbclust(iris.scaled, hcut, method='silhouette',hc_method='complete')

讨论:前面介绍的Elbow method得到的最佳值，需要我们手动去看，而Average silhouette method会直接提供一个最佳以供选择。而且K-means和PAM的推荐值是k=2，而层次聚类的推荐值是k=3。结合之前的Elbow method结果，设置k=3比较好。

2.3 Gap statistic method

Gap statistic method可以运用到任何聚类算法里面。该方法先比较不同k值聚类结果中组内变异量的总和（total within intracluster variation）。 利用统计学的假设检验来比较TSS值与那些随机分布的参考数据集之间是否显著差异。

计算过程：（看不懂）

根据不同的k值对实际数据进行聚类并计算$W_k$ 产生B个参考数据集（bootstrap法），按照不同的k值进行聚类，并计算Gap值：$Gap(k) = \frac{1}{B}\sum{b=1}^{B}log(Wk^*) - log(W_k)$ 让$\bar{w} = (\frac{1}{B}\sum{b}log(W{kb}^))$，计算标准差$sd(k) = \sqrt{\frac{1}{B}\sum_b(log(W_{kb}^)-\bar{w})^2}$和标准误$sk = sdk \times \sqrt{1+\frac{1}{B}}$ 选择满足$Gap(k) \ge Gap(k+1) - s_{k+1}$的最小k值

R语言里面的实现方法可以利用cluster包中的 clusGap()来计算

# clusGap(x, FUNcluster, K.max, B = 100, verbose = TRUE, ...)

对3种聚类方法进行测试：

library(cluster)
set.seed(123)
# 一般认为B=500就能得到一个比较好的结果，这里设为50以提高计算速度
gap_stat = clusGap(iris.scaled, FUN=kmeans, nstart=25, K.max=10, B=50)
print(gap_stat, method='firstmax')

# 绘图
plot(gap_stat, frame=FALSE, xlab='Number of cluster k')
abline(v=3, lty=2)

# 使用factoextra
fviz_gap_stat(gap_stat)

这里的最佳聚类数目是用firstmax方法（查看 ?cluster::maxSE）计算的，Tibshirani et al (2001)提出的方法可以参考下面脚本：

# Print
print(gap_stat, method = "Tibs2001SEmax")
# Plot
fviz_gap_stat(gap_stat, 
            maxSE = list(method = "Tibs2001SEmax"))
# Relaxed the gap test to be within two standard deviations
fviz_gap_stat(gap_stat, 
          maxSE = list(method = "Tibs2001SEmax", SE.factor = 2))

PAM和层次聚类的结果

# PAM聚类结果 
set.seed(123)
gap_stat <- clusGap(iris.scaled, FUN = pam, K.max = 10, B = 50)
# Plot gap statistic
fviz_gap_stat(gap_stat)

# Compute gap statistic
set.seed(123)
gap_stat <- clusGap(iris.scaled, FUN = hcut, K.max = 10, B = 50)
# Plot gap statistic
fviz_gap_stat(gap_stat)

原文参考：

  1. Determining the optimal number of clusters: 3 must known methods - Unsupervised Machine Learning
  2. k-means Clustering
  3. K-means another case
  4. 如何选聚类

#R #聚类
comments powered by Disqus