R语言从入门到精通：Day15（聚类分析）

2019-12-07 本文已影响0人 fc7587d03c83

上次推文的结尾预告了这一次的内容是关于聚类分析的，聚类分析广泛用于包括生物学在内的诸多领域。我们熟悉的，通过对DNA microarray data进行聚类分析来获得基因表达模式，从而帮助理解生物的正常发育以及导致许多疾病的根本原因。相信今天的这部分内容对大家会很有帮助。

------ 测试数据和代码见文末客服二维码

聚类分析是一种数据归约技术，旨在揭露一个数据集中观测值的子集。它可以把大量的观测值归约为若干个类。

这里的类被定义为若干个观测值组成的群组，群组内观测值的相似度比群间相似度高。这不是一个精确的定义，从而导致了各种聚类方法的出现。最常用的两种聚类方法是层次聚类(hierarchical agglomerative clustering)和划分聚类(partitioning clustering)。

在层次聚类中，每一个观测值自成一类，这些类每次两两合并，直到所有的类被聚成一类为止。在划分聚类中，首先指定类的个数K，然后观测值被随机分成K类，再重新形成聚合的类。对于划分聚类来说，最常用的算法是K均值(K-means)等。

在开始讨论这两类聚类方法之前，我们先要熟悉一下聚类方法都需要遵循的步骤。

**需提前按照的R包：cluster、NbClust、flexclust、fMultivar、ggplot2和rattle。

一个全面的聚类分析一般会包括以下11个典型步骤：

1.选择合适的变量；

2.缩放数据（最常用的方法是将每个变量标准化为均值=0和标准差=1的变量。其他的替代方法包括每个变量被其最大值相除或该变量减去它的平均值并除以变量的平均绝对偏差）；

3.寻找异常点；

4.计算距离；

5.选择聚类算法；

6.获得聚类方法；

7.确定类的数目（NbClust包的函数NbClust()可以解决这个问题）；

8.获得最终的聚类解决方案；

9.结果可视化；

10.解读类；

11.验证结果（包fpc、包clv和包clValid包含了评估聚类解的稳定性的函数）。

层次聚类

首先，对于层次聚类而言，起初每一个实例或观测值属于一类。每一次把两类聚成新的一类，直到所有的类聚成单个类为止，算法如下:

(1) 定义每个观测值(行或单元)为一类;

(2) 计算每类和其他各类的距离;

(3) 把距离最短的两类合并成一类，这样类的个数就减少一个;

(4) 重复步骤(2)和步骤(3)，直到包含所有观测值的类合并成单个的类为止。

不同的层次聚类方法主要区别在于步骤（2），即计算类之间的距离的方法不同。常见的层次聚类方法见表1。

表1，五种常见的层次聚类方法。

层次聚类方法可以用函数hclust()来实现，格式是hclust(d, method=)，其中参数d是通过 dist()函数产生的距离矩阵，并且参数method包括 "single"、"complete"、"average"、 "centroid"和"ward"（和表格里面的5种一一对应）。层次聚类的测试数据集主要来自于包flexclust中的数据集nutrient。

采用平均联动的聚类结果如图1。

图1：平均联动聚类结果

树状图应该从下往上读，它展示了这些条目如何被结合成类。每个观测值起初自成一类，然后相距最近的两类合并。合并继续进行下去，直到所有的观测值合并成一类。高度刻度代表了该高度类之间合并的判定值。但是这幅图并不能指出聚类的适当个数。包NbClust提供了众多的指数来确定在一个聚类分析里类的最佳数目。四个评判准则赞同聚类个数为2，四个判定准则赞同聚类个数为3，结果如图2。

图2：理想的聚类个数

我们尝试聚类个数为5（四个判定准则赞同，如图2）的情形。

图3：聚类个数为5

层次聚类

当需要嵌套聚类和有意义的层次结构时，层次聚类特别有用。在某种意义上分层算法是严苛的，一旦一个观测值被分配给一个类，它就不能在后面的过程中被重新分配。另外，层次聚类难以应用到有数百甚至数千观测值的大样本中。我们这里将讨论两种方法：K-means和基于中心点的划分(PAM)。

K-means

从概念上讲，K-means算法如下:

(1) 选择K个中心点(随机选择K行);

(2) 把每个数据点分配到离它最近的中心点;

(3) 重新计算每类中的点到该类中心点距离的平均值(也就说，得到长度为p的均值向量，这里的p是变量的个数);

(4) 分配每个数据到它最近的中心点;

(5) 重复步骤(3)和步骤(4)直到所有的观测值不再被分配或是达到最大的迭代次数(R把10次作为默认迭代次数)。

在R中K-means的函数格式是kmeans(x, centers)，这里参数x表示数值数据集(矩阵或数据框)，参数centers是要提取的聚类数目。

由于K均值聚类在开始要随机选择k个中心点，在每次调用函数时可能获得不同的方案。使用函数set.seed()可以保证结果是可复制的。此外，聚类方法对初始中心值的选择也很敏感。函数kmeans()有一个参数nstart尝试多种初始配置并输出最好的一个。例如，加上nstart=25 会生成25个初始配置。通常推荐使用这种方法。

不像层次聚类方法，K均值聚类要求你事先确定要提取的聚类个数。同样，包NbClust可以用来作为参考，同时我们也提供了一个自定义的函数wssplot来帮助你确定类的个数。

（划分聚类的测试数据集来自于包rattle的数据集wine，为了验证分类结果的准确性，我们选择先放弃第一个变量[类型]，进行聚类分析，再将结果和第一个变量对比，看看能否恢复已知的类型）

图4:函数wssplot的结果

图4表明从一类到三类下降得很快(之后下降得很慢)，建议选用聚类个数为三的解决方案。

图5:包NbClust的结果

图5也推荐聚为三类。

利用包flexclust中的函数randIndex()可以衡量两个分类之间的差异，结果为0.897。说明我们聚为3类的选择是很准确的。（代码联系文末客服，识别二维码领取）

基于中心点的划分

因为K-means聚类方法是基于均值的，所以它对异常值是敏感的。一个更稳健的方法是围绕中心点的划分(PAM)。与其用质心表示类，不如用一个最有代表性的观测值来表示(称为中心点)。K-means聚类一般使用欧几里得距离，而PAM可以使用任意的距离来计算。因此，PAM可以容纳混合数据类型，并且不仅限于连续变量（PAM算法和k-means聚类很类似，就不赘述了）。

包cluster中的函数pam()使用基于中心点的划分方法。格式是pam(x, k, metric="euclidean", stand=FALSE)，这里的参数x表示数据矩阵或数据框，参数k表示聚类的个数，参数metric表示使用的相似性/相异性的度量，而参数stand是一个逻辑值，表示是否有变量应该在计算该指标之前被标准化（这里测试数据集和k-means聚类方法一样）。

图6:PAM聚类方法的结果

但是，PAM在本例中的表现不如K均值，函数randIndex()结果下降为了0.699。

到这里，常见的两种聚类方法已经讨论完毕，但是聚类分析是一种旨在识别数据集子组的方法，并且在此方面十分擅长。事实上，它甚至能发现不存在的类。这里，代码中提供一个对完全随机数据进行PAM聚类的例子。图7展示了这个数据集的分布情况。其中并不存在聚类。但是函数wssplot()和包NbClust的结果都建议分类为2或3。利用PAM方法进行双聚类结果如图8。

图7：完全随机数据的分布

图8：PAM双聚类结果

很明显划分是人为的。实际上在这里没有真实的类。NbClust包中的立方聚类规则(Cubic Cluster Criteria，CCC)往往可以帮助我们揭示不存在的结构。当CCC的值为负并且对于两类或是更多的类递减时，就是典型的单峰分布（如图9）。

图9：CCC规则的示意图

当然，你也可以尝试多种聚类方法，如果结果都很类似，就可以确信你的聚类结果是准确的了。综合来说，聚类分析是一个宽泛的话题，而R有一些最全面的方案来实施现有的方法。想要了解更多，可以CRAN的聚类分析和有限混合模型部分，见如下链接。https://cran.r-project.org/web/views/Cluster.html

本期干货

- day15-聚类分析代码-

关注“科研猫”公众号，联系客服

胖雨小姐姐

折耳猫小姐姐

领取