划分聚类—k均值划分和围绕中心点划分及R操作

划分聚类—k均值划分（k-means）和围绕中心点划分（PAM）及R操作

在层次聚类中，无论“自底向上”的聚合策略（层次聚合），或 “自顶向下”的分拆策略（层次分划），所得聚类结果都是具有层级结构的。当需要嵌套聚类或有意义的层次结构时，层次聚类或许特别有用，在生物科学中这种情况很常见。

但在具有数百甚至数千观测值的大样本中，层次聚类通常表现不理想，这时不妨考虑非层次聚类（non-hierarchical clustering）的方法，顾名思义就是聚类结果无层级结构，对象关系比较简单便于解读。如划分聚类（partitioning clustering）是对一组对象进行简单分组的方法，首先指定期望获得的类的个数k，然后将观测值聚为k类，分类的依据是尽可能使组内对象之间的相似性比组间对象相似性更高（或距离更短）。

本篇简介两种常见的划分聚类方法，k均值划分（k-means）和围绕中心点划分（PAM），以及在R语言中的操作。

k均值划分（k-means）

k均值划分（k-means partitioning，k-means）从任意分配的k组数据开始，经迭代优化产生最终结果。 

k-means算法简述

对于给定数据集（a），通过k-means聚类的计算方法可概述如下。

（b）首先随机分配k个聚类中心，这里设定k=3；

（c）分配每个对象到离它最近的聚类中心，得到k个类；每个对象被分配到使下式得到最小值的那一类中：

x_ij表示第i个对象中第j个变量的值，`x_kj表示第k个类中第j个变量的均值，n是对象个数，p是变量个数。

（d）计算每个类的质心，将该质心作为新的聚类中心；

（e）分配每个对象到离它最近的聚类中心，得到k个类，分类标准同（c）；

（f）重复（d）（e）过程，直到对象分类不再变动或达到最大迭代次数为止。

可以看到，k-means是一种线性方法，使用数据局部结构构建聚类簇，把对象分成k组并使得对象到指定的聚类中心的平方总和为最小。它基于对象之间的欧几里得距离，通过不断迭代将总误差平方和（TESS，或称组内平方和）最小化，分组标准和Ward最小方差聚类相似。

k-means聚类在处理较大样本量的数据集时，通常比层次聚类更有优势。但是均值的使用意味着所有的变量必须是连续的，并且受异常值影响较大。对象聚类与初始指定k的数量密切相关，且对初始中心值的选择也很敏感。

和NMDS等迭代方法类似，如果到达最大迭代次数时仍未收敛，则可能会陷入局部最小值。为了尽可能避免该问题，通常运行多个同步的k-means并选择总体组内平方和最低的结果。

R语言执行k-means

在R中，可使用kmeans()执行k-means聚类。

示例数据和R代码的获取链接：

https://pan.baidu.com/s/1WEHBd2gSzA9Q3s99qhEBFg

一个简单示例

示例数据集GeneExpExample5000，记录了某项研究中获得的肾脏（Kidney）和肝脏（Liver）RNA样本的基因表达值信息，从中取了一部分数据作为示例数据集演示，共5000行（5000个基因的表达值）。

使用k-means的方法划分样本归类。

#示例数据集，该数据集实际上来自 DEGseq 包，详情可在加载 DEGseq 包后 ?GeneExpExample5000 查看
#已知的先验分组，Kidney_group=c(7, 9, 12, 15, 17, 18, 20)；Liver_group=c(8, 10, 11, 13, 14, 16, 19)
gene_express <- read.delim('GeneExpExample5000.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
reads_count <- t(gene_express[7:20])

#当变量量纲不同时，需要标准化数据使方差有意义；变量标准化还能降低极大值的影响
#本示例变量量纲一致，但如果觉得基因表达值的数量级相差过大可能导致较大离散度，可选标准化
#reads_count <- scale(reads_count)

#k-means 聚类，详情 ?kmeans
#由于 k-means 在开始时随机指定中心点，每次调用函数时可能获得不同的方案，使用 set.seed() 保证可重复性
#设定 k=2，并生成 25 种初始配置（初始中心点）并输出最好的一个
set.seed(123)
gene_kmeans <- kmeans(reads_count, centers = 2, nstart = 25)
summary(gene_kmeans)
gene_kmeans$cluster

#轮廓图评估
plot(silhouette(gene_kmeans$cluster, vegan::vegdist(reads_count, method = 'euclidean')))

显示k-means聚类方法可以区分两组样本，轮廓图表明分类准确。

对于最佳聚类数量的选择

上述示例是已知数据先验分组的情况。

但有时，我们可能并不知道样本来源，不知道聚为多少类是合适的。此时对于k值的选择，R中也给出了一些参考标准。

#随机生成模拟数据
set.seed(123)
rand <- rbind(matrix(rnorm(150, mean = 10, sd = 0.3), ncol = 5),
matrix(rnorm(150, mean = 3, sd = 0.2), ncol = 5),
matrix(rnorm(150, mean = 1, sd = 0.1), ncol = 5),
matrix(rnorm(150, mean = 6, sd = 0.3), ncol = 5),
matrix(rnorm(150, mean = 9, sd = 0.3), ncol = 5))

#NbClust 包的最佳聚类数量评估，详情 ?NbClust
library(NbClust)
nc <- NbClust(rand, distance = 'euclidean', min.nc = 2, max.nc = 10, method = 'kmeans', index = 'hubert')
plot(nc)

模拟生成了5组均值和方差不等的数据，NbClust的方法给出评估，聚类为5组是最佳的（在cluster=5的位置出现明显的拐点）。

#vegan 包的最佳聚类数量评估，详情 ?cascadeKM
library(vegan)
cc <- cascadeKM(rand, inf.gr = 2, sup.gr = 10, iter = 100, criterion = 'ssi')
plot(cc)

同样显示聚类为5组是最佳的（k=5时，变量的结构最清晰，ssi值最大）。

至于最终聚类是否和模拟生成的5组数据完全对应，运行k=5时的k-means对比观察即可，就不再演示了。

非欧式距离的应用

如果不期望使用欧式距离，或者只提供了非欧式距离的距离矩阵而缺少原始数据时，该方法可供参考。

例如期望通过对象间的Bray-curtis距离完成k-means。由于Bray-curtis距离并非欧氏距离，可以首先进行平方根转化，使其具有欧式几何属性，然后将转化后的距离应用主坐标分析（PCoA），并将所有PCoA轴作为k-means的输入，完成聚类。

#非欧式距离的应用，如 Bray-curtis 距离
#读取距离矩阵
dis_bray <- read.delim('bray_distance.txt', sep = '\t', row.names = 1, stringsAsFactors = FALSE, check.names = FALSE)

#平方根转化
dis <- (as.dist(dis_bray))^0.5

#PCoA
pcoa <- cmdscale(dis, k = (nrow(dis_bray) - 1), eig = TRUE)

#获取对象排序轴
site <- pcoa$point

#k-means，聚为 3 类
bray_kmeans <- kmeans(site, centers = 3, nstart = 25)
bray_kmeans$cluster

围绕中心点划分（PAM）

由于k-means聚类基于最小化组内欧式距离平方和，所以它对异常值敏感。相比之下，围绕中心点划分（Partitioning Around Mediods，PAM）是一个更稳健的方法。

PAM算法简述

PAM与k-means的过程大致一致，区别在于PAM用最具代表性的观测值来代替k-means中的质心，并且可以使用任意距离测度。

（1）随机选择k个观测对象作为聚类中心；

（2）计算每个聚类中心到每个对象的距离，可以为任意的距离测度；并分配每个对象到离它最近的聚类中心，得到k个类；

（3）计算每个聚类中心到该类中每个对象的距离的总和（D_sum）；

（4）在各类中，重新选择一个对象作为新的聚类中心，并计算新的聚类中心到每个对象的距离；并分配每个对象到离它最近的聚类中心，得到新的k个类；

（5）计算新的每个聚类中心到该类中每个对象的距离的总和（D’_sum）；

（6）如果D’_sum<D_sum，则使用新的聚类中心；

（7）重复（4）-（6）过程，直到对象分类不再变动为止。

可以看到，PAM的目的是使组内对象之间的相异性总和最小，这种相异性可以为任意的距离测度，不再局限于欧几里得距离。这个特征允许PAM容纳混合数据类型，不仅限于连续变量，因此PAM在方法上比k-means更灵活。

R语言执行PAM

在R中，可使用cluster包pam()执行PAM聚类。

以上文示例为例继续。

library(cluster)

#PAM 聚类，详情 ?pam
#可以输入原始数据，指定距离测度，同样以欧几里得距离为例
gene_pam <- pam(reads_count, k = 2, metric = 'euclidean')
summary(gene_pam)

#也可以首先计算距离矩阵，将距离矩阵作为输入，同样以欧几里得距离为例
dis_euc <- vegan::vegdist(reads_count, method = 'euclidean')
gene_pam <- pam(dis_euc, k = 2, diss = TRUE)
summary(gene_pam)

#可通过轮廓图评估上述分类是否可靠
plot(gene_pam)

显示PAM聚类方法可以区分两组样本，轮廓图表明分类准确。

友情链接