聚類樹和PCA等排序圖的組合繪製

聚類分析和排序分析（降維分析）都是用於探索多元數據結構的常用方法，二者的結果也可以結合在一起通過一張圖呈現，本篇展示一些常見的示例。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, '%'))