聚类树和PCA等排序图的组合绘制

聚类分析和排序分析（降维分析）都是用于探索多元数据结构的常用方法，二者的结果也可以结合在一起通过一张图呈现，本篇展示一些常见的示例。

示例文件、R脚本等的百度盘链接：

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

示例数据包含15个样本（对象），20个变量，下文对它执行聚类和降维，并作图展示。

排序图+聚类树

聚类树一般用以表示层次聚类的结果，排序图中也可添加聚类树展示对象间的层次结构。 

在前文简介主成分分析（PCA）的方法中，展示了一例在PCA排序图中添加聚类树的方法。

#vegan 包中的 PCA 方法
library(vegan)

dat <- read.delim('data.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
pca <- rda(dat, scale = TRUE)

pca1 <- round(100*pca$CA$eig[1]/sum(pca$CA$eig), 2)
pca2 <- round(100*pca$CA$eig[2]/sum(pca$CA$eig), 2)

#I 型标尺的 PCA 排序图
p <- ordiplot(pca, dis = 'site', type = 'n', choices = c(1, 2), scaling = 1,
xlab = paste('PCA1:', pca1, '%'), ylab = paste('PCA2:', pca2, '%'))
points(pca, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, cex = 1.2,
bg = c(rep('red', 5), rep('orange', 5), rep('green3', 5)), col = NA)

#添加聚类，如 UPGMA 聚类
tree <- hclust(dist(scale(dat), method = 'euclidean'), method = 'average')
plot(tree)

ordicluster(p, tree, col = 'gray')

排序图绘制以二维散点图，并将聚类树投影到二维平面中。

如果觉得二维平面难以观测层次结构，FactoMineR包提供了聚类树的三维可视化方案，同样以PCA为例展示。

#FactoMineR 包的三维聚类树
#参考链接：http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/117-hcpc-hierarchical-clustering-on-principal-components-essentials/
library(FactoMineR)

#PCA
res.pca <- PCA(dat, ncp = 3, scale.unit = TRUE, graph = FALSE)

#添加层次聚类，如 UPGMA 聚类
res.hcpc <- HCPC(res.pca, metric = 'euclidean', method = 'average', graph = FALSE)

#三维效果图
plot(res.hcpc, choice = '3D.map')

相比上述二维平面，层次结构更加清晰明了。

排序图+分类边界

非层次聚类方法，如k均值划分聚类、围绕中心点划分聚类、模糊c均值聚类等，一般无法以聚类树展示其结构。这种情况下，通常在排序图中添加分组椭圆或多边形，将同一类的对象“圈起来”，以反映分类信息。

当然上述的层次聚类结果也可以这样来表示。

以下是一些示例。 

例如前文模糊c均值聚类中，提到可以这样在排序图中标记分组。

#以 PCoA 为例降维，并在排序图中标注模糊 c 均值聚类结果

#计算距离，以欧几里得距离为例
dis_euc <- dist(scale(dat), method = 'euclidean')

#PCoA，事实上以欧几里得距离的 PCoA 等同于以原始数据的 PCA
pcoa <- cmdscale(dis_euc, k = (nrow(dat) - 1), eig = TRUE)
eig_prop <- 100*pcoa$eig/sum(pcoa$eig)
pcoa_site <- pcoa$point[ ,1:2]
plot(pcoa_site,
xlab = paste('PCoA axis1:', round(eig_prop[1], 2), '%'),
ylab = paste('PCoA axis2:', round(eig_prop[2], 2), '%'))

#模糊 c 均值聚类，聚为 3 类
library(cluster)

set.seed(123)
cm.fanny <- fanny(dis_euc, k = 3, diss = TRUE, maxit = 100)

#标注各对象最接近的聚类簇
for (i in 1:3) {
gg <- pcoa_site[cm.fanny$clustering == i, ]
hpts <- chull(gg)
hpts <- c(hpts, hpts[1])
lines(gg[hpts, ], col = 1+1)
}

#星图展示了各对象的成员值
stars(cm.fanny$membership, location = pcoa_site, draw.segments = TRUE,
add = TRUE, scale = FALSE, col.segments = 2:4, len = 0.5)

在这种软聚类方法中，星图展示了对象属于各聚类簇的成员值，面积越大代表成员值越高，即该对象在该组中的归属程度越大，并将各对象所属的最佳聚类簇用红线连接。数据集整体特征一目了然。

对于硬聚类的方法，由于每个对象只有一个确定的分类，因此这种分组边界的绘制就很简单，直接通过椭圆或多边形“圈起来”即可。

前文简介排序图绘制中，提到了一些标记对象分组的方法，同样也可用于标记已识别的分类结果。

#ggplot2 中的一些方法
library(ggplot2)

#以上文 PCA 的排序结果为例展示
#首先提取对象排序坐标，以 I 型标尺为例，以前两轴为例
pca_site <- data.frame(scores(pca, choices = 1:2, scaling = 1, display = 'site'))

#k-means 聚类，聚为 3 类
set.seed(123)
dat_kmeans <- kmeans(scale(dat), centers = 3, nstart = 25)
dat_kmeans$cluster
pca_site$cluster <- as.character(dat_kmeans$cluster)

#ggplot2 绘制二维散点图
p <- ggplot(data = pca_site, aes(x = PC1, y = PC2)) +
geom_point(aes(color = cluster)) +
scale_color_manual(values = c('red', 'blue', 'green3')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) +
labs(x = paste('PCA1:', pca1, '%'), y = paste('PCA2:', pca2, '%'))

#添加置信椭圆，可用于表示对象分类
p + stat_ellipse(aes(color = cluster), level = 0.95, linetype = 2, show.legend = FALSE)

p + stat_ellipse(aes(fill = cluster), geom = 'polygon', level = 0.95, alpha = 0.1, show.legend = FALSE) +
scale_fill_manual(values = c('red', 'purple', 'green'))

#将同类别对象连接在一起
source('geom_enterotype.r') #可获取自 https://pan.baidu.com/s/1KUA5owf4y13y1RJbff-l7g
p + geom_enterotype(aes(fill = cluster, color = cluster, label = cluster), show.legend = FALSE) +
scale_fill_manual(values = c('#ffa6a9', '#e8b4fd', '#c7ffc4'))

#多边形连接同类别对象的样式
library(plyr)

cluster_border <- ddply(pca_site, 'cluster', function(df) df[chull(df[[1]], df[[2]]), ])

p + geom_polygon(data = cluster_border, aes(fill = cluster), alpha = 0.1, show.legend = FALSE) +
scale_fill_manual(values = c('red', 'purple', 'green'))

上述提到的FactoMineR包也能用于绘制这类统计图，以上述添加层次聚类的PCA继续。

#上述已经展示了使用 FactoMineR 包可视化 PCA + 聚类树
#将聚类树更改为分组边界也是常见的可视化方案
library(factoextra)
library(FactoMineR)

fviz_dend(res.hcpc,
cex = 0.7, # Label size
palette = 'jco', # Color palette see ?ggpubr::ggpar
rect = TRUE, rect_fill = TRUE, # Add rectangle around groups
rect_border = 'jco', # Rectangle color
labels_track_height = 0.8) # Augment the room for labels

fviz_cluster(res.hcpc,
repel = TRUE, # Avoid label overlapping
show.clust.cent = TRUE, # Show cluster centers
palette = 'jco', # Color palette see ?ggpubr::ggpar
ggtheme = theme_minimal(),
main = 'Factor map')

还有三维的样式，以下展示一例，更多方法可参考三维排序图。

#以上文 PCA 的排序结果为例展示
#首先提取对象排序坐标，以 I 型标尺为例，以前 3 轴为例
pca_site <- data.frame(scores(pca, choices = 1:3, scaling = 1, display = 'site'))

pca1 <- round(100*pca$CA$eig[1]/sum(pca$CA$eig), 2)
pca2 <- round(100*pca$CA$eig[2]/sum(pca$CA$eig), 2)
pca3 <- round(100*pca$CA$eig[3]/sum(pca$CA$eig), 2)

#k-means 聚类，聚为 3 类
set.seed(123)
dat_kmeans <- kmeans(scale(dat), centers = 3, nstart = 25)
dat_kmeans$cluster
pca_site$cluster <- as.character(dat_kmeans$cluster)

#car+rgl 包，交互式三维图
library(car)
library(rgl)

scatter3d(pca_site$PC1, pca_site$PC2, pca_site$PC3, groups = as.factor(pca_site$cluster), surface = FALSE, ellipsoid = TRUE,
xlab = paste('PCA1:', pca1, '%'), ylab = paste('PCA2:', pca2, '%'), zlab = paste('PCA3:', pca3, '%'))

友情链接