层次聚合分类分析及其在R中的计算

层次聚合分类及其在R中的计算

层次聚类（hierarchical clustering）试图在不同层次对数据集进行划分，从而形成树形的聚类结构。数据集划分可采用“自底向上”的聚合策略，也可采用“自顶向下”的分拆策略。本篇首先简介“自底向上”的聚合策略。

层次聚合分类（hierarchical agglomerative classification），是一种有层次的、凝聚的数分值类方法，通过将单个样本聚集成更大的聚类群来创建分层嵌套的样本组。

层次聚合分类的结果一般以聚类树表示，它表示将样本或聚类群连接在一起的层次结构。在聚类树中，分支的高度代表了距离的远近，而节点周围分支的方向则可以是任意的。

在该方法中影响聚类结果的主要因素包括数据预处理方式、距离测度和聚类算法的选择。

某些原始数据可能不适合直接使用，例如受极端值、高权重的变量影响较大等，在分析前可选剔除离群点或进行合理的转化。数据预处理前后的聚类结果可能会相差明显。

选择恰当的距离测度对于层次聚合分类非常重要，但并非所有距离或相似度度量都可以用于聚类算法中。例如，Ward最小方差聚类方法要求距离测度必须是欧式距离属性。对于距离测度的简介，可参考前文。

聚类算法表示一组决定将样本分配给哪个组（聚类群）的规则。可以为该过程打个比方，在一个满是人（样本）的酒吧中，决定坐在哪一张桌子上（聚类群、组），有些人是我的好朋友（低距离，或高相似性），我倾向和他们坐在一起；有些人我比较陌生（高距离，或低相似性），我通常和他们离的较远。

接下来简介层次聚合分类及其在R语言中的计算方法。

层次聚合分类方法概述

首先简介几种常见的层次聚合分类模式。

单连接聚合聚类（单联动）

单连接聚合聚类（single linkage agglomerative clustering），或称最近邻体分类（nearest neighbour sorting），依据最短的成对距离（或最大的相似性）聚合对象（Lukaszewicz, 1951）。

以“坐在酒吧里的每个人”为例，都更倾向于和最近的朋友们坐在一桌（可理解为相似度最高或距离最短），而不考虑和陌生人凑在一桌（可理解为相似度较低或距离较大）。

首先计算相似度矩阵，然后根据对象对的组成相似度对对象进行排序。最相似的对将形成第一个簇，即样本212与样本214最相似归入第一簇，在S₂₀=0.600处合并；排第二位的最相似对将形成第二个簇，即样本431与样本432归入第二簇，在S₂₀=0.500处合并；以此类推，直到将所有对象加入簇中。

因此在单连接聚合聚类中，两个聚类群之间的距离定义为每个簇中两个对象之间的最短距离，这种计算规则使聚类很容易实现。

不难看出，单连接聚合聚类使一个对象很容易聚合到一个分类组中，因为一次连接足以导致融合，因此单连接聚合聚类也被称为“最亲密朋友”（closest friend）法。并且单连接聚类产生的梯度非常明显，即容易识别数据的梯度。

单连接聚类由于其是空间压缩的，这种算法使相异较大的对象很快合并在一起。因此其缺点是产生的分类组之间的差异通常不清晰，即从分区的角度来看比较难以解读，并且错分的可能性也较大。

完全连接聚合聚类（全联动）

与单连接聚合聚类相反，完全连接聚合聚类（Complete linkage agglomerative clustering），或称最远邻体分类（Furthest neighbour sorting），允许一个对象（或一个组）与另一个组聚合的依据是最远距离对，即一个组吸纳一个新成员的依据是比较该组内所有成员与备选新成员的最远距离（Lance and Williams, 1967）。在该模式下，两个组所有成员之间的距离都必须全部计算，然后再比较。

继续以“坐在酒吧里的每个人”举例，完全连接聚类可以理解为我加入了一个具有我不喜欢的人的桌子中，但是在我加入后该桌的关系也并不糟糕，因为其它每个桌子上都有我更不喜欢的人。

完全连接聚合聚类中，两个聚类群之间的距离定义为每个簇中两个对象之间的最长距离。

在完全连接聚类中，由于合并后的组与其它组的距离是原来距离的最大者，加大了并组与其它组间的距离，因此其是空间扩张的，对进一步合并起着阻碍作用。对于已知组，吸收一个新成员要求该组内所有成员都参与比较，组越大就越难再加入新成员。

因此与单连接聚类方法相比，完全连接聚类更倾向产生很多小的分离的组，并且分类组之间的差异也比较明显，其更适合寻找和识别数据的间断分布。

平均聚合聚类

单连接聚合聚类或完全连接聚合聚类均属于基于连接的聚类方法。

与上述方法相比，平均聚合聚类（average agglomerative clustering）是一类基于对象间平均相异性或聚类簇形心（centroid）的聚类方法，目的将对象加入到其平均距离最低的组。在“坐在酒吧里的每个人”例子中，平均聚合聚类模式可以理解为我综合衡量所有的桌中，对每张桌子上的人的平均好感度，然后加入具有最高平均好感度的那桌。

平均聚合聚类中，两个聚类群之间的距离定义为一个簇中每个对象到另一个簇中每个对象之间的平均距离。

平均聚合聚类包含4种主要类型，区别在于组的位置计算方式（中值或质心）和当计算融合距离时是否用每组包含的对象数量作为权重（Sneath and Sokal, 1973）。

下图形象展示了单连接聚类（a）、完全连接聚类（b）、中值聚类（c）和质心聚类（d）的区别（Greig-Smith 1983）。

对于中值法，一个对象加入一个组的依据是这个对象与该组每个成员之间的平均距离，两个组聚合的依据是一个组内所有成员与另一组内成员之间所有对象对的平均距离。中值法保持了合并后的组与其它对象（或组）间的距离，即空间保持的，因此通常比基于连接的聚类方法（单连接聚类或完全连接聚类）效果更好；缺点是非单调性（Gauch, 1982）。

质心法将两个组间的距离定义为两个组质心间的距离，组质心指组内所有对象间距离的均值。在单个对象或只有两个对象组成的组中，质心法与中值法完全相同，但当对象数大于3时，二者将不同。质心法的缺点在于，当两个组所含对象数目相差悬殊时，新合并的组将更接近较大的一组，使较小组的特征消失，而影响聚类结果（Greig-Smith, 1983）。

此外质心法还可能会导致聚类树翻转的问题，详见下文R语言运行部分。

Ward最小方差聚类

Ward最小方差聚类（Ward’s minimum variance clustering）是一种基于最小二乘法线性模型准则的聚类方法，分组的依据是使组内误差平方和（每个对象到其组质心的距离的平方和）最小化（Ward,1963）。

Ward聚类方法试图寻找一种解决方案，将两组对象（簇）合并在一起。一开始，每个对象形成一个独立的簇，误差平方和为零。在随后的每一步中，该方法将对象与对象、对象与簇或簇与簇进行合并，此时误差平方和增大，并最终找到具有最小误差平方和的组合，使总体误差平方和的增幅最低。

图示两种不同的Ward聚类计算过程。对于每个聚类簇k，计算单个对象（y_ij^(k)）到该簇质心（m_j^(k)）的距离的平方和（左式），或者取簇内所有对象对之间的平均平方距离（D_hi²，右式）。

由于对象到质心的距离是在欧几里得空间中计算的，因此在该方法中只能使用欧式距离测度（如欧几里得距离）。对于其它非欧式属性的距离测度的使用，可能会使聚类结果违背最小化组内方差的基本原则，使聚类结果不可靠。尽管如此，可以通过转化的方式实现它们的应用，例如Bray-curtis距离，可在平方根转化后获得欧式几何特征。

灵活聚类

Lance和Williams（1966，1967）提出了一种能够涵盖上述聚类方法的模型，称为灵活聚类（flexible clustering)，该方法中没有固定的聚合策略，通过指定不同的参数项（α_h、α_i、β、γ）实现不同的聚类效果。

式中，g、h、i代表了三个聚类群（组），S为组间相似性，n为各聚类群内的对象数量。

若使用距离测度而非相似性，则可使用下式计算，D代表了组间距离。

对于上述所提到的几种聚类方法，在灵活聚类中均可视为特殊形式。

可知，其中最主要的值是β，因此大多数情况下，灵活聚类也称为灵活beta连接聚类（flexible beta linkage clustering）。

β称为聚集强度系数（cluster-intersing coefficient），取值范围[-1,1)，当β=-1时，具有最低水平的聚类程度，各对象单独成簇；当β越接近于1，越具更高水平的聚类程度，或者说更多的对象倾向于聚为一簇。

图示基于Bray-curtis距离的灵活聚类，3种不同的结果分别对应于β=-0.95（左）、β=-0.25（中）和β=0.95（右）时的情形。

层次聚合分类在R语言中的计算

接下来展示上述聚类方法在R语言中的实现。

例如期望通过聚类识别不同组织的基因表达相似性，不同环境群落物种组成相似性等。

下文的示例数据，R代码等可见百度盘：

https://pan.baidu.com/s/1tM2iSX6-r0aWxZmkXD89cA

计算距离

计算样本（或分组）间的关联矩阵（或距离矩阵）是进行层次聚合分类分析的第一步，因此选择恰当的关联系数非常重要。

##示例 1，基因表达数据
dat_gene <- read.delim('gene_express.txt', row.names = 1)
dat_gene <- t(dat_gene)
 
#常用欧几里得距离表示各样本中的基因表达相异性，可由 vegan 包 vegdist() 计算距离测度
dis_euc <- vegan::vegdist(dat_gene, method = 'euclidean')
 
##示例 2，物种丰度数据
dat_otu <- read.delim('otu_table.txt', row.names = 1)
dat_otu <- t(dat_otu)
 
#群落数据的分析中，Bray-curtis 距离是常用的距离测度类型
dis_bray <- vegan::vegdist(dat_otu, method = 'bray')
 
#不过考虑到某些聚类方法只能以欧式距离作为输入
#因此可以通过平方根转化的形式，将 Bray-curtis 距离转化为具有欧式几何属性的距离测度
ade4::is.euclid(dis_bray)
dis_bray_euc <- dis_bray^0.5
ade4::is.euclid(dis_bray_euc)
 
##此外，如果提供了现有的已经计算好的某种距离矩阵，也可直接读取
#例如，外部的 Bray-curtis 距离矩阵
dis_bray <- as.dist(read.delim('bray_distance.txt', row.names = 1))
 
#平方根转化，使其具有欧式几何属性
dis_bray_euc <- dis_bray^0.5

运行层次聚合分类

以下展示上文提到的主要的层次聚合分类方法在R中的实现过程。

单连接聚合聚类&完全连接聚合聚类

可使用hclust()计算。

##基于连接的层次聚合分类，详情 ?hclust
#以物种丰度的 Bray-curtis 距离为例
 
#单连接聚合聚类，单联动
clust_single <- hclust(dis_bray, method = 'single')
summary(clust_single)
 
#完全连接聚合聚类，全联动
clust_complete <- hclust(dis_bray, method = 'complete')
summary(clust_complete)
 
#聚类树，默认效果
par(mfrow = c(1, 2))
plot(clust_single, main = 'single')
plot(clust_complete, main = 'complete')

平均聚合聚类

可使用hclust()计算。

##平均聚合聚类，详情 ?hclust
#以物种丰度的 Bray-curtis 距离为例
 
#使用中值的非权重成对组法，UPGMA
clust_average <- hclust(dis_bray, method = 'average')
summary(clust_average)
 
#使用质心的非权重成对组法，UPGMC
clust_centroid <- hclust(dis_bray, method = 'centroid')
summary(clust_centroid)
 
#使用中值的权重成对组法，WPGMA
clust_mcquitty <- hclust(dis_bray, method = 'mcquitty')
summary(clust_mcquitty)
 
#使用质心的权重成对组法，WPGMC
clust_median <- hclust(dis_bray, method = 'median')
summary(clust_median)
 
#聚类树，默认效果
par(mfrow = c(2, 2))
plot(clust_average, main = 'UPGMA')
plot(clust_centroid, main = 'UPGMC')
plot(clust_mcquitty, main = 'WPGMA')
plot(clust_median, main = 'WPGMC')

质心法的聚类方法有时会导致聚类树翻转（reversal），如上UPGMC或WPGMC所示，聚类树不再形成连续的嵌套分区，使分类结果难以解读。

对于该现象产生的成因，可简要描述如下。存在3个对象X₁、X₂和X₃（a），X₁和X₂首先归为一簇，因为它们之间具有最短的距离；尽管X₃-X₁或X₃-X₂的距离大于X₁-X₂的距离，但是X₃到X₁和X₂质心的距离小于X₁-X₂的距离，此时即产生聚类树翻转的情形（b）。

Legendre和Legendre（1998）在“Numerical Ecology”第342页前后对该现象作了详细描述，并建议用多分法（polychotomy）代替传统的二分法（dichotomy）解读这种图。

Ward最小方差聚类

可使用hclust()计算。

##Ward 最小方差聚类（包含两种不同的 Ward 聚类算法，详情 ?hclust）
#该方法中只能使用欧式几何属性的距离测度，因此更换为平方根转化后的 Bray-curtis 距离测度

#ward.D
clust_ward <- hclust(dis_bray_euc, method = 'ward.D')
summary(clust_ward)

#ward.D2
clust_ward2 <- hclust(dis_bray_euc, method = 'ward.D2')
summary(clust_ward2)

#聚类树，默认效果
par(mfrow = c(1, 2))
plot(clust_ward, main = 'ward.D')
plot(clust_ward2, main = 'ward.D2')

灵活聚类

可使用cluster包的agnes()计算。

##灵活聚类，详情 ?agnes
#也可在该函数中通过指定 method 执行上述的多种常规聚类，如 method='average' 即为 UPGMA
#当 method='flexible' 时，par.method 参数中的 4 个值分别对应了 αh、αi、γ、β 的设定值，例如
clust_flex <- cluster::agnes(dis_bray_euc, method = 'flexible', par.method = c(0.5, 0.5, 0, 0.5))
summary(clust_flex)

plot(clust_flex)

聚类评估

层次聚类是探索性分析，而非统计检验的过程。

对于如何比较和评估层次聚类结果，寻找可解读的聚类簇等，放在下篇介绍。

聚类树可视化调整

以下再简单展示一些聚类树的调整示例。

同样地，聚类树更详细的可视化方法将在后文使用一篇单独展示。

#plot() 绘制聚类树参数简要说明，以上述 UPGMA 结果为例
#默认
plot(clust_average)

#hang 参数（默认 0.1）可用于调整标签与树高的关系，当hang < 0时，强制标签从 0 开始并对齐，如下示例
par(mfrow = c(1, 2))
plot(clust_average, hang = 0.2)
plot(clust_average, hang = -1)
#main、sub、xlab、ylab 用于修改图片标题或坐标轴标题，如下示例
plot(clust_average, main = 'UPGMA\n(Bray-curtis distance)', sub = '', xlab = 'Sample', ylab = 'Height')
#col 可定义颜色，cex 可调整标签字体大小，lwd 可调整线宽，lty 可调整线型，如下示例
plot(clust_average, col = 'red', cex = 1, lwd = 1, lty = 2)
#设定选择观测的聚类群数量为 4
plot(clust_average, main = 'UPGMA\nBray-curtis distance', sub = '', xlab = 'Sample', ylab = 'Height')
rect.hclust(clust_average, k = 4, border = c('red', 'blue', 'green3', 'orange'))
#vegan 包函数 reorder() 可以实现对聚类树重排，详情 ?reorder
#尽可能使聚类树内对象的排列顺序和原始相异矩阵内对象的排列顺序一致
#不过效果似乎一般
library(vegan)

plot(reorder(clust_average, dis_bray), main = 'UPGMA\nBray-curtis distance', sub = 'after reorder', xlab = 'Sample', ylab = 'Height')
rect.hclust(reorder(clust_average, dis_bray), k = 4, border = c('red', 'blue', 'green3', 'orange'))

参考资料

https://www.davidzeleny.net/anadat-r/doku.php/en:hier-agglom

https://www.saedsayad.com/clustering_hierarchical.htm

DanielBorcard, FranoisGillet, PierreLegendre, et al. 2014. 数量生态学:R语言的应用（赖江山 译）. 高等教育出版社.

Gauch, H. G. Jr. 1982. Multivariate analysis in community ecology. Cambridge University Press, Cambridge. x + 298 pp.

Greig-Smith P. 1983. Quantitative Plant Ecology (3rd ed.). Oxford:Blackwell Scientific Publications.

Lance, G. N. & W. T. Williams. 1966. A generalized sorting strategy for computer classifications. Nature (Lond.) 212: 218.

Lance, G. N. & W. T. Williams. 1967. A general theory of classificatory sorting strategies. I. Hierarchical systems. Computer J. 9: 373-380.

Legendre P, Legendre L. 1998. Numerical Ecology. Second English edition. Developments in Environmental Modelling, 20, Elsevier.

Lukaszewicz, J. 1951. Sur la liaison et la division des points d'un ensemble fini. Colloq. Math. 2:282-285.

Sneath, P. H. A. & R. R. Sokal. 1973. Numerical taxonomy – The principles and practice of numerical classification. W. H. Freeman, San Francisco. xv + 573 pp.

Ward, J. H. 1963. Hierarchical grouping to optimize an objective function. J. Amer. Statist. Assoc. 58: 236-244.

友情链接