RDA、CCA的R2校正及约束轴的显著性检验概述

RDA、CCA的R²校正及约束轴的显著性检验

RDA排序涵括多元回归分析，关联程度的大小通过回归决定系数R²来反映，R²量化了线性模型中由环境变量所解释的物种组成的变差。RDA模型的解释变差反映了响应变量变化量的多少与解释变量有关，如果用比例表示，其值相当于多元回归的R²。然而，类似于多元回归的未校正R²，RDA的R²也是有偏差的，需要进行校正。

伴随着解释变量数量的增加或样方数的减少，无论解释变量是否与响应变量相关，R²一般也是会随之增加；由于随机相关的存在，解释变量的累加使得被解释方差表面上膨胀。解决方案是使用校正后的R²。此外，即使解释变量仅仅为随机生成，R²也是非零且为正（与此相比，其校正后的R²可能为0或者为负值），此时为了确定约束模型是否可信，通过蒙特卡洛置换检验（Monte Carlo permutation test）检验约束轴的显著性也是很必要的。

与此类似，CCA也需要校正R²，以及对约束轴执行检验。

校正R²

已知回归R²取决于数据集中的样方（也就是常说的样本数，以下统称为“样方”）数和解释变量的数量，样方数增加R²减小，预测变量数增加（即便这些变量仅为随机生成的）R²增大。如下图所示，即使对于随机变量也是如此，由此会导致我们错误地将实际上并不存在关联的解释变量与响应变量建立联系。因此若不进行校正，则依据R²计算所得约束轴解释率的值将不具备有效的参考价值。

图注：RDA解释变差（R²和校正后的R²）的比较。模拟生成群落数据集，得到梯度的样方个数以及解释变量（为随机变量）。由于解释变量只是随机生成的，我们不期望解释变量和响应变量产生关联。但是原始R²仍随着数据集中的样方数量的增加而减小，并随着解释变量数量的增加而增大；校正后的R²不受这两个因素的影响。

在RDA这种线性排序模型中（最小二乘多元回归），可以使用Ezekiel公式（Ezekiel，1930）获得校正后的R²：

其中n是样方数，m是解释变量的数量（更准确说，是模型的自由度）。只要m < n，即确保n – m – 1 > 0，就可以使用Ezekiel公式校正R²。一般而言，当m > n/2时，该公式可能过于保守。校正后的R²与样方数量和解释变量的数量无关。若校正后的R²趋近为零或为负值时，可以判断解释变量和响应变量无明显关联。

然而，在单峰排序模型中，Ezekiel公式返回的值被高估（并且R²对样方或解释变量数量的依赖性不能被完全消除）。Peres-Neto等（2006）提出依据置换的原理校正R²，对Ezekiel公式作变更，比较R²的观测值与使用相同数量的随机变量计算所得的期望R²（`R²_perm）的差异。需要注意，依据置换原理所得校正后的R²，每次结果会略有不同（如果置换次数设置得很高，每次结果间的差异将会小到可忽略的水平）。

图注：与原始R²相比，调整后的R²值可以达到零甚至为负值。若趋近为零，表明解释变量对响应变量的解释程度未能优于随机变量的解释程度，可以认为二者无明显关联；若为负值，表明解释变量对响应变量的解释程度甚至不及随机变量的解释程度。通常情况下，负值需被忽略不用于解释（这是很重要的，例如在变差分解中）。

备注：对于使用R包vegan中的方法获取的约束模型（RDA、CCA等），可使用RsquareAdj()获取校正后的R²。以RDA为例的示例：http://blog.sciencenet.cn/blog-3406804-1182489.html。

置换检验，检验约束轴的显著性

在经过R²校正后，约束模型的的合理性提升。接下来还需对约束模型以及各约束轴执行检验，当且仅当约束模型以及待选约束轴通过检验时，我们才具有充分的理由认为排序空间内变量间的关系是合理的。

由于生态学的数据普遍是非正态分布，因此传统的参数检验在生态学领域经常不适合，这也是目前大部分生态学数据分析方法尽可能使用置换检验的缘故。依据置换检验，将环境变量解释量的观测值与使用等数量随机变量所得的期望解释量进行比较，对零假设（零假设解释变量与响应变量无关）的真假做出判断，最终选择显著的解释变量参与解释。置换检验统计量称为伪F值（pseudo-F value），定义：

其中，n是样方数，m是约束轴的数量（或是模型自由度），SS(Ŷ)（被解释的方差）是拟合值矩阵的平方和，RSS（约束模型的残差平方和）是响应变量矩阵Y的平方和SS(Y)减去SS(Ŷ)。每个约束轴的显著性检验也是相同的原理。以检验第一约束轴为例，此时m = 1，第一约束轴的特征根即是F统计量公式的分子，分母为SS(Y)减去特征根后除以（n – m - 1）。第一轴之后的约束轴的置换检验稍复杂一些，必须以在它之前所有约束轴作为协变量重新进行偏RDA分析。如果在分析中不使用协变量，则伪F值可以用原始R²的值代替。

蒙特卡洛置换检验比较了检验统计量的观察值与在零假设下置换原始数据后所得检验统计量的期望分布，置换检验统计量伪F值（F_data）的观测值来自于原始解释变量的分析结果，零分布期望值（F_perm）为置换后解释变量的分析所得值。通过置换消除了解释变量与响应变量的联系，创建了解释变量与响应变量无关的零假设。置换需要重复多次以构建检验统计量能够与之比较的零分布。

图注：伪F值（F_data）的观察值与零分布期望值（F_perm）的比较。对于执行N次置换，将得到N个排序模型，也拥有N个F_perm值。在N个F_perm值的频数分布图中，比较观测回归方程的F_data的位置。如果F_data值比95%的F_perm值都大，即p < 0.05，此时我们可以拒绝零假设；相反如果F_data值未大于95%的F_perm值，即p > 0.05，此时我们不能拒绝零假设。

p值将通过以下公式计算：

其中n_x是置换检验统计量的零分布期望值高于观测值（F_perm > F_data）的置换次数，N是置换总次数。该公式还解释了为什么通常置换次数总设置以“9”结尾（例如99、199、499，而不使用100、200、500）：置换检验统计量的观测值也被认为是零分布的一部分，因此与期望值相加在一起（公式中的“+ 1”）。

置换检验中，可达到的最低p值取决于置换总次数（N）：

例如，当置换总次数N = 199时，可以获得的最低p值为“p_min = 1 /（199 + 1）= 0.005”，此时，将不能依据p < 0.001作出拒绝零假设判断。特别是在多重检验校正（multiple testing correction）中，例如RDA的前向选择或多个约束轴的检验，增加置换次数以确保校正后的最小p值低于选定的显著性阈值是很必要的。举例来说，499次置换能够获得的最低p值为“p_min = 1 /（499-1）= 0.002”；如果我们执行10个检验，例如对约束排序中的10个解释变量执行前向选择，并使用Bonferroni校正p值，则可获得的最小校正后p值为“p_{adj_min} = 1 /（499-1）* n_tests = 0.02”（此时n_tests = 10）。在这种情况下，若我们想根据p < 0.01作出拒绝零假设判断显然是不可能的，这个时候499次置换显得不足，需要增加置换的次数以便能确保在较低的p值下判断拒绝零假设的合理性。

检验时，首先需要对全模型执行检验，当且仅当全模型显著时才可依次对每一个约束轴执行检验，并最终确定应该在排序图中显示多少轴。否则，模型不成立或者不合理，有待思考原因并寻找其它的解决方案。对于每个轴约束轴的置换检验，推荐引入多重检验中的p值校正过程。而对于“额外获得”的残差非约束模型来讲，则无需根据置换的方法选择可能具有重要特征的非约束轴，必要时，选择使用约束模型中排序轴取舍的常用方法（例如断棍模型、Kaiser-Guttman准则等）确定残差轴更为合理。

备注：在R的vegan包中，约束模型的置换检验过程可使用anova()函数（很尴尬它就叫这个名字，R语言中方差分析的命令也是anova()……二者名称一样但用法完全不同，注意区分二者的用途）来实现；对于p值校正，将检验结果的p值输入至p.adjust()即可。以RDA为例的示例：http://blog.sciencenet.cn/blog-3406804-1182489.html。

参考资料

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

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

Ezekiel M. Methods of correlational analysis. Wiley, New York, 1930.

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.

友情链接