算法交流：聚类

╰半夏微凉° 2021-11-09 19:50 495阅读 0赞

“聚类分析是一种数据归约技术，旨在揭露一个数据集中观测值的子集。它可以把大量的观测值归约为若干个类。这里的类被定义为若干个观测值组成的群组，群组内观测值的相似度比群间相似度高。”——《R 语言实战》第二版

常用的两种聚类方法有：

1.  **层次聚类（hierachical clustering）**：每个数据点为一小类，两两通过树的方式合并，直到所有的数据点汇成一类。常用算法：
    
    1.  **单联动（single linkage）**：类 A 中的点与类 B 中的点间的最小距离。适合于细长的类。
    2.  **全联动（complete linkage ）**：类 A 中的点与类 B 中的点间的最大距离。适合于相似半径的紧凑类，对异常值敏感。
    3.  **平均联动（average linkage）**：类 A 中的点与类 B 中的点间的平均距离，也称为 UPGMA。适合于聚合方差小的类。
    4.  **质心（centroid）**：类 A 与类 B 的质心的距离。质心的定义是“类的变量均值向量”。对异常值不敏感，但表现可能稍弱。
    5.  **Ward 法（ward.D）**：两类之间的所有变量的方差分析平方和。适合于仅聚合少量值、类别数接近数据点数目的情况。
2.  **划分聚类（partitioning clustering）**：事先指定类数 KK，然后聚类。
    
    1.  **K均值（K-means）**：
    2.  **中心划分（Partitioning Around Medoids，即 PAM）**：

## 聚类步骤 ##

聚类是一个多步骤过程。典型的步骤有 11 步。

1.  **变量选取**：例如你需要对实验数据进行聚类，那么你需要仔细思考哪些变量会对聚类产生影响，而哪些变量是不需加入分析的。
2.  **缩放数据**：最常用的方法是标准化，将所有变量变为 ¯¯¯x=0,SE(x)=1x¯=0,SE(x)=1 的变量。
3.  **筛选异常**：筛选和删除异常数据对于某些聚类方法是很重要的，这可以借助 R 的 `outliers/mvoutlier` 包。或者，你可以换用一种受异常值干扰小的方法，比如中心划分聚类。
4.  **距离计算**：两个数据点间的距离度量有若干种，我们在[下一小节][Link 1]专门讨论。
5.  **选择聚类方法**：每个方法都有其优缺点，请仔细斟酌。
6.  **确定一种或多种聚类方法**
7.  **确定类数**：常用的方法是尝试使用不同的类数进行聚类，然后比较结果。R 中的 NbClust 包提供了一个拥有超过30个指标的 NbClust() 函数。
8.  **最终方案**
9.  **可视化**：层次聚类使用树状图；划分聚类使用可视化双变量聚类图。
10. **解释每个类**：通常会对每个类进行汇总统计（如果是连续型数据），或者返回类的众数/类别分布（如果含类别型数据）。
11. **验证**：聚类结果有意义吗？更换聚类方法能得到类似结果吗？R 中的 fpc, clv 与 clValid 包给出了评估函数。

## 距离计算：dist() 函数 ##

数据点之间的距离有多种度量方法。在 R 的 `dist()` 函数参数中，默认选项 `method=euclidean`。函数中内置的距离方法选项有：

*  **欧几里得距离（euclidean）**：L2L2 norm. 在拥有 nn 个变量的数据集中，数据点 ii 与 jj 的欧式距离是
    
    ![20190801173400667.png][]
 *  **最大距离（maximun）**：L∞L∞ norm. 两点之间的最大距离，即p→∞p→∞ 时的明科夫斯基距离。
 *  **曼哈顿距离（manhattan）**: L1L1 norm.
    
    ![20190801173419750.png][]
 *  **堪培拉距离（canberra）**：
    
    ![20190801173429889.png][]
 *  **二进制距离（binary）**：非 0 变量为 1，为 0 变量为 0.然后根据 0 的比例确定距离。
 *  **明科夫斯基距离（minkowski）**：LpLp norm.
    
    ![20190801173439257.png][]
    
    当 p=1p=1 时，即曼哈顿距离；p=2p=2 时，欧几里得距离；p→∞p→∞ 时，切比雪夫距离。

本文利用 [Iris data][]，数据内容是萼片、花瓣的长与宽。先读取数据：

datapath <- paste(getwd(), '/data/iris.data.csv', sep='')  # 我将其改成了 csv 格式
    iris.raw <- read.csv(datapath, head=F)
    head(iris.raw)
    
    # 去掉非数值的第 5 列
    iris <- iris.raw[,-c(5)]

<table> 
 <thead> 
  <tr> 
   <th>V1</th> 
   <th>V2</th> 
   <th>V3</th> 
   <th>V4</th> 
   <th>V5</th> 
  </tr> 
 </thead> 
 <tbody> 
  <tr> 
   <td>5.1</td> 
   <td>3.5</td> 
   <td>1.4</td> 
   <td>0.2</td> 
   <td>Iris-setosa</td> 
  </tr> 
  <tr> 
   <td>4.9</td> 
   <td>3.0</td> 
   <td>1.4</td> 
   <td>0.2</td> 
   <td>Iris-setosa</td> 
  </tr> 
  <tr> 
   <td>4.7</td> 
   <td>3.2</td> 
   <td>1.3</td> 
   <td>0.2</td> 
   <td>Iris-setosa</td> 
  </tr> 
  <tr> 
   <td>4.6</td> 
   <td>3.1</td> 
   <td>1.5</td> 
   <td>0.2</td> 
   <td>Iris-setosa</td> 
  </tr> 
  <tr> 
   <td>5.0</td> 
   <td>3.6</td> 
   <td>1.4</td> 
   <td>0.2</td> 
   <td>Iris-setosa</td> 
  </tr> 
  <tr> 
   <td>5.4</td> 
   <td>3.9</td> 
   <td>1.7</td> 
   <td>0.4</td> 
   <td>Iris-setosa</td> 
  </tr> 
 </tbody> 
</table>

### 算例：欧式距离 ###

R 内置的 `dist()` 函数默认使用欧式距离，以下与 `dist(iris, method='euclidean')` 等同。比如我们来计算 iris 的欧氏距离：

iris.e <- dist(iris)
    # 显示前 3 个数据点间的欧式距离。这是一个对角线全0的对称矩阵
    as.matrix(iris.e)[1:3, 1:3]

## 层次聚类算例 ##

层次聚类（HC）的逻辑是：**依次把距离最近的两类合并为一个新类，直至所有数据点合并为一个类。**

层次聚类的 R 函数是 `hclust(d, method=)` ，其中 d 通常是一个 `dist()` 函数的运算结果。

仍然使用上文的 Iris data 数据。

### 标准化 ###

尽管标准化不一定会用到，但是这是通常的手段之一。

iris.scaled <- scale(iris)
    head(iris)
    head(iris.scaled)

### 树状图与热力图 ###

层次聚类中的树状图是不可少的。一般使用欧式距离，聚类方法另行确定。

热力图不是必须的。

# 选择聚类方法
    dist_method <- "euclidean"
    cluster_method <- "ward.D"
    
    iris.e <- dist(iris.scaled, method=dist_method)
    iris.hc <- hclust(iris.e, method=cluster_method)
    plot(iris.hc, hang=-1, cex=.8, main="Hierachical Cluster Tree for Iris data")

![png][]

heatmap(as.matrix(iris.e),labRow = F, labCol = F)

![png][png 1]

### 确定聚类个数 ###

使用 NbClust 包。

library(NbClust)
    nc <- NbClust(iris.scaled, distance=dist_method, 
                  min.nc=2, max.nc=10, method=cluster_method)

*** : 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.

![png][png 2]

*** : 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:                                                
    * 9 proposed 2 as the best number of clusters 
    * 5 proposed 3 as the best number of clusters 
    * 4 proposed 5 as the best number of clusters 
    * 2 proposed 6 as the best number of clusters 
    * 1 proposed 8 as the best number of clusters 
    * 3 proposed 10 as the best number of clusters 
    
                       ***** Conclusion *****                            
     
    * According to the majority rule, the best number of clusters is  2 
     
     
    *******************************************************************

![png][png 3]

# 每个类数的投票数
    table(nc$Best.n[1,])
    barplot(table(nc$Best.n[1,]), xlab="Number of Clusters", ylab="Number of Supporting", 
            main="Determine the Number of Clustering")

0  2  3  5  6  8 10 
     2  9  5  4  2  1  3

![png][png 4]

### 完成聚类 ###

从上图我们可以确定聚类数量（两类），但是原数据指出应该是三类。以下 `cutree` 以及最后一个函数中使用3类。

*注：类别已确定是 3 种。在此前提下，ward.D 下的 HC 聚类效果最好，但是 NbClust 投票建议仍是分为 2 类。 complete 下的 HC 聚类投票建议是 3 类，但是效果反而不如 ward 法。*

cluster_num <- 3
    clusters <- cutree(iris.hc, k=cluster_num)
    table(clusters)  # 每类多少个值

clusters
     1  2  3 
    49 74 27

# 每类的各变量中位数
    aggregate(iris, by=list(cluster=clusters), median)

<table> 
 <thead> 
  <tr> 
   <th>cluster</th> 
   <th>V1</th> 
   <th>V2</th> 
   <th>V3</th> 
   <th>V4</th> 
  </tr> 
 </thead> 
 <tbody> 
  <tr> 
   <td>1</td> 
   <td>5.0</td> 
   <td>3.4</td> 
   <td>1.5</td> 
   <td>0.20</td> 
  </tr> 
  <tr> 
   <td>2</td> 
   <td>6.0</td> 
   <td>2.8</td> 
   <td>4.5</td> 
   <td>1.45</td> 
  </tr> 
  <tr> 
   <td>3</td> 
   <td>6.9</td> 
   <td>3.1</td> 
   <td>5.8</td> 
   <td>2.20</td> 
  </tr> 
 </tbody> 
</table>

# 画出矩形框
    plot(iris.hc, hang=-1, cex=.8, main="Hierachical Cluster Tree for Iris data")
    rect.hclust(iris.hc, k=cluster_num)

![png][png 5]

## 用 MDS 可视化结果 ##

使用多维缩放（Multidimensional Scaling）方法进行可视化。原数据的三个种类被标记为三种不同的点形状，聚类结果则以颜色显示。

可以看到setose品种聚类很成功，但有一些virginica品种的花被错误和virginica品种聚类到一起。

mds=cmdscale(iris.e,k=2,eig=T)
    x = mds$points[,1]
    y = mds$points[,2]
    
    library(ggplot2)
    p=ggplot(data.frame(x,y),aes(x,y))
    p+geom_point(size=3, alpha=0.8, aes(colour=factor(clusters),
                 shape=iris.raw[,5]))

![png][png 6]

### 附：全联动HC聚类图 ###

作为对比。全联动 NbClust 投票结果是3类，在此不再列出。

dist_method <- "euclidean"
    cluster_method <- "complete"
    
    iris.e <- dist(iris.scaled, method=dist_method)
    iris.hc <- hclust(iris.e, method=cluster_method)
    
    cluster_num <- 3
    clusters <- cutree(iris.hc, k=cluster_num)
    table(clusters)  # 每类多少个值
    
    # 画出矩形框
    plot(iris.hc, hang=-1, cex=.8, main="Hierachical Cluster Tree for Iris data")
    rect.hclust(iris.hc, k=cluster_num)

clusters
     1  2  3 
    49 24 77

![png][png 7]

mds=cmdscale(iris.e,k=2,eig=T)
    x = mds$points[,1]
    y = mds$points[,2]
    
    p=ggplot(data.frame(x,y),aes(x,y))
    p+geom_point(size=3, alpha=0.8, aes(colour=factor(clusters),
                 shape=iris.raw[,5]))

![png][png 8]

可以看出，setosa 聚类非常好； virginica 与 versicolor 的效果则是惨不忍睹。

*本文内容大量参考：*

1.  《R 语言实战》第二版 第16章。
2.  [此网页][Link 2]

[Link 1]: https://wklchris.github.io/R-clustering.html#%E8%B7%9D%E7%A6%BB%E8%AE%A1%E7%AE%97
[20190801173400667.png]: /images/20211109/8e541d6076714874b94c52aedb7d9658.png
[20190801173419750.png]: /images/20211109/c94dc65c11264837b48fc3267a97824b.png
[20190801173429889.png]: /images/20211109/329db51f322242899e15db20a12d9c5b.png
[20190801173439257.png]: /images/20211109/6e08da2386e2429bbabc26b69afc7b12.png
[Iris data]: https://archive.ics.uci.edu/ml/datasets/Iris
[png]: /images/20211109/8e293c1fd0d34073b2fa8ec8a643f65c.png
[png 1]: /images/20211109/3299adaebd2d4a5882f1f7e1f7287bd8.png
[png 2]: /images/20211109/6d7eb057fc5c4635a63b91c9fd523bdf.png
[png 3]: /images/20211109/dc2ceedaada9426c94ef257465e36486.png
[png 4]: /images/20211109/e1deff73a2db4381b93eae7e9050d825.png
[png 5]: /images/20211109/02d804917cef4e75a36d7bdff4e975e6.png
[png 6]: /images/20211109/3688453e31cd4e75961a612b0f022ecd.png
[png 7]: /images/20211109/0d548e48efe748a0804fce30a4a4afde.png
[png 8]: /images/20211109/6d560e96ee614ac39e88172bb1e8ff52.png
[Link 2]: http://xccds1977.blogspot.com/2012/01/r.html