CCA和RDA分析—那些你知道和不知道的

“知其然，必知其所以然”

"R实战"专题·第14篇

  编辑 | free傻孩子

  4456字 |12分钟阅读

本期推送内容

本节我们将为大家介绍两种非对称约束排序分析方法，分别是典范对应分析（CCA）冗余分析（RDA）。两者的计算原理我们就不再过多的赘述，网络上已经有很多相关介绍了。本节我们的重点主要是这两个排序分析能做什么？要注意什么？以及如何在R中绘图。（本节内容均为干货，建议收藏  转发）

1、如何选择RDA还是CCA

无论是RDA分析还是CCA分析，都需要两个数据框（data.frame）[已经对行名和列名进行了定义，行名为样方，列名为物种变量或者环境变量]，分别是环境因子数据和物种数据。

我们面临的第一问题是如何选择这两种排序分析，曾经有人建议使用vegan包里面的decorana函数来判断是选择RDA还是CCA。具体的，通过运行该函数，如果前4个轴的所有轴长度均小于3则选择RDA分析，如果4个轴中存在一个轴的长度大于4则选择CCA分析。

代码如下：

library(tidyverse)
library(vegan)
#load data
data("varespec")#species
data("varechem")#env
#选择RDA和CCA的依据
decorana(veg = varespec)

#Call:
#  decorana(veg = varespec)

#Detrended correspondence analysis with 26 segments.
#Rescaling of axes with 4 iterations.

#                DCA1 DCA2 DCA3 DCA4
#Eigenvalues 0.5235 0.3253 0.20010 0.19176
#Decorana values 0.5249 0.1572 0.09669 0.06075
#Axis lengths 2.8161 2.2054 1.54650 1.64864

很明显，我们看到所有的轴的长度均小于3，因此通过该方法预测，我们应该选择RDA分析。当然，如果你并不相信以上运算的结果，比较稳健的方法是同时计算RDA和CCA比较这两个分析中环境因子对物种分布的解释量大小，选择解释量大的那种方法即可。

2、RDA和CCA

test.rda<-rda(varespec,varechem)#RDA分析
test.rda#79.97%

test.cca<-cca(varespec,varechem)#RDA分析
test.cca#69.20%

通过同时做RDA和CCA我们发现，在RDA中环境因子对物种分布的解释量更高。

对结果的解读：

constrained（约束）指自变量（环境）矩阵能对因变量矩阵（物种）的整体解释量，如RDA分析中的79.97%和CCA分析中的69.20%。

unconstrained（非约束）指还剩下的没有被解释的部分，如RDA分析中的20.03%和CCA分析中的30.80%。

如何获得每一个轴的解释量：

1）使用每个轴对应的特征值（eigenvalues）除以对应的Total inertia

如：在RDA 分析中RDA1对物种分布的解释量 = 820.1/1825.6594 * 100% = 44.92%

在CCA分析中，CCA1对物种分布的解释量 = 0.4389/2.0832 *100% = 21.07%

2) 使用summary函数，代码和结果如下：

summary(test.rda)

注意：RDA和CCA中，每一轴对结果的解释量为“Eigenvalues, their contribution to the variance”（引擎值以及它们对方差的贡献）下的结果，而不是“Accumulated constrained eigenvalues”（累积约束引擎值）下的结果。后者，指的是每一个轴在整体解释量中的重要性。比如，RDA分析中整体解释量为79.97%，RDA1轴的解释量为44.92%，那么RDA1轴对整体解释量的重要性为：44.92/79.97 * 100%  = 56.17%

如何获得RDA和CCA分析的校正解释量和未校正解释量：

运行以下代码可以获得

RsquareAdj(test.rda)

我们看到RDA的校正解释量为48.80%；未校正解释量为79.97%

同理，CCA中校正解释量和非校正解释量也能用以上的方法算得（注：在后续的所有适用于RDA的分析同样都适用于CCA）。

校正RDA1和RDA2的解释量分别为：

RDA1: 0.4880 * 0.4492 = 0.2192，也就是21.92%

RDA2: 0.4880*0.2187 = 0.1067 ，也就是10.67%

上面这个获得每一轴校正解释量的方法是，（数量生态学）中给出的方法。但是，我觉得这种直接使用校正解释量乘以每个轴的解释量，获得每轴校正解释量的方法似乎并不符合逻辑。

我认为应该

1) 先计算校正值即：0.4880/0.7997 = 0.6102

2) 计算每一轴的校正解释量，即：

RDA1: 0.6102*0.4492 = 0.2741

RDA2: 0.6102*0.2187 = 0.1335

（如有错误之处，敬请指正）

如何获知，RDA分析中和CCA分析中整体模型是否显著以及每个变量的重要性

使用以下代码可以验证模型的显著性和每个变量的重要性：

permutest(test.rda,permu=999)

通过999次置换检验证明了在RDA中环境因子显著影响物种分布（p = 0.002）。

每个因子的重要性：

envfit(test.rda,varechem,permu=999)

结果表明，Mo对物种分布的影响不显著，可能是不重要的环境因子。

简单出图，对照以上结果：

plot(test.rda,display=c("si","bp"),scaling=3)#只显示样方数据位置

3、结果绘图

scrs<-scores(test.rda,display=c("sp","wa","lc","bp","cn"))
#a list
#计算箭头变化倍数并计算箭头的长度
multiplier <- vegan:::ordiArrowMul(scrs$biplot)
df_arrows<- scrs$biplot*multiplier
df_arrows %>%
data.frame() ->df_arrows2
df_arrows2

#提取样方坐标（一般我们都是通过样方来绘RDA的而不是物种）

sites <- scrs$sites %>%
as_tibble(rownames = "sample")
sites

这里sample下的数字就是我们的样方名，我们用的是vegan中自带的数据，所以没有样方名，只是一串数字，如果是自己的数据，那么这里应该是具体的样方名。

#分组
#我们根据环境变量varechem中的Humdepth对数据分成两个组，操作如下
varechem %>%
as_tibble(rownames = "sample") %>%
select(sample, Humdepth) %>%
mutate(group = if_else(Humdepth < 2.2,"deep","shallow")) ->group
group

#把分组加到样方数据（sites）中并最终绘图
sites %>%
left_join(group, by = "sample") %>%
ggplot(aes(x = RDA1, y=RDA2))+
geom_hline(aes(yintercept=0),colour="#d8d6d6",linetype=5)+
geom_vline(aes(xintercept=0),colour="#d8d6d6",linetype=5)+
geom_point(aes(color=group),size=3.5,shape=16) +
geom_segment(data=df_arrows2, aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
arrow = arrow(angle= 15,length = unit(0.25, "cm")),color="#000000",alpha=0.5)+#画箭头
geom_text(data=as.data.frame(df_arrows*1.1),aes(RDA1, RDA2, label = rownames(df_arrows2)),
color="#000000",alpha=0.7)+#箭头对应的信息
scale_color_manual(values = c("#005991","#910e00"))+#点的颜色
#scale_x_continuous(limits = c(-6,10))+
#scale_y_continuous(limits = c(-12,5))+
labs(x = "RDA1 (27.41%)", y = "RDA2 (13.35%)")+
theme_bw() +
theme(panel.grid=element_blank())

注意：如果不使用校正解释量也是可以的。那样的话这个RDA分析应该是这样的：

#把分组加到样方数据（sites）中并最终绘图
sites %>%
left_join(group, by = "sample") %>%
ggplot(aes(x = RDA1, y=RDA2))+
geom_hline(aes(yintercept=0),colour="#d8d6d6",linetype=5)+
geom_vline(aes(xintercept=0),colour="#d8d6d6",linetype=5)+
geom_point(aes(color=group),size=3.5,shape=16) +
geom_segment(data=df_arrows2, aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
arrow = arrow(angle= 15,length = unit(0.25, "cm")),color="#000000",alpha=0.5)+#画箭头
geom_text(data=as.data.frame(df_arrows*1.1),aes(RDA1, RDA2, label = rownames(df_arrows2)),
color="#000000",alpha=0.7)+#箭头对应的信息
scale_color_manual(values = c("#005991","#910e00"))+#点的颜色
#scale_x_continuous(limits = c(-6,10))+
#scale_y_continuous(limits = c(-12,5))+
labs(x = "RDA1 (44.92%)", y = "RDA2 (21.87%)")+
theme_bw() +
theme(panel.grid=element_blank())

创作不易：欢迎收藏、转发。