【原】【R分享|实战】方差分解分析(Variation partition analysis)及其显著性检验

科白君 2021-08-15

展开全文

‍‍‍‍

“ 世上所有的成功只是努力的另一个代名词。” --科白君

"R实战"专题·第13篇

编辑 | 科白维尼

3470字 | 8分钟阅读

本期推送内容

高分期刊中频频登场的方差分解分析(VPA)到底是啥？可能很多同学还没接触过该分析。目前，群落分析中常见的环境因子分析包括典范对应分析（canonical correspondence analysis）和冗余分析 (Redundancy analysis)，这两种类型分析都是基于降维的思想，将样本、物种、环境因子的信息映射到二维平面上，从而判断三者间的关系，可用于发现对群落结构有影响的环境变量。而VPA分析可以看做是CCA/RDA分析的一种升级。用几组解释变量（如环境、气候、土壤因子等数据）来共同解释一组响应变量的变化（如微生物数据），当我们需要某个解释变量所能够解释的方差变异程度信息（即某个环境因子对群落结构变化的贡献度）时，就可以采用VPA分析加以补充。

方差分解分析

方差分解的数学原理，可以参考该链接，有比较清晰且易理解的介绍https://mb3is./gustame/constrained-analyses/variation-partitioning。

从土壤微生物生态领域的角度：就是将每个解释变量（环境因子）进行独立运行CCA或RDA，获得每个解释变量（环境因子）对于响应变量（微生物群落）的方差变异的解释贡献度，之后通过将多组数据取交并集的方式获得每个解释变量（环境因子）的独立解释贡献度以及环境因子共同解释的贡献度。

今天目的是教会大家利用R语言进行VPA分析及其对应解释度的显著性检验。废话不多说，进入代码讲解~

利用R语言进行VPA分析

这里主要到"vegan"包中的varpart()函数并使用自带的数据集来进行VPA分析，并使用plot以图形式输出结果。温馨提示：关于VPA分析时所用的环境因子的分类是可以根据你自己选择来划分的，例如现有气候因子年均温、年降雨量；土壤理化因子有pH、SOC、TN、TP、TK、含水率、阳离子交换量。响应变量为土壤微生物丰度数据。因此，我将所有变量划分为两个解释变量：气候和土壤理化。

加载自带数据集 "mite"、"mite.env"、"mite.pcnm"。其中mite为物种丰度数据，mite.env和mite.pcnm为环境因子数据。

# intall.packages("vegan")
# 加载vegan包
library(vegan)
# 读取vegan包自带数据集
data(mite)
head(mite)
data("mite.env")
head(mite.env)
data("mite.pcnm")
head(mite.pcnm)

结果：

注意，数据格式的模仿：数据集的行为样本编号，列为物种或者环境变量。

接着，我们利用varpart()函数进行拟合模型。由于mite.pcnm的数据集太大了，我们就取了取了前5个作为第二个环境因子。

# 拟合模型
?varpart # 简单了解函数用法
fit <- varpart (mite, mite.env, mite.pcnm[,1:5], transfo="hel") 
fit

进行VPA分析时，需要记住：第一个数据物种数据，之后两个数据分别代表不同的环境变量，transfo表示对数据进行转换，hel为hellinger转换，避免“弓形效应”。

第一个结果：

为了让大家更好的理解VPA分析的结果，这里举例解读

如图所示微生物群落的方差总体解释量为Y，检测了两个环境因子数据X1、X2。VPA是将两个环境因子的数据各自独立运行RDA分析（RDA则为标准化后的解释变量对响应变量逐一进行多元回归分析，获得拟合值、残差值，最终整合成为拟合值矩阵，该矩阵进行PCA、然后排序所得结果），获得每个环境因子对于群落整体变差校正解释量，再运行两者共同存在时获得校正后解释度数据R方，具体得到如下结果：

[独立运行X1的RDA]：A+C部分解释贡献度占比

[独立运行X2的RDA]：B+C部分解释贡献度占比

[共同运行X1和X2的RDA]：A+B+C部分解释贡献度占比

最后可以得到：

[A]= [共同运行X1、X2的RDA]- [独立运行X2的RDA]

[B]= [共同运行X1、X2的RDA]- [独立运行X1的RDA]

[C]= [共同运行X1、X2的RDA]- [A]- [B]

[未被这两种环境因子解释到的残差(D)]=Y- [A]-[B]- [C]

第二个结果：

使用plot函数可视化得到的VPA分析结果：

plot(fit, bg = c("hotpink","skyblue")) # bg表示两个变量的背景颜色

结果：

通过维恩图的展示就可以清楚的看到每种环境因子对于总体变异的解释程度、共有解释程度以及残差。需要注意的事，图中共有的部分产生的原因在于环境因子数据对微生物解释存在着共线性而产生，如果环境因子完全相互独立理论上重叠部分=0；此外，如果解释贡献度出现负数，则说明这组环境因子数据对群落数据方差变化解释程度比使用随机变量的解释程度还要低，一般当做贡献度为0进行解释，建议在选择环境因子时减少共线性程度较高的环境因子以及贡献度为负数的环境因子数据，以保证结果准确。

VPA分析后续的显著性检验

实际上RDA和方差分解的结果是相同的（这里不做展示），可以发现环境X1与X2共同解释了26% 的群落结构变异，其中X1解释了 16% 的群落组成变异，气候解释了5% 的群落组成变异。然而对于二者是否显著，是可以做进一步检验。特殊情况，关于冗余分析（RDA）中环境因子共同解释部分出现负值的说明。我推荐赖江山老师对于该结果的深刻解读。http://blog.sciencenet.cn/blog-267448-1187530.html。

代码如下：

# 首先将两个环境因子合并
mite.total <- data.frame(mite.env,mite.pcnm[,1:5])

# 描述partial RDA的公式 
# X1
formula_X1 <- formula(mite ~ SubsDens+WatrCont+Substrate+Shrub+Topo+ Condition(V1)+Condition(V2)+Condition(V3)+Condition(V4)+Condition(V5))
# X2
formula_X2 <- formula(mite ~ Condition(SubsDens)+Condition(WatrCont)+Condition(Substrate)+Condition(Shrub)+Condition(Topo)+V1+V2+V3+V4+V5)
# X1和X2
formula_X1X2 <- formula(mite ~ SubsDens:WatrCont:Substrate:Shrub:Topo:V1:V2:V3:V4:V5)

# 利用partial RDA进行999次的置换检验，最后得出模型的显著性
anova(rda(formula_X1, data = mite.total))
anova(rda(formula_X2, data = mite.total))
anova(rda(formula_X1X2, data = mite.total))

其中：Condition是表示控制一部分变量，只探究非condition部分的变量，即控制mite.pcnm则只探究mite.env。反之，探究另一个环境因子。特别强调：最后一种情况是我理解上的共同解释部分的显著性检验，但是对错与否目前我也没找到正确答案。这里需要更多的统计学专业的同学给出建议来完善。

结果如下：