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不是单因素方差分析，只是名字一样，气人不？这里该语句用于逐因子检验，会给出每个环境因子的r²以及显著性P，但不能排除共线性的问题，即每个因子的r²和P值依赖于环境因子的顺序，排在前面的环境因子可能解释度最高且显著，后面的因子虽然不显著可能是因为被前面的解释过了。所以不能直接使用这部分结果。

vif.cca(rda_result) #计算方差膨胀系数，有些因子膨胀系数过大（达到几十或几百），需要排除掉一些因子，可进行前向选择保留必要因子

rda_tb_forward_r<-ordiR2step(rda(datarda~1,env,scale= FALSE), scope = formula(rda_result),  direction = 'forward', permutations = 999) ##前向选择，前向选择会保留一些因子，但总体解释的R²基本不变，但有时共线的是核心因子，保留的因子却是在逻辑上无关紧要的，很麻烦

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)） ##同样对执行完前向选择的结果进行检验，这里环境因子的R²和显著性可以使用

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) #层次分割，可给出所有因子的贡献以及显著性，个人感觉层次分割要比前向选择后显著性检验要简单一些。

转载本文请联系原作者获取授权，同时请注明本文来自牛国祥科学网博客。
链接地址：https://blog.sciencenet.cn/blog-3366717-1217929.html

上一篇：R语言使用ggplot2对冗余分析(RDA)等约束排序分析结果作图
下一篇：R语言PCA分析过程梳理及注意事项