聚类分析4—环境数据来解释 (数量生态学：R语言的应用-第四章)

2021-05-10 本文已影响0人 fafu生信小蘑菇

聚类分析4—环境数据来解释 (数量生态学：R语言的应用-第四章)

在这之前我们学习了聚类分析的基本概念、几种计算层次聚类的方法、进一步解读和比较层次聚类结果以及非层次聚类，这些聚类方法都是基于物种多度数据对样方进行分组，当然这些聚类方法也可以用于其他类型数据，特别是环境数据，所以本次就是介绍用环境数据来进行聚类分析。

本次的内容不多，主要分为两个部分：

用外部数据进行类型比较（方差分析途径）

双类型比较（列联表分析）

1.1 加载所需的包和数据

#加载包和数据
library(ade4)
library(adespatial)
library(vegan)
library(gclus)
library(cluster)
library(pvclust)
library(RColorBrewer)
library(labdsv)
library(rioja)
library(indicspecies)
library(mvpart)
library(MVPARTwrap)
library(dendextend)
library(vegclust)
library(colorspace)
library(agricolae)
library(picante)

#加载函数
source("drawmap.R")
source("drawmap3.R")
source("hcoplot.R")
source("test.a.R")
source("coldiss.R")
source("bartlett.perm.R")
source("boxplerk.R")
source("boxplert.R")

#从聚类结果获得二元差异矩阵的函数
grpdist <- function(X)
{
  require(cluster)
  gr <- as.data.frame(as.factor(X))
  distgr <- daisy(gr,"gower")
  distgr
}

#导入Doubs数据
load("Doubs.RData")
#剔除无物种数据的样方8
spe <- spe[-8,]
env <- env[-8,]
spa <- spa[-8,]
latlong <- latlong[-8,] #经纬度

1.2数据预处理

#先计算样方之间的弦距离矩阵
spe.norm <- decostand(spe,"normalize")
spe.ch <- vegdist(spe.norm,"euc")

# 计算Ward最小方差聚类
spe.ch.ward <- hclust(spe.ch, method="ward.D2")

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

# 设定聚类组的数量
k <- 4  
# 根据上面4个融合水平值图，可以观察到分4组水平在所有图里有小的跳跃
# 裁剪聚类树

spech.ward.g <- cutree(spe.ch.ward, k=k)
# 计算Ward聚类获得4组中各个物种的平均多度矩阵
groups <- as.factor(spech.ward.g)
spe.means <- matrix(0, ncol(spe), length(levels(groups)))
row.names(spe.means) <- colnames(spe)
for (i in 1:ncol(spe)) {
  spe.means[i, ] <- tapply(spe.norm[, i], spech.ward.g, mean)
}
# 平均物种多度矩阵作为初始点
startpoints <- t(spe.means)
# 基于初始点的k-均值划分
spe.kmeans2 <- kmeans(spe.norm, centers = startpoints)

# 最终分组的轮廓宽度值图
spech.ward.gk <- spe.kmeans2$cluster
dev.new(
  title = "轮廓-优化分区",
  noRStudioGD = TRUE
)
par(mfrow = c(1, 1))
k <- 4
sil <- silhouette(spech.ward.gk, spe.ch)
rownames(sil) <- row.names(spe)
plot(sil,
     main = "轮廓宽度值 - Ward & k均值",
     cex.names = 0.8,
     col = 2:(k + 1))

说明：这里数据处理方法，都是在聚类分析这一章学过的，不明白的同学，可以在回头去看看

这里用到的数据是非层次距离中得到的spech.ward.gk

2 用外部数据进行类型比较（方差分析途径）

我们之前学习的主要是内部的准则（例如轮廓法或其他聚类质量指数）都是仅仅依赖物种数据，还不足以选择最佳样方聚类结果。选择最终的聚类结果有时也需要基于生态学解释。生态学解释可视为样方聚类的外部验证。

利用外部独立的解释变量验证聚类结果（作为响应数据）可以用线性判别式分析（这个在后面我们会学习到）。

可以将样方分组当作因子对解释变量（当作响应变量值）进行方差分析，了解解释变量在各组间是否有显著差异。

下面的我们将学习用样方聚类簇为因子去对解释变量进行方差分析。

步骤：

检验某一环境变量是否符合方差分析假设（残差正态性和方差齐性）

用传统的（参数估计）单因素方差分析或非参数的Kruskal-Wallis检验解释变量在组间是否有显著差异

尽管在方差分析中，是将物种组成数据获得的聚类的分组结果作为解释变量，但是从生态学角度去分析，实际上是寻找环境因子对样方的分组的解释。

基于最优化的Ward聚类（分4组）的环境变量（为提高正态性进行某些简单的转化）箱线图

#绘制基于最优化的Ward聚类（分四组）的环境变量线箱线图
with(env, {
  dev.new(
    title = "定量环境变量箱线图",
    noRStudioGD = TRUE
  )
  par(mfrow = c(2, 2))
  boxplot(
    sqrt(ele) ~ spech.ward.gk,
    main = "海拔",
    las = 1,
    ylab = "sqrt(alt)",
    col = (1:k) + 1,
    varwidth = TRUE
  )
  boxplot(
    log(slo) ~ spech.ward.gk,
    main = "坡度",
    las = 1,
    ylab = "log(slo)",
    col = (1:k) + 1,
    varwidth = TRUE
  )
  boxplot(
    oxy ~ spech.ward.gk,
    main = "含氧量",
    las = 1,
    ylab = "oxy",
    col = (1:k) + 1,
    varwidth = TRUE
  )
  boxplot(
    sqrt(amm) ~ spech.ward.gk,
    main = "铵浓度",
    las = 1,
    ylab = "sqrt(amm)",
    col = (1:k) + 1,
    varwidth = TRUE
  )
})

最优化的Ward聚类（分四组）的环境变量线箱线图

# 方差分析假设检验
with(env, {
  # 残差的正态性
  shapiro.test(resid(aov(sqrt(ele) ~ as.factor(spech.ward.gk))))
  shapiro.test(resid(aov(log(slo) ~ as.factor(spech.ward.gk))))
  shapiro.test(resid(aov(oxy ~ as.factor(spech.ward.gk))))
  shapiro.test(resid(aov(sqrt(amm) ~ as.factor(spech.ward.gk))))
  
  #检验结果表明sqrt(alt)、log(pen)、oxy和sqrt(amm)的残差是正态分布。尝试为其他的环境变量寻找好的标准化
  
  # 方差齐性
  bartlett.test(sqrt(ele), as.factor(spech.ward.gk))
  bartlett.test(log(slo), as.factor(spech.ward.gk))
  bartlett.test(oxy, as.factor(spech.ward.gk))
  bartlett.test(sqrt(amm), as.factor(spech.ward.gk))
  #变量sqrt(alt)的方差不齐，所以参数检验的方差分析不适用
  
  # 可检验变量的方差分析A
  summary(aov(log(slo) ~ as.factor(spech.ward.gk)))
  summary(aov(oxy ~ as.factor(spech.ward.gk)))
  summary(aov(sqrt(amm) ~ as.factor(spech.ward.gk)))
  
  # 海拔Kruskal-Wallis 检验
  kruskal.test(ele ~ as.factor(spech.ward.gk))
  
})

注意：

在一系列分析的开头使用with（）来避免在每次分析中重复输入对象env的名称。这比使用attach（）和detach（）方便，因为当你的R控制台中有几个数据集，有些数据集碰巧有名字相同的变量时，attach（）可能会导致混淆。

Shapiro检验的零假设是变量正态分布；Bartlett检验的零假设是组间方差相等。因此，对于这两个检验，只有当p值大于显著性水平，即p>0.05时接受零假设，才能满足方差分析的假设。

参数Bartlett检验对偏离正态分布很敏感。对于非正态分布数据，可以使用bartlett.perm.R的函数来计算参数、置换和自助法（即替换的置换）Bartlett检验。

以下可以使用作者编写的通用函数，执行方差分析的多重比较和显示带有字母的环境变量分组后箱线图多重比较结果。不同字母表示组间有显著差异（按中位线递减顺序组）。

使用boxplert（）函数来执行方差分析和LSD多重比较检验

使用boxplerk（）函数执行Kruskal-Wallis 检验及其相应的多重比较（两者都使用“Holm”的p值校正）

dev.new(
  title = "ANOVA and Kruskal-Wallis 检验",
  noRStudioGD = TRUE
)
par(mfrow = c(2, 2))
#使用boxplert（）或boxplerk（）绘制事后检测的结果
#使用boxplert（）函数画矫正后的多重比较结果
with(env, {
  boxplerk(
    ele,
    spech.ward.gk,
    xlab = "",
    ylab = "ele",
    main = "海拔",
    bcol = (1:k) + 1,
    p.adj = "holm"
  )
  boxplert(
    log(slo),
    spech.ward.gk,
    xlab = "",
    ylab = "log(slo)",
    main = "坡度",
    bcol = (1:k) + 1,
    p.adj = "holm"
  )
  boxplert(
    oxy,
    spech.ward.gk,
    xlab = "",
    ylab = "oxy",
    main = "含氧量",
    bcol = (1:k) + 1,
    p.adj = "holm"
  )
  boxplert(
    sqrt(amm),
    spech.ward.gk,
    xlab = "",
    ylab = "sqrt(amm)",
    main = "铵浓度",
    bcol = (1:k) + 1,
    p.adj = "holm"
  )
})

image-20210509200652145

基于上面这些分析和图示，能确定这组鱼类群落的生态习性。

当然，我们也可以基于环境变量对样方进行聚类（类似获得生境类型的分组），然后通过指示种分析（以后会讲）检验不同生境内物种分布是否有差异。指示种分析过程中基于不同的生境类型物种需要逐个分析。因此，需要考虑多个物种指示种分析时会产生多重检验的统计学问题。

另外，作为替代方案，之后第6章会提出基于排序的多元方法，也可以直接描述和检验物种-生境关系。请期待。

3 双类型比较（列联表分析）

要是想直接比较分别基于物种数据和环境数据的样方聚类结果该怎么办呢？

直接生成一个含两种结果的列联表

用列联表Fishr精准检验比较两种样方聚类结果是否有显著差异

# 基于环境变量的样方聚类
env2 <- env[, -1]
env.de <- vegdist(scale(env2), "euc")
env.kmeans <- kmeans(env.de, centers = 4, nstart = 100)
env.kmeans.g <- env.kmeans$cluster

# 比较从物种和环境数据获得聚类4组的结果
table(spe.kmeans.g, env.kmeans.g)
#两种聚类结果是否相同？

比较从物种和环境数据获得聚类4组的结果

# 列联表 Fisher精确检验
fisher.test(table(spe.kmeans.g, env.kmeans.g))


# 用卡方检验分析两种聚类之间的差异 
chisq.test(table(spe.kmeans.g, env.kmeans.g))

#改用置换的方法进行卡方检验分析 
chisq.test(table(spe.kmeans.g, env.kmeans.g), 
                          simulate.p.value = TRUE)