|||
R语言RDA分析全过程梳理及注意事项
RDA分析不是简单的一行代码,有很多步骤,若舍去这些步骤,RDA分析就没有了意义。通过RDA分析,不仅需要给出哪些因子最重要,还要给出每个因子的解释率,最好有显著性检验。这样结果更明确,更加有说服力,这里梳理一下过程,避免以后分析有所遗漏。RDA分析详细步骤具体看刘尧科学网博客http://blog.sciencenet.cn/blog-3406804-1182489.html 有些语句也用了这里面的。
library(vegan)
library(permute)
library(lattice)
library(ggrepel) #作图使用此包,和ggplot2一起用,使文字尽可能不重叠
library (ggplot2)
datarda<-read.table("clipboard",header=T) # 物种数据
#datarda-hel<- decostand(datarda, method = 'hellinger') # 物种数据有时需要转化,特别是有0出现的情况,基于转化的RDA分析tb-RDA
env<-read.table("clipboard",header=T) #环境数据,有时环境数据需要数据转换以满足数据分析要求,但很多文章里土壤pH是不做转换的
rda_result<-rda(datarda~.,env,scale=FALSE)
summary(rda_result)
r2<-RsquareAdj(rda_result)
rda_noadj<-r2$r.squared #原始R2,不能直接使用
rda_adj <- r2$adj.r.squared #校正后的R2乘以RDA1和2的解释度,才是作图用的实际解释度
anova.cca(rda_result2,step=1000) #全模型检验,有时文章需要提供
anova.cca(rda_result2,by='axis',step=1000) # 逐轴检验,看每条轴是否显著,原则上不显著的轴要舍去
#anova(rda_result, by = "term", permutations=999) #这里的anova不是单因素方差分析,只是名字一样,气人不?这里该语句用于逐因子检验,会给出每个环境因子的r2以及显著性P,但不能排除共线性的问题,即每个因子的r2和P值依赖于环境因子的顺序,排在前面的环境因子可能解释度最高且显著,后面的因子虽然不显著可能是因为被前面的解释过了。所以不能直接使用这部分结果。
vif.cca(rda_result) #计算方差膨胀系数,有些因子膨胀系数过大(达到几十或几百),需要排除掉一些因子,可进行前向选择保留必要因子
rda_tb_forward_r<-ordiR2step(rda(datarda~1,env,scale= FALSE), scope = formula(rda_result), direction = 'forward', permutations = 999) ##前向选择,前向选择会保留一些因子,但总体解释的R2基本不变,但有时共线的是核心因子,保留的因子却是在逻辑上无关紧要的,很麻烦
plot(rda_tb_forward_r, scaling = 1, main = '前向选择后,I 型标尺', display = c('wa', 'cn'))
rda_tb_forward_vp <- varpart(phylum_hel, env['pH'], env[c('DOC', 'SOM', 'AP', 'AK', 'NH4')])##前向选择后对结果执行方差分解,可以将环境因子分成不同组,如该语句将pH归为一组,DOC', 'SOM', 'AP', 'AK', 'NH4 归为一组,结果可以使用韦恩图表示出来。
plot(rda_tb_forward_vp, digits = 2, Xnames = c('TC', 'forward_env'), bg = c('blue', 'red'))##韦恩图
anova(rda(phylum_hel,env['pH'],env[c('DOC','SOM','AP','AK','NH4')]),permutations=999)) ##同样对执行完前向选择的结果进行检验,这里环境因子的R2和显著性可以使用
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')) #图
之前说了,因为共线性问题,环境因子对解释变量是否有影响以及每个因子的贡献不易判断,所以可以使用层次分割来分析每个因子的贡献以及显著性。当然增强回归树模型同样可以做出每个因子贡献,但只能适用于解释变量为一个的情况。赖江山老师开发的新包rdacca.hp可以结合RDA或CCA分析进行层次分割。具体看http://blog.sciencenet.cn/blog-267448-1200984.html
library(devtools)
library(hier.part)
library(scales)
library(rdacca.hp)
rdacca.hp (datarda,env,pieplot = "tv",type="RDA") # 每个解释变量占总变化量的比例
rdacca.hp (datarda,env,pieplot = "tev",type="RDA") #每个解释变量占总解释变化量的比例
permu.hp (datarda,env, type="RDA", permutations=999) #层次分割,可给出所有因子的贡献以及显著性,个人感觉层次分割要比前向选择后显著性检验要简单一些。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 07:09
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社