生态学|R语言-基于距离的RDA分析（dbRDA）--两种不同方法
生态学|R语言-基于距离的RDA分析（dbRDA）–两种不同方法
这两天帮一个同学做dbRDA分析，在此之前我并没有接触过这个方法，去了解了一下发现原理非常简单。
普通RDA是获取群落数据与环境因子之间的线性关系区别就在于，RDA是基于群落多度做的，dbRDA是先基于Bray-Curtis距离（或其他距离）进行主坐标分析（Principal Coordinates Analysis，PCoA），随后将特征值大于0的坐标轴与环境因子进行RDA分析，即dbRDA。
RDA侧重于研究群落结构与环境变量之间的线性关系；
dbRDA侧重于解释环境变量对群落结构变异的贡献，
知道原理了之后，我们就可以在R中一步步进行。当然vegan包中的dbrda函数可以一键执行。
下面就将两种方法都列出来，并使用envfit函数进行置换检验（permutation test），给出dbRDA的画图代码。
方法一：已知原理的分步骤进行，暂且封之为“爱学习版dbRDA”
animal为群落数据，格式为行为样地，列为物种，值为多度；
env为环境数据，格式为行是样地，列是环境变量
#读入物种数据
animal<-read.csv("abundance.csv",header = TRUE, row.names = 1)
#读取环境数据
env<-read.csv("habitat.csv", header = TRUE, row.names = 1)
获得PCoA特征值大于0的坐标轴
#计算Bray-curtis 距离，也可以在method里选其他的
beta_div <- vegdist(animal, method = "bray")
#基于距离进行PCoA分析
pcoa_result <- cmdscale(beta_div, k = nrow(animal)-1, eig = TRUE, add = TRUE)#add=TRUE是校正负的特征值，K这里选择了行数-1
#提取主坐标轴信息
pcoa_scores <- as.data.frame(pcoa_result$points)#特征值>0的坐标轴
将pcoa结果与环境变量一起进行rda分析（这一步之前也可以先筛一下环境变量）
db_rda <- rda(pcoa_scores, env, scale = FALSE)
summary(db_rda)
一起看一下结果，环境变量总共解释了85.04%的群落结构变异。
Eigenvalues, and their contribution to the variance 部分的解释率是指RDA轴对总体（total）的解释率，因此在一些同学的分析中会出现这部分数值非常低的情况。
Accumulated constrained eigenvalues部分的解释率是指RDA轴对于约束惯性(Constrained)部分的解释率。
画图的时候可以注意一下到底标注哪一个解释率。
方法2：dbrda函数简化版
这个方法省去了自己进行PCoA的步骤，我们还是从头开始
#读入物种数据
animal<-read.csv("abundance.csv",header = TRUE, row.names = 1)
#读取环境数据
env<-read.csv("habitat.csv", header = TRUE, row.names = 1)
#计算Bray-curtis 距离，也可以在method里选其他的
beta_div <- vegdist(animal, method = "bray")
进行rda分析
db_rda<- dbrda( beta_div ~ ., env, dist="bray")
summary(db_rda)
~~这里我对比了两种方法的结果，有微微的不同，但不知道问题出在哪里，如果有知道的同学可以评论区告诉我，十分感谢！
置换检验
对dbRDA结果进行检验，评估每个环境变量与排序轴的相关性显著性，最终得到环境变量对群落结构的解释程度和方向。
# 置换检验（评估环境因子的整体显著性）
RDA.perm=permutest(db_rda,permu=999)
RDA.perm
# 检查每个环境因子与群落变化的相关性
RDA.env=envfit(db_rda,env,permu=999)
RDA.env
最后提取绘图用的变量
# 提取样本和环境的得分
site_scores <- scores(db_rda, display = "sites") # 样本点坐标
soil_scores <- scores(db_rda, display = "bp") # 环境因子箭头坐标
# 转换为数据框
df_sites <- as.data.frame(site_scores)
df_soil <- as.data.frame(soil_scores)
# 添加箭头缩放因子
multiplier <-1.1 #记得调整
df_soil <- df_soil * multiplier
# 绘图
a<-ggplot() +
geom_point(data = df_sites, aes(x = dbRDA1, y = dbRDA2,
color =animal$group), alpha = 0.7, size = 4) +
#group是分组变量，可以自行导入
scale_color_brewer(palette = "Set1") +
#添加环境因子箭头
geom_segment(data = df_soil,
aes(x = 0, y = 0, xend = dbRDA1, yend = dbRDA2),
arrow = arrow(length = unit(0.2, "cm")),
color = "darkgreen", alpha = 0.8, linewidth = 0.8) +
geom_text(data = df_soil,
aes(x = dbRDA1 * 1.1, y = dbRDA2 * 1.1, label = rownames(df_soil)),
color = "darkgreen", size = 4, check_overlap = TRUE) +
# 添加坐标轴标签并显示解释率
labs(x = paste0("RDA1 (", round(db_rda$CCA$eig[1]/sum(db_rda$CCA$eig)*100,2), "%)"),
y = paste0("RDA2 (", round(db_rda$CCA$eig[2]/sum(db_rda$CCA$eig)*100,2), "%)"),
title = "") +
# 美化主题
theme_bw() +
theme(panel.grid = element_blank()) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey80") +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey80")
a
