数量生态_db-RDA之前世今生

1. 从多因子设计实验说起

包含假设检验的设计实验（structured experiment）对生态学至关重要，这种实验通常包含复杂的、多因子设计。而方差分析等统计方法提供了复杂模型中单个处理因子效应的假设检验。以两因子为例，如果存在显著交互作用，这意味着一个因子A的效应在另一个因子B的不同水平下是不一样的。但是，各因子的交互作用信息不能通过单因子实验得到。如果是多因子实验的话，情况更为复杂。

然而，一些生态学假设即是探讨整个群落组成（多变量数据）对环境因子的响应。参数化的多变量方差分析（MANOVA）在生态学中的应用存在诸多限制。首先，大多数生态学数据很难满足正态性和方差齐性。群落组成的数据矩阵中由稀有种造成的0值会破坏数据正态性。其次，传统的MANOVA受制于有限的变量个数，样本数一般要大于变量数。对于少量变量，可以通过特定数量的重复来解决。而生态学数据中，物种数远远超过重复数是很常见的，尤其是微生物群落组成数据。近年来，逐渐出现一些非参数化的多变量检验模型，如Mantel’s test，ANOSIM等等。这些模型依然采用方差分析的方法，但是采用置换检验，因此不需要满足多元正态性假设。研究人员可以灵活选择适合数据的距离矩阵，而不仅仅是采用经典MANOVA所要求的欧式距离。

尽管这些非参数模型在生态学中很有吸引力和实用性，但是它们却有着不能在ANOVA框架下检验因子间多变量交互作用的缺陷。特别是，这些方法不能够在结构化模型下区分多变量变异的不同部分。原因有二：其一，距离矩阵不是度量的（metric），因此线性模型不能直接使用。 其二，其对置换检验的依赖。多变量交互项不能够通过距离矩阵的置换来检验。这是因为置换检验的零假设不能用重复之间的可交换性来表达（翻译能力有限，原文The reason for this is that the null hypothesis cannot be articulated in terms of the exchangeability of the original replicates）。那么，如何采用基于距离的方法处理多变量数据中的因子交互项？distance-based RDA应运而生。

综上，db-RDA主要解决一下三个问题：

采用主坐标分析将非度量（nonmetric）和半度量（semi-metric）距离矩阵转换到欧式空间，这样可以使用线性ANOVA
基于多元回归的线性模型，将多变量的RDA统计数与单变量ANOVA的F统计数对应起来
通过置换过程，采用RDA检验多变量实验的因子交互项

db-RDA过程包含样本间距离矩阵计算、PCoA、RDA等过程，其中预测变量矩阵X包含一个ANOVA模型中的dummy variables，响应变量矩阵Y包含主坐标轴。

Legendre P.

2. 拆解db-RDA

2.1 主坐标分析

首先来回顾下PCoA。Gower 于1966年提出PCoA分析，也即是度量多维标度。PCoA采用样本间任意类型的距离矩阵，并形成每个样本对应的笛卡尔坐标，全维度的主坐标空间的坐标保留了样本间原始距离关系。PCoA的主要过程简述如下：

首先对距离矩阵D进行转换，获得矩阵A， $\alpha_{ij}=-1/2 d_{ij}^2$
将矩阵A中心化，获得矩阵 $\Delta_1$ ，其元素 $\delta_{ij} = \alpha_{ij} - \overline{a_i} - \overline{a_j} + \overline{\alpha}$ (分别为矩阵A行均值，矩阵A列均值，矩阵A均值)
计算矩阵 $\Delta_1$ 的特征值和特征向量
获得主坐标轴，用对应的特征值平方根对特征向量进行scale

对于n个样本来讲，我们最多需要n-1个主成分轴来在欧式空间展示所有的样本点。然而，由于0特征值的存在，主成分轴的个数小于等于n-1。样本重复数、响应变量数和距离系数的选择都会影响主成分轴的数量。对于欧式距离来讲，如果响应变量矩阵Y的行数大于列数，最大主坐标数量即是Y的变量数，并且主坐标即是主成分。对其他距离系数来讲，PCoA分析会产生更多的主坐标轴。

对于度量性的距离来讲（欧式距离、卡方距离等），PCoA分析的主坐标轴将会保留样本间原有的距离关系D。对于非度量或半度量性的距离（Bray-Curtis距离）来讲，PCoA分析只会保留原有距离矩阵的欧式空间部分，而其余部分呈现为负值，不存在实际的主坐标轴（imaginary axis）。这些负值对应的变异，不能在欧式空间中体现。如果仅对应正特征值的主坐标用来进行RDA分析，那么RDA统计数将会是预测变量X解释响应变量Y变异的有偏估计。ape包的pcoa()函数默认选取正值的PCoA轴，在有限空间的排序。

# 以Iris数据为例，展示不同距离矩阵的特征值区间
Dbray = vegdist(iris[,1:4], method="bray")
DEuc = vegdist(iris[,1:4], method="euclidean")
Dchi = vegdist(iris[,1:4], method="chisq")
Resbray <- pcoa(Dbray)
Reseuc <- pcoa(DEuc)
Reschi <- pcoa(Dchi)
a = Resbray$values
a1 = Reseuc$values
a2 = Reschi$values

> range(a$Eigenvalues)  # 负值
[1] -0.05937494  2.34727592
> range(a1$Eigenvalues)
[1]   3.551429 630.008014
> range(a2$Eigenvalues)
[1]  0.09966624 10.72547109

2.2 负值校正

针对非欧式空间部分变异，Gower和Legendre总结了两种不同的校正方法。

常数c₁校正，将原始距离矩阵进行 $d_{ij}^\prime = (d_{ij}^2+2c_1)^{0.5}$ 转换得到新的距离矩阵。常数c₁是矩 $\Delta_1$ 最大负特征值的绝对值。这种校正来源于Lingoes的早期工作，函数pcoa()中设置correction="lingoes"。常数c₁是产生欧式坐标的最小值，任何大于c₁的值都可以消除负特征值，使整个体系满足欧式空间。校正后，新的PCoA分析最多包含n-2个正特征值，至少2个空特征值，不包含负值。
常数c₂校正，通过 $d_{ij}^\prime = d_{ij}+c_2$ 转换得到新的距离矩阵，常数c₂是以下非对称阵的最大特征值。
$\begin{bmatrix} 0 & 2\Delta_1 \\ -I & -4\Delta_2 \end{bmatrix}$ 此中， $\Delta_2$ 是原始距离矩阵D通过 $-1/2d_{ij}$ 转换并标准化得到。常数c₂是产生欧式坐标的最小值，任何大于c₂的值都可以消除负特征值，使整个体系满足欧式空间。

上述两种方法中，矩阵A和D对角线的元素均没有发生变化（均为0）。通过校正之后，主坐标分析形成了样本点在欧式空间的全排布。这些主坐标轴可以利RDA模型进行线性模型分析。值得注意的是，某些情况下，PCoA的负特征值可以通过其他类型的数据转换消除。例如，对于Bray-Curtis矩阵来说，原始距离矩阵通过平方根转换可以消除PCoA分析产生的负特征值。尽管这个还没有得到数学证明，但是Bray-Curtis矩阵的二进制形式——Sørensen距离经过平方根转换确实可以实现。

2.3 “四个问题”

尽管这些校正过程经过了数学验证，但是实践生态学家依然心存疑惑：

1）常数c₁和c₂到底有多大？
2）距离矩阵经过校正发生了怎样变化？
3）校正过程是否对F统计数的构建和假设检验产生影响？
4）到底该用哪种校正？

2.3.1 问题1

针对第一个问题，我们依据标准对数正太分布来模拟一个物种丰度的群落组成，物种权重来源于一个均匀分布。对响应变量矩阵计算Bray-Curtis距离，值的区间为[0,1]。校正系数的值随着样本数和物种数的比值增加而增加。在物种极少的生态系统中，校正系数最大。

问题1-模拟结果

2.3.2 问题2

针对第二问题，直接比较原始Bray-Curtis距离矩阵和不同校正下的欧式矩阵，构建线性回归模型。如下图所示，m₀为未校正结果，；m₁是Lingoes校正，整体上呈现凹曲线，值区间呈现收缩；m₂是Cailliez校正；m₃是对Bray-Curtis矩阵进行平方根转换后进行PCoA分析，值区间也呈现收缩，并且曲线方向与m₁相反（如下图）。

问题2-比较

2.3.3 问题3和4

针对第三个问题，依然采用同一套模拟数据，看不同校正方法下的F统计数变化。不同的校正方法下，所有的F统计数都收缩，越来越靠近1（见下图问题34-1）。这意味着，不管采用什么校正方法，样本间的相对距离变得不再明显。对于方法1，F统计数变化是单调的（见下图问题34-2），校正后的概率与原始矩阵的是一致的。因此，Type I error和powers在校正前后不变（见下图问题34-3，问题34-4）。

相反，方法2和3的F统计数不是单调的。F统计数参考值相比收缩的置换分布显得大了，因此，P值变小了，这导致了膨胀的Type I error和人为造成的powers增加。这些问题在方法3中更严重。因此，这两个方法均不适合于置换检验和F统计数构建（详见附录A、B、C）。因此，校正方法1是唯一合适的负特征值处理办法，并且这个过程对F统计数的可能性不会产生影响

问题34-1 F统计数收缩

问题34-2 统计数比较

问题34-3 表1

问题34-3 表2

2.4 构建多变量分析的 F^#statistic

RDA分析最初由Rao提出，由Van den Wollenberg发扬光大。将redundancy看作一个量，Gittins这样描述：一个测量域（measurement domain）由其他域（other domain）的线性组合所预测的总变异的比例。因此，单词redundancy即是解释的变异。RDA常用于生态学中的直接排序分析，即响应变量矩阵Y与预测变量矩阵X之间的关系。CNOCO提供了RDA分析过程。然而，RDA也可以用于多变量假设检验（multivariate hypothesis-testing）。此种情况下，“other domain”包含了一系列描述实验处理因子的哑变量（dummy variables）。为了理解通过RDA构建多变量统计数，我们有必要以多元回归的方式理解方差分析原理。

Rao C. R.

单变量方差分析的F统计数与RDA分析的多变量 $F^\#$ 统计数对应关系陈述如下。这种对应关系的理解可以参考一下两点：1）单变量方差分析和哑变量回归的等同性；2）RDA算法本质上是一系列多元回归的组成（参见https://www.jianshu.com/p/814bfadbefd7）。RDA中典范特征值的和等于矩阵Y被矩阵X解释的方差。通过计算过程，这些特征值可以转换为单位量，对应矩阵Y的总平方和、矩阵Y的总变异或者矩阵Y变异的一部分。RDA进行假设检验的 $F^\#$ 统计数公式为：
$F^\# =\frac{trace/q}{RSS/(N-q-1)}$ 式中，trace为典范特征值总和，RSS为残差平方和（sum of all unconstrained eigenvalues − trace），q为矩阵X的变量数。

在RDA分析应用于典范梯度排序的背景下，该分析的目的是探索矩阵Y与矩阵X的关系。其零假设是物种组成变异与环境因子的线性组合之间没有显著关系。而在探索交互项效应的背景下，矩阵X并不包含环境变量数据，而是设计实验中因子水平对应的哑变量。因此，因子A不存在显著效应的零假设可以通过 $F^\#$ 统计数以类似单变量方差分析的形式来进行假设检验。Verdonschot and ter Braak也将之称为“stacked”F statistic。我们将之称为“pseudo”F statistic，它包含1）施加于所有物种上的因子平方和（除以自由度）作为分子，2）施加于所有物种上的残差平方和（除以自由度）作为分母（如下表）。

单变量与多变量统计数对应关系

2.5 置换检验

采用RDA进行多变量假设检验，必须注意分母均方。每一项的F统计数构建都要与对应的零假设一致。虽然置换检验对数据分布没有要求，但是实验设计和零假设都会影响置换过程。Manly认为无交互作用的零假设与主效应没关系，因此可以直接对所有重复进行置换即可。而Edgington则不认同。ter Braak提出了一个方法来解决这个问题，即在拟合协变量和矩阵X的变量之后，对矩阵Y变异的残差部分进行置换。此种方法的优势在于：1）交互项可以被检测；2）预测变量矩阵X和协变量矩阵Xc的结构在置换过程中没有变化；3）残差置换的解释力更强。CANOCO默认采用此种方法进行置换检验。模拟结果表明ter Braak的置换方法在特定显著水平下保持了Type I error，对几乎等同于多元回归的备择假设具有解释力。此种方法被ter Braak称为“permutations under the alternative hypothesis”和“permutations under the full model”。

3. 生态案例

澳洲新南威尔士的海岸潮间带生物群落包含牡蛎（oysters）、藤壶（barnacles）、多毛类动物（polychaetes）、藻类（algae）和其他无脊椎动物（other invertebrates）。研究人员设计了实验，验证腹足类食草动物（gastropod grazers）对无脊椎动物和藻类的群落聚集的影响。这个研究主要为了揭示db-RDA两个主要用途：1）检验食草动物在不同时间下的多变量效应一致性；2）食草动物的间接效应。

整个实验设计包含两个因子：grazing和time。Grazing因子设定了三个水平：open，caged，caged control。Time因子包含三个时间段：93年10月-94年4月，94年1月到7月，94年10月到95年1月。该实验旨在验证grazers对物种聚集的影响，并且探讨这种影响随时间的变化。换句话说，我们期望检验多变量交互项的显著性：Grazers×Time。

此处的多变量分析中，响应变量矩阵包含25个变量，其中18个为无脊椎动物，7个是大型藻类。原始数据经过转换（y^0.25），计算Bray距离矩阵，并进行NMDS排序。根据排序结果，我们发现open和cage control并没有在排序空间区分开来，表明cages本身并没有影响到群落组成。而捕食者去除的caged处理明显改变了群落组成。同时，time的效应也非常显著。ANOSIM检验结果也表明grazers和time的显著效应。尽管ANOSIM没有检验交互作用，但是排序图上也可以看出端倪，time1时期的群落组成差异（cages，open）明显低于time2和time3。我们采用db-RDA计算多变量交互项。矩阵X为表征交互项的哑变量，矩阵Xc为表征主效应的哑变量，矩阵Y为表征从Bray距离矩阵得到PCoA轴，校正方法采用lingoes。db-RDA检验出了显著的交互项：F^#=5 1.358, P=5 0.001, 999 permutations。由于grazers偏爱取食大型藻类，这些结果是可以预测的。因此，我们去除藻类植物，再进行分析，最终结果依然类似。从生态意义上来讲，这意味着，grazers效应本质在时间上是不一致的，即是群落组成的变化并不仅仅因为其对藻类的影响。而这一点在Anderson and Underwood当时的文中并没有体现出来。

案例1

4. 总结

总之，db-RDA存在诸多优势：

1）灵活选择距离矩阵类别
2）PCoA将距离信息映射到欧式空间中，进而应用于线性模型
3）负特征值校正不影响后续置换检验
4）采用多元回归的方差分析形式，RDA过程可以揭示因子及其交互效应。
5）构建类似单变量分析的多元统计数进行假设检验
6）置换检验不需要数据满足多变量正态性，对RDA模型的变量数也没有要求
7）残差置换
8）检测交互项显著性

而在db-RDA的实际应用中，需要着重注意多变量统计数的构建以及适应复杂实验设计的置换过程，同时，db-RDA对于不平衡数据、组内异质性比较敏感。

上述内容译自Legendre, P. and Anderson, M.J. (1999), DISTANCE-BASED REDUNDANCY ANALYSIS: TESTING MULTISPECIES RESPONSES IN MULTIFACTORIAL ECOLOGICAL EXPERIMENTS. Ecological Monographs, 69: 1-24. 各位看官如有兴趣，可下载原文。之后，在2001年，Anderson, M.J.及其合作者Brian H. McArdle又提出了另外一种方法，直接对距离矩阵进行分解，而不进行校正和特征值分析，即使这种距离矩阵是半定量的（后面填坑）。

依据此文，db-RDA最初的应用是探索多因子实验的交互项效应。随着技术发展，在vegan包中，dbrda()或者capscale()所执行的db-RDA逐渐演变成了一种重要的群落排序方法（个人见解）。二者对负特征值的处理不同：函数capscale()直接将虚轴舍弃，仅保留正值；函数dbrda()是直接对距离矩阵的总变异进行分解（即Brian H. McArdl文章的方法？？后面填坑）。

（未完待续）