精华内容
下载资源
问答
  • 可以应用于数据的多元分析,降维以及冗余分析,文件内包含多个程序代码,数据文件请根据自身自行导入,这里不提供相应数据文件,该代码偏向初学者,若有不足,欢迎各位指正
  • RDA冗余分析及方差分解,对于解决环境主要影响因子作用明显。
  • R、冗余分析(RDA)、ggplot2、置信椭圆

    万次阅读 多人点赞 2018-10-12 11:28:47
    R、冗余分析(RDA)、ggplot2  在生态环境领域中,冗余分析(RDA)是我们常用的分析方法,分析目的为“解释变量”对“响应变量”的影响情况。类似RDA的方法,还有CCA。这里以RDA为例→数据处理、分析过后,我们...

    R、冗余分析(RDA)、ggplot2、置信椭圆

     在生态环境领域中(实际中,其他专业也用到),冗余分析(RDA)是我们常用的分析方法,分析目的为“解释变量”对“响应变量”的影响情况。类似RDA的方法,还有CCA。这里以RDA为例→数据处理、分析过后,我们需要对结果进行可视化,R语言ggplt2程序包无疑是可视化神器,然而,怎样利用ggplot2对RDA结果进行可视化,需要我们对RDA结果进行了解,提取需要展示的元素。许多论文中在排序分析的图中添加了置信椭圆,ggplot2自带函数目前支持每个分组样本量≥4时画椭圆,这里我们使用gglayer包的geom_ord_ellipse函数给样方添加置信椭圆(支持样本量≥3)。

     library(vegan)
     library(plyr)
      library(gglayer)
     library(ggplot2)
     library(ggrepel)
    fc=read.csv("D:\\R\\factor.csv",header = T,row.names = 1)#读取解释变量数据
    sp=read.csv("D:\\R\\sp.csv",header = T,row.names = 1)#读取响应变量数据
    spp=decostand(sp,method = "hellinger")#对响应变量做转化
    fcc=log10(fc)#对解释变量做转化
    uu=rda(spp~.,fcc)#RDA分析
    ii=summary(uu)  #查看分析结果
    sp=as.data.frame(ii$species[,1:2])*2#可根据出图结果,对画图数据做一定的放大或缩小,下同
    st=as.data.frame(ii$sites[,1:2])
    yz=as.data.frame(ii$biplot[,1:2])
    grp=as.data.frame(c(rep("a",3),rep("b",3),rep("c",4),rep("d",6)))#根据样方类型分组,“a”有3个样本,“b”有3个样本……,共16个。注意样本的顺序和个数!
    colnames(grp)="group"
    ggplot() +
      geom_point(data = st,aes(RDA1,RDA2,shape=grp$group,fill=grp$group),size=4)+
      scale_shape_manual(values = c(21:25))+
      geom_ord_ellipse(aes(st$RDA1,st$RDA2,color=grp$group,group=grp$group),###注意,是在这里添加椭圆
                       ellipse_pro = 0.68,linetype=3,size=1)+###注意,是在这里添加椭圆
      geom_segment(data = sp,aes(x = 0, y = 0, xend = RDA1, yend = RDA2), 
                   arrow = arrow(angle=22.5,length = unit(0.35,"cm"),
                                 type = "closed"),linetype=1, size=0.6,colour = "red")+
      geom_text_repel(data = sp,aes(RDA1,RDA2,label=row.names(sp)))+
      geom_segment(data = yz,aes(x = 0, y = 0, xend = RDA1, yend = RDA2), 
                   arrow = arrow(angle=22.5,length = unit(0.35,"cm"),
                                 type = "closed"),linetype=1, size=0.6,colour = "blue")+
      geom_text_repel(data = yz,aes(RDA1,RDA2,label=row.names(yz)))+
      labs(x=paste("RDA 1 (", format(100 *ii$cont[[1]][2,1], digits=4), "%)", sep=""),
           y=paste("RDA 2 (", format(100 *ii$cont[[1]][2,2], digits=4), "%)", sep=""))+
      geom_hline(yintercept=0,linetype=3,size=1) + 
      geom_vline(xintercept=0,linetype=3,size=1)+
      guides(shape=guide_legend(title=NULL),color=guide_legend(title=NULL),
             fill=guide_legend(title=NULL))+
      theme_bw()+theme(panel.grid=element_blank())

    在这里插入图片描述

      RDA探索完毕,上面的图比较骚气,但,不能用,因为需要检验模型是否显著、解释变量是否存在共线性,在论文写作中应删掉某些图层,使图更加美观、更简洁。这里只是简单的举个例,在进行约束排序分析之前,我们先要检查数据,是符合RDA还是CCA,网上有很多例子,这里不再赘述。当RDA或者CCA不能很好的解决我们的问题时,我们需要结合其它分析方法,如gbm、RF能求出相对重要性等,也可利用SEM去探索。原则上,不同的数据需要用不同的方法探索,选择一种最理想的结果。
     刚接触R的朋友,可能会因为R的“难”入门而选择较容易的CANOCO软件,当然后者也是生态环境领域的数据分析神器,但是,当我们想要个性化分析、个性化出图时,目前后者无法满足,R,应是首选。实验设计、实验过程、数据处理及分析、绘图、写作等是科研的必然过程,为促进相互进步、资源共享,我们创建了学术交流QQ群:335774366。欢迎有兴趣的朋友加入→指导。
    声明:以上内容仅为作者个人理解,有不对的地方,欢迎指正。

    展开全文
  • 微生物生态数据分析——冗余分析

    千次阅读 2018-10-25 17:13:27
    微生物生态数据分析——冗余分析 sa=read.table("NRRDA.csv",header=T,row.names=1,sep=",") env=read.table("NRenv.csv&...

    微生物生态数据分析——冗余分析

    sa=read.table("NRRDA.csv",header=T,row.names=1,sep=",")
    env=read.table("NRenv.csv",header=T,sep=",",row.names=1)
    name <- read.table("name.csv", header=F,sep=",",colClasses=c("character","character"))#用作产生图例,对样本进行分类
    library(vegan)
    decorana(sa)#DCA函数,用于决定选择进行RDA或者CCA
    #如果DCA排序前4个轴中最大值超过4,选择单峰模型排序更合适。如果是小于3,则选择线性模型更好(Lepx & Smilauer 2003)。如果介于3-4之间,单峰模型和线性模型都可行
    sam=rda(sa,env,center=T,scale=T)
    plot(sam,scaling=3)#plot(gts.rda,display=c("sp","bp"),scaling=3)
    # display=c("sp","bp") 表示显示物种与环境因子。如果显示样方与环境因子,可以表示为display=c("si","bp"),如果物种、样方和环境因子同时显示,可以设定display=c("sp","si","bp")。
    library(ggplot2)
    new=sam$CCA
    samples<-data.frame(sample=row.names(new$u),RDA1=new$u[,1],RDA2=new$u[,2]) 
    name=name$V2#分类文件
    species<-data.frame(spece=row.names(new$v),RDA1=new$v[,1],RDA2=new$v[,2])  
    envi<-data.frame(en=row.names(new$biplot),RDA1=new$biplot[,1],RDA2=new$biplot[,2]) 
    rda1 =round(sam$CCA$eig[1]/sum(sam$CCA$eig)*100,2) #第一轴标签,显示解释度
    rda2 =round(sam$CCA$eig[2]/sum(sam$CCA$eig)*100,2) #第二轴标签,显示解释度
    line_x = c(0,envi[1,2],0,envi[2,2],0,envi[3,2])  #行数与ling_g数量一致,envi[line_g,2]
    line_x  
    line_y = c(0,envi[1,3],0,envi[2,3],0,envi[3,3])  
    line_y  
    line_g = c("grazing","grazing","Tot bio","Tot bio","ND","ND")  
    line_g  
    line_data = data.frame(x=line_x,y=line_y,group=line_g)  
    line_data
    #开始重绘RDA图
    #填充样本数据,分别以RDA1,RDA2为X,Y轴,不同样本以颜色区分
    ggplot(data=samples,aes(RDA1,RDA2)) + geom_point(aes(color=sample),size=2)#生成图例的数据使用的是sample,不是samples!,不然会提示错误: Aesthetics must be either length 1 or the same as the data (6): colour, x, y
    #geom_point(aes(color=name),size=2
    #填充微生物物种数据,不同物种以图形区分,seq增加形状数量
     + geom_point(data=species,aes(shape=spece),size=2) + scale_shape_manual(values=seq(0,15))+
    #填充环境因子数据,直接展示
    geom_text(data=envi,aes(label=en),color="blue") + 
    #添加0刻度纵横线
    geom_hline(yintercept=0) + geom_vline(xintercept=0)+ 
    #添加原点指向环境因子的带箭头的直线,size可以调节直线粗细
    geom_line(data=line_data,aes(x=x,y=y,group=group),color="green") + 
    geom_segment(data=line_data,aes(x=x,y=y,xend = line_data[,1], yend = line_data[,2],group=group),color="black",size=1.5,arrow=arrow(angle=35, length=unit(0.3, "cm"))) +
    #添加横纵坐标轴标签
    labs(title="RDA Plot",x=paste("RDA1 ",rda1," %"),y=paste("RDA2 ",rda2," %")) +
    #标题字体格式设置
    #theme(text=element_text(family="serif"))+
    #去除背景颜色及多余网格线
    theme_bw() + theme(panel.grid=element_blank()) 
    #大功告成,保存为矢量图等等
    ggsave("NRRDA2.PDF")      
    #RDA更详细的分析,
    summary(sam)
    #检验环境因子相关显著性(Monte Carlo permutation test)
    permutest(sam,permu=999) # permu=999是表示置换循环的次数
    ef=envfit(sam,env,permu=999)#每个环境因子显著性检验
    ef
    

    看到有很多同学问导入的数据格式问题,我重新写了一份详细的排序分析文章,发在了公众号“科学研究进行时”上,文章名称:R绘图-RDA排序分析,里面有数据格式的截图,大家可以扫描二维码关注查看。
    在这里插入图片描述

    展开全文
  • R语言实现冗余分析(RDA)完整代码

    万次阅读 多人点赞 2020-05-21 21:07:16
    话不多说,先贴代码。 需要数据和文档、数据私信我! library(vegan) ##读取数据 #读入物种数据,细菌门水平丰度表(OTU 水平数据量太大,后续的置换检验和变量选择过程很费时间,不方便做示例演示) ...

    话不多说,先贴代码。

    文末有数据的获取地址

    library(vegan)
    
    ##读取数据
    #读入物种数据,细菌门水平丰度表(OTU 水平数据量太大,后续的置换检验和变量选择过程很费时间,不方便做示例演示)
    phylum <- read.delim('phylum_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
    
    #读取环境数据
    env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
    
    ##RDA
    #调用格式 1
    #rda(Y, X, W)
    #或者 2
    #rda(Y~var1+var2+var3+factorA+var2*var3+Condition(var4))
    
    #直接使用原始数据,不做转化。对于群落物种组成数据来讲(因为通常包含很多 0 值),不是很推荐
    rda_result <- rda(phylum~., env, scale = FALSE)
    
    ##tb-RDA
    #物种数据 Hellinger 预转化(处理包含很多 0 值的群落物种数据时,推荐使用)
    phylum_hel <- decostand(phylum, method = 'hellinger')
    
    #使用全部的环境数据
    rda_tb <- rda(phylum_hel~., env, scale = FALSE)
    
    #若只关注局部环境数据,除了在原始表格中修改变量个数外,还可直接在 rda() 中指定
    #例如只考虑 pH、TC、TN、AP、AK 这 5 种环境变量
    #rda_tb <- rda(phylum_hel~pH+TC+TN+AP+AK, env, scale = FALSE)
    
    #默认情况下,rda(phylum_hel~., env),不包含环境变量间的交互作用
    #若想关注某两个环境变量间的交互作用,例如添加 TC 和 TN 的交互作用项
    #rda_tb <- rda(phylum_hel~pH+TC+DOC+SOM+TN+NO3+NH4+AP+AK+TC*TN, env, scale = FALSE)
    
    ##偏 RDA
    #例如控制土壤 pH 影响后(pH 作为协变量),观测其它环境因素的影响;物种数据 Hellinger 预转化
    rda_part <- rda(phylum_hel~TC+DOC+SOM+TN+NO3+NH4+AP+AK+Condition(pH), data = env, scale = FALSE)
    
    ##db-RDA(以 Bray-curtis 距离为例,注意这里直接使用了原始的物种丰度矩阵)
    #计算距离
    dis_bray <- vegdist(phylum, method = 'bray')
    
    #PCoA 排序
    pcoa <- cmdscale(dis_bray, k = nrow(phylum) - 1, eig = TRUE, add = TRUE)
    
    #提取 PCoA 样方得分(坐标)
    pcoa_site <- pcoa$point
    
    #db-RDA,使用全部的环境数据
    rda_db <- rda(pcoa_site, env, scale = FALSE)
    
    #或者,capscale() 提供了直接运行的方法
    rda_db <- capscale(phylum~., env, distance = 'bray', add = TRUE)
    
    #若基于欧氏距离,则和常规 RDA 的结果一致
    #物种数据 Hellinger 转化后,分别执行 tb-RDA 与使用欧氏距离的 db-RDA
    rda_tb_test <- rda(phylum_hel~., env)
    rda_db_test <- capscale(phylum_hel~., env, distance = 'euclidean')
    
    par(mfrow = c(1, 2))
    plot(rda_tb_test, scaling = 1)
    plot(rda_db_test, scaling = 1)
    
    ##非线性 RDA *
    #本文不做介绍
    #可参见 “DanielBorcard, FranoisGillet, PierreLegendre, et al. 数量生态学:R语言的应用(赖江山 译). 高等教育出版社, 2014.” 170-174 页内容
    
    ##RDA 结果解读,以下以 tb-RDA 结果为例
    #查看统计结果信息,以 I 型标尺为例
    rda_tb.scaling1 <- summary(rda_tb, scaling = 1)
    rda_tb.scaling1
    
    #作图查看排序结果,三序图,包含 I 型标尺和 II 型标尺
    par(mfrow = c(1, 2))
    plot(rda_tb, scaling = 1, main = 'I 型标尺', display = c('wa', 'sp', 'cn'))
    rda_sp.scaling1 <- scores(rda_tb, choices = 1:2, scaling = 1, display = 'sp')
    arrows(0, 0, rda_sp.scaling1[ ,1], rda_sp.scaling1[ ,2], length =  0, lty = 1, col = 'red')
    plot(rda_tb, scaling = 2, main = 'II 型标尺', display = c('wa', 'sp', 'cn'))
    rda_sp.scaling2 <- scores(rda_tb, choices = 1:2, scaling = 2, display = 'sp')
    arrows(0, 0, rda_sp.scaling2[ ,1], rda_sp.scaling2[ ,2], length =  0, lty = 1, col = 'red')
    
    #隐藏物种信息,以 I 型标尺为例展示双序图,并查看分别使用物种加权计算的样方坐标以及拟合的样方坐标的差异
    par(mfrow = c(1, 2))
    plot(rda_tb, scaling = 1, main = 'I 型标尺,加权', display = c('wa', 'cn'))
    plot(rda_tb, scaling = 1, main = 'I 型标尺,拟合', display = c('lc', 'cn'))
    
    ##RDA 结果提取
    #scores() 提取排序得分(坐标),以 I 型标尺为例,分别提取前两个约束轴中的样方(排序对象)、物种(响应变量)、环境因子(解释变量)的排序坐标
    rda_site.scaling1 <- scores(rda_tb, choices = 1:2, scaling = 1, display = 'wa')	#使用物种加权和计算的样方坐标
    rda_sp.scaling1 <- scores(rda_tb, choices = 1:2, scaling = 1, display = 'sp')
    rda_env.scaling1 <- scores(rda_tb, choices = 1:2, scaling = 1, display = 'cn')
    #若需要输出在本地,以样方坐标为例,以 csv 格式为例
    write.csv(data.frame(rda_site.scaling1), 'rda_site.scaling1.csv')
    
    #coef() 提取典范系数
    rda_coef <- coef(rda_tb)
    
    #其它提取方式,首先查看结果包含的所有信息
    names(rda_tb.scaling1)
    #然后提取相关内容,例如我们想提取前两轴的“样方得分”,可如此做
    rda_site.scaling1 <- rda_tb.scaling1$sites[ ,1:2]
    #该结果和上述“scores(rda_tb, choices = 1:2, scaling = 1, display = 'wa')”的结果是一致的
    #同理,可如此输出在本地
    write.csv(data.frame(rda_site.scaling1), 'rda_site.scaling1.csv')
    
    ##R2 校正
    #RsquareAdj() 提取 R2
    r2 <- RsquareAdj(rda_tb)
    rda_noadj <- r2$r.squared	#原始 R2
    rda_adj <- r2$adj.r.squared	#校正后的 R2
    
    #关于约束轴的解释量,应当在 R2 校正后手动计算
    
    ##置换检验
    #所有约束轴的置换检验,以 999 次为例
    rda_tb_test <- anova(rda_tb, permutations = 999)
    #或者使用
    rda_tb_test <- anova.cca(rda_tb, step = 1000)
    
    #各约束轴逐一检验,以 999 次为例
    rda_tb_test_axis <- anova(rda_tb, by = 'axis', permutations = 999)
    #或者使用
    rda_tb_test_axis <- anova.cca(rda_tb, by = 'axis', step = 1000)
    
    #p 值校正(Bonferroni 为例)
    rda_tb_test_axis$`Pr(>F)` <- p.adjust(rda_tb_test_axis$`Pr(>F)`, method = 'bonferroni')
    
    ##断棍模型和 Kaiser-Guttman 准则确定残差轴
    #提取残差特征值
    pca_eig <- rda_tb$CA$eig
    
    #Kaiser-Guttman 准则
    pca_eig[pca_eig > mean(pca_eig)]
    
    #断棍模型
    n <- length(pca_eig)
    bsm <- data.frame(j=seq(1:n), p = 0)
    bsm$p[1] <- 1/n
    for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i))
    bsm$p <- 100*bsm$p/n
    bsm
    
    # 绘制每轴的特征根和方差百分比 
    par(mfrow = c(2, 1))
    barplot(pca_eig, main = '特征根', col = 'bisque', las = 2)
    abline(h = mean(pca_eig), col = 'red')
    legend('topright', '平均特征根', lwd = 1, col = 2, bty = 'n')
    barplot(t(cbind(100 * pca_eig/sum(pca_eig), bsm$p[n:1])), beside = TRUE, main = '% 变差', col = c('bisque', 2), las = 2)
    legend('topright', c('% 特征根', '断棍模型'), pch = 15, col = c('bisque', 2), bty = 'n')
    
    ##变量选择
    #计算方差膨胀因子
    vif.cca(rda_tb)
    
    #vegan 包 ordistep() 前向选择,基于 999 次置换检验
    rda_tb_forward_p <- ordistep(rda(phylum_hel~1, env, scale = FALSE), scope = formula(rda_tb), direction = 'forward', permutations = 999)
    
    #vegan 包 ordiR2step() 前向选择,基于 999 次置换检验
    rda_tb_forward_r <- ordiR2step(rda(phylum_hel~1, env, scale = FALSE), scope = formula(rda_tb), R2scope = rda_adj, direction = 'forward', permutations = 999)
    
    #以 rda_tb 和 rda_tb_forward_r 为例,简要绘制双序图比较变量选择前后结果
    par(mfrow = c(1, 2))
    plot(rda_tb, scaling = 1, main = '原始模型,I 型标尺', display = c('wa', 'cn'))
    plot(rda_tb_forward_r, scaling = 1, main = '前向选择后,I 型标尺', display = c('wa', 'cn'))
    
    #细节部分查看
    summary(rda_tb_forward_r, scaling = 1)
    
    #比较选择前后校正后 R2 的差异
    RsquareAdj(rda_tb)$adj.r.squared
    RsquareAdj(rda_tb_forward_r)$adj.r.squared
    
    #packfor 包 forward.sel() 前向选择
    library(packfor)
    forward.sel(phylum_hel, env, adjR2thresh = rda_adj)
    
    ##变差分解 varpart(),以前向选择后的简约模型 rda_tb_forward_r 为例(包含 6 个环境解释变量)
    #以两组环境变量为例,运行变差分解
    rda_tb_forward_vp <- varpart(phylum_hel, env['pH'], env[c('DOC', 'SOM', 'AP', 'AK', 'NH4')])
    rda_tb_forward_vp
    
    plot(rda_tb_forward_vp, digits = 2, Xnames = c('pH', 'CNPK'), bg = c('blue', 'red'))
    
    #查看前向选择中被剔除的环境变量“TC”,与这 6 个被保留的环境变量之间解释变差的“共享程度”
    rda_tb_forward_vp <- varpart(phylum_hel, env['TC'], env[c('pH', 'DOC', 'SOM', 'AP', 'AK', 'NH4')])
    plot(rda_tb_forward_vp, digits = 2, Xnames = c('TC', 'forward_env'), bg = c('blue', 'red'))
    
    #解释变差的置换检验,以 pH 所能解释的全部变差为例;999 次置换
    anova(rda(phylum_hel, env['pH']), permutations = 999)
    #若考虑 pH 单独解释的变差部分,需将其它变量作为协变量;999 次置换
    anova(rda(phylum_hel, env['pH'], env[c('DOC', 'SOM', 'AP', 'AK', 'NH4')]), permutations = 999)
    
    ##plot() 作图示例,以前向选择后的简约模型 rda_tb_forward_r 为例,展示前两轴,II 型标尺,双序图,默认使用物种加权和计算的样方坐标
    pdf('rda_test1.pdf', width = 5, height = 5)
    plot(rda_tb_forward_r, choices = c(1, 2), scaling = 2, type = 'n')
    text(rda_tb_forward_r, choices = c(1, 2), scaling = 2, dis = 'cn', col = 'blue', cex = 0.8)
    points(rda_tb_forward_r, choices = c(1, 2), scaling = 2, pch = 21, bg = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), col = NA, cex = 1.2)
    dev.off()
    
    ##ggplot2 作图,以前向选择后的简约模型 rda_tb_forward_r 为例,展示前两轴,II 型标尺,双序图,默认使用物种加权和计算的样方坐标
    #提取样方和环境因子排序坐标,前两轴,II 型标尺
    rda_tb_forward_r.scaling2 <- summary(rda_tb_forward_r, scaling = 2)
    rda_tb_forward_r.site <- data.frame(rda_tb_forward_r.scaling2$sites)[1:2]
    rda_tb_forward_r.env <- data.frame(rda_tb_forward_r.scaling2$biplot)[1:2]
    
    #读取样本分组数据(附件“group.txt”)
    group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
    
    #合并样本分组信息,构建 ggplot2 作图数据集
    rda_tb_forward_r.site$sample <- rownames(rda_tb_forward_r.site)
    rda_tb_forward_r.site <- merge(rda_tb_forward_r.site, group, by = 'sample')
    
    rda_tb_forward_r.env$sample <- NA
    rda_tb_forward_r.env$group <- rownames(rda_tb_forward_r.env)
    
    #ggplot2 作图
    library(ggplot2)
    
    p <- ggplot(rda_tb_forward_r.site, aes(RDA1, RDA2)) +
    geom_point(aes(color = group)) +
    scale_color_manual(values = c('red', 'orange', 'green3')) +
    theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.title = element_blank(), legend.key = element_rect(fill = 'transparent')) + 
    labs(x = 'RDA1 (42.91%)', y = 'RDA2 (9.80%)') +
    geom_vline(xintercept = 0, color = 'gray', size = 0.5) + 
    geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
    geom_segment(data = rda_tb_forward_r.env, aes(x = 0,y = 0, xend = RDA1,yend = RDA2), arrow = arrow(length = unit(0.1, 'cm')), size = 0.3, color = 'blue') +
    geom_text(data = rda_tb_forward_r.env, aes(RDA1 * 1.1, RDA2 * 1.1, label = group), color = 'blue', size = 3)
    
    ggsave('rda_test2.pdf', p, width = 5, height = 4)
    ggsave('rda_test2.png', p, width = 5, height = 4)
    

    evn说明:

    ENV    描述
    pH    土壤 pH 值
    TC    土壤总碳(Total Carbon)含量
    DOC    土壤溶解性有机碳(Dissolved Organic Carbon)含量
    SOM    土壤有机质(Soil Organic Matter)含量
    TN    土壤总氮(Total Nitrogen)含量
    NO3    土壤硝态氮(Nitrate Nitrogen)含量
    NH4    土壤氨态氮(Ammonia Nitrogen)含量
    AP    土壤速效磷(Available Phosphorous)含量
    AK    土壤速效钾(Available Potassium)含量
    数据获取地址:

    链接:https://pan.baidu.com/s/1Mt5Sr5mA7PQ7XYrHJHy2pw 
    提取码:qruo 
     

    展开全文
  • R语言绘图实战:RDA冗余分析

    万次阅读 多人点赞 2017-05-13 10:36:24
    #如果小于3.0, RDA的结果会更合理(基于线性模型,冗余分析) #以RDA分析为例 sp0 (sp ~ 1, se) sp0 plot(sp0) #加入所有环境变量排序,RDA分析 sp1 (sp ~ ., se) sp1 plot(sp1) #到这里RDA图已经出来了,很多文章...
    #载入vegan包
    library(vegan)
    #读取“样本-物种”文件
    sp <- read.table(file=file.choose(),sep="\t",header=T,row.names=1)
    sp
    #读取“样本-环境因子”文件
    se <- read.table(file=file.choose(),sep="\t",header=T,row.names=1)
    se
    #选择用RDA还是CCA分析?先用“样本-物种”文件做DCA分析!
    decorana(sp) 
    #根据看分析结果中Axis Lengths的第一轴的大小
    #如果大于4.0,就应选CCA(基于单峰模型,典范对应分析)
    #如果在3.0-4.0之间,选RDA和CCA均可
    #如果小于3.0, RDA的结果会更合理(基于线性模型,冗余分析)
    #以RDA分析为例
    sp0 <- rda(sp ~ 1, se)  
    sp0
    plot(sp0)
    #加入所有环境变量排序,RDA分析
    sp1 <- rda(sp ~ ., se)  
    sp1
    plot(sp1)
    #到这里RDA图已经出来了,很多文章也都直接放这张图,但如果想追求更美的话,还可以找ggplot2包借个衣服,包装自己
    #准备作图数据:提取RDA分析结果的数据,作为新图形元素
    new<-sp1$CCA
    new
    #提取并转换“样本”数据
    samples<-data.frame(sample=row.names(new$u),RDA1=new$u[,1],RDA2=new$u[,2])
    samples
    #提取并转换“物种”数据
    species<-data.frame(spece=row.names(new$v),RDA1=new$v[,1],RDA2=new$v[,2])
    species
    #提取并转换“环境因子”数据
    envi<-data.frame(en=row.names(new$biplot),RDA1=new$biplot[,1],RDA2=new$biplot[,2])
    envi
    #构建环境因子直线坐标
    line_x = c(0,envi[1,2],0,envi[2,2],0,envi[3,2],0,envi[4,2],0,envi[5,2],0,envi[6,2])
    line_x
    line_y = c(0,envi[1,3],0,envi[2,3],0,envi[3,3],0,envi[4,3],0,envi[5,3],0,envi[6,3])
    line_y
    line_g = c("pH","pH","T","T","S2","S2","NH4","NH4","NO2","NO2","Fe2","Fe2")
    line_g
    line_data = data.frame(x=line_x,y=line_y,group=line_g)
    line_data
    #载入ggplot2包
    library(ggplot2)
    #开始重绘RDA图
    #填充样本数据,分别以RDA1,RDA2为X,Y轴,不同样本以颜色区分
    ggplot(data=samples,aes(RDA1,RDA2)) + geom_point(aes(color=sample),size=2) +
    #填充微生物物种数据,不同物种以图形区分
     geom_point(data=species,aes(shape=spece),size=2) + 
    #填充环境因子数据,直接展示
     geom_text(data=envi,aes(label=en),color="blue") +
    #添加0刻度纵横线
     geom_hline(yintercept=0) + geom_vline(xintercept=0)+
    #添加原点指向环境因子的直线
     geom_line(data=line_data,aes(x=x,y=y,group=group),color="green") +
    #去除背景颜色及多余网格线
     theme_bw() + theme(panel.grid=element_blank())
    #大功告成,保存为矢量图等等
    ggsave("RDA2.PDF")

    展开全文
  • 欢迎转载~ 工具Github链接地址 ... 工具功能 ...第一张图表格是预制体的使用资源情况,第二张图是分析预制体的依赖资源具体内容,主要是图片MD5作为key值存于字典中,字典的值为{图集Tag,SpriteAtla...
  • 传统BIRA结构存在多次地址比较的问题,严重影响了存储器的修复速度与读写性能。为了解决这一问题,提出了基于布鲁姆过滤器的BIRA技术。新型BIRA结构在传统结构的基础上增加了一个布鲁姆过滤器,通过减少地址比较次数...
  • 通过分析冗余并联机器人的驱动支链,建立动力学模型,得到其质量矩阵,
  • Bundle冗余分析

    2017-11-21 16:08:47
    http://blog.csdn.net/akof1314/article/details/78141789
  • 参考分析过程和结论总结部分——矩阵表(重要)、相关系数表、冗余分析表.pdf
  • 面向网络结构的语义推理冗余分析,王国栋,原野,合理构建本体,减少推理产生的冗余信息可有效提高语义网推理效率。本文以农业AGROVOC本体为例,首先将语义网中的概念及关系抽象成网络
  • 关于微生物群落功能冗余的思考,鲁海燕,徐艾诗,土壤是微生物的大本营,它承载了大部分生命的基因多样性。微生物群落在各种生态进程中具有重要作用,然而对于微生物多样性与执行
  • #资源达人分享计划#
  • canoco中文简明教程

    2018-04-14 16:45:17
    Canoco for Windows 是新一代的 CANOCO 软件,是生态学应用软件中用于约束与非约 束排序的最流行工具。Canoco for Windows 整合了排序以及回归和排列方法学,...● PrCoord:对特定数据集进行主坐标分析以及冗余分析
  • 常见的一种特征筛选手段,可以从大量变量中筛选特征变量实现保留变量与目标之间的最大相关性,而彼此间的重复信息最小
  • 冗余代码检测与分析

    2017-10-09 18:03:00
    \本文要点\\代码冗余的原因多种多样,从未使用的变量到未完成的变更,再到废弃的代码;\\t冗余代码会产生一系列的影响,包括源代码臃肿、可靠性及可维护性降低。在某些情况下,死代码也会影响性能;\\t为了检测冗余...
  •  在一个桥接的局域网里,为了增强可靠性,必然要建立一个冗余的路径,网段会用冗余的网桥连接。但是,在一个透明桥桥接的网络里,存在冗余的路径就能建立一个桥回路,桥回路对于一个局域网是致命的。它会带来如下...
  • 数据库冗余数据分析

    千次阅读 2014-04-06 22:55:59
    数据库物理层面的冗余指数据库存储的硬件资源的冗余,逻辑结构层面的冗余是指包括表、记录、字段、属性值以及索引、数据字典中的冗余,由于数据库逻辑实现的基础是各种硬件资源,所以物理层面的冗余影响数据库逻辑...
  • UPS的冗余性技术分析

    2020-08-19 21:53:25
    作为一种电子设备,UPS本身没有容错能力。传统在线式UPS系统虽然实现了蓄能供电及旁路转换过程的不间断供电,但随着电力需求标准的提高,用户渴望得到更为安全的UPS系统...
  • CANOCO5.0教程

    2018-11-30 13:05:48
    该文档是CANOCA5.0教程资源,不过是英文教程,目前市场上较完整中文教程还不多见。
  • Unity AssetBundle 冗余检测与资源分析

    千次阅读 2017-09-30 13:10:15
    检测冗余可以在未打包前对将要打包的资源做分析,但是这无法完全保证打包之后的 AssetBundle 完全无冗余,一是分析时无法保证正确无冗余,二是引用的内置资源无法剔除冗余,所以对打包之后的 AssetBundle 包进行检测...
  • 对如何实现安全制动并联冗余的回油通道进行了设计和分析,采用每一回油通道均并联2个阀位监测电磁换向阀,即使有一个电磁换向阀出现换向故障,另一个还能实现本次安全制动,并监测到换向故障的电磁换向阀,能报警并闭锁...
  • CANOCO-CCA分析简明教程

    2017-12-14 01:39:57
    CCA分析简明教程,详细介绍了如何用CANOCO进行CCA分析
  • 主要介绍了mysql重复索引与冗余索引,简单说明了重复索引与冗余索引的概念、应用场景并结合实例形式分析了mysql重复索引与冗余索引相关操作技巧,需要的朋友可以参考下
  • canoco5+教程

    2019-03-03 11:21:36
    生 态 学 软 件 canoco 5 运行INSTALL.EXE安装,使用以下试用许可证: Name: Any Comany: Any License code: c599999 根据系统位数选择Runasdate_x64或者Runasdate_x86,复制到某个软件存放位置 ...
  • 在一个桥接的局域网里,为了增强可靠性,必然要建立一个冗余的路径,网段会用冗余的网桥连接。但是,在一个透明桥桥接的网络里,存在冗余的路径就能建立一个桥回路,桥回路对于一个局域网是致命的。
  • 第 2卷 第 3期 华北科技学院学报 2005年 9月 循环冗余校验算法分析和实现 杜杏菁 ,刘春梅 (华北科技学院计算机系 ,北京东燕郊 101601) 摘 要 :在网络中传输报文时 ,噪声干扰或传输中断等因素往往使接收端收到的报文...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 174,092
精华内容 69,636
关键字:

冗余分析