精华内容
下载资源
问答
  • RDA分析
    千次阅读
    2020-12-23 16:09:37

    RDA或者CCA[1]是基于对应分析发展而来的一种排序方法,将对应分析与多元回归分析相结合,每一步计算均与环境因子进行回归,又称多元直接梯度分析。此分析是主要用来反映菌群与环境因子之间关系。RDA是基于线性模型,CCA是基于单峰模型。分析可以检测环境因子、样本、菌群三者之间的关系或者两两之间的关系。

    1. RDA或CCA模型的选择原则:先用Species Table数据做DCA分析,看分析结果中Axis lengths的**轴的大小,如果大于4.0,就应该选CCA,如果3.0-4.0之间,选RDA和CCA均可,如果小于3.0,RDA的结果要好于CCA。

    2.通过bioenv函数判断环境因子与样本群落分布差异的**Pearson相关系数,通过**相关系数得到环境因子子集。

    3.将样本物种分布表与环境因子或环境因子子集分别做CCA或者RDA分析。

    4.通过类似于ANOVA的permutest分析来判断CCA或者RDA分析的显著性。

    输入:

    1、Species Table文件,由分析模块"Summary the representation of taxonomic groups"生成。

    2、样本环境因子信息文件,示例:

    !!注:输入的环境因子个数必须大于等于2,并且小于样本个数。

    Sample     Temperature     Salinity     DO   pH

    1410X       20.30        18.10        8.39 8.80

    1501X       11.00        21.10        11.27        8.12

    1504X       14.90        19.30        5.96 8.04

    1507X       27.10        17.40        7.24 8.29

    3、如果表示环境因子的射线太短或太长,通过设置如下所示的参数进行等比例放大或缩小。

    输出:

    decorana.txt:判断用CCA分析还是RDA分析的文件;

    bioenv.txt:bioenv分析判断结果;

    rda.txt:记录环境因子(Constrained)与非环境因子(Unconstrained)对样本群落分布的影响分别所占的比例;

    permutest.rda.txt:permutest显著性检验结果;

    rda.sites.xls.txt:各样本在各个排序轴上的值;

    rda.biplot.xls.txt:各环境因子在各排序轴的计算值;

    envfit.rda.txt:各环境因子对排序结果的相关性系数及显著性检验值;

    rda.cont.xls.txt:各排序轴的特征值及解释度,一般取前两个排序轴作图;

    图1. Multiple samples RDA/CCA analysis

    图2. Multiple samples RDA/CCA analysis

    图中数字表示样本名,不同颜色或形状表示不同环境或条件下的样本组;箭头表示环境因子;图2中蓝色下三角表示不同的细菌类型;物种与环境因子之间的夹角代表物种与环境因子间的正、负相关关系(锐角:正相关;钝角:负相关;直角:无相关性);由不同的样本向各环境因子做垂线,投影点越相近说明样本间该环境因子属性值越相似,即环境因子对样本的影响程度相当。

    分析模块引用了R语言(v2.12.1)vegan包(v2.0-1)中的RDA分析。

    相关文献如下所示:Sheik CS, Mitchell TW, Rizvi FZ, Rehman Y, Faisal M, et al. (2012)Exposure of Soil Microbial Communities to Chromium and Arsenic Alters Their Diversity and Structure. PLoS ONE 7(6): e40059. doi:10.1371/journal.pone.0040059.

    Legendre, P. and Legendre, L. (1998) Numerical Ecology. 2nd English ed. Elsevier.

    McCune, B. (1997) Influence of noisy environmental data on canonical correspondence analysis. Ecology 78, 2617-2623.

    Palmer, M. W. (1993) Putting things in even better order: The advantages of canonical correspondence analysis. Ecology 74,2215-2230.

    Ter Braak, C. J. F. (1986) Canonical Correspondence Analysis: a new eigenvector technique for multivariate direct gradient analysis. Ecology 67, 1167-1179.

    更多相关内容
  • 一般微生物测序实验会包含三个生物学重复,然后获得有重复的OTU table和环境因子的数据。许多文献在对OTUs进行冗余分析时都包含了...下面我们以基于距离的冗余分析(db-RDA)作为例子,尝试分析一下,看看两者的区别。
  • 冗余分析(redundancy analysis, RDA)自己之前也听过,好像是生态学研究中用的比较多,主要是用来探索环境和一些样本指标之间的关系。最近自己在看一些群体遗传相关的内容,发现RDA也可以用在群体遗传方面,比如这个...

    冗余分析(redundancy analysis, RDA)自己之前也听过,好像是生态学研究中用的比较多,主要是用来探索环境和一些样本指标之间的关系。最近自己在看一些群体遗传相关的内容,发现RDA也可以用在群体遗传方面

    ,比如这个参考链接 https://popgen.nescent.org/2018-03-27_RDA_GEA.html 就介绍了这个分析,主要研究内容自己还没有看明白:大体好像是利用芯片技术测了一些狼的基因型,同时采集了狼生活地点的环境数据,利用RDA同时分析基因型数据和环境数据。这个看的还有些模棱两可,还需要仔细看看。这个链接对应的两篇论文

    找资料的时候还找到了另外一篇论文

    3975b289b280

    image.png

    3975b289b280

    image.png

    3975b289b280

    image.png

    今天的推文重复一下这个论文里的冗余分析的代码

    首先是读入数据

    sim1.csv这个数据集1:14列是环境数据,后面都是基因型数据

    geno

    env

    geno[1:6,1:6]

    head(env)

    对基因型数据进行过滤

    这里又涉及到了最小等位基因频率这个概念

    MAF

    frequencies

    maf MAF & frequencies < (1-MAF))

    geno

    接下来就是RDA分析了

    library(vegan)

    RDA

    library(ggplot2)

    p1

    geom_line(aes(x=c(1:length(RDA$CCA$eig)), y=as.vector(RDA$CCA$eig)), linetype="dotted", size = 1.5, color="darkgrey") +

    geom_point(aes(x=c(1:length(RDA$CCA$eig)), y=as.vector(RDA$CCA$eig)), size = 3, color="darkgrey") +

    scale_x_discrete(name = "Ordination axes", limits=c(1:9)) +

    ylab("Inertia") +

    theme_bw()

    #library(robustbase)

    #install.packages("robust")

    # library("robust")

    # library(qvalue)

    rdadapt

    {

    loadings

    resscale

    resmaha

    lambda

    reschi2test

    qval

    q.values_rdadapt

    return(data.frame(p.values=reschi2test, q.values=q.values_rdadapt))

    }

    res_rdadapt

    p2

    geom_point(aes(x=c(1:length(res_rdadapt[,1])), y=-log10(res_rdadapt[,1])), col = "gray83") +

    geom_point(aes(x=c(1:length(res_rdadapt[,1]))[which(res_rdadapt[,2] < 0.1)], y=-log10(res_rdadapt[which(res_rdadapt[,2] < 0.1),1])), col = "orange") +

    xlab("SNPs") + ylab("-log10(p.values)") +

    theme_bw()

    which(res_rdadapt[,2] < 0.1)

    p3

    geom_point(aes(x=RDA$CCA$v[,1], y=RDA$CCA$v[,2]), col = "gray86") +

    geom_point(aes(x=RDA$CCA$v[which(res_rdadapt[,2] < 0.1),1], y=RDA$CCA$v[which(res_rdadapt[,2] < 0.1),2]), col = "orange") +

    geom_segment(aes(xend=RDA$CCA$biplot[,1]/10, yend=RDA$CCA$biplot[,2]/10, x=0, y=0), colour="black", size=0.5, linetype=1, arrow=arrow(length = unit(0.02, "npc"))) +

    geom_text(aes(x=1.2*RDA$CCA$biplot[,1]/10, y=1.2*RDA$CCA$biplot[,2]/10, label = colnames(env[,2:11]))) +

    xlab("RDA 1") + ylab("RDA 2") +

    theme_bw() +

    theme(legend.position="none")

    library(patchwork)

    p1/(p2+p3)

    3975b289b280

    image.png

    代码能够跑完,但是具体是什么意思还不太明白,仔细看看这些论文

    欢迎大家关注我的公众号

    小明的数据分析笔记本

    展开全文
  • 在生态学统计分析中用的PCA、CCA、RDA均能通过canoco实现
  • 群落分析的冗余分析RDA)概述

    万次阅读 2021-01-14 01:24:23
    约束排序之冗余分析(RDA)概述前篇先后简介了主成分分析(PCA)、对应分析(CA)、主坐标分析(PCoA)以及非度量多维尺度分析(NMDS)。这些排序方法均属于非约束排序,只涉及一个数据矩阵,并在低维空间中尽可能呈现原始的...

    约束排序之冗余分析(RDA)概述

    前篇先后简介了主成分分析(PCA)、对应分析(CA)、主坐标分析(PCoA)以及非度量多维尺度分析(NMDS)。这些排序方法均属于非约束排序,只涉及一个数据矩阵,并在低维空间中尽可能呈现原始的数据结构。非约束排序方法中不存在解释变量(对于物种多度数据而言,解释变量通常指代环境因素),尽管可以通过相关性或多元回归的方式被动添加至排序空间中。与此相比,约束排序则可以从排序开始直接加入解释变量进行运算,它涉及两个数据矩阵,响应变量矩阵以及解释变量矩阵。本篇继续以群落分析为例,对约束排序方法之一的冗余排序(RDA)作个简述。

    RDA的基本方法描述

    冗余分析(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算法可以简要总结如下(详细过程可参阅Legendre和Legendre(1998)“Numerical Ecology”,579-584页的内容)。其中矩阵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)

    偏冗余分析(Partial canonical ordination,偏RDA)相当于多元偏线性回归分析(Davies和Tso,1982),在实际应用中同样广泛。

    与偏线性回归相似,解释变量被分为两组:X和W,其中X表示模型中将被考虑的解释变量,W表示对响应变量Y的影响被控制的协变量。若变量矩阵W对响应变量Y的影响已知,在这种情况下,我们常期望在控制W对Y的影响的前提下,将关注点集中在另一组变量矩阵X对响应变量Y的影响上。例如,可以以气候变量X作为解释变量,土壤因子变量W作为协变量,对群落物种数据进行RDA分析。这样分析的目的是在控制土壤因子影响后,展示单独能够被气候变量线性模型解释的物种格局。

    在解释变量的前向选择过程中,偏RDA应用广泛。

    基于距离的冗余分析(db-RDA)

    尽管tb-RDA的应用拓展了RDA的适用范围,但无论常规的RDA或tb-RDA,样方或物种的降维过程实质上均以欧氏距离为基础。有时我们可能期望关注非欧式距离样方或物种关系的RDA。

    Legendre和Anderson(1999)提出的基于距离的冗余分析(Distance-based redundancy analysis,db-RDA)用于解决这个问题,并且证明RDA能够以方差分析方式分析由用于选择的任何距离矩阵。db-RDA将主坐标分析()计算的样方得分矩阵应用在RDA中,其好处是可以基于任意一种距离测度(例如Bray-curtis距离等,而不再仅仅局限于欧氏距离)进行RDA排序,因此db-RDA在生态学统计分析中被广泛使用。当然,若是我们直接计算样方群落间的欧式距离矩阵并将其输入至db-RDA中计算时,也将会得到和常规的RDA一致的结果。

    db-RDA首先基于物种多度数据(根据实际需求,选择使用原始物种数据或经过某种转化后的物种数据)计算相异矩阵(如Bray-curtis距离矩阵),作为PCoA的输入,之后将所有PCoA排序轴上的样方得分矩阵用于执行RDA,而不再使用原始的物种数据以及解释变量(环境因子数据)直接作为RDA的输入。由于在PCoA中可能会产生负特征值,必要时需要引入一些有效的校正方法。

    尽管物种信息在相异矩阵的计算过程中丢失,但主坐标矩阵依然可以视为表征数据总变差的距离矩阵,因此db-RDA结果反映了解释变量对从整个响应数据中得出的样方相似性(相异矩阵)的影响。若期望补充物种得分,方法与PCoA中过程一致,物种得分可以通过与它们所在样方得分的多度加权平均与PCoA轴建立关联而投影到最终的排序图中,用以表明响应变量(即物种信息)对PCoA排序的贡献程度。

    非线性关系的冗余分析

    RDA是通过多元线性回归分析后获得的拟合值矩阵的PCA分析,因此,很多用于多元回归的技术都可以在RDA中使用。上面所有的RDA模型都仅使用一阶的解释变量。然而,物种分布对环境梯度通常是单峰响应(如下图所示),即物种通常有最适合的生态区域。因此,对于这种情况一阶线性模型可能不太适用。

    图注:单一物种沿环境梯度的单峰响应曲线。D为该物种的最适环境梯度,E为较适环境,A、B、C则为位于该物种生态位外(对于该物种来讲此时环境很极端,难以生存或增长)。

    绘制所有响应变量与单个解释变量之间的散点图检验非线性关系非常繁琐,识别和拟合单峰响应最简单模式是在一阶函数基础上加入二阶解释变量(即二次项)并运行变量的前向选择,以保留某些一阶或二阶变量。然而,含有二阶解释变量的RDA结果很难解读,所以只有充分理由认为是非线性关系时才能使用非线性的RDA。三阶的解释变量也可以用于拟合单峰响应关系,但需要高偏态分布的响应变量。原始的多项式通常高度相关,因此建议使用正交的多项式。

    需注意,所有的非线性RDA仅使用非转化的响应变量。

    RDA结果解读

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

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

    变差(方差)解释程度

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

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

    RDA的约束轴

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

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

    RDA的非约束轴

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

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

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

    排序轴特征值和解释率

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

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

    RDA排序图

    维度选择

    对于排序对象、解释变量以及响应变量的相互关系,最终通过排序图直观呈现。一般而言,我们仅选择前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排序图中,环境向量的长度表示样方物种的分布与该环境因子相关性的大小;向量与约束轴夹角的大小表示环境因子与约束轴相关性的大小,夹角小说明关系密切,若正交则不相关。

    图注:二维散点图展示某RDA结果的前两个约束轴。其中A图为双序图,展示样方和环境变量;B图为三序图,展示样方、物种和环境变量。

    PCA中有I、II型标尺的选择,RDA中也是如此。

    无论RDA双序图或三序图,均需要同时展示对象和变量,但没有同时可视化对象和变量的最优化方法。通常根据关注侧重点的不同,考虑使用不同的标度(或称尺度、标尺)展示排序图,以重点反映关键的信息。得分集作为RDA输出的典型特征,将根据使用的标度而变化。一般有两种标尺模式,不同模式的排序图有不同的解读方式。如果对排序样方之间的距离更感兴趣,或者大多数解释变量为因子类型,则考虑I型标尺;如果对变量之间的相关关系更感兴趣,则考虑II型标尺。

    在两种标尺的RDA排序图内,对于双序图,由于未展示物种信息故无需考虑样方和物种关系的解读方式;对于三序图,样方和物种箭头的解读同PCA,将样方对象点垂直投影到物种变量箭头上位置表示该物种在该样方内数值在所有样方内的排序位置。但是,解释变量的箭头或因子变量的形心的解读另有规则。

    以下是对RDA I型标尺和II型标尺的简要描述(详细过程可参阅Legendre和Legendre(1998)“Numerical Ecology”,585-587页的内容)。

    RDA的I型标尺

    I型标尺关注的是对象之间的关系。矩阵U的特征向量(代表了约束轴中的响应变量得分,参见开头“RDA算法简述”,下述X、Y、Ŷ等同述),被标准化为单位长度。X空间内的样方得分由公式Z=ŶU获得,这些向量具有相同的方差λk。Y空间内的样方得分由公式F=YU获得,这些向量的方差通常显著高于λk,因为Y同时涵括了拟合值和残差值因此其总方差高于Ŷ(仅为拟合值)的方差。通过双序图表示矩阵Z和U,或F和U时,特征向量与样方得分矩阵的结果较好地还原了原始矩阵:ZU’=Ŷ或FU’=Y。对于定量解释变量x与样方得分拟合值的相关性,乘以(λk/Y的总方差)1/2后在各排序轴中呈现,其中λk为约束轴k的特征值。在这种标尺中,各约束轴中的样方得分具有不同的方差。

    I型标尺的RDA排序图主要展示如下信息:

    (1)样方点垂直投影到响应变量或定量解释变量的箭头或延长线上,投影点近似于该样方内该响应变量或解释变量的数值沿着变量的位置;

    (2)响应变量与解释变量箭头之间的夹角反映了它们之间的相关性(但响应变量之间的夹角无此含义);

    (3)定性解释变量的形心与响应变量(物种)箭头之间的解读如同样方点与响应变量之间的解读(因为定性解释变量的形心也是一组样方的形心);

    (4)定性解释变量的形心之间或形心与样方点之间的距离近似他们之间的欧式距离。

    RDA的II型标尺

    II型标尺中,每个特征向量被标准化为特征根的平方根,关注的是变量之间的关系。

    通过UɅ1/2转化,将矩阵U的特征向量(代表了约束轴中的响应变量得分,参见开头“RDA算法简述”,下述X、Y、Ŷ等同述)标准化为(λk)1/2。X空间内的样方得分通过Z=ŶU获得后,再经ZɅ-1/2转化缩放为单位方差。同样地,Y空间内的样方得分通过F=YU获得后,再经FɅ-1/2转化缩放,并且所得向量的方差通常显著高于由X所得向量的方差,原因同上述“I型标尺”中的描述。通过双序图表示矩阵Z和U,或F和U时,特征向量与样方得分矩阵的结果较好地还原了原始矩阵:ZU’=Ŷ或FU’=Y。定量解释变量x与样方得分拟合值的相关性直接呈现在各排序轴中。

    II型标尺的RDA排序图主要展示如下信息:

    (1)将样方点垂直投影到响应变量或定量解释变量的箭头或延长线上,投影点位置近似于该响应变量或解释变量在该样方内的数值;

    (2)响应变量与解释变量箭头之间的夹角反映它们之间的相关性,响应变量之间和解释变量之间也同样解读;

    (3)定性解释变量的形心与响应变量箭头之间的解读如同样方点与响应变量之间的解读,可以通过投影点判断;

    (4)定性解释变量的形心之间或形心与样方之间的距离不再近似欧式距离。

    看完描述是不是有点懵?接下来通过两个图来简要说明下对RDA排序图的解读吧。

    (1)a图,存在排序样方(样本)I和II,解释变量(环境变量)1,探究I或II与1的关系时,将I或II垂直投影在1的向量(箭头)上,根据交叉点的位置判断变量1在I或II中的值。交叉点越靠近该变量向量的正方向,则表明所对应的样方中,该变量的数值越大。例如,假设变量1为土壤碳含量,样方I投影在1的正方向,样方II投影在1的负方向上(图中红色虚线反向延长线部分),两个交叉点相比较,I与1的交叉点更位于1延伸方向,因此可知I中的土壤碳含量要比II中的土壤碳含量要高。

    若1为响应变量(物种变量),观察方法同样适用。例如变量1为物种species1,同样据此可判断物种species1在I中的丰度高于在II中的丰度。

    注:无论I型标尺或II型标尺,均可据此判断变量在样方中的相对数值大小。

    (2)b图,根据向量(箭头)夹角判断变量间的相关性。∠a接近90°,即接近正交,表明变量1和2之间的相关性很小,二者相互之间几乎不存在影响。∠b小于90°,夹角为锐角,表明变量2和3之间存在正相关;锐角角度越小,则正相关性越大。∠c大于90°,夹角为钝角,表明变量3和4之间存在负相关;钝角角度越大,则负相关性越大。

    注:对于I型标尺,仅可据此观测解释变量与响应变量间的相关性;对于II型标尺,既可以据此观测解释变量与响应变量间的相关性,也可以观测解释变量之间、或响应变量之间的相关性。

    (3)对于因子类型的解释变量5(定性变量,非数值型变量),在图中以点表示而非以向量表示,探究变量5与其它变量间的相关性时需要根据投影判断。例如,变量5垂直投影在变量4的正方向,表明与变量4存在正相关;投影在变量2的负方向,表明与变量2存在负相关;相关性的大小,可以通过垂线交叉点与原点(0,0)的距离来表示。

    注:对于I型标尺,仅能据此观测定性解释变量与响应变量间的相关性;对于II型标尺,既可以据此观测定性解释变量与响应变量间的相关性,也可以观测其与定量解释变量之间的相关性。

    (4)若为I型标尺,还可根据图中样方点之间的距离判断样方群落之间的相似性。两个样方距离越近,则群落相似性越大;反之越低。

    (5)此外,还可通过比较解释变量(环境变量)向量在约束轴上投影的相对长度,判断环境变量对群落特征的贡献度。例如在图b中,将变量2和变量1均投影至RDA2轴,此时变量2的投影长度相对更长,表明变量2比变量1对RDA2轴形成的贡献更大。

    解释变量向量与约束轴夹角的大小同样具有意义,表示解释变量与约束轴相关性的大小,夹角小说明关系密切,若正交则不相关。例如在图b中,变量2的向量与RDA2轴的夹角比与RDA1轴的夹角更小,表明变量2与RDA2的关联程度比与RDA1的关联程度要高,即相较之下变量2更贡献于RDA2轴。

    注:无论I型标尺或II型标尺,均可据此判断。

    坐标位置的特殊调整

    此外,在绘制排序图时,通常会考虑将排序对象或变量的得分(即坐标)乘以一个常量,即对排序图中的点或向量坐标同比例放大或缩小一定的数值后展示,以产生易于解读的生态排序结果(可参见Legendre和Legendre,1998,404页)。这种做法可能与严格的数学定义相违背,但是对于解答生物学问题通常是没有影响的。

    也可参考前述博文的末尾所述。

    其它内容

    本篇关于RDA的基础部分,大致就先介绍这些吧,对于初步理解绝对够用了。

    RDA作为多元统计分析中的常见方法,其使用范围广泛,因此目前可用于运行RDA的软件也有很多可供选择。例如,生态学统计分析的常用R包vegan中,提供了关于RDA模型计算、I型或II型标尺选择、R2校正、置换检验、变量选择、变差分解等一系列的方法,因此下一篇将简介R包vegan的RDA分析。

    参考资料

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

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

    Davies P T, Tso K S. Procedures for Reduced-Rank Regression. Journal of the Royal Statistical Society. Series C (Applied Statistics), 1982, 31(3):244-255.

    GUSTA ME Blog:https://mb3is.megx.net/gustame/constrained-analyses/rda

    Legendre P, Anderson M J . Distance-Based Redundancy Analysis: Testing Multispecies Responses in Multifactorial Ecological Experiments. Ecological Monographs, 1999, 69(1):1-24.

    Legendre P, Gallagher E D . Ecologically Meaningful Transformations for Ordination of Species Data. Oecologia, 2001, 129(2):271-280.

    Legendre P, Legendre L. Numerical Ecology. Second English edition. Developments in Environmental Modelling, 1998, 20, Elsevier

    Mccune B P, Grace J B. Analysis of Ecological Communities. Journal of Experimental Marine Biology and Ecology, 2002, 289(2).

    Rao C R. The use and interpretation of principal component analysis in applied research. Sankhya, A, 1964, 26(4):329-358.

    展开全文
  • 数量生态学笔记||冗余分析RDA

    千次阅读 2021-01-30 07:01:21
    上一节数量生态学笔记||冗余分析(RDA)概述中,我们回顾了RDA的计算过程,不管这个过程我们有没有理解透彻,我希望你能知道的是:RDA是响应变量矩阵与解释变量之间多元多重线性回归的拟合值矩阵的PCA分析。...

    de1bd2531883

    上一节数量生态学笔记||冗余分析(RDA)概述中,我们回顾了RDA的计算过程,不管这个过程我们有没有理解透彻,我希望你能知道的是:RDA是响应变量矩阵与解释变量之间多元多重线性回归的拟合值矩阵的PCA分析。本节我们就是具体来看一个RDA的分析案例,来看看里面的参数以及结果的解读。

    # CHAPTER 6 - CANONICAL ORDINATION

    # ********************************

    # 载入所需程序包

    library(ade4)

    library(vegan)

    library(packfor)# 可在http://r-forge.r-project.org/R/?group_id=195下载,但是好像在R 3.5.1上加载不了,所以这篇我用R3.4来做的。packfor已经不用,函数都搬到adespatial

    # 如果是MacOS X系统,packfor程序包内forward.sel函数的运行需要加载

    # gfortran程序包。用户必须从"cran.r-project.org"网站内选择"MacOS X",

    # 然后选择"tools"安装gfortran程序包。

    rm(list = ls())

    setwd("D:\\Users\\Administrator\\Desktop\\RStudio\\数量生态学\\DATA")

    library(MASS)

    library(ellipse)

    library(FactoMineR)

    # 附加函数

    source("evplot.R")

    source("hcoplot.R")

    # 导入CSV数据文件

    spe

    env

    spa

    # 删除没有数据的样方8

    spe

    env

    spa

    # 提取环境变量das(离源头距离)以备用

    das

    # 从环境变量矩阵剔除das变量

    env

    # 将slope变量(pen)转化为因子(定性)变量

    pen2

    pen2[env$pen <= quantile(env$pen)[4]] = "steep"

    pen2[env$pen <= quantile(env$pen)[3]] = "moderate"

    pen2[env$pen <= quantile(env$pen)[2]] = "low"

    pen2

    table(pen2)

    # 生成一个含定性坡度变量的环境变量数据框env2

    env2

    env2$pen

    # 将所有解释变量分为两个解释变量子集

    # 地形变量(上下游梯度)子集

    envtopo

    names(envtopo)

    #水体化学属性变量子集

    envchem

    names(envchem)

    # 物种数据Hellinger转化

    spe.hel

    使用vegan包运行RDA

    vegan包运行RDA有两种不同的模式。第一种是简单模式,直接输入用逗号隔开的数据矩阵对象到rda()函数:

    math?formula=simpleRDA%20%3D%20rda(Y%2CX%2CW)

    式中

    math?formula=Y为响应变量矩阵,

    math?formula=X为解释变量矩阵,

    math?formula=W为偏RDA分析需要的协变量矩阵。

    此公式有一个缺点:

    math?formula=Y%2CW不能有因子变量(定性变量)。如果有因子变量,建议使用第二种模式:

    math?formula=formulaRDA%3C-rda(Y%20%5Csim%20var1%20%2B%20factorA%20%2B%20var2*var3%2BCondition(var4)%2Cdata%3DXWdata%20)

    式中,

    math?formula=Y为响应变量矩阵。解释变量矩阵包括定量变量(var1)、因子变量(factorA)以及变量2和变量3的交互作用项,协变量(var4)被放到Condition()里。所用的数据都放在XWdata的数据框里。

    这个公式与lm()函数以及其他回归函数一样,左边是响应变量,右边是解释变量。

    # 基于Hellinger 转化的鱼类数据RDA,解释变量为对象env2包括的环境变量

    # 关注省略模式的公式

    spe.rda

    summary(spe.rda) # 2型标尺(默认)

    #这里使用一些默认的选项,即 scale=FALSE(基于协方差矩阵的RDA)和#scaling=2

    RDA结果的摘录:

    RDA formula :

    Call:

    rda(formula = spe.hel ~ alt + pen + deb + pH +

    dur + pho + nit + amm + oxy + dbo, data = env2)

    方差分解(Partitioning of variance):总方差被划分为约束和非约束两部分。约束部分表示响应变量

    math?formula=Y矩阵的总方差能被解释变量解释的部分,如果用比例来表示,其值相当于多元回归的

    math?formula=R%5E2。在RDA中,这个解释比例值也称作双多元冗余统计。然而,类似多元回归的未校正的

    math?formula=R%5E2,RDA的

    math?formula=R%5E2是有偏差的,需要进一步校正。

    Partitioning of variance:

    Inertia Proportion

    Total 0.5025 1.0000

    Constrained 0.3654 0.7271

    Unconstrained 0.1371 0.2729

    特征根以及对方差的贡献率(Eigenvalues, and their contribution to the variance ):当前这个RDA分析产生了12个典范轴(特征根用RDA1 至RDA12表示)和16个非约束轴(特征根用PC1至PC16表示)。输出结果不仅包含每轴特征根同时也给出累积方差解释率(约束轴)或承载轴(非约束轴),最终的累计值必定是1.12 个典范轴累积解释率也代表响应变量总方差能够被解释变量解释的部分。

    两个特征根的重要区别:典范特征根RDAx是响应变量总方差能够被解释变量解释的部分,而残差特征根RCx响应变量总方差能够被残差轴解释的部分,与RDA无关。

    Eigenvalues, and their contribution to the variance

    Importance of components:

    RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10 RDA11 RDA12 PC1 PC2 PC3

    Eigenvalue 0.2281 0.0537 0.03212 0.02321 0.008707 0.007218 0.004862 0.002919 0.002141 0.001160 0.0009134 0.0003406 0.04580 0.02814 0.01529

    Proportion Explained 0.4539 0.1069 0.06393 0.04618 0.017328 0.014363 0.009676 0.005809 0.004260 0.002308 0.0018176 0.0006778 0.09115 0.05600 0.03042

    Cumulative Proportion 0.4539 0.5607 0.62467 0.67085 0.688176 0.702539 0.712215 0.718025 0.722284 0.724592 0.7264100 0.7270878 0.81824 0.87424 0.90466

    PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16

    Eigenvalue 0.01399 0.009841 0.007676 0.004206 0.003308 0.002761 0.002016 0.001752 0.0009851 0.0005921 0.0004674 0.0002127 0.0001004

    Proportion Explained 0.02784 0.019584 0.015276 0.008371 0.006583 0.005495 0.004013 0.003486 0.0019604 0.0011783 0.0009301 0.0004233 0.0001998

    Cumulative Proportion 0.93250 0.952084 0.967360 0.975731 0.982314 0.987809 0.991822 0.995308 0.9972684 0.9984468 0.9993768 0.9998002 1.0000000

    累积约束特征根(Accumulated constrained eigenvalues)表示在本轴以及前面所有轴的典范轴所能解释的方差占全部解释方差的比例累积。

    Accumulated constrained eigenvalues

    Importance of components:

    RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10 RDA11 RDA12

    Eigenvalue 0.2281 0.0537 0.03212 0.02321 0.008707 0.007218 0.004862 0.002919 0.002141 0.001160 0.0009134 0.0003406

    Proportion Explained 0.6242 0.1470 0.08792 0.06352 0.023832 0.019755 0.013308 0.007990 0.005859 0.003174 0.0024999 0.0009322

    Cumulative Proportion 0.6242 0.7712 0.85913 0.92265 0.946483 0.966237 0.979545 0.987535 0.993394 0.996568 0.9990678 1.0000000

    物种得分(Species scores)双序图和三序图内代表响应变量的箭头的顶点坐标。与PCA相同,坐标依赖标尺Scaling的选择。

    Scaling 2 for species and site scores

    * Species are scaled proportional to eigenvalues

    * Sites are unscaled: weighted dispersion equal on all dimensions

    * General scaling constant of scores: 1.93676

    Species scores

    RDA1 RDA2 RDA3 RDA4 RDA5 RDA6

    CHA 0.13383 0.11623 -0.238180 0.018611 0.043221 -0.029737

    TRU 0.64238 0.06648 0.123713 0.181572 -0.009691 0.029793

    VAI 0.47475 0.07015 -0.010218 -0.115369 -0.045317 -0.030033

    LOC 0.36260 0.06972 0.041240 -0.190586 -0.046881 0.006448

    OMB 0.13079 0.10709 -0.239224 0.043603 0.065881 0.003458

    BLA 0.06587 0.12475 -0.216900 -0.004157 0.021793 -0.004195

    HOT -0.17417 0.06778 -0.008426 -0.016419 -0.079730 0.044706

    TOX -0.12683 0.16052 -0.035733 -0.016087 -0.089768 -0.001880

    VAN -0.07963 0.04200 0.007636 -0.059179 -0.033596 -0.121440

    CHE -0.10903 -0.17552 -0.090099 -0.168373 0.019444 0.008745

    BAR -0.18528 0.21154 -0.073087 -0.006879 -0.012995 0.040484

    SPI -0.16064 0.15513 -0.014309 -0.002488 -0.060810 0.011045

    GOU -0.20537 0.02484 -0.007973 -0.017742 -0.049137 -0.096231

    BRO -0.10734 0.02848 0.090055 0.012324 0.075184 -0.057088

    PER -0.09164 0.10506 0.070393 -0.057443 0.013870 -0.009906

    BOU -0.20907 0.16002 0.025500 0.012078 -0.011477 0.022035

    PSO -0.22799 0.11121 0.018800 -0.009474 -0.027431 0.024517

    ROT -0.16098 0.01348 0.041628 0.032398 0.054117 -0.094582

    CAR -0.17384 0.14901 0.022262 0.009534 0.004991 -0.005396

    TAN -0.14025 0.10632 0.078290 -0.122627 0.054162 0.031256

    BCO -0.18594 0.12222 0.053881 0.026170 0.044015 0.014577

    PCH -0.14630 0.08894 0.061880 0.034763 0.083530 0.004396

    GRE -0.30881 0.01606 0.039366 0.029254 -0.011141 -0.052412

    GAR -0.31982 -0.16601 -0.018225 -0.115454 0.054341 0.064772

    BBO -0.23897 0.09090 0.051627 0.010224 0.007004 0.036497

    ABL -0.43215 -0.22639 -0.108190 0.138807 -0.083920 0.008460

    ANG -0.19442 0.14149 0.033659 0.017387 0.008110 0.017638

    样方得分(Site scores (weighted sums of species scores))物种得分的加权和:使用响应变量矩阵

    math?formula=Y计算获得的样方坐标。

    Site scores (weighted sums of species scores)

    RDA1 RDA2 RDA3 RDA4 RDA5 RDA6

    1 0.40151 -0.154306 0.55539 1.600773 0.191866 0.916893

    2 0.53523 -0.025084 0.43389 0.294615 -0.518456 0.458860

    3 0.49430 -0.014605 0.49409 0.169038 -0.246166 0.163421

    4 0.33452 0.001173 0.51626 -0.321009 0.088716 -0.219837

    5 0.02794 -0.194357 0.44612 -0.559210 0.853768 -1.115654

    6 0.24422 -0.130778 0.41372 -0.696264 0.181514 -0.273473

    7 0.46590 -0.125982 0.31674 -0.137834 -0.548635 -0.061703

    9 0.03662 -0.605060 -0.07022 -1.260916 0.669108 1.164986

    10 0.31381 -0.198654 0.10764 -0.635139 -0.741448 -0.990236

    11 0.48116 -0.039598 -0.37851 0.181924 0.221494 0.254511

    12 0.49162 0.014263 -0.37983 0.163103 0.223730 0.324672

    13 0.49848 0.212367 -0.67408 0.518823 0.400091 0.221622

    14 0.38202 0.229538 -0.75771 0.223651 0.515712 -0.139740

    15 0.28739 0.218713 -0.71887 -0.210821 0.176392 -0.553185

    16 0.09129 0.400192 -0.34443 -0.376097 -0.366868 -0.575230

    17 -0.05306 0.423994 -0.41009 -0.188492 -0.726152 0.151876

    18 -0.14185 0.385926 -0.36814 -0.217143 -0.644298 -0.001052

    19 -0.28204 0.275528 -0.01877 -0.371457 -0.691725 -0.062230

    20 -0.39683 0.209468 0.11547 -0.177972 -0.387121 0.048690

    21 -0.42851 0.278256 0.22010 -0.005993 -0.027083 -0.042209

    22 -0.46553 0.251819 0.22784 0.040192 0.152965 0.032185

    23 -0.28123 -1.145599 -0.50543 0.300015 -0.004403 1.157206

    24 -0.40893 -0.752909 -0.26785 0.428851 -0.645606 0.643084

    25 -0.35205 -0.770380 -0.12186 0.459170 0.078892 -1.725973

    26 -0.46923 0.093958 0.23058 0.107865 0.310432 0.132556

    27 -0.47021 0.213521 0.24887 0.084219 0.331685 0.125439

    28 -0.47266 0.233922 0.27053 0.105776 0.381436 0.093719

    29 -0.37457 0.393260 0.10423 0.202115 0.282621 0.021834

    30 -0.48932 0.321417 0.31431 0.278218 0.487541 -0.151031

    样方约束——解释变量的线性组合(Site constraints (linear combinations of constraining variables)):使用解释变量矩阵

    math?formula=X计算获得的样方坐标,是拟合的(fitted)样方坐标。

    Site constraints (linear combinations of constraining variables)

    RDA1 RDA2 RDA3 RDA4 RDA5 RDA6

    1 0.55135 0.002395 0.47774 0.626878 -0.210700 0.31511

    2 0.29737 0.105715 0.64862 0.261161 -0.057741 0.09322

    3 0.36834 -0.185376 0.59788 0.324212 -0.002385 0.31090

    4 0.44348 -0.066086 0.33260 -0.344463 -0.279591 -0.37079

    5 0.27004 -0.366721 0.17992 -0.453274 0.716614 -0.06545

    6 0.37148 -0.255624 0.40847 0.217259 0.023374 0.34360

    7 0.53874 -0.179999 0.06845 -0.424896 0.024884 -0.33454

    9 -0.04438 -0.362632 0.12371 -1.180662 0.348188 0.26352

    10 0.16289 -0.154212 0.22252 -0.241425 -0.573048 -0.02867

    11 0.29912 0.176150 -0.08233 0.003924 0.164806 -0.44603

    12 0.36843 0.197492 -0.41095 0.300566 -0.053922 -0.01139

    13 0.42626 0.190107 -0.59764 0.100988 0.118714 -0.21022

    14 0.34373 0.134362 -0.80378 0.063879 0.665797 0.48126

    15 0.21385 0.237182 -0.56341 -0.001099 -0.028564 0.01655

    16 0.03056 0.352192 -0.12110 -0.202316 0.058413 -0.43542

    17 -0.10499 0.178587 -0.26925 0.046988 -0.608314 0.21237

    18 -0.11204 0.221631 -0.24024 -0.302957 -0.251346 -0.01448

    19 -0.05479 0.311860 -0.30701 0.010366 -0.481829 -0.12855

    20 -0.25684 0.303770 -0.06768 -0.036587 -0.562578 0.13698

    21 -0.39177 0.196355 0.01877 -0.281086 -0.383524 0.39310

    22 -0.21361 0.180414 0.01066 0.074301 -0.036849 -0.02429

    23 -0.21654 -1.016853 -0.57298 0.548175 0.182594 0.51443

    24 -0.52578 -0.645438 -0.11182 -0.240149 -0.512492 0.32420

    25 -0.38886 -0.867381 -0.08079 0.482839 -0.106743 -1.21305

    26 -0.48440 0.031510 0.14065 -0.114545 0.425712 -0.17989

    27 -0.61221 0.138191 0.32316 -0.015795 0.232397 0.38288

    28 -0.46921 0.459843 0.22002 0.078870 0.278747 -0.36504

    29 -0.38450 0.344487 0.20622 0.353008 0.504693 0.17747

    30 -0.42572 0.338080 0.24960 0.345839 0.404692 -0.13777

    解释变量双序图得分(Biplot scores for constraining variables):排序图内解释解释变量箭头的坐标,按照下面的过程获得:运行解释变量与拟合的样方坐标之间的相关分析,然后将所有相关系数转化为双序图内坐标。所有的变量包括

    math?formula=k个水平的因子口可以有自己的坐标对因子变量在排序轴的坐标,用各个因子的形心表示更合适。

    Biplot scores for constraining variables

    RDA1 RDA2 RDA3 RDA4 RDA5 RDA6

    alt 0.8239 -0.203257 0.46604 -0.16936 0.003229 0.10407

    penmoderate -0.3592 -0.008729 -0.21727 -0.18278 0.157934 0.50094

    pensteep 0.2791 0.156027 -0.06876 0.01927 0.176390 -0.15469

    penvery_steep 0.6129 -0.148496 0.45379 0.03618 -0.191046 -0.04715

    deb -0.7770 0.254952 -0.17470 0.30995 0.227580 -0.11938

    pH 0.1023 0.178431 -0.30131 0.03959 0.298584 0.04848

    dur -0.5722 0.044963 -0.56040 -0.14813 0.275617 -0.24527

    pho -0.4930 -0.650488 -0.19868 0.29286 0.055893 -0.39345

    nit -0.7755 -0.203836 -0.23285 0.19744 -0.078849 -0.35073

    amm -0.4744 -0.687577 -0.16648 0.28470 -0.051768 -0.33852

    oxy 0.7632 0.575528 -0.16425 0.08026 -0.136143 0.13748

    dbo -0.5171 -0.791805 -0.15652 0.22064 0.075568 -0.09105

    因子解释变量形心(Centroids for factor constraints):因子变量各个水平形心点的坐标,即每个水平所用标识为一的样方的形心。

    Centroids for factor constraints

    RDA1 RDA2 RDA3 RDA4 RDA5 RDA6

    penlow -0.2800 0.005530 -0.09025 0.07614 -0.07860 -0.18390

    penmoderate -0.2093 -0.005086 -0.12660 -0.10650 0.09203 0.29189

    pensteep 0.1965 0.109867 -0.04842 0.01357 0.12420 -0.10892

    penvery_steep 0.3908 -0.094679 0.28933 0.02307 -0.12181 -0.03006

    在rda()函数中大家感兴趣的典范特征系数(即每个解释变量与每个典范轴之间的回归系数),可用coef()函数获得:

    #如何从rda()输出结果中获得典范系数

    coef(spe.rda)

    alt 0.0004482548 7.559499e-05 0.0005205825 0.0003883845 0.001857206 -6.313946e-05 -0.001355362 0.001120849 -0.0002530083 0.001189659

    penmoderate -0.0123961693 -1.660194e-02 0.0160069104 -0.0278562534 0.276128846 1.310695e-01 -0.022918427 0.018830063 -0.3113354204 -0.278006278

    pensteep 0.0478390238 4.893302e-02 0.1022577908 0.1347997771 0.393812929 -1.795824e-01 0.046319741 0.123642821 0.0963820046 -0.447099975

    penvery_steep 0.0180005587 -5.691933e-02 0.2322637617 0.1002359565 0.036814635 -1.742060e-01 0.517299284 0.067564014 -0.2262450630 -0.590022907

    deb -0.0014063766 4.440084e-03 0.0089298486 0.0164715901 0.013318449 2.705757e-03 -0.002419805 0.010394632 -0.0006430624 0.004427139

    pH 0.0188227657 -3.163592e-02 -0.0482021704 0.1141647498 0.412886847 1.091403e-01 0.139806409 -0.436295510 -0.0215841003 -0.904063764

    dur 0.0025580344 -1.955496e-03 -0.0065901935 -0.0093556696 0.005228707 -6.098098e-03 0.002195518 0.010536248 -0.0006844877 0.003353110

    pho 0.1031541920 4.583584e-02 -0.1000153096 -0.1050243435 0.422991234 -3.694132e-01 0.035874664 -0.701043138 -0.2865315085 0.245917738

    nit -0.0123824749 1.041485e-01 0.0625396719 0.0774218297 0.234401221 -3.541252e-02 -0.240428544 0.128403162 -0.0686364968 0.113090041

    amm -0.1084411839 -4.407786e-01 0.0057247742 0.0538542716 -1.812468883 4.798631e-03 0.281937862 1.068013480 0.3084215704 -1.217501183

    oxy 0.0686692124 1.980446e-02 -0.0894153251 0.1200884061 0.032052566 3.880445e-02 -0.058026043 0.061374900 -0.0196444146 0.089881042

    dbo 0.0108322463 -2.696114e-02 -0.0253225230 0.0745175780 0.067082880 9.276548e-02 -0.019719504 0.047865971 0.0359365102 0.065416035

    RDA11 RDA12

    alt 0.0006826822 0.0009471677

    penmoderate 0.0398080898 -0.2974896027

    pensteep 0.2445035939 -0.3475535793

    penvery_steep 0.2457103975 -0.1848717482

    deb -0.0022565029 0.0064858596

    pH 0.0696045266 0.5756301035

    dur 0.0007758175 0.0062181193

    pho -0.0015544897 -0.6309008260

    nit 0.3983543655 0.0942246573

    amm -1.5964107514 0.8979015239

    oxy 0.0627415675 0.0258937429

    dbo 0.1113928484 0.0403158432

    提取。解读和绘制vegan包输出的RDA结果

    校正

    math?formula=R%5E2

    math?formula=R%5E2_%7Badj%7D%3D1-%5Cfrac%7Bn-1%7D%7Bn-m-1%7D%5Cleft(%201-R%5E2%20%5Cright)

    # 提取校正R2

    # **********

    # 从rda的结果中提取未校正R2

    (R2

    [1] 0.7270878

    # 从rda的结果中提取校正R2

    (R2adj

    [1] 0.5224036

    #可以看出,校正R2总是小于R2。校正R2作为被解释方差比例的无偏估计,后#面的变差分解部分所用的也是校正R2。

    # RDA三序图

    现在绘制RDA的排序图。如果一张排序图中有三个实体:样方、响应变量、解释变量,这种排序图称为三序图(triplot)为了区分响应变量和解释变量,定量解释变量用箭头表示,响应变量用不带箭头的线表示。

    # RDA三序图

    # **********

    # 1型标尺:距离三序图

    plot(spe.rda, scaling=1, main="RDA三序图:spe.hel~env2 - 1型标尺- 加权和样方坐标")

    #此排序图同时显示所有的元素:样方、物种、定量解释变量(用箭头表示)

    #和因子变量的形心。为了与定量解释变量区分,物种用不带箭头的线表示。

    spe.sc

    arrows(0, 0, spe.sc[, 1], spe.sc[, 2], length=0, lty=1, col="red")

    de1bd2531883

    plot(spe.rda, main="RDA三序图:spe.hel~env2 - 2型标尺- 加权和样方坐标")

    spe2.sc

    arrows(0, 0, spe2.sc[, 1], spe2.sc[, 2], length=0, lty=1, col="red")

    de1bd2531883

    # 样方坐标是环境因子线性组合

    # 1型标尺

    plot(spe.rda, scaling=1, display=c("sp", "lc", "cn"),

    main="RDA三序图:spe.hel~env2 - 1型标尺- 拟合的样方坐标")

    arrows(0, 0, spe.sc[, 1], spe.sc[, 2], length=0, lty=1, col="red")

    de1bd2531883

    # 2型标尺

    plot(spe.rda, display=c("sp", "lc", "cn"),

    main="RDA三序图:spe.hel~env2 - 2型标尺- 拟合的样方坐标")

    arrows(0, 0, spe2.sc[,1], spe2.sc[,2], length=0, lty=1, col="red")

    de1bd2531883

    RDA 结果的置换检验

    # RDA所有轴置换检验

    anova.cca(spe.rda, step=1000)

    Permutation test for rda under reduced model

    Permutation: free

    Number of permutations: 999

    Model: rda(formula = spe.hel ~ alt + pen + deb + pH + dur + pho + nit + amm + oxy + dbo, data = env2)

    Df Variance F Pr(>F)

    Model 12 0.36537 3.5522 0.001 ***

    Residual 16 0.13714

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    # 每个典范轴逐一检验

    anova.cca(spe.rda, by="axis", step=1000)

    Permutation test for rda under reduced model

    Forward tests for axes

    Permutation: free

    Number of permutations: 999

    Model: rda(formula = spe.hel ~ alt + pen + deb + pH + dur + pho + nit + amm + oxy + dbo, data = env2)

    Df Variance F Pr(>F)

    RDA1 1 0.228081 26.6098 0.001 ***

    RDA2 1 0.053696 6.2646 0.003 **

    RDA3 1 0.032124 3.7478 0.361

    RDA4 1 0.023207 2.7075 0.763

    RDA5 1 0.008707 1.0159 1.000

    RDA6 1 0.007218 0.8421 1.000

    RDA7 1 0.004862 0.5673 1.000

    RDA8 1 0.002919 0.3406 1.000

    RDA9 1 0.002141 0.2497 1.000

    RDA10 1 0.001160 0.1353 1.000

    RDA11 1 0.000913 0.1066 1.000

    RDA12 1 0.000341 0.0397 1.000

    Residual 16 0.137141

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    #每个典范轴的检验只能输入由公式模式获得的rda结果。有多少个轴结果是

    #显著的?

    # 使用Kaiser-Guttman准则确定残差轴

    spe.rda$CA$eig[spe.rda$CA$eig > mean(spe.rda$CA$eig)]

    PC1 PC2 PC3 PC4 PC5

    0.045802781 0.028143080 0.015288209 0.013987518 0.009841239

    #很明显,还有一部分有意思的变差尚未被目前所用的这套环境变量解释。

    偏RDA分析

    偏RDA:固定地形变量影响后,水体化学属性的效应

    # 简单模式:X和W可以是分离的定量变量表格

    spechem.physio

    spechem.physio

    Call: rda(X = spe.hel, Y = envchem, Z = envtopo)

    Inertia Proportion Rank

    Total 0.5025 1.0000

    Conditional 0.2087 0.4153 3

    Constrained 0.1602 0.3189 7

    Unconstrained 0.1336 0.2659 18

    Inertia is variance

    Eigenvalues for constrained axes:

    RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7

    0.09137 0.04591 0.00928 0.00625 0.00387 0.00214 0.00142

    Eigenvalues for unconstrained axes:

    PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8

    0.04642 0.02072 0.01746 0.01325 0.00974 0.00588 0.00512 0.00400

    (Showed only 8 of all 18 unconstrained eigenvalues)

    summary(spechem.physio)

    # 公式模式;X和W必须在同一数据框内

    class(env)

    spechem.physio2

    + Condition(alt + pen + deb), data=env)

    spechem.physio2

    Call: rda(formula = spe.hel ~ pH + dur + pho + nit + amm + oxy + dbo + Condition(alt + pen + deb), data = env)

    Inertia Proportion Rank

    Total 0.5025 1.0000

    Conditional 0.2087 0.4153 3

    Constrained 0.1602 0.3189 7

    Unconstrained 0.1336 0.2659 18

    Inertia is variance

    Eigenvalues for constrained axes:

    RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7

    0.09137 0.04591 0.00928 0.00625 0.00387 0.00214 0.00142

    Eigenvalues for unconstrained axes:

    PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8

    0.04642 0.02072 0.01746 0.01325 0.00974 0.00588 0.00512 0.00400

    (Showed only 8 of all 18 unconstrained eigenvalues)

    #上面两个分析的结果完全相同。

    偏RDA检验(使用公式模式获得的RDA结果,以便检验每个轴)

    anova.cca(spechem.physio2, step=1000)

    Permutation test for rda under reduced model

    Permutation: free

    Number of permutations: 999

    Model: rda(formula = spe.hel ~ pH + dur + pho + nit + amm + oxy + dbo + Condition(alt + pen + deb), data = env)

    Df Variance F Pr(>F)

    Model 7 0.16024 3.0842 0.001 ***

    Residual 18 0.13360

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    anova.cca(spechem.physio2, step=1000, by="axis")

    Permutation test for rda under reduced model

    Forward tests for axes

    Permutation: free

    Number of permutations: 999

    Model: rda(formula = spe.hel ~ pH + dur + pho + nit + amm + oxy + dbo + Condition(alt + pen + deb), data = env)

    Df Variance F Pr(>F)

    RDA1 1 0.091373 12.3108 0.001 ***

    RDA2 1 0.045907 6.1851 0.010 **

    RDA3 1 0.009277 1.2499 0.964

    RDA4 1 0.006251 0.8422 0.993

    RDA5 1 0.003866 0.5209 0.996

    RDA6 1 0.002142 0.2886 1.000

    RDA7 1 0.001425 0.1920 0.997

    Residual 18 0.133599

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    # 偏RDA三序图(使用拟合值的样方坐标)

    # 1型标尺

    plot(spechem.physio, scaling=1, display=c("sp", "lc", "cn"),

    main="RDA三序图:spe.hel~chem ︳Tope- 1型标尺- 拟合的样方坐标")

    spe3.sc

    arrows(0, 0, spe3.sc[, 1], spe3.sc[, 2], length=0, lty=1, col="red")

    de1bd2531883

    # 2型标尺

    plot(spechem.physio, display=c("sp", "lc", "cn"),

    main="RDA三序图:spe.hel~chem ︳Tope- 2型标尺- 拟合的样方坐标")

    spe4.sc

    arrows(0, 0, spe4.sc[,1], spe4.sc[,2], length=0, lty=1, col="red")

    de1bd2531883

    解释变量向前选择

    每个变量的共线性程度可以用变量的方差膨胀因子(variance inflation factor,VIF)度量,VIF是衡量一个变量的回归系数的方差由共线性引起的膨胀比例。如果VIF值超过20,表示共线性很严重。实际上,VIF超过10则可能会有共线性问题,需要处理。

    # 两个RDA结果的变量方差膨胀因子(VIF)

    # ********************************************

    # 本章第一个RDA结果:包括所有环境因子变量

    vif.cca(spe.rda)

    alt penmoderate pensteep penvery_steep deb pH dur pho nit

    20.397021 2.085921 2.987679 4.594983 6.684187 1.363575 3.049937 30.614913 18.953092

    amm oxy dbo

    40.114746 6.854703 17.980889

    vif.cca(spechem.physio) # 偏RDA

    alt pen deb pH dur pho nit amm oxy dbo

    16.188416 1.873974 6.711460 1.205235 3.268620 25.363359 16.080319 30.685907 6.904214 17.782727

    # 使用双终止准则(Blanchet等,2008a)前向选择解释变量

    # 1.包含所有解释变量的RDA全模型

    spe.rda.all

    Call: rda(formula = spe.hel ~ alt + pen + deb + pH + dur + pho + nit + amm + oxy + dbo, data = env)

    Inertia Proportion Rank

    Total 0.5025 1.0000

    Constrained 0.3689 0.7341 10

    Unconstrained 0.1336 0.2659 18

    Inertia is variance

    Eigenvalues for constrained axes:

    RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7 RDA8 RDA9 RDA10

    0.22803 0.05442 0.03382 0.03008 0.00749 0.00566 0.00443 0.00281 0.00138 0.00079

    Eigenvalues for unconstrained axes:

    PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8

    0.04642 0.02072 0.01746 0.01325 0.00974 0.00588 0.00512 0.00400

    (Showed only 8 of all 18 unconstrained eigenvalues)

    # 2.全模型校正R2

    (R2a.all

    [1] 0.5864353

    # 3.使用packfors 包内forward.sel()函数选择变量

    # library(packfor) #如果尚未载入packfors包,需要运行这一步

    forward.sel(spe.hel, env, adjR2thresh=R2a.all)

    Testing variable 1

    Testing variable 2

    Testing variable 3

    Testing variable 4

    Procedure stopped (adjR2thresh criteria) adjR2cum = 0.594764 with 4 variables (superior to 0.586435)

    variables order R2 R2Cum AdjR2Cum F pval

    1 alt 1 0.32806549 0.3280655 0.3031790 13.182488 0.001

    2 oxy 9 0.16402853 0.4920940 0.4530243 8.396715 0.001

    3 dbo 10 0.09733024 0.5894243 0.5401552 5.926448 0.001

    4 pen 2 0.06323025 0.6526545 0.5947636 4.368924 0.007

    #注意,正如这个函数的提示信息所示,选择最后一个变量其实违背了

    #adjR2thresh终止准则,所以pen变量最终不应该在被选变量列表中。

    # 使用vegan包内ordistep()函数前向选择解释变量。该函数可以对因子变量进# 行选择,也可以运行解释变量的逐步选择和后向选择。

    step.forward

    direction="forward", pstep = 1000)

    Start: spe.hel ~ 1

    Df AIC F Pr(>F)

    + alt 1 -28.504 13.1825 0.005 **

    + oxy 1 -27.420 11.7086 0.005 **

    + deb 1 -26.872 10.9840 0.005 **

    + nit 1 -26.716 10.7795 0.005 **

    + dbo 1 -23.172 6.4340 0.005 **

    + dur 1 -22.499 5.6673 0.005 **

    + pho 1 -22.189 5.3200 0.005 **

    + amm 1 -22.047 5.1620 0.005 **

    + pen 1 -20.155 3.1305 0.005 **

    + pH 1 -17.489 0.4839 0.815

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Step: spe.hel ~ alt

    Df AIC F Pr(>F)

    + oxy 1 -34.620 8.3967 0.005 **

    + dbo 1 -32.103 5.5373 0.005 **

    + amm 1 -30.777 4.1281 0.010 **

    + pho 1 -30.560 3.9032 0.010 **

    + nit 1 -29.451 2.7810 0.035 *

    + pen 1 -29.049 2.3847 0.040 *

    + deb 1 -27.972 1.3504 0.170

    + dur 1 -27.954 1.3332 0.290

    + pH 1 -27.426 0.8403 0.515

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Step: spe.hel ~ alt + oxy

    Df AIC F Pr(>F)

    + dbo 1 -38.789 5.9264 0.005 **

    + pho 1 -37.052 4.1280 0.010 **

    + amm 1 -36.527 3.6055 0.015 *

    + pen 1 -36.399 3.4797 0.015 *

    + dur 1 -34.740 1.8964 0.095 .

    + deb 1 -34.388 1.5714 0.125

    + nit 1 -33.474 0.7474 0.620

    + pH 1 -33.035 0.3605 0.950

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Step: spe.hel ~ alt + oxy + dbo

    Df AIC F Pr(>F)

    + pen 1 -41.639 4.3689 0.015 *

    + dur 1 -39.394 2.2555 0.025 *

    + deb 1 -38.436 1.4019 0.190

    + pho 1 -37.789 0.8420 0.455

    + nit 1 -37.577 0.6611 0.655

    + amm 1 -37.583 0.6656 0.700

    + pH 1 -37.316 0.4399 0.910

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Step: spe.hel ~ alt + oxy + dbo + pen

    Df AIC F Pr(>F)

    + dur 1 -41.502 1.5255 0.130

    + deb 1 -41.058 1.1528 0.380

    + pho 1 -40.822 0.9570 0.470

    + amm 1 -40.587 0.7641 0.630

    + nit 1 -40.560 0.7418 0.635

    + pH 1 -40.375 0.5912 0.760

    > RsquareAdj(rda(spe.hel ~ alt, data=env))$adj.r.squared

    [1] 0.303179

    > RsquareAdj(rda(spe.hel ~ alt+oxy, data=env))$adj.r.squared

    [1] 0.4530243

    > RsquareAdj(rda(spe.hel ~ alt+oxy+dbo, data=env))$adj.r.squared

    [1] 0.5401552

    > RsquareAdj(rda(spe.hel ~ alt+oxy+dbo+pen, data=env))$adj.r.squared

    [1] 0.5947636

    > #简约的RDA分析

    > # **************

    > spe.rda.pars

    > spe.rda.pars

    Call: rda(formula = spe.hel ~ alt + oxy + dbo, data = env)

    Inertia Proportion Rank

    Total 0.5025 1.0000

    Constrained 0.2962 0.5894 3

    Unconstrained 0.2063 0.4106 25

    Inertia is variance

    Eigenvalues for constrained axes:

    RDA1 RDA2 RDA3

    0.21802 0.05088 0.02729

    Eigenvalues for unconstrained axes:

    PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8

    0.05835 0.04357 0.02568 0.01740 0.01451 0.01230 0.00787 0.00646

    (Showed only 8 of all 25 unconstrained eigenvalues)

    > anova.cca(spe.rda.pars, step=1000)

    Permutation test for rda under reduced model

    Permutation: free

    Number of permutations: 999

    Model: rda(formula = spe.hel ~ alt + oxy + dbo, data = env)

    Df Variance F Pr(>F)

    Model 3 0.29619 11.963 0.001 ***

    Residual 25 0.20632

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    > anova.cca(spe.rda.pars, step=1000, by="axis")

    Permutation test for rda under reduced model

    Forward tests for axes

    Permutation: free

    Number of permutations: 999

    Model: rda(formula = spe.hel ~ alt + oxy + dbo, data = env)

    Df Variance F Pr(>F)

    RDA1 1 0.218022 26.4181 0.001 ***

    RDA2 1 0.050879 6.1651 0.001 ***

    RDA3 1 0.027291 3.3069 0.002 **

    Residual 25 0.206319

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    > vif.cca(spe.rda.pars)

    alt oxy dbo

    1.223386 3.544201 3.402515

    > (R2a.pars

    [1] 0.5401552

    # 1型标尺

    plot(spe.rda.pars, scaling=1, display=c("sp", "lc", "cn"),

    main="RDA三序图:spe.hel~alt+oxy+dbo - 1型标尺 - 拟合的样方坐标")

    spe4.sc = scores(spe.rda.pars, choices=1:2, scaling=1, display="sp")

    arrows(0, 0, spe4.sc[, 1], spe4.sc[, 2], length=0, lty=1, col="red")

    de1bd2531883

    # 2型标尺

    plot(spe.rda.pars, display=c("sp", "lc", "cn"),

    main="RDA三序图:spe.hel~alt+oxy+dbo - 2型标尺 - 拟合的样方坐标")

    spe5.sc = scores(spe.rda.pars, choices=1:2, display="sp")

    arrows(0, 0, spe5.sc[,1], spe5.sc[,2], length=0, lty=1, col="red")

    #如果第三典范轴也显著,可以选择绘制轴1和轴3、轴2和轴3的三序图。

    de1bd2531883

    变差分解

    # 变差分解说明图

    par(mfrow=c(1,3))

    showvarparts(2) # 两组解释变量

    showvarparts(3) #三组解释变量

    showvarparts(4) #四组解释变量

    # 1.带所有环境变量的变差分解

    spe.part.all

    spe.part.all

    Partition of variance in RDA

    Call: varpart(Y = spe.hel, X = envchem, envtopo)

    Explanatory tables:

    X1: envchem

    X2: envtopo

    No. of explanatory tables: 2

    Total variation (SS): 14.07

    Variance: 0.50251

    No. of observations: 29

    Partition table:

    Df R.squared Adj.R.squared Testable

    [a+b] = X1 7 0.60579 0.47439 TRUE

    [b+c] = X2 3 0.41526 0.34509 TRUE

    [a+b+c] = X1+X2 10 0.73414 0.58644 TRUE

    Individual fractions

    [a] = X1|X2 7 0.24135 TRUE

    [b] 0 0.23304 FALSE

    [c] = X2|X1 3 0.11205 TRUE

    [d] = Residuals 0.41356 FALSE

    ---

    Use function ‘rda’ to test significance of fractions of interest

    plot(spe.part.all, digits=2)

    #这些图内校正R2是正确的数字,但是韦恩图圆圈大小相同,未与R2的大小成比例。

    > # 2.分别对两组环境变量进行前向选择

    > spe.chem

    > R2a.all.chem

    > forward.sel(spe.hel, envchem, adjR2thresh=R2a.all.chem, nperm=9999)

    Testing variable 1

    Testing variable 2

    Testing variable 3

    Testing variable 4

    Procedure stopped (adjR2thresh criteria) adjR2cum = 0.487961 with 4 variables (superior to 0.474388)

    variables order R2 R2Cum AdjR2Cum F pval

    1 oxy 6 0.30247973 0.3024797 0.2766456 11.708553 0.0001

    2 dbo 7 0.09015052 0.3926303 0.3459095 3.859122 0.0043

    3 nit 4 0.11522115 0.5078514 0.4487936 5.852965 0.0005

    4 amm 5 0.05325801 0.5611094 0.4879610 2.912325 0.0083

    > spe.topo

    > R2a.all.topo

    > forward.sel(spe.hel, envtopo, adjR2thresh=R2a.all.topo, nperm=9999)

    Testing variable 1

    Testing variable 2

    Testing variable 3

    Procedure stopped (alpha criteria): pvalue for variable 3 is 0.228900 (superior to 0.050000)

    variables order R2 R2Cum AdjR2Cum F pval

    1 alt 1 0.32806549 0.3280655 0.3031790 13.182488 0.0001

    2 pen 2 0.05645105 0.3845165 0.3371717 2.384674 0.0469

    > # 解释变量简约组合(基于变量选择的结果)

    > names(env)

    [1] "alt" "pen" "deb" "pH" "dur" "pho" "nit" "amm" "oxy" "dbo"

    > envchem.pars

    > envtopo.pars

    > # 变差分解

    > (spe.part

    Partition of variance in RDA

    Call: varpart(Y = spe.hel, X = envchem.pars, envtopo.pars)

    Explanatory tables:

    X1: envchem.pars

    X2: envtopo.pars

    No. of explanatory tables: 2

    Total variation (SS): 14.07

    Variance: 0.50251

    No. of observations: 29

    Partition table:

    Df R.squared Adj.R.squared Testable

    [a+b] = X1 3 0.50785 0.44879 TRUE

    [b+c] = X2 2 0.38452 0.33717 TRUE

    [a+b+c] = X1+X2 5 0.66351 0.59036 TRUE

    Individual fractions

    [a] = X1|X2 3 0.25318 TRUE

    [b] 0 0.19561 FALSE

    [c] = X2|X1 2 0.14156 TRUE

    [d] = Residuals 0.40964 FALSE

    ---

    Use function ‘rda’ to test significance of fractions of interest

    > plot(spe.part, digits=2)

    > # 所有可测部分的置换检验

    > # [a+b]部分的检验

    > anova.cca(rda(spe.hel, envchem.pars), step=1000)

    Permutation test for rda under reduced model

    Permutation: free

    Number of permutations: 999

    Model: rda(X = spe.hel, Y = envchem.pars)

    Df Variance F Pr(>F)

    Model 3 0.25520 8.5992 0.001 ***

    Residual 25 0.24731

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    > # [b+c]部分的检验

    > anova.cca(rda(spe.hel, envtopo.pars), step=1000)

    Permutation test for rda under reduced model

    Permutation: free

    Number of permutations: 999

    Model: rda(X = spe.hel, Y = envtopo.pars)

    Df Variance F Pr(>F)

    Model 2 0.19322 8.1216 0.001 ***

    Residual 26 0.30929

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    > # [a+b+c]部分的检验

    > env.pars

    > anova.cca(rda(spe.hel, env.pars), step=1000)

    Permutation test for rda under reduced model

    Permutation: free

    Number of permutations: 999

    Model: rda(X = spe.hel, Y = env.pars)

    Df Variance F Pr(>F)

    Model 5 0.33342 9.0704 0.001 ***

    Residual 23 0.16909

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    > # [a]部分的检验

    > anova.cca(rda(spe.hel, envchem.pars, envtopo.pars), step=1000)

    Permutation test for rda under reduced model

    Permutation: free

    Number of permutations: 999

    Model: rda(X = spe.hel, Y = envchem.pars, Z = envtopo.pars)

    Df Variance F Pr(>F)

    Model 3 0.14020 6.3565 0.001 ***

    Residual 23 0.16909

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    > # [c]部分的检验

    > anova.cca(rda(spe.hel, envtopo.pars, envchem.pars), step=1000)

    Permutation test for rda under reduced model

    Permutation: free

    Number of permutations: 999

    Model: rda(X = spe.hel, Y = envtopo.pars, Z = envchem.pars)

    Df Variance F Pr(>F)

    Model 2 0.078219 5.3197 0.001 ***

    Residual 23 0.169091

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    > #各个部分置换检验有不显著的吗?

    > # 3.无变量nit(硝酸盐浓度)的变差分解

    > envchem.pars2

    > (spe.part2

    Partition of variance in RDA

    Call: varpart(Y = spe.hel, X = envchem.pars2, envtopo.pars)

    Explanatory tables:

    X1: envchem.pars2

    X2: envtopo.pars

    No. of explanatory tables: 2

    Total variation (SS): 14.07

    Variance: 0.50251

    No. of observations: 29

    Partition table:

    Df R.squared Adj.R.squared Testable

    [a+b] = X1 2 0.39263 0.34591 TRUE

    [b+c] = X2 2 0.38452 0.33717 TRUE

    [a+b+c] = X1+X2 4 0.65265 0.59476 TRUE

    Individual fractions

    [a] = X1|X2 2 0.25759 TRUE

    [b] 0 0.08832 FALSE

    [c] = X2|X1 2 0.24885 TRUE

    [d] = Residuals 0.40524 FALSE

    ---

    Use function ‘rda’ to test significance of fractions of interest

    > plot(spe.part2, digits=2)

    de1bd2531883

    RDA 作为多元方差分析(MANOVA)的工具

    # 基于RDA的双因素MANOVA

    # **************************

    # 生成代表海拔的因子变量(3个水平,每个水平含9个样方)

    alt.fac

    # 生成近似模拟pH值的因子变量

    pH.fac

    2, 1, 2, 3, 2, 1, 1, 3, 3))

    # 两个因子是否平衡?

    table(alt.fac, pH.fac)

    table(alt.fac, pH.fac)

    # 用Helmert对照法编码因子和它们的交互作用项

    alt.pH.helm

    contrasts=list(alt.fac="contr.helmert", pH.fac="contr.helmert"))

    alt.pH.helm

    (Intercept) alt.fac1 alt.fac2 pH.fac1 pH.fac2 alt.fac1:pH.fac1 alt.fac2:pH.fac1 alt.fac1:pH.fac2 alt.fac2:pH.fac2

    1 1 -1 -1 -1 -1 1 1 1 1

    2 1 -1 -1 1 -1 -1 -1 1 1

    3 1 -1 -1 0 2 0 0 -2 -2

    4 1 -1 -1 1 -1 -1 -1 1 1

    5 1 -1 -1 0 2 0 0 -2 -2

    6 1 -1 -1 -1 -1 1 1 1 1

    7 1 -1 -1 0 2 0 0 -2 -2

    8 1 -1 -1 1 -1 -1 -1 1 1

    9 1 -1 -1 -1 -1 1 1 1 1

    10 1 1 -1 1 -1 1 -1 -1 1

    11 1 1 -1 -1 -1 -1 1 -1 1

    12 1 1 -1 0 2 0 0 2 -2

    13 1 1 -1 0 2 0 0 2 -2

    14 1 1 -1 1 -1 1 -1 -1 1

    15 1 1 -1 -1 -1 -1 1 -1 1

    16 1 1 -1 -1 -1 -1 1 -1 1

    17 1 1 -1 1 -1 1 -1 -1 1

    18 1 1 -1 0 2 0 0 2 -2

    19 1 0 2 1 -1 0 2 0 -2

    20 1 0 2 -1 -1 0 -2 0 -2

    21 1 0 2 1 -1 0 2 0 -2

    22 1 0 2 0 2 0 0 0 4

    23 1 0 2 1 -1 0 2 0 -2

    24 1 0 2 -1 -1 0 -2 0 -2

    25 1 0 2 -1 -1 0 -2 0 -2

    26 1 0 2 0 2 0 0 0 4

    27 1 0 2 0 2 0 0 0 4

    attr(,"assign")

    [1] 0 1 1 2 2 3 3 3 3

    attr(,"contrasts")

    attr(,"contrasts")$alt.fac

    [1] "contr.helmert"

    attr(,"contrasts")$pH.fac

    [1] "contr.helmert"

    #检查当前对照法产生的表格,哪一列代表海拔因子、pH值和交互作用项?

    # 检查Helmert对照表属性1:每个变量的和为1

    apply(alt.pH.helm[, 2:9], 2, sum)

    # 检查Helmert对照表属性2:变量之间不相关

    cor(alt.pH.helm[, 2:9])

    alt.fac1 alt.fac2 pH.fac1 pH.fac2 alt.fac1:pH.fac1 alt.fac2:pH.fac1 alt.fac1:pH.fac2 alt.fac2:pH.fac2

    alt.fac1 1 0 0 0 0 0 0 0

    alt.fac2 0 1 0 0 0 0 0 0

    pH.fac1 0 0 1 0 0 0 0 0

    pH.fac2 0 0 0 1 0 0 0 0

    alt.fac1:pH.fac1 0 0 0 0 1 0 0 0

    alt.fac2:pH.fac1 0 0 0 0 0 1 0 0

    alt.fac1:pH.fac2 0 0 0 0 0 0 1 0

    alt.fac2:pH.fac2 0 0 0 0 0 0 0 1

    # 使用函数betadisper()(vegan包)(Marti Anderson检验)验证组内协方差矩阵# 的齐性

    spe.hel.d1

    # 海拔因子

    (spe.hel.alt.MHV

    Homogeneity of multivariate dispersions

    Call: betadisper(d = spe.hel.d1, group = alt.fac)

    No. of Positive Eigenvalues: 26

    No. of Negative Eigenvalues: 0

    Average distance to median:

    1 2 3

    0.5208 0.5175 0.3881

    Eigenvalues for PCoA axes:

    PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8

    6.5329 1.7407 1.2269 1.0591 0.6117 0.4683 0.3987 0.3207

    anova(spe.hel.alt.MHV)

    Analysis of Variance Table

    Response: Distances

    Df Sum Sq Mean Sq F value Pr(>F)

    Groups 2 0.1032 0.051602 0.8763 0.4292

    Residuals 24 1.4133 0.058889

    permutest(spe.hel.alt.MHV) # 置换检验

    Permutation test for homogeneity of multivariate dispersions

    Permutation: free

    Number of permutations: 999

    Response: Distances

    Df Sum Sq Mean Sq F N.Perm Pr(>F)

    Groups 2 0.1032 0.051602 0.8763 999 0.439

    Residuals 24 1.4133 0.058889

    > # pH值因子

    > (spe.hel.pH.MHV

    Homogeneity of multivariate dispersions

    Call: betadisper(d = spe.hel.d1, group = pH.fac)

    No. of Positive Eigenvalues: 26

    No. of Negative Eigenvalues: 0

    Average distance to median:

    1 2 3

    0.6658 0.6716 0.7019

    Eigenvalues for PCoA axes:

    PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8

    6.5329 1.7407 1.2269 1.0591 0.6117 0.4683 0.3987 0.3207

    > anova(spe.hel.pH.MHV)

    Analysis of Variance Table

    Response: Distances

    Df Sum Sq Mean Sq F value Pr(>F)

    Groups 2 0.00676 0.0033802 0.1587 0.8542

    Residuals 24 0.51124 0.0213018

    > permutest(spe.hel.pH.MHV) #置换检验

    Permutation test for homogeneity of multivariate dispersions

    Permutation: free

    Number of permutations: 999

    Response: Distances

    Df Sum Sq Mean Sq F N.Perm Pr(>F)

    Groups 2 0.00676 0.0033802 0.1587 999 0.855

    Residuals 24 0.51124 0.0213018

    > #组内协方差齐性,可以继续分析。

    # 首先检验交互作用项。海拔因子和pH因子构成协变量矩阵

    interaction.rda

    anova(interaction.rda, step=1000, perm.max=1000)

    #交互作用是否显著?显著的交互作用表示一个因子的影响依赖另一个因子

    #的水平,这将妨害主因子变量的分析。

    # 检验海拔因子的效应,此时pH值因子和交互作用项作为协变量矩阵

    factor.alt.rda

    anova(factor.alt.rda, step=1000, perm.max=1000, strata=pH.fac)

    #海拔因子影响是否显著?

    #检验pH值因子的效应,此时海拔值因子和交互作用项作为协变量矩阵

    factor.pH.rda

    alt.pH.helm[, c(2:3, 6:9)])

    anova(factor.pH.rda, step=1000, perm.max=1000, strata=alt.fac)

    #pH值影响是否显著?

    # RDA和显著影响的海拔因子三序图

    alt.rda.out

    plot(alt.rda.out, scaling=1, display=c("sp", "wa", "cn"),

    main="Multivariate ANOVA, factor altitude - scaling 1 - wa scores")

    spe.manova.sc

    arrows(0, 0, spe.manova.sc[, 1], spe.manova.sc[, 2], length=0, col="red")

    de1bd2531883

    基于距离的RDA分析

    > # 基于距离的RDA分析(db-RDA)

    > # ****************************

    > # 1.分步计算

    > spe.bray

    > spe.pcoa

    > spe.scores

    > # 交互作用的检验。从协变量矩阵获得Helmert对照编码海拔因子和pH值因子

    > interact.dbrda

    > anova(interact.dbrda, step=1000, perm.max=1000)

    Permutation test for rda under reduced model

    Permutation: free

    Number of permutations: 999

    Model: rda(X = spe.scores[1:27, ], Y = alt.pH.helm[, 6:9], Z = alt.pH.helm[, 2:5])

    Df Variance F Pr(>F)

    Model 4 0.021857 0.441 1

    Residual 18 0.223016

    > #交互作用是否显著?如果没有,可以继续检验主因子的效应(此处未显示)

    > # 2.直接使用vegan包内capscale()函数运行。只能以模型界面运行。响应变量

    > #可以是原始数据矩阵。

    > interact.dbrda2

    > anova(interact.dbrda2, step=1000, perm.max=1000)

    Permutation test for capscale under reduced model

    Permutation: free

    Number of permutations: 999

    Model: capscale(formula = spe[1:27, ] ~ alt.fac * pH.fac + Condition(alt.fac + pH.fac), distance = "bray", add = TRUE)

    Df SumOfSqs F Pr(>F)

    Model 4 0.4667 0.4811 1

    Residual 18 4.3658

    > # 或者响应变量可以是相异矩阵。

    > interact.dbrda3

    > anova(interact.dbrda3, step=1000, perm.max=1000)

    Permutation test for capscale under reduced model

    Permutation: free

    Number of permutations: 999

    Model: capscale(formula = spe.bray ~ alt.fac * pH.fac + Condition(alt.fac + pH.fac), add = TRUE)

    Df SumOfSqs F Pr(>F)

    Model 4 0.4667 0.4811 1

    Residual 18 4.3658

    非线性关系的RDA分析

    > # 二阶解释变量的RDA

    > # *******************

    > # 生成das和das正交二阶项(由poly()函数获得)矩阵

    > das.df

    > colnames(das.df)

    > # 验证两个变量是否显著

    > forward.sel(spe, das.df)

    Testing variable 1

    Testing variable 2

    variables order R2 R2Cum AdjR2Cum F pval

    1 das 1 0.44777219 0.4477722 0.4273193 21.892865 0.001

    2 das2 2 0.07870749 0.5264797 0.4900550 4.321662 0.005

    > # RDA和置换检验

    > spe.das.rda

    > anova(spe.das.rda)

    Permutation test for rda under reduced model

    Permutation: free

    Number of permutations: 999

    Model: rda(formula = spe ~ das + das2, data = as.data.frame(das.df))

    Df Variance F Pr(>F)

    Model 2 36.304 14.454 0.001 ***

    Residual 26 32.652

    ---

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    > # 三序图(拟合的样方坐标,2型标尺)

    > plot(spe.das.rda, scaling=2, display=c("sp", "lc", "cn"),

    + main="RDA三序图:spe~das + das2 - 2型标尺 - 拟合的样方坐标")

    > spe6.sc

    > arrows(0, 0, spe6.sc[, 1], spe6.sc[, 2], length=0, lty=1, col="red")

    de1bd2531883

    # 4种鱼类的分布地图

    # ******************

    par(mfrow=c(2, 2))

    plot(spa$x, spa$y, asp=1, col="brown", cex=spe$TRU,

    xlab="x (km)", ylab="y (km)", main="褐鳟")

    lines(spa$x, spa$y, col="light blue")

    plot(spa$x, spa$y, asp=1, col="brown", cex=spe$OMB,

    xlab="x (km)", ylab="y (km)", main="鳟鱼")

    lines(spa$x, spa$y, col="light blue")

    plot(spa$x, spa$y, asp=1, col="brown", cex=spe$ABL,

    xlab="x (km)", ylab="y (km)", main="欧鮊鱼")

    lines(spa$x, spa$y, col="light blue")

    plot(spa$x, spa$y, asp=1, col="brown", cex=spe$TAN,

    xlab="x (km)", ylab="y (km)", main="鲤鱼")

    lines(spa$x, spa$y, col="light blue")

    de1bd2531883

    自写代码角

    为了能够正确自写RDA分析代码,有必要参考Legendre和Legendre

    (1998)第11.1节相关内容。下面是计算步骤(基于协方差矩阵的RDA)

    1.计算中心化的物种数据矩阵与标准化解释变量矩阵的多元线性回归,获得拟合值矩阵;

    2.计算拟合值矩阵的PCA;

    3.计算两类样方坐标;

    4.结果输出。

    下面代码解释部分使用的公式编码与Legendre和Legendre(1998)一致。

    下面的代码集中在RDA约束部分,目的是使读者对RDA数学过程感兴趣,而不是最优化程序。

    myRDA

    # 1.数据的准备

    # *************

    Y.mat

    Yc

    X.mat

    Xcr

    # 2.多元线性回归的计算

    # *********************

    # 回归系数矩阵 (eq. 11.4)

    B

    # 拟合值矩阵(eq. 11.5)

    Yhat

    # 残差矩阵

    Yres

    #维度

    n

    p

    m

    # 3. 拟合值PCA分析

    # ******************

    # 协方差矩阵 (eq. 11.7)

    S

    # 特征根分解

    eigenS

    # 多少个典范轴?

    kc 0.00000001))

    # 典范轴特征根

    ev

    # 矩阵Yc(中心化)的总方差(惯量)

    trace = sum(diag(cov(Yc)))

    # 正交特征向量(响应变量的贡献,1型标尺)

    U

    row.names(U)

    # 样方坐标(vegan包内'wa' 坐标,1型标尺eq. 11.12)

    F

    row.names(F)

    # 样方约束(vegan包内'lc' 坐标,2型标尺eq. 11.13)

    Z

    row.names(Z)

    # 典范系数 (eq. 11.14)

    CC

    row.names(CC)

    # 解释变量

    # 物种-环境相关

    corXZ

    # 权重矩阵的诊断

    D

    # 解释变量双序图坐标

    coordX

    coordX2

    row.names(coordX)

    row.names(coordX2)

    # 相对特征根平方根转化(为2型标尺)

    U2

    row.names(U2)

    F2

    row.names(F2)

    Z2

    row.names(Z2)

    # 未校正R2

    R2

    # 校正R2

    R2adj

    # 4.残差的PCA

    # *******************

    # 与第5章相同,写自己的代码,可以从这里开始...

    # eigenSres

    # evr

    # 5.输出Output

    result

    names(result)

    "Species_sc1", "wa_sc1", "lc_sc1", "Biplot_sc1", "Species_sc2",

    "wa_sc2", "lc_sc2", "Biplot_sc2")

    result

    }

    #将此函数应用到Doubs鱼类数据和环境数据的RDA分析

    doubs.myRDA

    summary(doubs.myRDA)

    Length Class Mode

    Total_variance 1 -none- numeric

    R2 1 -none- numeric

    R2a 1 -none- numeric

    Can_ev 10 -none- numeric

    Can_coeff 100 -none- numeric

    Species_sc1 270 -none- numeric

    wa_sc1 290 -none- numeric

    lc_sc1 290 -none- numeric

    Biplot_sc1 100 -none- numeric

    Species_sc2 270 -none- numeric

    wa_sc2 290 -none- numeric

    lc_sc2 290 -none- numeric

    Biplot_sc2 100 -none- numeric

    展开全文
  • ​Q1:什么是RDA分析RDA分析(Redundancy analysis),即冗余分析,对比主成分分析可以发现,其实冗余分析就是约束化的主成分分析。冗余分析(redundancy analysis, RDA)或者典范对应分析(canonical ...
  • 原来你是这样的排序分析杨慧宏基因组2017-12-16微生态相关文献中都经常出现。这些分析成图相似,且都是通过样本点之间的距离反映样本间菌群结构的相似性和差异...PCA、PCoA和NMDS分析属于非约束性排序分析,而RDA/CC...
  • R语言RDA分析结果环境因子丢失

    千次阅读 2020-12-30 23:32:48
    用R语言做RDA分析,分析环境因子与物种的相关性。有27的样本,31个物种,14个环境因子。但是最后输出的环境因子得分表“Biplot scores for constraining variables”只始终只显示13个环境因子是什么情况?!如下图:...
  • env=log10(env) rda_analysis<-rda(phylum1~.,env) result(rda_analysis) result sp=as.data.frame(result$species[,1:2])*5###提取相应变量坐标,乘以5是使图美观,不影响分析 st=as.data.frame(result$sites[,1:2]...
  • 使用CANOCO进行CCA或RDA教程分析.ppt

    千次阅读 2021-01-17 18:51:15
    如何用CANOCO 4.5 进行CCA/RDA分析 郑有用 由于制作者水平有限,本篇不涉及蒙特卡洛检验。进一步学习应参考《Canoco for Windows 4.5 中文简明教程》 与《Multivariate Analysis of Ecological Data Using Canoco》...
  • R语言实现冗余分析RDA)完整代码

    万次阅读 多人点赞 2020-05-21 21:07:16
    #读入物种数据,细菌门水平丰度表(OTU 水平数据量太大,后续的置换检验和变量选择过程很费时间,不方便示例演示) phylum <- read.delim('phylum_table.txt', row.names = 1, sep = '\t', stringsAsFactors = ...
  • 在进行微生物多样性分析时,大家一定会想微生物群落与环境因子之间的关联性分析,其中RDA分析便能帮助我们解析微生物群落与环境因子之间的相关性问题。RDA分析即冗余分析,是环境因子约束化的PCA分析,可以将样本...
  • oracle诊断工具-RDA使用

    2020-12-23 16:09:42
    收藏[@more@]RDA是Remote Diagnostic Agent 的简称,是oracle用来收集、分析数据库的工具,运行该工具不会改变系统的任何参数,RDA收集的相关数据非常全面,可以简化我们日常监控、分析数据库的工作,Oracle Support...
  • 与PCA、PCoA、CA及NMDS等非约束排序方法不同的是:RDA分析涉及到两个矩阵间的降维排序分析,受RDA排序算法创始人的影响,RDA分析一般需要提供名为“群落物种”和“环境变量”矩阵数据。RDA的优势是结合了主成分分析...
  • 冗余分析(Redundancy analysis,RDA)是一种回归分析结合主成分分析的排序方法,也是多响应变量(multi-response)回归分析的拓展。在群落分析中常使用RDA,将物种多度的变化分解为与环境变量相关的变差(variation...
  • Remote Diagnostic Agent (RDA) 是一个工程师用Perl语言编写的命令行诊断工具,RDA提供统一的诊断工具支持包和预防的解决方法。...RDA不会对系统任何的修改,它只为Oracle支持收集有用的数据,如果需要可...
  • 在线画图工具-CCA与RDA分析

    万次阅读 2017-04-27 09:44:21
    RDA分析(Redundancy analysis),即冗余分析,对比主成分分析可以发现,其实冗余分析就是约束化的主成分分析。 RDA或CCA的选择问题:RDA是基于线性模型,CCA是基于单峰模型。一般我们会选择CCA来直接梯度分析。但是...
  • 我想用RDA作非生态类的分析,只有一组数据,这种情况相当于生态类的样方只有一种,可以作RDA分析
  • # wascores()函数以多度加权平均方式将物种被动投影到样方的PCoA图中,获得物种坐标信息 spe.wa (pcoa$points, spe) spe.wa # 提取特征根 eig = pcoa$eig eig # 绘制PCoA排序图,其实排序分析做完之后,都是使用...
  • R语言绘图实战:RDA冗余分析

    万次阅读 多人点赞 2017-05-13 10:36:24
    #载入vegan包 library(vegan) #读取“样本-物种”文件 sp (file=file.choose(),sep="\t",header=T,row.names=1) sp #读取“样本-环境因子”文件 se (file=file.choose(),sep="\t",header=T,row.names=1) ...#选择用RDA
  • RDA进行微生物环境因子分析

    万次阅读 2018-11-02 09:48:36
    在进行微生物多样性分析时,大家一定会α,Β多样性分析。α多样品通俗来讲就是样本内的物种多样性。Β多样性是指在地区尺度上,物种组成沿着某个梯度方向从一个群落到另一个群落的变化率。即沿着某一环境梯度,...
  • RDA分析(Redundancy analysis),即冗余分析,对比主成分分析可以发现,其实冗余分析就是约束化的主成分分析。 RDA或CCA的选择问题:RDA是基于线性模型,CCA是基于单峰模型。一般我们会选择CCA来直接梯度分析。但是...
  • 使用SPSS软件进行数据分析一个严谨的实验,在获得实验数据的那一刻,数据的处理与分析,就显得尤为重要,今天先聊聊主成分分析,以城市分析为例子,大家延伸思考吧,文档已通过论证。想跟师兄交流实验数据问题,加我...
  • 前言在进行微生物多样性分析时,大家一定会α,β多样性分析。通俗来讲,α多样性就是样本内的物种多样性。β多样性是指在地区尺度上,物种组成沿着某个梯度方向从一个群落到另一个群落的变化率。即...
  • 前言在进行微生物多样性分析时,大家一定会α,β多样性分析。通俗来讲,α多样性就是样本内的物种多样性。β多样性是指在地区尺度上,物种组成沿着某个梯度方向从一个群落到另一个群落的变化率。即...
  • 原标题:欧易/鹿明生物云平台:点点鼠标,轻松完成RDA/CCA分析编者按常规做RDA、CCA分析一般使用CANOCO软件或者R语言的Vegan包去做分析,前者太贵,后者复杂,肿么办?欧易/鹿明云平台小工具别担心,云平台小工具来...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,230
精华内容 492
关键字:

RDA分析怎么做