R语言层次聚类算法及可视化

Posted 2021-04-08 PLANTOMIX

点击上方 蓝字  关注我们

---

关于层次聚类

层次聚类（Hierarchical Clustering）是一种和 K- 均值聚类类似的聚类算法。和 K-均值聚类不同的是，层次聚类不需要预先确定一个层次值；除此之外 ，层次聚类比 K-均值聚类更有优势的一点在于层次聚类的结果可以用树状图（dendrogram）来展示。

准备工作

本节需要用到的 R 包：

# 加载包
library(dplyr)       # 数据预处理
library(ggplot2)     # 数据可视化
library(cluster)     # 实现聚类算法
library(factoextra)  # 聚类结果可视化

构建测试数据集（数据来自 R 包 AmesHousing）：

# 数据预处理
ames_scale <- AmesHousing::make_ames() %>%
        select_if(is.numeric) %>%  # 选择数据类型为数字
        select(-Sale_Price) %>%    # 去除特定的数据
        mutate_all(as.double) %>%  # 强制转换类型
        scale()                    # 中心化标准化数据

层次聚类算法

层次聚类算法可以分成两种：

1. 聚集聚类（Agglomerative clustering），通常叫 AGNES，理解成从个体到整体的聚类方式；
2. 分裂层次聚类（Divisive hierarchical clustering），通常叫 DIANA，理解成从整体到个体的扩散。

从图上可以看出两者是相反的过程。

和 K- 均值类似，也是先计算观察变量之间的距离（如欧氏距离、曼哈顿距离等）。其中欧氏距离是用得最多的一种。在层次聚类中一个重要的问题是：如何衡量两个变量（clusters）之间的相异性？常见的方法有如下几种：

1. Maximum or complete linkage clustering：计算 cluster1 和 cluster2 中两两元素之间的相异性，选择最大值作为距离 cluster1 和 cluster2 的相异值，这种算法更倾向于产生更加紧凑的 cluster；
2. Minimum or single linkage clustering：分别计算 cluster1 和 cluster2 内的元素的相异性，求均值作为 cluster1 和 cluster2 的相异值，这种算法通常生成一个更长的树状图；
3. Centroid linkage clustering：计算 cluster1 和 cluster2 的 “重心” 的相异值；
4. Ward’s minimum variance method：最小化集群内的总方差，每一步都将簇间距离最小的一对簇进行合并。

还有多其他的聚类算法，但是上面这 4 种是最常见使用最广泛的。

聚类算法在 R 语言中的实现

Agglomerative hierarchical clustering

实现 Agglomerative hierarchical clustering 的函数是 hclust。首先需要用函数 dist 计算距离，然后选择不同的方法进行聚类分析 (i.e. "complete", "average", "single", or "ward.D").。

# 设置随机数种子
set.seed(123)

# 计算距离
d <- dist(ames_scale, method = "euclidean")

# 聚类
par(mfrow = c(2,2))

# for 循环批量计算聚类并绘图
for (i in c("complete","average","single","ward.D")) {
        hc1 <- hclust(d, method = i)
        
        plot(hc1, cex = 0.6, hang = -1,main = i)
}

由于样本有点多，图并不是很好看，但是可以看出差异。

除了 hclust 这个函数，还可以使用 agnes 这个函数，这个函数可以计算出聚类强度值（agglomerative coefficient ，AC）。AC 越接近 1 说明聚类越平衡。需要注意的是，AC 会随着观测值的增加而变大，所以 AC 值不能用来比较不同大小的数据集。我们来看看上面 4 中聚类算法的 AC 值分别是多少：

# 设置随机数种子
set.seed(123)
agnes.ac = NULL
# 批量计算AC值
# 方法选择
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

# 构建函数提取AC
ac <- function(x) {
        agnes(ames_scale, method = x)$ac
}

# 查看AC值
purrr::map_dbl(m, ac)

# 结果
  average    single  complete      ward 
0.9139303 0.8712890 0.9267750 0.9766577 

Divisive hierarchical clustering

这种聚类算法的函数是 diana，没有方法可以进行选择。

# 计算 divisive hierarchical clustering
hc4 <- diana(ames_scale)

# Divise coefficient
hc4$dc

## [1] 0.9191094

确定最佳聚类

虽然层次聚类能够展示所有样本之间的关系，但我们还是可以像 K- 均值聚类那样看看聚成多少类比较合适。

# Plot cluster results
p1 <- fviz_nbclust(ames_scale, FUN = hcut, method = "wss", 
                   k.max = 10) +
  ggtitle("(A) Elbow method")
p2 <- fviz_nbclust(ames_scale, FUN = hcut, method = "silhouette", 
                   k.max = 10) +
  ggtitle("(B) Silhouette method")
p3 <- fviz_nbclust(ames_scale, FUN = hcut, method = "gap_stat", 
                   k.max = 10) +
  ggtitle("(C) Gap statistic")

# Display plots side by side
gridExtra::grid.arrange(p1, p2, p3, nrow = 1)

树状图基本可视化

# 树状图可视化
# 选择前100个样本作为示例
d <- dist(ames_scale[1:50,], method = "euclidean")
hc <- hclust(d, method = 'average')

# 直接绘图
plot(hc)

让所有样品名称排在同样的高度：

plot(hc, hang = -1, cex = 0.8)

plot.dendrogram() 函数

使用 plot.dendrogram() 函数对树状图进行展示需要先将聚类的结果转换成 dendrogram：

# 格式转换
hcd <- as.dendrogram(hc)

绘制两种不同形状的图：

# 两种不同形状
par(mfrow = c(1,2))

# Default plot
plot(hcd, type = "rectangle", 
     ylab = "Height",main = 'rectangle')
# Triangle plot
plot(hcd, type = "triangle", 
     ylab = "Height",main = 'triangle')

dev.off()

使用 ggtree 进行可视化

Y 叔的 ggtree 可以说是树状图（进化树）可视化的神级包了：

# ggtree
hc.ggtree = as.phylo(hc)

ggtree(hc.ggtree, layout = 'fan')+
        geom_tiplab()+
        geom_point(aes(shape=isTip, color=isTip), size=2)

iTOL 进行可视化

iTOL（https://itol.embl.de/）是一个树状图可视化的神器，聚类的树状图也可以用 iTOL 进行可视化。

# 先将聚类结果转换成可用的数据
write.tree(as.phylo(hc), file = 'test.nwk')

然后将导出的文件导入到 iTOL 网站就可以了，写点配置文件就可以了。配置文件的写法可以参考：

http://ynaulx.cn/2020/01/09/iTOL%E4%BF%AE%E9%A5%B0%E8%BF%9B%E5%8C%96%E6%A0%91/#more