扩增子分析视频课程推荐：https://study.omicsclass.com/index

典范冗余分析(Redundancy analysis, RDA)

多个解释变量如何解释多个响应变量---用RDA分析

分析步骤

    ①Y, X矩阵都是中心化后的响应/解释变量矩阵

    ②Y矩阵中每个响应变量与所有解释变量的多元回归，获得每个响应值变量的拟合值矩阵Y^和残差矩阵。

③将拟合值矩阵进行PCA分析，将产生一个典范特征根向量和特征向量矩阵U。

④使用矩阵U计算两套样方排序得分，一套使用中心化后的原始数据Y，获得在原始变量Y空间内的样方排序坐标，即计算YU，在vegan里面称为“样方得分”；另一套使用拟合值矩阵Y^获得在X空间内的样方排序坐标，即计算Y^U，vegan里面称为“样方约束”。

⑤将第二步获得的残差矩阵进行PCA分析获得残差非约束排序。（不属于RDA内容了）

 

RDA的排序轴实际上是解释变量的线性组合，目的是寻找能最大程度解释响应变量矩阵变差的一系列解释变量的线性组合，是被解释变量约束的排序。

 

(1)R包vegan初步运行RDA

RDA包含标准全因子RDA：

tb-RDA(Transformation-basedredundancy analysis，基于转化的冗余分析)

偏RDA(Partial canonical ordination，偏冗余分析)

db-RDA(Distance-basedredundancy analysis，基于距离的冗余分析)(CAP)

非线性RDA等多种模式。

实际的数据分析中，需根据情况选择合适的分析模式。

标准RDA分析

vegan包中运行RDA的命令为rda()， rda()的调用格式为：

 

#调用格式 1

rda(Y, X, W)

#或者 2

rda(Y~var1+var2+var3+factorA+var2*var3+Condition(var4))

 

第1种调用格式，Y为响应变量矩阵，X为解释变量矩阵，W为偏RDA分析需要输入的协变量矩阵。当不考虑协变量时，可直接输入响应变量矩阵Y和解释变量矩阵X，即“rda(Y, X)”。这种写法虽然简单，但要求解释变量矩阵或协变量矩阵中不能含有因子变量（定性变量）。

第2种写法。其中Y为响应变量矩阵，var1、var2、var3为定量解释变量，factorA为因子变量，var2*var3为变量2和变量3的交互作用项，var4为协变量。与上述第1种写法相比，不仅能够识别因子变量，还能够对变量间的交互作用加以考虑。

 

#使用全部的环境数据

rda_result <- rda(Y矩阵~., X矩阵的名称, scale = FALSE)

 

式中scale = FALSE，意为使用拟合值的协方差矩阵运算RDA；若为TRUE，则基于相关矩阵运算RDA。

 

tb-RDA

通常情况下，对于生态学中的群落物种数据，不建议直接用于RDA排序。特别是当物种丰度表中出现很多“0”值的情况下，应当首先进行Hellinger转化后，再输入至RDA中。这种首先对原始数据作一定的转化后在执行RDA的方法称为tb-RDA（基于转化的冗余分析）。

因此我们首先对示例物种丰度数据进行Hellinger转化，并使用转化后的数据执行RDA。

 

##tb-RDA

#物种数据 Hellinger 预转化（处理包含很多 0 值的群落物种数据时，推荐使用）

Y_trans <- decostand(原始Y矩阵, method = 'hellinger')

 

#使用全部的环境数据

rda_tb <- rda(Y_trans ~., X矩阵的名称, scale = FALSE)

 

#参考下文中对 RDA 结果解读的说明，可分别使用 plot(rda_result) 和 plot(rda_tb) 查看比较 Hellinger 转化前后二者的差异

 

rda()命令中，使用Hellinger转化后的物种数据作为响应变量，原始环境因子数据（数据框env）作为解释变量。式中scale = FALSE，意为使用拟合值的协方差矩阵运算RDA；若为TRUE，则基于相关矩阵运算RDA。

 

偏RDA

在控制变量矩阵W对响应变量Y的影响的前提下，将关注点集中在另一组变量矩阵X对响应变量Y的影响上，应用到RDA中即为偏RDA。

例如对于示例数据，我们期望控制土壤pH影响后（pH作为协变量），观测其它环境因素的影响。

 

##偏 RDA

#控制土壤 pH 影响后（pH 作为协变量），观测其它环境因素的影响；物种数据 Hellinger 预转化

rda_part <- rda(Y_trans ~部分X矩阵+Condition(pH), data = X矩阵的名称, scale = FALSE)

 

rda()命令中，使用Hellinger转化后的Y矩阵作为响应变量，原始X矩阵数据（作为解释变量。由于加入了协变量，此时“~”后不再直接使用“.”表示全部的环境因子，而是通过指定的方式，将所考虑的解释变量的名称以“+”连接，并将协变量（示例为pH）添加至Condition()里。式中scale同上。

db-RDA（CAP）

常规的RDA基于欧氏距离，但有时我们可能需要用到基于非欧式距离的RDA，就可以使用db-RDA（基于距离的冗余分析）来解决。db-RDA首先基于原始物种数据计算相异矩阵，作为PCoA（主坐标分析）的输入，之后将所有PCoA排序轴上的样方得分矩阵用于执行RDA。db-RDA将PCoA计算的样方得分矩阵应用在RDA中，其好处是可以基于任意一种距离测度（例如Bray-curtis距离、Unifrac距离等，而不再仅仅局限于欧氏距离）进行RDA排序，极大拓宽了RDA的适用范围，因此db-RDA在生态学统计分析中使用广泛。

例如此处我们期望基于样方群落间的Bray-curtis距离，执行db-RDA，以探索环境因子对群落物种组成的解释程度。根据上述对db-RDA计算原理的简要描述，我们可以通过以下过程来实现。

 

##db-RDA

#计算距离（以 Bray-curtis 距离为例，注意这里直接使用了原始的物种丰度矩阵）

dis_bray<- vegdist(原始Y矩阵, method = 'bray')

 

#PCoA 排序

pcoa<- cmdscale(dis_bray, k = nrow(原始Y矩阵) - 1, eig = TRUE, add = TRUE)

 

#提取 PCoA 样方得分（坐标）

pcoa_site<- pcoa$point

 

#db-RDA，使用全部的环境数据

rda_db<- rda(pcoa_site, X矩阵名称, scale = FALSE)

 

以上命令中，首先使用vegdist()计算样方群落间的Bray-curtis距离，由于不再基于欧式距离，故可以直接使用原始物种丰度表，无需Hellinger预转化。然后使用cmdscale()基于获得的Bray-curtis距离矩阵运行PCoA，其中add = TRUE表示为距离矩阵加一个常数，避免产生负的特征根，也称为Cailliez校正。（PCoA运算完毕后，提取样方排序得分（赋值给矩阵pcoa_site，主坐标矩阵视为表征数据总变差的距离矩阵）输入至rda()中作为响应变量，原始环境因子数据作为解释变量，执行RDA。

或者，vegan包中capscale()提供了直接运行db-RDA的方法，一步完成无需分解步骤。同样地，基于样方群落间的Bray-curtis距离，对示例数据执行db-RDA。

 

#或者，capscale() 提供了直接运行的方法

rda_db<- capscale(原始Y矩阵., X矩阵名称, distance = 'bray', add = TRUE)

 

在capscale()中输入原始物种组成矩阵作为响应变量，原始环境因子数据作为解释变量，“~.”表示将所有的环境因子作为解释变量带入计算；distance参数用于指定计算所依据的距离，此处使用Bray-curtis距离；add = TRUE表示为距离矩阵加一个常数，避免产生负的特征根，也称为Cailliez校正。

若是当我们在db-RDA中使用欧式距离的话，会得到和常规的RDA一样的结果。下面我们使用上述Hellinger转化后的物种数据，分别运行RDA（由于物种数据预转化了，其实应当称为tb-RDA）与使用欧氏距离的db-RDA。

非线性关系的RDA

当响应变量与解释变量明显为非线性关系（例如单峰关系，当然，若响应变量与解释变量存在明显的单峰响应模式，可使用CCA来解决），常规的线性RDA模型不适合时，可考虑使用这种方法来解决。

非线性RDA只能使用非转化的响应变量。

参考赖江山老师译作“数量生态学:R语言的应用”170-174页内容

往上可以看到四种RDA的出图结果，可以进行对比

 

(2)变量选择（前向选择）

约束排序中，有时解释变量太多，需要想办法减少解释变量的数量。一般来讲，减少解释变量的数量可能有两个并不冲突的原因：第一是寻求简约的模型，利于我们对模型的解读；第二是当解释变量过多时会导致模型混乱，例如有些解释变量之间可能存在较强的线性相关，即共线性问题，可能会造成回归系数不稳定。

方差膨胀因子

每个变量的共线性程度可以使用变量的方差膨胀因子（Variance Inflation Factor，VIF）度量，VIF是衡量一个解释变量的回归系数的方差由共线性引起的膨胀比例。如果VIF值超过20，表示共线性很严重。实际上，VIF值超过10则可能就会有共线性的问题，需要处理。

vegan包中计算约束模型中变量VIF值的命令为vif.cca()。以上文tb-RDA结果为例计算VIF。

 

#计算方差膨胀因子

vif.cca(rda_tb)

 

可以看出有几个变量的VIF值特别大（超过10表示存在共线性，超过20表示共线性很严重），所以有必要剔除一些变量。

多元回归变量筛选通常有三种模式：前向选择（forwardselection）、后向选择（backward selection）以及逐步选择（stepwise selection，也称双向选择，forward-backward selection）。其中前向选择在RDA分析中最为常用，以下将以前向选择为例简介RDA模型中的变量选择方法。

R2的校正

与多元回归的R2一样，RDA的R2是有偏的，首先，随着解释变量的增加，无论它们是否与响应变量有关，R2都会增大；另一方面，由于随机相关的存在。同样也会使得解释方差膨胀。

普通回归分析中的校正公式同样也适用于RDA中的校正：

n是样方数量，m是解释变量的数量，只要保证n-m-1>0就可以（侧面印证样本采集要多）。适用vegan包中RsquareAdj()函数获得。

r2 <- RsquareAdj(rda_tb)

 

此例校正后的R2为0.5928396.

对于我们的示例数据的tb-RDA结果，尽管在R2校正后约束模型的总体解释量有所下降（由0.7338下降至0.5928），但看起来似乎仍很可观。特别是的前两个约束轴的解释量仍然非常好：RDA1和RDA2的解释量总和为“0.43177(0.72835*0.5928)+0.09571(0.161446*0.5928)=0.52748”，即承载了响应变量矩阵52.75%的方差。

 

前向选择

在R中，解释变量的前向选择常使用vegan包中ordistep()函数、ordiR2step()函数，或者packfor包中的forward.sel()函数来完成。不同的命令执行前向选择所依据的原理是不同的。

（1）依次分别运行每个解释变量与响应变量的RDA排序。

（2）选择“最好”的解释变量。对于ordistep()，根据置换检验中F统计量的显著性水平（p值）选择最好的变量，即最小的p值。如果遇到p值相等的情况，拥有最小AIC信息准则（Akaike information criterion，AIC）的变量应该入选。最好的变量也就是最显著的变量。

（3）寻找模型中第二个、第三个、第四个……解释变量。上一步选取了只含有一个最好变量的模型，现在再加入某一个剩下的解释变量。对于ordistep()，选择偏RDA置换检验最小p值的模型，即表示新加入的变量对模型的贡献最显著。如果遇到p值相等的情况，具有最小AIC值的模型应该入选。

（4）这个过程一直持续，直到无显著的解释变量为止。即如果加入新变量的偏RDA置换检验显著性p值大于或等于α（例如，以α=0.05为准），选择过程即被终止。

 

ordiR2step()

ordiR2step()的前向选择原理与ordistep()类似，但引入全模型R2adj作为第二个终止原则（即增添了Borcard终止准则）：如果当前所选模型的R2adj达到或超过全模型的R2adj，或备选变量的偏RDA置换检验p值不显著，则变量选择将停止。

与上文一致，以下仍然以示例tb-RDA模型为例做演示。和ordistep()用法一致，除了前向选择，ordiR2step()也可以执行后向选择和逐步选择，并且可以直接将rda()的输出作为输入。

 

result1_forward<-ordiR2step(rda(tydata~1, xdata, scale = FALSE), scope = formula(result1),R2scope = adjr2$adj.r.squared, direction = 'forward', permutations = 999)

 

#简要绘制双序图比较变量选择前后结果

plot(results1,scaling = 1, main = '原始模型，I 型标尺', display = c('wa', 'cn'))

plot(result1_forward,scaling = 1, main = '前向选择后，I 型标尺', display = c('wa', 'cn'))

 

#细节部分查看

summary(rda_tb_forward_r,scaling = 1)

 

#比较选择前后校正后 R2 的差异

RsquareAdj(results1)$adj.r.squared

RsquareAdj(results1_forward)$adj.r.squared

 

ordiR2step()中，rda(phylum_hel~1,env)意为从env（环境解释变量数据集矩阵）中的第一个变量开始执行选择；rda_adj即为RDA全模型的校正R2，在上文中已经通过RsquareAdj()获得，参数scope = formula(rda_tb)即指定全模型的校正R2作为判断的依据；direction = 'forward'，前向选择；permutations = 999，999次置换。

 

全模型中共计9个环境变量，经过ordistep()前向选择后保留了其中的6个。尽管解释变量的数量少了1/3，但是我们发现模型的总解释量几乎没有改变，即变量选择后在提供了简约的模型的同时，并没有牺牲解释能力。

然后再根据排序图，或者summary统计，解读模型，观察比较变量选择前后的差异；并校正模型R2以及检验约束轴，必要时增添p值校正，最终选择合适的约束轴用于解释生态学现象等。

 

变量选择中的注意事项

尽管变量选择策略具有明显的优势，但一定切记不应盲目依赖自动选择程序在回归模型中选择相关的环境变量，因为可能得到生态学上不相关的模型，或者被忽视的其它变量组合（当前保留变量集中的部分解释变量+某些被剔除的解释变量）在某种程度上也可能比当前模型更具代表性（Legendre和Legendre，1998）。即仅通过统计学手段得到的最优变量集合可能并非代表了最具生态学意义的模型。

 

(3)RDA结果的初步查看、解读和信息提取

RDA结果总概

以上述tb-RDA排序结果为例，对vegan所得RDA结果的主要结构内容进行简要的说明。

 

#查看统计结果信息，以 I 型标尺为例

rda_tb.scaling1<- summary(rda_tb, scaling = 1)

rda_tb.scaling1

 

summary()展示结果时，根据所重要要关注的内容，可选根据I型标尺或II型标尺。如果对排序样方之间的距离更感兴趣，或者大多数解释变量为因子类型，则考虑I型标尺；如果对变量之间的相关关系更感兴趣，则考虑II型标尺。

关于I型标尺和II型标尺。这里我们展示了按I型标尺缩放后排序结果（参数scaling = 1；同理，scaling = 2展示II型标尺；若使用scaling = 0，则还原最原始的未经任何缩放的结果，在生态学数据的RDA分析中似乎不常使用。

 

summary(rda分析结果)主要有8个部分：

   1) partitioning of variance. 总方差被分为约束型和非约束型。约束型部分表示响应变量Y能被解释变量X解释的部分。相当于多元回归中未被矫正的R2，见下图。

   2) eigenvalues, their contribution to thevariance. 给出了每个典范轴的特征跟，同时也给出了累积方差解释率(约束轴)和方差承载率(非约束轴)，最终累积一定是1. 典范轴的累积解释率等于前面的constrained proportion.

   3) accumulated constrained eigenvalues. 表示每个典范轴占所能解释方差的比例。比如该例中RDA1占constrained proportion的0.73786(总和为1).

   4) scaling 1/2 for species and site scores.

--------------------------------------------------------下面的解释涉及到作图方面----------------------------------------------

   5) species scores. 双序图或三序图中代表响应变量的箭头的顶点坐标。依赖scaling的选择。是display中的“sp”。

   6) site scores-weighted sums of speciesscores. 使用响应变量矩阵Y计算获得的样方坐标。是display中的“wa”。

①提取第一列(RDA1)，第二列(RDA2)数值

②增加一列group，即把样地编号写进去

   7) site constraints: linear combinations ofconstraining variables. 使用解释变量矩阵X计算获得的样方坐标，是拟合的样方坐标。是display中的“lc”。

   8) biplot scores for constraining variables. 排序图内解释变量箭头的坐标：运行解释变量与拟合的样方坐标之间的相关分析，然后将所得相关系数转化为双序图内坐标。所有的变量包括K个水平的因子都有自己的坐标。

③提取第一列(RDA1)，第二列(RDA2)数值

 

RDA结果的典范特征系数（每个解释变量与每个典范轴之间的回归系数）没显示，用coef()函数调取。

 

ggplot第一步将site scores-weighted sums of species scores 中RDA1和RDA2设为x轴和y轴。

 

PCx与RDA模型无关

RDA1=0.73786表示该轴的解释量占约束模型总解释量的0.72835(总和为1)；而约束模型可以解释总响应变量的0.6781(总和为1).

RDA的排序图

用plot()绘制简单的三序图：三种实体，样方，响应变量，解释变量。可使用?plot.cca()查看详细说明。

 

plot(results4, scaling = 1, main="RDAplot", display=c("wa","sp","cn"))

results4_sca1<-scores(results4,choices=1:2, scaling=1,display="sp")

arrows(0,0, results4_sca1[,1], results4_sca1[,2],length=0,lty=1,col="red")

 

      display函数表示选择展示什么。Wa表示物种加权和计算的样方坐标/lc表示拟合值(解释变量的线性组合)计算的样方坐标，cn代表约束成分(解释变量)。

      scores()函数的用法，用于提取物种、样方