R语言 RDA分析(去冗余物种)

Python014

R语言 RDA分析(去冗余物种),第1张

也做了挺多次RDA分析,自己现在小结一下RDA分析流程:

就我个人而言,虚线前面都是不太经历的步骤,我一般不会主动删去样品的环境信息,因为我接触的菌群这块本来就没有什么多余的环境信息-_-||,所以我的重点放在怎么去除多余OTU或菌群上面。

一般而言,我首先会做一次差异分析,挑选有差异的OTU或菌群进行展示(phyloseq推荐使用DESeq2和edgeR,详见 Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible ),这里不是重点不在赘述。

但是差异OTU或菌群还有可能太多,RDA呈现出来密密麻麻的,调也得调好久,最后还是好不美观。

偶然间,发现envfit不仅可以评估环境因子的显著性,也可以评估物种的相关性和显著性,这为我们进一步去取冗余物种提供了条件,值得记录下来学习。

示例:

方法/步骤

1.录入原始数据。如图所示,原始数据一般采用excel表格来录入,第一列为决策单元序列,比如公司、行业等;后续各列依次是产出和投入变量,切忌产出变量一定要在投入变量前面。

2.分析效率情况。如图所示,将原始数据的格式进行统一调整之后,导入deap分析软件中,设定好相应的程序和命令后,即可运行出数据分析的结果。其中firm是公司序号,crste是技术效率,vrste是纯技术效率,scale是规模效率,最后一列是规模报酬的状态,irs是规模报酬递增,drs是规模报酬递减,-是规模报酬不变。

3

2.分析冗余情况。如图所示,DEA数据分析结果会分别给出投入、产出的冗余量,其中产出冗余数值是表示产出少了多少,而投入冗余则是表示投入多了多少。

4

4.分析参考单元。如图所示,peers表示的是可以作为效率改进参照的公司序号。由结果可见,5和13的决策单元的效率值为一,这样其他公司以此作为参照,对投入产出量进行调整,便可实现DEA有效。

冗余分析(RDA)和基于转化的冗余分析(tb-RDA)

Rao(1964)首次提出冗余分析(Redundancy analysis,RDA),从概念上讲,RDA是响应变量矩阵与解释变量矩阵之间多元多重线性回归的拟合值矩阵的PCA分析,也是多响应变量(multi-response)回归分析的拓展。在群落分析中常使用RDA,将物种多度的变化分解为与环境变量相关的变差(variation;或称方差,variance,因为RDA中变差=方差;由约束/典范轴承载),用以探索群落物种组成受环境变量约束的关系。

包含很多零值的物种多度数据在执行多元回归或其它基于欧式距离的分析方法之前必须被转化,Legendre和Gallagher(2001)提出的基于转化的RDA(Transformation-based redundancy analysis,tb-RDA)用于解决这个问题。tb-RDA在分析前首先对原始数据做一定的转化(例如Hellinger预转化包含很多零值的群落物种数据),并使用转化后的数据执行RDA。即除了第一步增添了数据转化外,其余过程均和常规的RDA相同,只是在原始数据本身做了改动,RDA算法本质未变。

RDA算法可以简要总结如下。其中矩阵Y是标准化的响应变量矩阵,X矩阵是标准化的解释变量矩阵。RDA中通常使用标准化后的解释变量,因为在很多情况下解释变量具有不同的量纲,解释变量标准化的意义在于使典范系数的绝对值(即模型的回归系数)能够度量解释变量对约束轴的贡献,解释变量的标准化不会改变回归的拟合值和约束排序的结果。在群落分析中,响应变量矩阵一般即为物种多度数据,解释变量矩阵即为环境变量数据。

(1)先将矩阵Y中的每个响应变量分别与矩阵X中的所有解释变量进行多元回归,通过回归模型获得每个响应变量的拟合值(fitted values,即在回归线上对应的值)以及残差(residuals,响应变量的观测值和拟合值之间的差值),最终得到包含所有响应变量拟合值及残差的拟合值矩阵Ŷ以及残差矩阵Yres)。

(2)对拟合值矩阵Ŷ运行PCA,得到典范特征向量(eigenvectors)矩阵U。使用矩阵U计算两套样方排序得分(坐标):一套使用中心化的原始数据矩阵Y获得在原始变量Y空间内的样方排序坐标(即计算YU,所获得的坐标称为“样方得分”,即物种得分的加权和);另一套使用拟合值矩阵Ŷ获得在解释变量X空间内的样方排序坐标(即计算ŶU,所获得的坐标称为“样方约束”,即约束变量的线性组合)。

(3)一般来讲,RDA过程执行到上步就算完成了。但一般情况下我们会同时对残差矩阵Yres运行PCA,获得残差非约束排序。非约束轴即代表了解释变量未能对响应变量作出解释的部分,严格地来说不属于RDA的范畴,但能够帮助我们获取更多信息。

Zelený博士使用仅包含一个解释变量(环境变量)的数据形象化地展示了RDA过程(原文: https://www.davidzeleny.net/anadat-r/doku.php/en:rda_cca )。

(1)执行物种spe1与环境变量env1的线性回归(由于此处示例中仅存在一个环境解释变量,故此回归为一元线性回归;当存在多解释变量时,即为多元线性回归),将回归模型拟合的物种丰度值存储在拟合值矩阵,物种丰度的残差存储在残差矩阵。见下图1中所示的过程。

(2)如此对物种组成矩阵中的所有物种重复相同的操作,最终获得包含所有物种丰度拟合值及残差的两个矩阵。见下图2中所示的两个矩阵。(1)(2)过程即形象化地展示了RDA中的回归细节部分。

(3)回归过程执行完毕后,使用PCA,在拟合值矩阵中提取约束的排序轴,并在残差值矩阵中提取非约束的轴。见下图2中所示的过程,在该示例中,由于仅有一个解释变量(环境变量env1),因此仅得到一个约束的排序轴(排序图中的垂直轴是第一个非约束轴)。

RDA排序结果产生的约束轴的数量为min[p, m, n - 1];如果同时获得非约束排序结果(即PCA),则非约束轴数量为min[p, n - 1]。其中,p为响应变量数量;m为定量解释变量数量以及定性解释变量(因子变量)的因子水平的自由度(即该变量因子水平数减1);n为排序对象数量。

从上述解释计算RDA的步骤中即可看出,RDA的排序轴实际上是解释变量的线性组合(即线性模型拟合值的排序)。换句话说,RDA的目的是寻找能最大程度解释响应变量矩阵变差的一些列的解释变量的线性组合,因此RDA是被解释变量约束的排序。约束排序与非约束排序的区别很明显:约束排序过程中解释变量矩阵控制排序轴的权重(特征根)、正交性和方向。在RDA中,排序轴解释或模拟(从统计意义上讲)依赖矩阵(响应变量)的变差,并可以检验响应变量矩阵Y与解释变量矩阵X的线性相关显著性;非约束排序PCA分析则不存在这种情况。尽管在非约束模型中,可以通过在排序后被动地加入解释变量以达到解释排序轴的目的( 详见前文 ),但此举与约束排序相比具有本质区别。

在群落分析中,对于非约束排序模型(如 PCA ),我们感兴趣的信息主要是排序图中样方和物种变量得分的相对位置、部分排序轴的相对重要性(根据特征值判断)以及排序轴的生态解释等;而对于约束排序模型(如RDA),我们通常更关注环境变量对物种组成的影响(即环境变量所能解释的变差,以及 解释程度的显著性 )、哪些环境变量对于群落结构的解释更为重要( 变量选择 )以及获知各变量或变量集解释的变差(变差分解)等。这些相关的延伸(也很重要)内容不在本文中介绍,若有需要可点击对应链接阅读。

通常情况下我们在执行RDA时(如使用R语言vegan包的rda()函数运行RDA),能够同时获得约束轴(即解释变量能够解释的部分,以约束轴呈现)和非约束轴(即解释变量未能解释的部分,多元回归的残差部分,该部分以非约束轴呈现)两部分信息,原始响应变量矩阵的总变差为约束轴解释变差和非约束轴解释变差的加和。

同前述非约束排序 PCA ,在RDA概念中,变差=方差。

约束模型解释变差反映了响应变量变化量的多少与解释变量有关,如果用比例表示,其值相当于多元回归的R2,这个解释比例值也称作双多元冗余统计(Bimultivariate redundancy statistic)。然而,类似于多元回归的未校正R2,RDA的R2也是有偏差的,需要进行校正。

同时,并非每一个约束轴都是合理有效的,还需依据置换的原理检验各约束轴的显著性,对约束轴进行取舍( 详见前文 )。因此,与非约束模型 PCA 等的非约束轴等不同,RDA约束轴的评判方法比较严格,若约束轴未通过检验,则不应被选择。(PCA只是探索性分析方法,非约束轴的选择并无严格的标准;RDA已经涉及了统计检验的过程,显著性通过p值衡量)

如上所述,也就是约束轴未能解释的,多元回归的残差部分,额外以非约束PCA轴作为呈现。对RDA进行解读时,最好同时结合约束轴和非约束轴中的信息,尽管非约束部分严格来讲不属于RDA范畴,但很多情形中仍具参考价值。

群落分析中,常通过RDA描述环境变量(解释变量)解释样方物种组成(响应变量)的差异。如果约束轴解释的变差大于非约束轴解释的变差,表明响应数据的大部分变化量均可通过解释变量作出解释,群落物种组成分布真实地由给定环境因子所影响(对于RDA结果,即二者呈现出较好的线性梯度);如果约束轴解释变差低于非约束轴解释变差,或者约束轴解释变差仅占总变差的较小比例,此时应谨慎对待,因为模型并未显示出给定环境因子能够对群落物种的组成作出有效的解释。

约束模型解释量偏低的原因可能是还有重要的解释变量尚未考虑,或是解释变量之间存在交互作用,或者归因于实际群落中物种和环境的复杂关系通常很难仅通过简单的模型有效描述出等(例如常规的RDA基于一阶线性模型,但物种和环境的关系多数情况下并非一阶线性关系,这种情况下,物种分布可能并非不受这些环境因子的约束,仅仅归因于简单的一阶线性模型无法有效描述其关系)。

RDA中每一个约束轴的特征值(eigenvalue)与特征值总和(约束轴和非约束轴特征值总和)的比例即为该轴的解释率。所有约束轴解释率总和即R2。因此,对于合理的RDA模型来讲,选定轴(通常选取特征值最高的前2-3轴用来观测)的解释量不能太低。

少数情况下,残差之间的排序或相关性(非约束轴)可能比具有良好特征的约束轴更具生物学意义。如上所述,对于RDA的残差,即额外以PCA轴的形式呈现。如果有必要,通过观测非约束空间中的样方和物种的相对位置可以帮助解读这些残差的特征。

维度选择

对于排序对象、解释变量以及响应变量的相互关系,最终通过排序图直观呈现。一般而言,我们仅选择前2-3个特征值较高(且显著)的约束轴用于观测(并尝试对其做出解释),并表示为二维/三维散点图的样式,少数情况下也会根据实际情况选择特定的排序轴(例如第二轴的趋势不明显,第三轴反而明显,因此跳过了第二轴,使用二维点图对第一、三轴可视化并做出解释;有时也会选择使用两个二维点图,分别展示并解释第一、二和三、四轴等)。有一点需要切记,就是不要试图解释太多的轴,太多的生态维度反而意义不大,正如McCune和Grace(2002)所说:“Very few ecologists have dared to venture into the uncertain waters of four or more dimensions”。

排序图表示

以R语言vegan包分析群落数据为例,排序对象(样方)、响应变量(物种)以及解释变量(环境变量)在各个约束轴中的排序结果分别报告为样方得分(site scores)、物种得分(species scores)以及解释变量得分(explanatory variable scores),投影到排序图中即表示为坐标轴上对应位置处的坐标。根据是否展示物种向量,排序图可分为双序图(仅展示样方和环境变量二者关系)和三序图(展示样方、物种及环境变量三者关系)。

RDA排序图中,样方直接在对应坐标处绘制为点。物种变量则呈现为向量,由原点(0,0)起始,指向物种得分的对应坐标处,向量的方向表示了该物种丰度增加的方向。解释变量得分(explanatory variable scores)同样以向量的形式表示在RDA排序图中,环境向量的长度表示样方物种的分布与该环境因子相关性的大小;向量与约束轴夹角的大小表示环境因子与约束轴相关性的大小,夹角小说明关系密切,若正交则不相关。