RDA、db-RDA（CAP）及CCA的变差分解

常规RDA、db-RDA（CAP）及CCA的变差分解

在常规RDA、db-RDA（CAP）以及CCA等约束排序中，属于不同类别的两组或两组以上解释变量共同解释一组响应变量的现象非常常见，如果这些解释变量之间存在某些联系，那么解释变差（variation）将由密不可分的多个解释变量所共享。

当存在两个或更多解释变量时，我们可能对它们中的每一个所能解释的变差产生兴趣。有时不仅仅控制一组变量后分析另一组变量的偏分析，而且需要量化两组或多组变量单独及共同解释的变差。处理这个问题的方法即为变差分解（variation partitioning），将总变差划分为由每个变量（或变量集）所能独自解释的变差部分以及共同解释的变差部分。在生态多元数据分析中，Borcard等（1992）首先提出了变差分解的概念及分解过程，Peres-Neto等（2006）提出使用校正R²促进变差分解的使用。

变差分解概述

以两组解释变量为例阐述，下图Venn图概括了响应变量的总变差组成。

当同时存在两组解释变量A和B时，图中响应变量Y的总变差可以分解为以下几个部分：

[a+b]变量A能够解释的总变差（变量A的边际效应，marginal effect，即当约束模型中仅存在变量A这一个解释变量时，约束模型所能够被解释的总变差）；

[b+c]变量B能够解释的总变差；

[a]仅能被变量A所单独解释的变差（变量A的局部效应，partial effect，即将变量B作为协变量处理时，变量A所能解释的变差部分）；

[c]仅能被变量B所单独解释的变差；

[b]变量A和B共同解释的变差（由于变量A和B具有相关性，无法归其所属）；

[d]变量A和B均无法解释的变差，即残差。

常规RDA的变差分解

在RDA中响应变量的总变差=总方差，且早期变差分解也仅能在RDA中使用，因此RDA的变差分解也常称为方差分解（variance partitioning）。

以上述两组解释变量的变差分解为例，若量化其中 [a]、[b]、[c]、[d]四部分的变差，必须运行3次RDA分析。变差分解的概念步骤如下。

（1）如果有必要，执行变量选择（例如对所有解释变量执行前向选择），只保留显著的变量。

（2）执行RDA分析，获得各解释变量所能解释的变差。对应于上图，即：

单独以A作为解释变量进行Y的RDA分析，可以获得[a+b]部分的值；

单独以B作为解释变量进行Y的RDA分析，可以获得[b+c]部分的值；

将A和B一起作为解释变量进行Y的RDA分析，可以获得[a+b+c]部分的值。

（3）R²校正也是不可或缺的。

如果每组解释变量（变量集）包含相同数量的变量（例如，两个土壤和两个气候变量），则组间解释的变差是可以直接比较的；若变量数量不一致，则需要在比较前校正R²。同时由于随机相关的存在，R²一般会伴随解释变量数量的增加或样方数的减少而增大，解释变量的累加会造成被解释变差表面上膨胀。考虑到这些原因，我们通常会计算上述3个RDA分析的校正后R²。

（4）通过减法计算获得各部分校正后的变差：

[a]_adj = [a+b+c]_adj - [b+c]_adj；

[c]_adj = [a+b+c]_adj - [a+b]_adj；

[b]_adj = [a+b]_adj - [a]_adj = [b+c]_adj - [c]_adj；

[d]_adj = 1 - [a+b+c]_adj。

上述3个RDA分析都可以应用置换检验。

同上文所述，[a]为变量A的局部效应（partial effect），即将变量B作为协变量处理时，变量A所能解释的变差部分；[c]的理解同理。因此对于[a]和[c]部分的计算和检验也可以通过偏RDA分析进行。由于目前没有公式可以直接计算偏RDA的校正R²，但上面描述的减法过程可以避免这个问题，Peres-Neto等（2006）也是用这种方法获得各部分变差的无偏估计。

同上文所述，[b]由变量A、B所共享，归因于变量A和B存在某种相关性，无法明确[b]的所属，因此[b]部分不能直接估计和检验。同时可知，若变量A和B互相独立（不相关、正交），那么二者将不存在共同解释的部分，即此时[b]将为0。

在某些情况下，获得的校正R²可能是负值，表明解释变量对响应变量的解释程度甚至不及随机变量的解释程度，该解释变量的存在意义不大。在对分析结果作生态解释时，负的校正R²可以忽略（当作0看待）。

在R中，RDA的变差分解可使用vegan包中varpart()函数来实现。

接前文RDA的数据作展示。

##常规 RDA 的变差分解
#数据的网盘链接：https://pan.baidu.com/s/18xTqAZPa2Al2sDQpS1ysFA

#读取物种数据，细菌门水平丰度表
phylum <- read.delim('phylum_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#推荐在对纯物种多度数据执行 PCA、RDA 等线性排序方法前，执行 Hellinger 转化
phylum_hel <- decostand(phylum, method = 'hellinger')

#读取环境变量
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

library(vegan)

#运行变差分解，详情 ?vegdist
#以两组（pH 为一组，5 种碳氮磷钾指标为另一组）环境变量为例
rda_vp <- varpart(phylum_hel, env['pH'], env[c('DOC', 'SOM', 'AP', 'AK', 'NH4')])
rda_vp

plot(rda_vp, digits = 2, Xnames = c('pH', 'CNPK'), bg = c('blue', 'red'))

这里选择了6个环境变量（前文RDA分析中，前向选择后保留了这6种环境变量，此处拿来展示），根据环境变量的属性，将它们划分为两组：一组为土壤化学属性，pH；另一组则包括了土壤碳氮磷钾含量，DOC、SOM、AP、AK、NH4。对这两组环境变量执行变差分解后，结果如下。

其中，pH即对应于上图中的“a+b”部分，解释变差约13.2%，其中单独解释部分为“a”，解释变差约1%；同理，DOC、SOM、AP、AK、NH4这5个的合集即对应于上图中的“b+c”部分，解释变差约57.1%，其中单独解释部分为“c”，解释变差约44.9%；两组环境变量共同解释部分即为“b”部分，解释变差约12.2%；两组环境变量均不能解释的部分为“d”，约41.9%。

#再以四组为例，碳、氮、磷、钾指标间比较
rda_vp <- varpart(phylum_hel, env[c('DOC', 'SOM')], env['NH4'], env['AP'], env['AK'])
rda_vp

plot(rda_vp, digits = 2, Xnames = c('C', 'N', 'P', 'K'), bg = c('blue', 'red', 'green', 'yellow'))

在前文RDA中，我们发现在经过变量的前向选择后，“TC”这个环境变量被剔除了。通常由于有些解释变量之间可能存在较强的线性相关，即共线性问题，可能会造成回归系数不稳定，所以前向选择程序一般会将这些变量剔除。那么，对于“TC”这个环境变量来讲，是否是因为这个原因呢？

#结合前文 RDA 中的示例分析
#查看前向选择中被剔除的环境变量“TC”，与这 6 个被保留的环境变量之间解释变差的“共享程度”
rda_vp <- varpart(phylum_hel, env['TC'], env[c('pH', 'DOC', 'SOM', 'AP', 'AK', 'NH4')])
plot(rda_vp, digits = 2, Xnames = c('TC', 'forward_env'), bg = c('blue', 'red'))

事实的确如此，TC所解释的变差，均与这6个被保留的（前向选择后所保留的）环境变量之间产生交集，即存在非常强的共线性问题。而对于TC本身来讲，它所能单独解释的变差部分为0，所以前向选择程序自动将它剔除了。

还可以通过置换检验，检验各解释变差部分的显著性。例如，对pH所解释变差部分的置换检验如下所示。

##置换检验确定解释变差的显著性
#需结合 RDA 执行式，vegan 中为 rda()

#解释变差的置换检验，以 pH 所能解释的全部变差为例；999 次置换，详情 ?anova.cca
anova.cca(rda(phylum_hel, env['pH']), permutations = 999)

#若考虑 pH 单独解释的变差部分，需将其它变量作为协变量；999 次置换
anova.cca(rda(phylum_hel, env['pH'], env[c('DOC', 'SOM', 'AP', 'AK', 'NH4')]), permutations = 999)

注意观察，存在不显著的解释变差部分吗？对于这部分，应当如何处理呢？

db-RDA（CAP）的变差分解

db-RDA（CAP）也是在常规RDA算法的基础上衍生而来的，对于db-RDA的变差分解，仍然可以通过R包vegan的varpart()函数执行。

接前文db-RDA的数据作展示。

##db-RDA 的变差分解
#数据的网盘链接：https://pan.baidu.com/s/1itLEWnV_sLpKohHleMwggA

#读入物种数据，以细菌 OTU 水平丰度表为例
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- data.frame(t(otu))

#读取环境数据
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

#计算样方距离，以 Bray-curtis 距离为例，详情 ?vegdist
dis_bray <- vegdist(otu, method = 'bray')
#或者直接使用已准备好的距离矩阵，网盘中同样提供了已计算好的 Bray-curtis 距离矩阵
dis_bray <- as.dist(read.delim('bray_distance.txt', row.names = 1, sep = '\t', check.names = FALSE))

library(vegan)

#以两组环境变量（DOC、AP+AK）为例，运行变差分解，详情 ?varpart
#输入 varpart() 的为已经计算好的 Bray-curtis 距离测度
#参数 add=T 意为校正 PCoA 的负特征根，这个最好和先前计算 db-RDA 的执行式保持一致（即二者要么都校正，要么都不校正）
db_rda_vp <- varpart(dis_bray, env['DOC'], env[c('AP', 'AK')], add = TRUE)
db_rda_vp

plot(db_rda_vp, digits = 2, Xnames = c('DOC', 'AP+AK'), bg = c('blue', 'red'))

对于该示例，DOC（“a+b”部分）解释了总变差的14.5%，其中单独（“a”部分）解释的比例为4.2%；AP+AK单独（“c”部分）解释了9.4%，和DOC共同（“b”部分）解释了10.3%；DOC+AP+AK未能（“d”部分）解释的变差占比76.0%。

##置换检验确定解释变差的显著性
#需结合 db-RDA 执行式，vegan 中提供了 capscale()、dbrda() 等函数可直接执行 db-RDA
#下式均以 capscale() 为例

#解释变差的置换检验，以 DOC 所能解释的全部变差为例；999 次置换，详情 ?anova.cca
anova.cca(capscale(dis_bray~DOC, env, add = TRUE), permutations = 999)

#若考虑 DOC 单独解释的变差部分，需将其它变量作为协变量；999 次置换
anova.cca(capscale(dis_bray~DOC+Condition(AP+AK), env, add = TRUE), permutations = 999)

CCA的变差分解

由于在CCA中，响应变量的总变差通过惯量（而非方差）表征，因此“方差分解”这个名称最好不要在CCA中使用，还是使用“变差分解”为好。

早期由于计算方法的原因，CCA的变差分解难以实现。例如在RDA中常用的Ezekiel校正法在CCA中不适用；且RDA变量选择中的常用方法，如Borcard终止准则（Blanchet等，2008）也不能很好地用于CCA。这些阻碍了CCA分析中正确估计解释变量能够解释的变差比例（即无法准确计算校正后R²），也过度夸大了I类错误，导致了无法对CCA模型执行变差分解。

后来随着技术的更新，一系列的方法被引入。同样以R包vegan为例，RsquareAdj()已能够对CCA的R²实现校正（早期是不能的，仅能用于提取CCA的R²但不能校正），ordiR2step()也已将Borcard终止准则引入至CCA的变量前向选择中，类似地，varpart()函数也可用于执行CCA的变差分解。

接前文CCA的数据作展示。

##CCA 的变差分解
#数据的网盘链接：https://pan.baidu.com/s/1UWlvIUWrnl4gyoEhpCGVtw

#读入物种数据，以细菌 OTU 水平丰度表为例
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- data.frame(t(otu))

#读取环境数据
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

library(vegan)

#以两组环境变量（DOC、AP+AK）为例，运行变差分解，详情 ?varpart
#参数 chisquare = TRUE，执行 CCA 的变差分解；默认情况下 chisquare = FALSE，即执行 RDA 的变差分解
#注：如果 varpart() 不支持 CCA，请更新 R 版本（如 R3.6 的 vegan）
cca_vp <- varpart(otu, env['DOC'], env[c('AP', 'AK')], chisquare = TRUE)
cca_vp

plot(cca_vp, digits = 2, Xnames = c('DOC', 'AP+AK'), bg = c('blue', 'red'))

对于该示例，DOC（“a+b”部分）解释了总变差的12.4%，其中单独（“a”部分）解释的比例为3.8%；AP+AK单独（“c”部分）解释了8.9%，和DOC共同（“b”部分）解释了8.6%；DOC+AP+AK未能（“d”部分）解释的变差占比78.7%。

##置换检验确定解释变差的显著性
#需结合 CCA 执行式，vegan 中为 cca()

#解释变差的置换检验，以 DOC 所能解释的全部变差为例；999 次置换，详情 ?anova.cca
anova.cca(cca(otu~DOC, env), permutations = 999)

#若考虑 DOC 单独解释的变差部分，需将其它变量作为协变量；999 次置换
anova.cca(cca(otu~DOC+Condition(AP+AK), env), permutations = 999)

注意事项

最后，还需注意几点。

变差分解的一个主要目的是计算共同解释的变差，但在解读时必须非常谨慎，因为共同部分很难确定与哪组变量有关。

变差分解共同解释部分（对应于前图“b”部分）与方差分析中因子的交互作用（interaction）不同。交互作用是双因素方差分析中，一个因子不同水平与其它因子不同水平之间的协同作用。变差分解是由于解释变量共线性引起解释变差的重叠，不是协同作用的结果。换句话说，如果两个（组）变量相互独立（不相关，正交），则共同解释部分“[b]”将为0。

参考资料

DanielBorcard, FranoisGillet, PierreLegendre, et al. 数量生态学:R语言的应用（赖江山 译）. 高等教育出版社, 2014.

Blanchet F G , Legendre P , Borcard D . Forward selection of explanatory variables. Ecology, 2008, 89(9):2623-2632.

Borcard D , Drapeau L P . Partialling out the Spatial Component of Ecological Variation. Ecology, 1992, 73(3):1045-1055.

David Zelený博士：https://www.davidzeleny.net/anadat-r/doku.php/en:varpart

Peres-Neto P R , Legendre P , Dray, Stéphane, et al. Variation partitioning of species data matrices: estimation and comparison of fractions. Ecology, 2006, 87(10):2614-2625.

友情链接