数量生态学：R语言的应用—第四章聚类分析2-解读和比较层次聚类结

作者: fafu生信小蘑菇 | 来源:发表于2021-05-07 23:36 被阅读0次

数量生态学：R语言的应用—第四章聚类分析2

今天来给大家介绍一下如何解读和比较层次聚类结果

由于今天内容比较多，所以我们先来看看思维导图，看看大致有什么

解读和比较层次聚类结果

这里强调一下，聚类分析是一种探索性分析，而非统计检验。影响聚类结果的因素包括聚类方法本身和用于聚类分析的关联系数的选择。因此，选择与分析目标一致的方法非常重要。由函数hclust（）产生的对象包含很多聚类分析结果信息，也是绘制聚类树的依据。

在R里键入summary（聚类结果对象名）可以获得聚类结果相关数量信息列表。
summary（）的信息有助于解读和比较聚类的结果。

1. 同表型相关

一个聚类树内两个对象之间的同表型距离是两个对象在同一组分类水平内的距离。
任意两个对象，在聚类树上从一个对象向上走，到达与另外一个对象交汇节点向下走，势必会到达第二个对象：交汇节点所在的层次水平即是两个对象同表型距离。
同表型矩阵是所有对象对的同表型距离矩阵。
同表型相关是原始的距离矩阵和同表型距离矩阵之间的Pearson 相关系数。
具有最高的同表型相关系数的聚类方法可视为原始矩阵最好的聚类模型。

但是，同表型矩阵由原始距离矩阵推导而来，两组距离矩阵不独立，因此不能检验同表型相关的显著性。此外，同表型相关强烈依赖于数据的聚类方法的选择。

我们可以通过stats程序包内cophenetic函数计算同表型矩阵和同表型相关。

#同表型相关

# 单连接聚类同表型相关
spe.ch.single.coph <- cophenetic(spe.ch.single)
cor(spe.ch, spe.ch.single.coph)

# 完全连接聚类同表型相关
spe.ch.comp.coph <- cophenetic(spe.ch.complete)
cor(spe.ch, spe.ch.comp.coph)

# 平均聚类同表型相关
spe.ch.UPGMA.coph <- cophenetic(spe.ch.UPGMA)
cor(spe.ch, spe.ch.UPGMA.coph)

# Ward聚类同表型相关
spe.ch.ward.coph <- cophenetic(spe.ch.ward)
cor(spe.ch, spe.ch.ward.coph)

#哪个聚类树保持与原始的弦距离矩阵最接近的关系？同表型相关也可以用spearman秩相关或Kendall秩相关表示
cor(spe.ch, spe.ch.ward.coph, method="spearman")

同表型相关

为了更好的展示一个距离矩阵与通过不同聚类方法得到的同表型矩阵之间的相关性，可以绘制原始距离对阵同表型距离的Shepard图。

# Shepard图
par(mfrow=c(2,2))
plot(spe.ch, 
     spe.ch.single.coph, 
     xlab="弦距离", 
     ylab="同表型距离", 
     asp=1, 
     xlim=c(0,sqrt(2)),
     ylim=c(0,sqrt(2)), 
     main=c("单连接",paste("同表型相关系数=",
                        round(cor(spe.ch, spe.ch.single.coph),3))))
abline(0,1)
lines(lowess(spe.ch, spe.ch.single.coph), col="red") #添加LOWESS趋势平滑线

plot(spe.ch, 
     spe.ch.comp.coph, 
     xlab="弦距离", 
     ylab="同表型距离", 
     asp=1, 
     xlim=c(0,sqrt(2)), 
     ylim=c(0,sqrt(2)),
     main=c("完全连接", paste("同表型相关系数=",
                          round(cor(spe.ch, spe.ch.comp.coph),3))))
abline(0,1)
lines(lowess(spe.ch, spe.ch.comp.coph), col="red")#添加LOWESS趋势平滑线

plot(spe.ch, 
     spe.ch.UPGMA.coph, 
     xlab="弦距离", 
     ylab="同表型距离", 
     asp=1, 
     xlim=c(0,sqrt(2)), 
     ylim=c(0,sqrt(2)), 
     main=c("UPGMA", paste("同表型相关系数=",
                           round(cor(spe.ch, spe.ch.UPGMA.coph),3))))
abline(0,1)
lines(lowess(spe.ch, spe.ch.UPGMA.coph), col="red")#添加LOWESS趋势平滑线

plot(spe.ch, 
     spe.ch.ward.coph, 
     xlab="弦距离", 
     ylab="同表型距离", 
     asp=1, 
     xlim=c(0,sqrt(2)), 
     ylim=c(0,max(spe.ch.ward$height)),
     main=c("Ward聚类", paste("同表型相关系数=",
                            round(cor(spe.ch, spe.ch.ward.coph),3))))
abline(0,1)
lines(lowess(spe.ch, spe.ch.ward.coph), col="red")#添加LOWESS趋势平滑线

#哪个聚类方法产生的同表型距离（纵坐标）跟原始的距离（横坐标）线性关系最好？

比较弦距离（物种数据）与四种聚类方法同表型距离的SHEPARD图

另一个比较聚类结果的指标是Gower距离，它等于原始距离与同表型距离之间差值的平方和。

具有最小Gower距离的聚类方法也可视为原始矩阵最好的聚类模型。

但同表型相关系数和Gower距离分析结果并不总是一致。

# Gower距离
(gow.dist.single <- sum((spe.ch-spe.ch.single.coph)^2))
(gow.dist.comp <- sum((spe.ch-spe.ch.comp.coph)^2))
(gow.dist.UPGMA <- sum((spe.ch-spe.ch.UPGMA.coph)^2))
(gow.dist.ward <- sum((spe.ch-spe.ch.ward.coph)^2))
#这些代码加括号可以直接在产生对象后立刻在屏幕上展示出来

Gower距离

2. 寻找可解读的聚类簇

为了解读和比较聚类的结果，通常需要寻找可解读的聚类簇，这意味着需要决定聚类树裁剪到哪个水平。在某些情况下，在裁剪后的聚类树上添加其他信息，或将裁剪后的聚类结果作为附加信息添加到其他分析结果中，都OK。

2.1 融合水平值图

聚类树的融合水平值是聚类树中两个分支融合处的相异性的数值。绘制融合水平值的变化图有利于定义裁剪的水平。

# 融合水平值图
par(mfrow=c(2,2))
summary(spe.ch.complete)  # 总结聚类分析的结果
# 绘制完全连接聚类融合水平值图
plot(
  spe.ch.complete$height, 
  nrow(spe):2, 
  type="S", 
  main="融合水平-弦距离-完全连接", 
  ylab="k (聚类簇数量)",
  xlab="h (节点高度)", 
  col="grey")
# 函数text()的使用，在图上直接标注聚类簇的数字标识
text(spe.ch.complete$height, 
     nrow(spe):2, 
     nrow(spe):2, 
     col="red", 
     cex=0.8)

# 绘制UPGMA聚类融合水平值图
plot(spe.ch.UPGMA$height, 
     nrow(spe):2, 
     type="S", 
     main="融合水平-弦距离-UPGMA", 
     ylab="k （聚类簇数量）", 
     xlab="h (节点高度）", 
     col="grey")
text(spe.ch.UPGMA$height, 
     nrow(spe):2, 
     nrow(spe):2, 
     col="red", 
     cex=0.8)

#绘制Ward聚类融合水平值图
plot(spe.ch.ward$height, 
     nrow(spe):2, 
     type="S", 
     main="融合水平-弦距离-Ward", 
     ylab="k （聚类簇数量）", 
     xlab="h (节点高度）",
     col="grey")
text(spe.ch.ward$height, 
     nrow(spe):2, 
     nrow(spe):2, 
     col="red", 
     cex=0.8)

#绘制beta模糊聚类融合水平值图(beta=-0.25)
plot(spe.ch.beta2$height, 
     nrow(spe):2, 
     type="S", 
     main="融合水平-弦距离-Beta灵活聚类", 
     ylab="k （聚类簇数量）", 
     xlab="h （节点高度）",
     col="grey")
text(spe.ch.beta2$height, 
     nrow(spe):2, 
     nrow(spe):2, 
     col="red", 
     cex=0.8)

4种聚类树的融合水平值图

每幅图从右向左看，距离比较长的水平线是被建议可能分组的水平，例如，完全连接聚类中9组和Ward聚类中的4组。

上面4个图看起来差异很大。记住，这些解决方案中任何一个都不是绝对正确，每个方案都可以为数据分组提供独特见解

使用cutree()函数设定分类组的数量，并使用列联表比较分类差异。

# 设定聚类组的数量
k <- 4  
# 根据上面4个融合水平值图，可以观察到分4组水平在所有图里有小的跳跃
# 裁剪聚类树
spech.single.g <- cutree(spe.ch.single, k=k)
spech.complete.g <- cutree(spe.ch.complete, k=k)
spech.UPGMA.g <- cutree(spe.ch.UPGMA, k=k)
spech.ward.g <- cutree(spe.ch.ward, k=k)
spech.beta.g <- cutree(spe.ch.beat2, k=k)
# 通过列联表比较分类结果
# 单连接vs完全连接
table(spech.single.g, spech.complete.g)

# 单连接vs UPGMA
table(spech.single.g, spech.UPGMA.g)

#单连接vs Ward
table(spech.single.g, spech.ward.g)

# 完全连接vs UPGMA
table(spech.complete.g, spech.UPGMA.g)

# 完全连接vs Ward
table(spech.complete.g, spech.ward.g)

# UPGMA vs Ward
table(spech.UPGMA.g, spech.ward.g)

#beta-灵活聚类 vs Ward
table(spech.beta.g, spech.ward.g)

image-20210506215707465

image-20210506215804039

如果两个聚类的结果完全一样，那么这个列联表每行和每列只有一个非零数字，其他应该为0。此处并没有出现这种情况。如何解读这些列联表呢？例如，单连接聚类第2组含有26个样方，这些样方在Ward聚类中被分散到4个组里。

2.2 比较两个聚类树以突出共同的子树

可以通过聚类树的对比来确定多个算法共同的分组。我们可以使用Dendextend包中tangelgram()函数。

#聚类结果原始“hclust”类对象必须首先转化为"dendrogram"
#类的对象
class(spe.ch.ward)  #[1]"hclust"
dend1 <- as.dendrogram(spe.ch.ward)
class(dend1)  #[1]"dendrogram"
dend2 <- as.dendrogram(spe.ch.complete)
dend12 <- dendlist(dend1,dend2)
tanglegram(
  untangle(dend12),
  sort=FALSE,
  common_subtrees_color_branches=TRUE,
  main_left="Ward聚类",
  main_right="完全连接"
)

两个聚类树的成对比较

在上面这个图中样方的位置尽可能的重新排列，为了让两棵聚类树更好的匹配。彩色的图用来突出共同的聚类簇，黑色的线所连着样方点表示在两棵树上不同的位置。

2.3 多尺度自助重采样

非约束聚类分析属于不以验证先验假设为目标的数据探索方法，然而，自然变异导致样本的变异，分类的结果很可能反映这种现象，因此评估分类的不确定性（或其对应的稳健性）是合理的。

常用的是自助抽样法，其原理是从数据集中随机抽取子集数据，然后进行这些子集数据的聚类分析。

Pvclust程序包提供绘制带有自助抽样P值聚类树，所有AU-p值都用红色字体标注，不太准确的 BP值以绿色标注。高AU值（例如大于0.95或更大）的聚类簇表示受到数据的高度支持。

注意函数Pvclust（）所使用的对象必须是我们通常的数据结构（即行是变量）的转置矩阵。

#计算聚类树中每个聚类簇(节点)的P值
spech.pv <-
  pvclust(t(spe.norm),
          method.hclust = "ward.D2",
          method.dist = "euc",
          parallel=TRUE)

# 绘制带P值的聚类树
par(mfrow = c(1, 1))
plot(spech.pv,
     main="带AU/BP值（%）的聚类树")


#此函数为给定对象绘制一个带有p值的树状图属于pvclust类。AU p值（默认为红色打印）
#是“近似无偏”p值的缩写，即通过多尺度引导重采样计算。BP值（打印默认为绿色）是“引导概率”值，
#它不如AU值作为p值准确。我们可以认为具有高AU值（例如95%）的团簇（边缘）是强烈的
#凸显高au值的聚类簇
pvrect(spech.pv, alpha = 0.95, pv = "au")
lines(spech.pv)
pvrect(spech.pv, alpha = 0.91, border = 4)

image-20210507153341443

image-20210507153444858

参数parallel=TRUE设置是通过平行运算提高运算速度，这时必须要求parallel程序包参与运作

用红色框括起来并用下划线表示的聚类簇具有显著的AU-p值（p≥0.95)，用蓝色框括起来的簇具有较小的P值（本例中P≥0.91）。此分析有助于突显最稳健的样方组。但请注意，对于相同的数据结果可能因每次运行的情况不同，因为这是一个基于随机的抽样过程。

一旦我们选择了一种聚类树并评估聚类树的稳健性，接下任务就是利用三种方法来确定合适的聚类簇的数量，分别是轮廓宽度值，聚类比较和诊断物种。

2.4 平均轮廓宽度图

轮廓宽度是指一个对象与所属聚类簇归属程度的测度，是一个对象与同组内其他对象的平均距离和该对象与最领近聚类簇内所有对象的平均距离比较。轮廓宽度范围从-1到1。

使用cluster程序包内的silhouette()函数绘制轮廓宽度图。该函数的帮助文件提供了轮廓宽度的正式定义，

简单的说轮廓宽度值越大，对象的聚类就越好，负值意味着该对象有可能被错分到当前聚类簇。

平均轮廓宽度值可以用于度量分组的质量，在每一种分类水平中，计算能够衡量对象与其所在组的关系紧密程度的轮廓宽度值，然后选择最大平均轮廓宽度值所对应的分组数量为最优的结果。

#平均轮廓宽度图
# 选择和重新命名聚类树（"hclust"对象）
hc <- spe.ch.ward
# hc <- spe.ch.beta2
# hc <- spe.ch.complete
par(mfrow = c(1, 1))
#平均轮廓宽度值(Rousseeuw质量指数)
Si <- numeric(nrow(spe))
for (k in 2:(nrow(spe) - 1))
{
  sil <- silhouette(cutree(hc, k = k), spe.ch)
  Si[k] <- summary(sil)$avg.width
}
k.best <- which.max(Si)
plot(
  1:nrow(spe),
  Si,
  type = "h",
  main = "轮廓宽度-最佳聚类簇数量",
  xlab = "k (组数)",
  ylab = "平均轮廓宽度"
)
axis(
  1,
  k.best,
  paste("最优", k.best, sep = "\n"),
  col = "red",
  font = 2,
  col.axis = "red"
)
points(k.best,
       max(Si),
       pch = 16,
       col = "red",
       cex = 1.5
)

image-20210507153547667

上图显示分组数K=2~29时所有轮廓宽度值图，在这个标准下，最大平均轮廓宽度值为最优分组结果，这里的最佳分组数量是2组，4组次之，再次6组，但从生态学角度分析，4组和6组可能是更有意义的。

2.5 距离矩阵和代表分组的二元矩阵的比较

这种方法是计算原始距离与代表不同分类水平的二元矩阵（从聚类数计算获得）之间的相关性，然后选择最高的相关系数所对应的分类水平作为最优的分组方案。

想法是选择两者之间距离相关性最高的分组水平，但统计学检验是不可能的，因为这两个矩阵之间并不独立（代表聚类树的二元矩阵实际上是从原始的距离矩阵推导而来的）。

为了计算代表不同分类水平下表示样方之间关系的二元相异矩阵，需要使用作者编写的grpdist（）函数计算向量化分类组。

# 通过矩阵相关选择最佳聚类簇数量

kt <- data.frame(k = 1:nrow(spe), r = 0)
for (i in 2:(nrow(spe) - 1)) 
{
  gr <- cutree(hc, i)
  distgr <- grpdist(gr)
  mt <- cor(spe.ch, distgr, method = "pearson")
  kt[i, 2] <- mt
}
k.best <- which.max(kt$r)
plot(
  kt$k,
  kt$r,
  type = "h",
  main = "矩阵相关-最优聚类簇数",
  xlab = "k (组数)",
  ylab = "Pearson相关系数"
)
axis(
  1,
  k.best,
  paste("最优", k.best, sep = "\n"),
  col = "red",
  font = 2,
  col.axis = "red"
)
points(k.best,
       max(kt$r),
       pch = 16,
       col = "red",
       cex = 1.5)

原始距离矩阵和代表不同分组情况的二元矩阵之间相关系数柱状图

这个矩阵图显示从3到6个组的时候可以获得两个矩阵相关系数的峰值

2.6 物种保真度分析

评估分组质量的另一个内部标准是基于物种保真度分析。基本思想是保留能够最大程度被诊断物种（也称为”指示种“、”典型种，”特征种“、或”差异种“表征的聚类簇，诊断物种即是在一组样方中多度相对更多且更均匀的物种。

具体而言最佳，分组结果将是最大化：ⅰ.指标值的总和 ⅱ.具有显著指示物种的聚类簇的比例。

IndVal指数集成了特异性和保真度

#通过指示种选择最优聚类簇的数量
# 分析IndVal, Dufrene-Legendre; package: labdsv
IndVal <- numeric(nrow(spe))
ng <- numeric(nrow(spe))
for (k in 2:(nrow(spe) - 1))
{
  iva <- indval(spe, cutree(hc, k = k), numitr = 1000)
  gr <- factor(iva$maxcls[iva$pval <= 0.05])
  ng[k] <- length(levels(gr)) / k
  iv <- iva$indcls[iva$pval <= 0.05]
  IndVal[k] <- sum(iv)
}
k.best <- which.max(IndVal[ng == 1]) + 1
col3 <- rep(1, nrow(spe))
col3[ng == 1] <- 3

dev.new(
  title = "IndVal-based search for optimal number of clusters",
  width = 12,
  height = 6,
  noRStudioGD = TRUE
)
par(mfrow = c(1, 2))
plot(
  1:nrow(spe),
  IndVal,
  type = "h",
  main = "IndVal-最优聚类簇数",
  xlab = "k (组数)",
  ylab = "IndVal 和",
  col = col3
)
axis(
  1,
  k.best,
  paste("最优", k.best, sep = "\n"),
  col = "red",
  font = 2,
  col.axis = "red"
)
points(
  which.max(IndVal),
  max(IndVal),
  pch = 16,
  col = "red",
  cex = 1.5
)
text(28, 15.7, "a", cex = 1.8)

plot(
  1:nrow(spe),
  ng,
  type = "h",
  xlab = "k (组数)",
  ylab = "比例",
  main = "具有显著指示物种的聚类簇比例",
  col = col3
)
axis(1,
     k.best,
     paste("最优", k.best, sep = "\n"),
     col = "red",
     font = 2,
     col.axis = "red")
points(k.best,
       max(ng),
       pch = 16,
       col = "red",
       cex = 1.5)
text(28, 0.98, "b", cex = 1.8)

image-20210507160354915

上面这些代码，除了搜索满足两个标准的分组方案之外，还用绿色凸显具有显著标指标物种的分组方案（对象ng）。对于数据集结果表明分两组或三组是最优方案。

条形图显示物种指标值（IndVal）的总和：（a）和各种分组情况具有显著指示物种的聚类簇比例（b）所有具有显著指标物种的分组方案均以虚线凸显。

2.7 最终分组的轮廓图

在我们这个案例中，基于轮廓值、矩阵相关和指示值分析给出的结果并不相同，从2组到6组都有。但是呢，分组4不算太多，也不算少是比较合适的分组结果。可以选择分组是作为最终的分组结果。

# 选择聚类簇的数量
k <- 4
# 最终分组的轮廓图
spech.ward.g <- cutree(spe.ch.ward, k = k)
sil <- silhouette(spech.ward.g, spe.ch)
rownames(sil) <- row.names(spe)
plot(
  sil,
  main = "轮廓图-弦距离-Ward聚类",
  xlab="轮廓宽度值S",
  cex.names = 0.8,
  col = 2:(k + 1),
  nmax = 100
)
#组1和组3最连贯，同时组2可能含有被错分的对象。

分组4时每个样方的轮廓宽度图

2.8 利用绘图工具修饰的最终聚类树

使用上面生成最终分4组的聚类树，利用绘图工具修饰的最终聚类树。

# 重新排列从函数hclust()获得的聚类树
spe.chwo <- reorder.hclust(spe.ch.ward, spe.ch)

# 绘制重排后带组标识的聚类数
plot(
  spe.chwo,
  hang = -1,
  xlab = "4 组",
  sub = "",
  ylab = "高度",
  main = "弦距离 - Ward 聚类(重排)",
  labels = cutree(spe.chwo, k = k)
)
rect.hclust(spe.chwo, k = k)

image-20210507162711855

#用组颜色绘制最终的树状图（RGBCMY…）

#使用附加hcoplot（）函数的快速方法：
#用法：
#hcoplot（tree=hclust.object，
#diss=相异矩阵，
#lab=对象标签（默认为空），
#k=nb.簇，
#title=paste（“重新排序的树状图从”，deparse（tree$call），
#sep=“\n”））
hcoplot(spe.ch.ward, spe.ch, k = 4)

image-20210507163817646

函数reorder.hclust（）的作用是重新排列函数从hclust（）获得的聚类树，使聚类树内对象的排列顺序与原始相异矩阵内对象的排列顺序尽可能一致。重排不影响聚类树的结构。函数reorder.hclust（）也可以直接给距离矩阵。

当使用rect.hclust()函数为聚类簇加框时，分类组数量（参数k=）的设定可以变为融合水平值的设定（参数h=），在自编函数hsoplot（）中也可以做同样的转换。

另一个函数identify.hclust（）允许通过人机交互的方式确定聚类树，将鼠标停留在任何位置并点击左键就可以选择与该位置相应的分类水平，这个操作也可以自动提取所在位置分类簇的所含对象的列表。

参数hang=-1是设定聚类树的分支终端从0的位置开始，并且聚类簇的标识将在这个设定值之下。

dendextend程序包提供各种函数来提高聚类数的表现形式。首先必须将“hclust”对象转化为“dendrogram”对象。可以聚类树不同的分支用不同的颜色和线条区分（基于最终分区）

# 将“hclust”对象转化为“dendrogram”对象
dend <- as.dendrogram(spe.chwo)
#dendextend程序包的函数对不同的分组进行标色
dend %>% set("branches_k_color", k = k) %>% plot

# 使用标准色显示聚类簇
clusters <- cutree(dend, k)[order.dendrogram(dend)]
dend %>% set("branches_k_color", 
             k = k, value = unique(clusters) + 1) %>% plot
# 增加颜色条带
colored_bars(clusters + 1,
             y_shift = -0.5,
             rowLabels = paste(k, "clusters"))

最终选择4组时带不同颜色分组的聚类树

2.9 聚类结果空间分布图

#这个绘制图版大小
dev.new(title = "Four Ward clusters on river",
        width = 9,
        noRStudioGD = TRUE)
# 4个Ward聚类簇在Doub河的分布情况
plot(spa, asp=1, type="n", main="4个Ward聚类组", 
     xlab="x坐标 (km)", ylab="y坐标 (km)")
lines(spa, col="light blue")
text(70, 10, "上游", cex=1.2, col="red")
text(15, 115, "下游", cex=1.2, col="red")
# 添加4组分类信息
grw <- spech.ward.g
k <- length(levels(factor(grw)))
for (i in 1:k) {
  points(spa[grw==i,1], spa[grw==i,2], pch=i+20, cex=3, col=i+1, bg=i+1)
}
text(spa, row.names(spa), cex=0.8, col="white", font=2)
legend("bottomright", paste("组",1:k), pch=(1:k)+20, col=2:(k+1), 
       pt.bg=2:(k+1), pt.cex=1.5, bty="n")
#请比较当前生成的地图与第2章生成的4种鱼类物种分布地图。

#第二种方法
#4个Ward聚类簇在Doub河的分布
#（see Chapter2）
drawmap(xy = spa,
        clusters = spech.ward.g,
        main = "4个Ward聚类簇在Doub河的分布")

image-20210507172957681

image-20210507173101004

2.10 热图和排序的群落表

那么如何用彩色的方阵表达聚类树，颜色的深浅代表样方之间的相似度

# 用聚类结果重排距离矩阵的热图
dend <- as.dendrogram(spe.chwo)
heatmap(
  as.matrix(spe.ch),
  Rowv = dend,
  symm = TRUE,
  margin = c(3, 3)
)

image-20210507180900289

提示：请注意，hclust类对象已转换为dendrogram类对象。只有dendro-gram才允许树状图的高级图形操作

这里说明一下，这个图的颜色好像和《数量生态学：R语言的应用》第二版的附图不一样，试了好久，有知道的小伙伴欢迎交流。

最后，直接探索聚类簇内的物种信息非常有用，即根据分组的情况重新排列原始数据表。

为了避免较大的值影响图形展示效果，vegan程序包里编写了vegemtie（）函数，专门利用外部的信息重新排列和展示样方一物种矩阵表格。如果没有提供物种的顺序，物种将按照基于样方得分的加权平均进行排列。

请注意，vegemtie（）函数所使用的多度数据必须是1位的计数数值（例如从0到9）。
任何两位或更多位数值都会导致vegemite（）停止工作。这里我们的多度数据以0~5的标尺表示，因此没有问题。其他数据可能必须重新编码。

vegemite（）内部含有专门对植被数据进行重新编码的函数covercale（）。

# 重排群落表格
# 物种按照在样方得分加权平均进行排列
#点代表缺失
or <- vegemite(spe, spe.chwo)

image-20210507181734947

# 基于聚类树的双排列群落表格的热图
heatmap(
  t(spe[rev(or$species)]),
  Rowv = NA,
  Colv = dend,
  col = c("white", brewer.pal(5, "Greens")),
  scale = "none",
  margin = c(4, 4),
  ylab = "物种 (样方的加权平均)",
  xlab = "样方"
)