RDA 综合教程
pmgh
2024-12-30
1.RDA基本理论
基本概念
冗余分析（redundancy analysis, RDA）是响应变量矩阵与解释变量矩阵之间多元多重线性回归的拟合值矩阵的PCA分析，也是多响应变量（multi-response）回归分析的拓展。简单一点来说，RDA是通过线性回归分析结合主成分分析的排序方法，目的是寻找能最大程度解释响应变量矩阵变差的一些列的解释变量的线性组合，也就是环境对于样本的影响，因此RDA是被解释变量约束的排序。
方差分解（Partitioning of variance）：总方差被划分为约束和非约束两部分。约束部分表示响应变量Y矩阵的总方差能被解释变量解释的部分，如果用比例来表示，其值相当于多元回归的R2。在RDA中，这个解释比例值也称作双多元冗余统计。然而，类似多元回归的未校正的R2，RDA的R^2是有偏差的，需要进一步校正。
特征根以及对方差的贡献率(Eigenvalues, and their contribution to the variance ):结果中的RDA即为典范轴，PC为非约束轴，输出结果不仅包含每轴特征根同时也给出累积方差解释率（约束轴）或承载轴（非约束轴），最终的累计值必定是1，所有典范轴累积解释率也代表响应变量总方差能够被解释变量解释的部分。两个特征根的重要区别：典范特征根RDAx是响应变量总方差能够被解释变量解释的部分，而残差特征根RCx响应变量总方差能够被残差轴解释的部分，与RDA无关。
累积约束特征根（Accumulated constrained eigenvalues）表示在本轴以及前面所有轴的典范轴所能解释的方差占全部解释方差的比例累积。
物种得分（Species scores）双序图和三序图内代表响应变量的箭头的顶点坐标。与PCA相同，坐标依赖标尺Scaling的选择。
样方得分（Site scores (weighted sums of species scores)）物种得分的加权和：使用响应变量矩阵Y计算获得的样方坐标。
样方约束——解释变量的线性组合（Site constraints (linear combinations of constraining variables)）：使用解释变量矩阵X计算获得的样方坐标，是拟合的（fitted）样方坐标。
解释变量双序图得分（Biplot scores for constraining variables）：排序图内解释解释变量箭头的坐标，按照下面的过程获得：运行解释变量与拟合的样方坐标之间的相关分析，然后将所有相关系数转化为双序图内坐标。所有的变量包括k个水平的因子口可以有自己的坐标对因子变量在排序轴的坐标，用各个因子的形心表示更合适。
因子解释变量形心（Centroids for factor constraints）：因子变量各个水平形心点的坐标，即每个水平所用标识为一的样方的形心。
典范系数（canonical coefficients）：即每个解释变量与每个典范轴之间的回归系数。几个重要的作用：（1）解释变量对约束轴的贡献：典范系数表示每个解释变量对每个约束轴的贡献大小和方向。较大的绝对值表示该解释变量对该轴的贡献较大，正负符号表示正相关或负相关。（2）模型解释：通过查看典范系数，可以了解哪些解释变量对响应变量的变异有显著影响。这有助于解释环境变量如何影响物种组成或其他响应变量。（3）变量选择：典范系数可以帮助识别哪些解释变量在模型中最为重要，从而进行变量选择和模型优化。（4）生态解释：在生态学研究中，典范系数可以帮助解释环境因子如何影响生态系统的结构和功能。例如，某个环境因子的高典范系数可能表明它对物种分布有重要影响。
方差膨胀因子（variance inflation factor,VIF）：是衡量一个变量的回归系数的方差由共线性引起的膨胀比例。如果VIF值超过20，表示共线性很严重。实际上，VIF超过10则可能会有共线性问题，需要处理。
2.问题准备
研究背景
对三种环境中的土壤进行了采样（环境A、B、C，各环境下采集9个土壤样本），结合二代测序观测了土壤细菌群落物种组成，并对多种土壤理化指标进行了测量，相应结果包括如下几个数据：
- RDA-01_1_phylum_table.txt(csv):
通过16S测序所得的细菌类群丰度表格，列名为34个细菌“门水平”，行为各门在27个样点中的丰度
- RDA-01_2_env_table.txt(csv):
27个样点的土壤理化性质数据(RDA-01_4_env说明)，列名含义如下
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）含量
- RDA-01_3_group.txt:
采样点划分数据
问题假设
这里将“细菌门”简化为物种，注意这是由研究目的决定的，因为OTU水平的数量决定的，划分越细，数据量越大
将“理化性质”简化为环境因子
研究目的
判断物种组成与分布对环境因子的响应关系，同时判断哪些环境因子的影响更为显著，由此采用RDA进行分析
3.执行RDA
说明
rda函数调用格式有两种
rea(Y, X, W)或rda(Y~var1+var2+factorA+var2*var3+Condition(var4))。对第1种，Y为响应变量矩阵，X为解释变量矩阵，W为协变量矩阵。这种写法虽然简单，但要求解释变量矩阵或协变量矩阵中不能含有因子变量（定性变量）（当仅存在响应变量矩阵Y时，则执行PCA排序）。第2种写法（推荐）中，Y为响应变量矩阵，var1、var2、var3为定量解释变量，factorA为因子变量，var2*var3为变量2和变量3的交互作用项，var4为协变量。scale=FALSE，使用拟合值的协方差矩阵运算RDA；若TRUE，则基于相关矩阵运算RDA。
“~.”表示将数据框env中所有的列变量作为解释变量带入RDA排序（在不考虑协变量以及变量间的交互作用项的情况下，可以直接使用），能够自动识别定量解释变量以及定性解释变量类型。
tb-RDA（基于转化的冗余分析）采用Hellinger转化的数据进行RDA分析，适用于物种丰度表中出现很多“0”值的情况（一般情况下都用这个）
偏RDA，用于在控制变量矩阵W对响应变量Y的影响的前提下，另一组变量矩阵X对响应变量Y的影响上。如本例中，控制土壤pH影响后（pH作为协变量），观测其它环境因素的影响。
常规RDA基于欧氏距离，db-RDA基于其他距离进行冗余分析。db-RDA首先基于原始物种数据计算相异矩阵，作为PCoA（主坐标分析）的输入，之后将所有PCoA排序轴上的样方得分矩阵用于执行RDA，好处是可以基于任意一种距离测度（例如Bray-curtis距离、Unifrac距离等，而不局限于欧氏距离）进行RDA排序，极大拓宽了RDA的适用范围，因此db-RDA在生态学统计分析中使用广泛。
使用cmdscale()基于获得的Bray-curtis距离矩阵运行PCoA，其中add=TRUE表示为距离矩阵加一个常数，避免产生负的特征根，也称为Cailliez校正。当然，若我们已经存在了某个距离矩阵文件（例如已经提前计算好的Bray-curtis距离矩阵文件），也可将该文件直接导入至R中并运行PCoA排序，方法与上述一致。PCoA运算完毕后，提取样方排序得分（赋值给矩阵pcoa_site，主坐标矩阵视为表征数据总变差的距离矩阵）输入至rda()中作为响应变量，执行RDA。
非线性RDA：当响应变量与解释变量明显为非线性关系（例如单峰关系，当然，若响应变量与解释变量存在明显的单峰响应模式，可使用CCA来解决），常规的线性RDA模型不适合时，可考虑此方法。
代码
# 加载包
library(vegan)
library(ggplot2)
library(magrittr)
library(rdacca.hp) #变差分解
# 读取数据
phylum <- read.delim('./#data/RDA-01_1_phylum_table.txt', row.names = 1) #物种数据
env <- read.delim('./#data/RDA-01_2_env_table.txt', row.names = 1) #环境数据
# 判断适用模型
# 根据看分析结果中Axis Lengths的第一轴的大小
# 如果大于4.0,就应选CA/CCA（基于单峰模型，典范对应分析）
# 如果在3.0-4.0之间，选PCA/RDA或者CA/CCA均可
# 如果小于3.0, CA/RDA的结果会更合理（基于线性模型，冗余分析）
decorana(phylum) ##
## Call:
## decorana(veg = phylum)
##
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
## Total inertia (scaled Chi-square): 0.1097
##
## DCA1 DCA2 DCA3 DCA4
## Eigenvalues 0.06344 0.011844 0.009860 0.007282
## Additive Eigenvalues 0.06344 0.011477 0.009837 0.007368
## Decorana values 0.06345 0.009089 0.003835 0.002045
## Axis lengths 0.67326 0.319171 0.275812 0.361713# 计算RDA
# 1.基础RDA
rda_result0 <- rda(phylum~., env, scale = FALSE) #全部环境变量RDA分析
rda_result1 <- rda(phylum~pH+TC+TN, env, scale = FALSE) #指定环境变量
rda_result2 <- rda(phylum~pH+TC+TC*TN, env, scale = FALSE) #指定交互变量
rda_result3 <- rda(phylum~TC+Condition(pH), env, scale = FALSE) #指定协变量（控制PH后看影响）
# 2.tb-RDA（加入数据转化）
phylum_hel <- decostand(phylum, method = 'hellinger') #物种数据Hellinger预转化
rda_tb <- rda(phylum_hel~., env, scale = FALSE) #使用全部数据进行
# 也有用log1p转化数据的，也是为了去除0值，但是比较复杂
# 3.偏RDA（加入控制变量）
rda_part <- rda(phylum_hel~TC+DOC+SOM+TN+Condition(pH), data = env, scale = FALSE)
# 4.db-RDA（基于非欧距离）
dis_bray <- vegdist(phylum, method = 'bray') #计算Bray-curtis距离，无需Hellinger预转化
# dis_bray <- as.dist(read.delim('bray.txt', row.names = 1, sep = '\t'))#可读取已有距离
pcoa <- cmdscale(dis_bray, k = nrow(phylum)-1, eig = TRUE, add = TRUE) #运行PCoA并Cailliez校正
pcoa_site <- pcoa$point #提取 PCoA 样方得分（坐标）
rda_db <- rda(pcoa_site, env, scale = FALSE) #db-RDA，使用全部的环境数据
# 4.db-RDA（基于距离的冗余分析）（vegan内置方法）
rda_db <- capscale(phylum~., env, distance = 'bray', add = TRUE) #效果同上
# 5 非线性RDA
# 这里数据不适合，详情见赖江山老师译作“数量生态学:R语言的应用”170-174页内容。4.结果分析
说明
- summary()展示结果，可选参数scaling表示缩放: scaling = 0：不缩放 scaling = 1：适用于样本之间的距离更感兴趣的情况，通常用于解释变量为因子的情况。 scaling = 2：适用于变量之间的相关关系更感兴趣的情况，通常用于解释变量为连续变量的情况。
代码
# 1.对于RDA、tb-RDA、db-RDA等
rda_tb <- rda(phylum_hel~., env, scale = FALSE)
rda_tb.scaling1 <- summary(rda_tb, scaling = 1)
# rda_tb.scaling1 #基本结果
# rda_tb.scaling1$species #物种得分
# rda_tb.scaling1$sites #样方得分
# rda_tb.scaling1$constraints #样方约束得分，拟合值
# rda_tb.scaling1$biplot #解释变量得分
# 对偏RDA，有以下不同
phylum_hel <- decostand(phylum, method = 'hellinger') #物种数据Hellinger预转化
rda_part <- rda(phylum_hel~TC+DOC+SOM+TN+Condition(pH), data = env, scale = FALSE)
rda_part.scaling1 <- summary(rda_part, scaling = 1)
# rda_part.scaling1结果解读
- rda_tb.scaling1:
- Call: 所用RDA命令
- Partitioning of variance: 方差分解，下面分别为约束轴解释的部分，以及非约束轴解释的部分
- RDA为约束轴，RDA的Eigenvalue和为0.01971，Proportion和为0.7338
- PCA为非约束轴，PCA的的Eigenvalue和为0.00715，Proportion和为0.2662
- Accumulated constrained eigenvalues，这里RDA的Proportion和Cumulative Proportion为相对的，解释方差占约束轴能够解释的总方差的比例
- rda_tb.scaling1$species，物种得分，即响应变量（物种数据）在约束轴中的排序坐标
- rda_tb.scaling1$sites，物种得分，即响应变量（物种数据）在各约束轴中的排序坐标
- rda_tb.scaling1$constraints，样方约束，解释变量的线性组合，使用解释变量矩阵获得的排序对象（样本）在各约束轴中的排序坐标，是你和的样方坐标
- rda_tb.scaling1$biplot，解释变量得分，即解释变量（环境因子）在各约束轴中的排序坐标
- centroids of factor constraints，当存在定性变量时，结果中会出现，是因子变量各水平形心点的坐标，本例中不存在。
- 对于偏RDA，方差分解多了一项“conditioned”，表示已经被协变量（此处为PH）解释的部分，而约束部分则表示为剔除协变量外，能被解释变量解释的部分
结果（tbRDA）表明，该RDA模型承载总方差的比例为73.38%，其中RDA1和RDA2的解释量总和为53.445%+11.847%=65.292%，即承载了响应变量矩阵65.29%的方差（因为scale=FALSE了）。虽然结果很好，但是该结果未经校正和检验，目前不可靠。
5.信息提取
代码
# 信息提取函数主要有以下几个
# 1. scores，提取排序得分，即坐标，display参数与plot中的相同
# scores(rda, choices=1:2, scaling=2, display='wa') #样点，等于"sites"，等于rda$sites
# scores(rda, choices=1:2, scaling=2, display='sp') #响应变量，等于"species"，等于rda$species
# scores(rda, choices=1:2, scaling=2, display='bp') #解释变量，等于rda$biplot
# scores(rda, choices=1:2, scaling=2, display='cn') #约束轴分数，空的
# scores(rda, choices=1:2, scaling=2, display='lc') #线性组合分数
# scores(rda, choices=1:2, scaling=2, display='reg') #回归分数
# scores(rda, choices=1:2, scaling=2, display='all') #全部
# 2. coef()，提取典范系数
# coef(rda_tb)
# 3. 常规方法
# 注意根据?summary.cca，scaling默认是"species"，即2
# rda_tb.scaling1 <- summary(rda_tb, scaling = 1)
# rda_tb.scaling1$sites[ ,1:2] #等于scores(rda_tb, choices = 1:2, scaling = 1, display = 'wa')6.R2和p值校正、置换检验
说明
约束模型解释的总变差（即响应变量矩阵的总方差能够被解释变量解释的部分，如果用比例表示，即等同于多元回归的R2，也就是约束轴的总解释量）以及各约束轴解释量等，不建议直接在原始结果中获取。原因是初始所得的RDA模型需要经过校正和检验后才可使用，否则可能会带来较大的偏差
代码
# RsquareAdj() 提取 R2
r2 <- RsquareAdj(rda_tb)
rda_noadj <- r2$r.squared #原始 R2
rda_adj <- r2$adj.r.squared #校正后的 R2结果解释
rda_tb.scaling1的结果如下 Partitioning of variance: Inertia Proportion Total 0.02686 1.0000 Constrained 0.01971 0.7338 Unconstrained 0.00715 0.2662
使用调整后R2乘以对应方差，即为修正后的方差： + 总方差：0.5928×0.02686＝0.01592 + 解释方差占约束轴能够解释的总方差的比例：0.5928×0.72835=0.43177，0.5928×0.16145=0.09571
因此，尽管在R2校正后约束模型的总体解释量有所下降（由0.7338下降至0.5928），但看起来似乎仍很可观。特别是的前两个约束轴的解释量仍然非常好：RDA1和RDA2的解释量总和为“0.43177+0.09571=0.52748”，即承载了响应变量矩阵52.75%的方差，达到了一半以上。
说明
至此RDA的结果仍然不可靠，还需要建议约束轴的显著性。
首先需对全模型执行置换检验查看显著性，仅当全模型通过检验时，才能依次对每一个约束轴执行检验判断是否能够保留。否则，RDA模型则从全局上来讲就不合理，本次的RDA结果不能再被使用，需要考虑更改分析方案（更换排序模型，或者再尝试对数据进行有效的转化后使用等）。
vegan包中执行置换检验的命令是anova()，注意这个方法全称是anova.cca()，可使用?anova.cca()查看帮助：
+ anova()中，用permutations参数指定置换次数 +
anova.cca()中，用step参数设定产生参照分布的样本总数，step数值减1即为实际的置换次数
+
默认情况下，对所有约束轴的置换检验，可以轴逐一检验用by = 'axis'
+ 此外，anova()还可以通过by = ’
margin’检验协变量等，详情见?anova.cca()
在进行多重比较时（例如在多个统计测试中），不进行校正会增加出现假阳性结果的概率。P值校正方法（如Bonferroni校正）可以帮助控制这种错误率，确保结果的可靠性。
代码
# 置换检验
# 所有约束轴的置换检验，以 999 次为例，以下两种方式相同
rda_tb_test <- anova(rda_tb, permutations = 999)
rda_tb_test <- anova.cca(rda_tb, step = 1000)
# rda_tb_test
# 各约束轴逐一检验，以 999 次为例，以下两种方式相同
rda_tb_test_axis <- anova(rda_tb, by = 'axis', permutations = 999)
rda_tb_test_axis <- anova.cca(rda_tb, by = 'axis', step = 1000)
# rda_tb_test_axis
# p 值校正（Bonferroni 为例）
rda_tb_test_axis$`Pr(>F)` <- p.adjust(rda_tb_test_axis$`Pr(>F)`, method = 'bonferroni')
# rda_tb_test_axis结果解释
置换检验的结果，根据rda_tb_test的p显示，模型整体是显著的（Pr(>F)为0.001）。但是根据rda_tb_test_axis仅第一轴（Pr(>F)为0.001）和第二轴（Pr(>F)为0.004）是显著的，后面几个不显著需要被剔除。
使用rda_tb_test_axis校正后，第一轴（Pr(>F)为0.005）和第二轴（Pr(>F)为0.020）是显著的
根据结果，仅能提取前两个约束轴的信息进行更进一步的分析。但实际上，根据校正后的R2，RDA1和RDA2的解释量总和为“0.43177+0.09571=0.52748”（即承载了响应变量矩阵52.75%的方差，超过一半；再结合置换检验结果，即便不考虑其它的约束轴，约束模型所能解释的方差占比也较高），表明给定的环境变量对群落物种组成分布具有相当高的解释度。事实上，即便其他约束轴通过检验，一般也不将它们列为考虑的范畴，原因是它们的解释量很低，强制加入模型解释物种组成成因可能会带来混乱的结果。我们更期望的还是具有较高解释量的约束轴能够被保留，以及样方或物种的排序分布具有明显趋势的模型不要被舍弃等。
注意，非约束模型中排序轴取舍的常用方法不适用于约束模型。例如断棍模型、凯撒格特曼规则（Kaiser-Guttman rule）等，我们可能在PCA分析中常用它们作为排序轴的选择标准，但它们在RDA约束轴的取舍中意义不大。对于约束轴的取舍，还需通过置换检验获得。然而，对于RDA的非约束残差轴，断棍模型等仍然适用。
7.非约束残差轴重要性确定
说明
当RDA模型的承载方差比例较低，存在相当一部分比例的残差未参与解释时，可能需要审视非约束残差轴。
此处以tb_rda为例，提取非约束残差轴的特征值，并分别依据断棍模型和Kaiser-Guttman准则确定残差轴。
代码
# 断棍模型和 Kaiser-Guttman 准则确定残差轴
pca_eig <- rda_tb$CA$eig #提取残差特征值
# pca_eig[pca_eig > mean(pca_eig)] #Kaiser-Guttman 准则
#断棍模型
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
# 绘制每轴的特征根和方差百分比
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')
结果解释
结果表明，综合考虑，大约3个残差非约束轴是比较重要的。若有必要，我们可能需要绘制前2-3个非约束轴的排序图，查看为何还有相当一部分方差未能解释。
8.变量选择（前向选择）
说明
约束排序中，有时解释变量太多，需要想办法减少解释变量的数量。一般来讲，减少解释变量的数量可能有两个并不冲突的原因：第一是寻求简约的模型，利于我们对模型的解读；第二是当解释变量过多时会导致模型混乱，例如有些解释变量之间可能存在较强的线性相关，即共线性问题，可能会造成回归系数不稳定。处理步骤如下：
- 使用方差膨胀因子（Variance Inflation Factor,
VIF）度量多重共线性，若VIF>10，则有共线性问题，若VIF>20则共线性很严重。vegan包用
vif.cca()计算VIF。 - 使用用vegan包中ordistep()函数、ordiR2step()函数，或者packfor包中的forward.sel()函数前向选择。 3.变量选择中的注意事项：尽管变量选择策略具有明显的优势，但一定切记不应盲目依赖自动选择程序在回归模型中选择相关的环境变量，因为可能得到生态学上不相关的模型，或者被忽视的其它变量组合（当前保留变量集中的部分解释变量+某些被剔除的解释变量）在某种程度上也可能比当前模型更具代表性（Legendre和Legendre，1998）。即仅通过统计学手段得到的最优变量集合可能并非代表了最具生态学意义的模型。
1.ordistep()函数
说明
ordistep()函数，执行前向选择（或后向，具体见?ordistep()），过程大致如下：
- 依次分别运行每个解释变量与响应变量的RDA排序。
- 选择“最好”的解释变量。对于ordistep()，根据置换检验中F统计量的显著性水平（p值）选择最好的变量，即最小的p值。如果遇到p值相等的情况，拥有最小AIC信息准则（Akaike information criterion）的变量应该入选。最好的变量也就是最显著的变量。
- 寻找模型中第二个、第三个、第四个……解释变量。上一步选取了只含有一个最好变量的模型，现在再加入某一个剩下的解释变量。对于ordistep()，选择偏RDA置换检验最小p值的模型，即表示新加入的变量对模型的贡献最显著。如果遇到p值相等的情况，具有最小AIC值的模型应该入选。
- 这个过程一直持续，直到无显著的解释变量为止。即如果加入新变量的偏RDA置换检验显著性p值大于或等于α（例如，以α=0.05为准），选择过程即被终止。
但是ordistep()的选择标准过于宽松，有时会选择显著但不包含任何变量的模型（夸大I类错误），或选择包括过多变量的模型（夸大被解释方差的量）。为了解决这个问题，可将Borcard终止准则（Blanchet等，2008）引入作为前向选择的第二个终止原则，加强版的前向选择程序，即ordiR2step()和forward.sel()。
在使用ordistep()时，若是发现结果似乎异常，例如变量选择后的校正R2高于了全模型的校正R2，则推荐使用ordiR2step()或forward.sel()执行变量选择过程（其实也推荐在一开始就直接使用这种加强版的前向选择程序ordiR2step()或forward.sel()，忽略ordistep()）。
代码
# 计算方差膨胀因子
# vif.cca(rda_tb) #结果显示TC、DOC、SOM、TN、AP、AK都超过10
# vegan 包 ordistep() 前向选择
# 其中，phylum_hel~1，意为从env（环境解释变量数据集矩阵）中的第一个变量开始执行选择
# direction = 'forward'，前向选择
# permutations = 999，999次置换。
rda_tb_forward_p <- ordistep(rda(phylum_hel~1, env, scale = FALSE), scope = formula(rda_tb), direction = 'forward', permutations = 999)
# 细节部分查看
# summary(rda_tb_forward_p, scaling = 1)
# 比较选择前后校正后 R2 的差异
# RsquareAdj(rda_tb)$adj.r.squared
# RsquareAdj(rda_tb_forward_p)$adj.r.squared
# 简要绘制双序图比较变量选择前后结果
par(mfrow = c(1, 2))
plot(rda_tb, scaling = 1, main = '原始模型，I 型标尺', display = c('wa', 'cn'))
plot(rda_tb_forward_p, scaling = 1, main = '前向选择后，I 型标尺', display = c('wa', 'cn'))
结果解释
全模型中共计9个环境变量，经过ordistep()前向选择后保留了其中的3个，总解释量从0.5928396变为0.461555几乎没有牺牲解释能力。
2.ordiR2step()函数
说明
ordiR2step()的前向选择原理与ordistep()类似，但引入全模型R2adj作为第二个终止原则（即增添了Borcard终止准则）：如果当前所选模型的R2adj达到或超过全模型的R2adj，或备选变量的偏RDA置换检验p值不显著，则变量选择将停止。
和ordistep()用法一致，除了前向选择，ordiR2step()也可以执行后向选择和逐步选择，并且可以直接将rda()的输出作为输入（此外，ordiR2step()也同样可用于CCA模型，vegan包cca()的输出）。详细请使用?ordiR2step()。
代码
# vegan 包 ordiR2step() 前向选择，其中：
# phylum_hel~1意为从env（环境解释变量数据集矩阵）中的第一个变量开始执行选择
# rda_adj即为RDA全模型的校正R2，在上文中已经通过RsquareAdj()获得
# scope = formula(rda_tb)即指定全模型的校正R2作为判断的依据
# direction = 'forward'，前向选择
# permutations = 999，999次置换。
rda_tb_forward_r <- ordiR2step(rda(phylum_hel~1, env, scale = FALSE), scope = formula(rda_tb), R2scope = rda_adj, direction = 'forward', permutations = 999)
# 细节部分查看
# summary(rda_tb_forward_r, scaling = 1)
# 比较选择前后校正后 R2 的差异
# RsquareAdj(rda_tb)$adj.r.squared
# RsquareAdj(rda_tb_forward_r)$adj.r.squared
# 简要绘制双序图比较变量选择前后结果
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'))
结果解释
全模型中共计9个环境变量，经过ordiR2step()前向选择后保留了其中的6个，总解释量从0.5928396变为0.5814848几乎没有牺牲解释能力。
3.forward.sel()函数
说明
由于forward.sel()和ordiR2step()的选择原理相似，都是引入了Borcard终止准则后的加强版的选择程序，所以一般情况下二者前向选择的结果是一致的（尽管具有不同的参数设置逻辑）。 但是与ordiR2step()相比，forward.sel()目前仅适用于RDA（改良后的ordiR2step()既可用于RDA，vegan包rda()的输出；也可用于CCA，vegan包cca()的输出），且仅可用前向选择（ordiR2step()同时支持前向、后向、逐步选择），所以看起来似乎并没有ordiR2step()好用
packfor包安装失败
9.变差分解（方差分解）
说明
生态学研究中，属于不同类别的两组或两组以上解释变量共同解释一组响应变量时，分析每一个解释变量所单独能解释的变差（或者两组或多组变量共同解释的变差），需要用到变差分解（variation partitioning，或称方差分解，variance partitioning）。
以两组解释变量A和B为例，响应变量Y的总变差可以分解为以下几个部分：
- 变量A能够解释的总变差（A的部分，变量A的边际效应，marginal_effect，即当约束模型中仅存在A时能够被解释的总变差）；
- 变量B能够解释的总变差；
- 仅能被变量A所单独解释的变差（A排除B的部分，即变量A的局部效应，partial_effect，即将变量B作为协变量处理时，变量A所能解释的变差部分）；
- 仅能被变量B所单独解释的变差；
- 变量A和B共同解释的变差（AB相交部分，由于变量A和B具有相关性，无法归其所属）；
- 变量A和B均无法解释的变差（整体排除A、B的部分），即残差。
vegan包内的varpart()命令可用于执行变差分解过程。这里我们以上述经过ordiR2step()前向选择后所得简约tb-RDA模型“rda_tb_forward_r”为例（包含6个环境解释变量），简要展示varpart()的用法。对于此命令的详细信息请使用?varpart()查看帮助文档。
变差分解的一个主要目的是计算共同解释的变差，但在解读时必须非常谨慎，因为共同部分很难确定与哪组变量有关。
变差分解共同解释部分（相交部分）与方差分析中因子的交互作用（interaction）不同。交互作用是双因素方差分析中，一个因子不同水平与其它因子不同水平之间的协同作用。变差分解是由于解释变量共线性引起解释变差的重叠，不是协同作用的结果。换句话说，如果两个（组）变量相互独立（不相关，正交），则共同解释部分将为0。
代码
# 以两组环境变量为例，运行变差分解
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'))
结果解释
上文在前向选择后所得简约RDA模型“rda_tb_forward”中，保留了6个环境变量。我们根据环境因子的属性，将它们划分为两组：一组为土壤化学属性，pH；另一组则包括了土壤碳氮磷钾含量，DOC、SOM、AP、AK、NH4。对这两组环境变量执行变差分解后，结果如下。
其中，pH为0.010+0.122，DOC、SOM、AP、AK、NH4为0.122+0.449，这两组总共为0.581，等于ordiR2step() 前向选择模型的R2.adj= 0.5814848
使用变差分解查看多重共线性
对未经过变量前向选择的原始RDA模型执行变差分解，通过比较变量选择前后所得模型中各组环境变量解释变差的差异，被剔除的环境变量所解释变差的大小及其和其它变量间交叉区域的大小等，以据此加深对通过变量选择构建简约模型的必要性的理解。
例如，我们发现在经过变量的前向选择后，“TC”这个环境变量被剔除了。如前所述，由于有些解释变量之间可能存在较强的线性相关，即共线性问题，可能会造成回归系数不稳定，所以前向选择程序一般会将这些变量剔除。那么，对于“TC”这个环境变量来讲，是否是因为这个原因呢？下面我们通过变差分解的方式查看。
# 查看前向选择中被剔除的环境变量“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'))
结果表明，TC所解释的变差均与这6个被保留的（前向选择后所保留的）环境变量之间产生交集，即存在非常强的共线性问题。而对于TC本身来讲，它所能单独解释的变差部分居然是0。所以，前向选择程序自动将它剔除了。
但是有一点需要特别注意，正如上文提到的，尽管变量选择策略具有明显的优势，但自动选择程序所选择的变量一般来讲仅供参考，因为仅通过统计学手段得到的最优变量集合可能并非代表了最具生态学意义的模型。
使用置换检验判断显著性
还可以通过置换检验，检验各解释变差部分的显著性。例如，对pH所解释变差部分的置换检验如下所示（置换检验命令见上文所述anova()）。
# 解释变差的置换检验，以 pH 所能解释的全部变差为例做999 次置换
anova(rda(phylum_hel, env['pH']), permutations = 999) #结果p=0.06显著## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(X = phylum_hel, Y = env["pH"])
## Df Variance F Pr(>F)
## Model 1 0.0044511 4.9667 0.016 *
## Residual 25 0.0224048
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 若考虑 pH 单独解释的变差部分（排除相交部分），需将其它变量作为协变量；999 次置换
anova(rda(phylum_hel, env['pH'], env[-1]), permutations = 999) #结果p=0.39不显著## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(X = phylum_hel, Y = env["pH"], Z = env[-1])
## Df Variance F Pr(>F)
## Model 1 0.0004210 1.001 0.395
## Residual 17 0.0071496rdacca.hp包变差分解
rdacca.hp(phylum_hel, env, method = "RDA", type = "adjR2", var.part = F) #如果var.part=T，则展示联合解释## $Method_Type
## [1] "RDA" "adjR2"
##
## $Total_explained_variation
## [1] 0.593
##
## $Hier.part
## Unique Average.share Individual I.perc(%)
## pH 0.0000 0.0979 0.0979 16.51
## TC -0.0108 0.0122 0.0014 0.24
## DOC 0.1484 -0.0253 0.1231 20.76
## SOM 0.0296 -0.0116 0.0180 3.04
## TN -0.0128 0.0134 0.0006 0.10
## NO3 0.0367 0.0040 0.0407 6.86
## NH4 0.1310 -0.0588 0.0722 12.18
## AP 0.0782 -0.0137 0.0645 10.88
## AK 0.2999 -0.1254 0.1745 29.43
##
## attr(,"class")
## [1] "rdaccahp"10.RDA排序图
说明
在作图前，首先需要确认以下信息。
- 根据置换检验，已知前两个约束轴是显著的，则可以使用；
- 根据“RsquareAdj(rda_tb_forward_r)$adj.r.squared”可知简约模型的校正后R2为“0.5814848”，然后再根据“summary(rda_tb_forward_r)”的结果可知约束轴RDA1和RDA2的承载方差分别占约束模型总方差的“0.73786”和“0.168480”，那么由此推断约束轴RDA1和RDA2的解释量分别约为“0.5814848×0.73786 = 42.91%”和“0.5814848×0.168480 = 9.80%”。解释量可观；
- 根据所关注的重点，确定使用哪种标尺。例如这次想重点观测变量间的相关性，那么就以II型标尺为例展示；
- 使用“plot(rda_tb_forward_r, scaling = 2, display = c(‘wa’, ‘cn’))”查看简约作图结果，排序空间内的样方分布是否区分明显，环境变量对排序空间的贡献是否达到了预期等。
1.ordiplot绘图
代码
# ordiplot是vegan包自带绘图函数
ordiplot(rda_tb)
2.plot绘图
说明
plot()命令中，display参数用于定义在排序图中展示哪些信息，具体见
?plot.cca()。其中：- sp：物种，
- wa：使用物种加权和计算的样方坐标，
- lc：以拟合值（解释变量的线性组合）计算样方的坐标，
- cn：约束成分（即解释变量）。可知对于样方坐标有两种可选展示方式，
scores()用于提取物种、样方或环境变量的排序坐标，同样以display参数定义。
vegan默认绘图中，排序向量的位置与summary有所不同，原因是做了同比例缩放，以便于展示和分析
代码1（基础绘图）
# 对于RDA（更多画图说明，查看plot.cca）
rda_tb <- rda(phylum_hel~., env, scale = FALSE)
rda_tb.scaling1 <- summary(rda_tb, scaling = 1)
par(mfrow = c(1, 2))
# 图1，1型标尺
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')
# 图2，2型标尺
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')
代码2（简要信息绘图）
# 隐藏物种信息，以 I 型标尺为例展示双序图，并查看分别使用物种加权计算的样方坐标（wa）以及拟合的样方坐标（lc）的差异
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'))
代码3（添加点）
# 以前向选择后的简约模型 rda_tb_forward_r 为例，展示前两轴，II 型标尺，双序图，默认使用物种加权和计算的样方坐标
# rda_tb_forward_r <- ordiR2step(rda(phylum_hel~1, env, scale = FALSE), scope = formula(rda_tb), R2scope = rda_adj, direction = 'forward', permutations = 999) #建立前向选择模型
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)
3.ggplot2绘图
代码1（教程自带）
# 使用教程自带数据，代码略有修改
rda.summ <- summary(rda_tb_forward_r, scaling = 2)
plot.sit <- rda.summ$sites %>% data.frame() %>% '['(1:2) #提取样本
plot.env <- rda.summ$biplot %>% data.frame() %>% '['(1:2) #提取环境因子
group <- read.delim('./#data/RDA-01_3_group.txt', sep='\t', stringsAsFactors=F) #读取样本分组数据
plot.sit$sample <- rownames(plot.sit) #添加样本名
plot.sit <- merge(plot.sit, group, by = 'sample') #合并分组信息
plot.env$index <- rownames(plot.env) #添加环境因子信息
rda.contri <- rda.summ$cont$importance[2,1:2]*100 #提取贡献率
labs.x <- sprintf("RDA1 (%.2f%%)", rda.contri[1])
labs.y <- sprintf("RDA1 (%.2f%%)", rda.contri[2])
#画图
ggplot(plot.sit, aes(RDA1, RDA2)) +
geom_point(aes(color = group), size=3) +
geom_segment(data = plot.env, aes(x = 0,y = 0, xend = RDA1, yend = RDA2),
arrow = arrow(length = unit(0.2, 'cm')), linewidth = 0.6, color = 'blue') +
geom_text(data = plot.env, aes(RDA1*1.1, RDA2*1.1, label=index), color='blue', size=3)+
scale_color_manual(values = c('red', 'orange', 'green3')) +
geom_vline(xintercept = 0, color = 'gray', linewidth = 0.5) +
geom_hline(yintercept = 0, color = 'gray', linewidth = 0.5) +
labs(x = labs.x, y = labs.y) +
theme_bw()+
theme(aspect.ratio = 1/1,panel.grid = element_blank()) 
代码2（论文RDA绘图代码）
# 加载包
library(data.table)
library(ggplot2)
library(magrittr)
library(vegan)
library(ggrepel)
library(boot)
# 读取数据，这里前几个Al2O3不一样，是因为三个月的数据，先矫正再平均的
rda_tributary <- fread("./#data/RDA-04_1_CSP_tributary.csv")
# rda_tributary列名
#1 id #2 label #3 Subbasin #4 sampling_site #5 dem_min
#6 dem_max #7 dem_dif #8 Car_percent #9 Sli_percent #10 esmean
#11 mPre_year #12 mPre_month #13 mTem_year #14 mTem_month #15 mEt0_year
#16 mEt0_month #17 AI_year
#18 K2O #19 Na2O #20 CaO #21 MgO
#22 all_k #23 all_na #24 all_ca #25 all_mg #26 SCP
# 生成数据
col_name <- names(rda_tributary)
col_id <- c(1,2)
col_y <- c(22,23,24,25) #4个被解释变量，CSP钙镁纳钾
col_x <- c(7,8,9,10,17,13) #6个解释变量：demdif/car/sli/es/ai/tem
col_y_n <- length(col_y)
col_x_n <- length(col_x)
y.colname <- c("K2O","Na2O","CaO","MgO")
x.colname <- c("Elevation Drop", "Carbonate rock", "Siliciclastic rock",
"Soil Erosion", "Aridity Index", "Annual temperature")
dax <- rda_tributary[, .SD, .SDcols = col_x] #年气温、干旱指数
day <- rda_tributary[, .SD, .SDcols = col_y] #CSP
dan <- rda_tributary$Subbasin
# bootstrap抽样
# 定义一个函数，用于返回引导样本的均值
bootstrap_function <- function(data, indices) {
sample_data <- data[indices, ]
return(as.vector(t(sample_data)))
}
dat <- cbind(day,dax,dan)
results <- boot(data = dat, statistic = bootstrap_function, R = 50)
datbs <- matrix(results$t[,1:ncol(dat)], ncol = ncol(dat), byrow = F)
dan2 <- datbs[,col_y_n+col_x_n+1] %>% as.data.frame() %>% set_colnames("name") #提取样点名
datbs.n <- datbs[,-(col_y_n+col_x_n+1)] %>% as.numeric() %>% matrix(nrow=nrow(datbs)) #转化为数值
day2 <- datbs.n[,1:col_y_n] %>% as.data.frame() %>% set_colnames(y.colname)
dax2 <- datbs.n[,(col_y_n+1):(col_y_n+col_x_n)] %>% as.data.frame() %>% set_colnames(x.colname)
# 计算RDA，有hellinger转换，则属于tb-RDA
day2.h <- decostand(day2,method = "hellinger") # Hellinger转化，适用于物种0值较多情况，
rda.r <- rda(day2.h~.,dax2)
rda.s <- summary(rda.r, scaling = 2)
# 绘图数据
xz <- as.data.frame(rda.s$biplot[,1:2]) #解释变量
yz <- as.data.frame(rda.s$species[,1:2]) #响应变量
st <- as.data.frame(rda.s$site[,1:2]) #样点
# 文字标记位置调整
st_text <- st
st_text$l <- dan2$name
st_text <- unique(st_text)
# RDA轴贡献标记
# 注意不要用下面这个，会导致末位0不显示
# labs(x=paste("RDA 1 (", format(100 *rda.s$cont[[1]][2,1], digits=4), "%)", sep=""),
# y=paste("RDA 2 (", format(100 *rda.s$cont[[1]][2,2], digits=4), "%)", sep=""))
rda.contri <- rda.s$cont$importance[2,1:2]*100 #提取贡献率
labs.x <- sprintf("RDA1 (%.2f%%)", rda.contri[1])
labs.y <- sprintf("RDA1 (%.2f%%)", rda.contri[2])
# 有无文字画两次图，因为文字有可能重叠，没有文字的图，需手动ps文字
ggplot() +
geom_point(data = st,aes(RDA1,RDA2),color='gray',size=4)+
geom_segment(data = yz,aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
arrow = arrow(angle=15,length = unit(0.35,"cm"), type = "closed"),
linetype=1, linewidth=0.6, colour="red")+
geom_segment(data = xz,aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
arrow = arrow(angle=15,length = unit(0.35,"cm"), type = "closed"),
linetype=1, linewidth=0.6,colour = "blue")+
geom_text_repel(data = st_text,aes(RDA1,RDA2,label=l),color='gray',size=3)+
geom_text_repel(data = yz, aes(RDA1, RDA2), colour="red", label=names(day), parse = T)+
geom_text_repel(data = xz,aes(RDA1,RDA2,label=row.names(xz)),colour = "blue")+
geom_hline(yintercept=0,linetype=3,linewidth=1)+
geom_vline(xintercept=0,linetype=3,linewidth=1)+
# scale_x_continuous(limits = c(-0.62, 0.62), breaks = round(seq(-0.6,0.6,by=0.2),1))+
# scale_y_continuous(limits = c(-0.82, 0.82), breaks = seq(-0.8,0.8,by=0.2))+
labs(x=labs.x, y=labs.y)+
theme_bw()+
theme(panel.grid=element_blank(),
text = element_text(family = "serif"))
11.完整案例
说明
重点展示代码和图，次要展示结果和解释
代码
# 加载包
library(vegan)
library(ggplot2)
library(ggrepel)
library(gridExtra)
### 1.数据准备
grp <- read.table('./#data/RDA-05_1_grp.txt', row.names = 1, header = T, sep = '\t')
otu <- read.table("./#data/RDA-05_2_OTU.txt",header=T,row.names = 1, sep = '\t') #读取因变量矩阵
env <- grp[,-c(1:2)]
otu <- otu[,rownames(grp)] %>% .[rowSums(.) !=0,] %>% t() #去除均为0的行，对列排序，转置
otu.helli <- decostand(otu, method = "hellinger") #标准化
### 2. 模型选择
decorana(otu.helli) #执行DCA去趋势分析：结果都小于3，适合做RDA
### 3. 前向选择构建模型
fit.0 <- rda (otu.helli ~ 1, data = env) #因变量矩阵和指示变量逐个进行RDA分析
fit.all <- rda (otu.helli ~ ., data = env) #因变量矩阵与所有环境变量进行RDA分析
fit <- ordiR2step (fit.0, scope = formula(fit.all)) #执行前向选择
### 4. 模型检验
total <- summary(fit) #查看RDA分析结果
vif.cca(fit) #方差膨胀因子：结果没有大于10，无线性相关
summary(eigenvals(fit)) #查看每个排序轴的特征值和能够解释的变异
anova.cca(fit) #查看RDA模型显著性：结果0.001，显著
anova.cca(fit, by="term", parallel=2) #查看解释变量的显著性：结果前两个AN为0.001，AP为0.005，显著
anova.cca(fit, by="axis") #查看排序轴的显著性：结果前两轴RDA1为0.001，RDA2为0.001，显著
RsquareAdj(fit) #查看R2：结果R2为0.4625488，调整R2为0.3132568
coef(fit) #查看典范系数
### 5. 各个环境因子对整体变异的贡献率 #####
# 当前模型（经过前向选择后）包含的环境因子
env.fac <- attr(fit$terms, "term.labels")
# 当前模型的解释能力
expl.tot <- total$tot.chi #总体方差
expl.con <- total$constr.chi #变量所能解释的方差
expl.unc <- total$unconst.chi #未解释方差
# 当前模型约束轴内，各环境变量的解释能力
expl.con.fac <- c()
for (i in env.fac){
di <- setdiff(env.fac, i)
trda <- rda(otu.helli, env[i], env[di]) #将其他因子作为控制变量输入RDA
expl.con.fac <- c(expl.con.fac, trda$CCA$eig) #结果为环境因子单独解释的特征根
}
expl.con.fac <- c(expl.con.fac, expl.con-sum(expl.con.fac)) #再添一项，为joint部分
names(expl.con.fac) <- c(env.fac, "Joint")
expl.con.fac <- expl.con.fac/expl.con
# 生成画图数据
# 总贡献
prob.com <- c(expl.con, expl.unc)/expl.tot
df.tot <- data.table(class = c("Constrained", "Unconstrained"), prob = prob.com)
df.tot[, lab_y:=cumsum(df.tot$prob) - 0.5*df.tot$prob]
df.tot[, lab_t:=sprintf("%s (%.2f%%)",class,prob*100)]
# 约束贡献
expl.con.fac <- sort(expl.con.fac)
df.con <- data.table(class=factor(names(expl.con.fac), levels = names(expl.con.fac)), prob=expl.con.fac)
df.con[, lab_y:=cumsum(df.con$prob) - 0.5*df.con$prob]
df.con[, lab_t:=sprintf("%s (%.2f%%)",class,prob*100)]
p1 <- ggplot(df.tot, aes(x = 1, y = prob, fill = class))+
geom_bar(stat = "identity", position = position_stack(reverse = T), width = 1, show.legend = F)+
geom_text(aes(x = 1.5, y = lab_y, fill = NULL, label = lab_t), color = "black", size = 2.5)+
coord_polar(theta = "y", direction = -1, start = 0)+
labs(x = NULL, y = NULL, title = "Variance explained by each factor")+
theme_bw()+
theme(aspect.ratio = 1/1,
axis.ticks = element_blank(),
axis.text = element_blank(),
panel.grid = element_blank(),
panel.border = element_blank(),
plot.title = element_text(hjust = 0.5))
p2 <- ggplot(df.con, aes(x = 1, y = prob, fill = class))+
geom_bar(stat = "identity", position = position_stack(reverse = T), width = 1, show.legend = F)+
geom_text(aes(x = 1.9, y = lab_y, fill = NULL, label = lab_t), color = "black", size = 3)+
coord_polar(theta = "y", direction = -1, start = 0)+
labs(x = NULL, y = NULL, title = "Variance explained by each factor")+
theme_bw()+
theme(aspect.ratio = 1/1,
axis.ticks = element_blank(),
axis.text = element_blank(),
panel.grid = element_blank(),
panel.border = element_blank(),
plot.title = element_text(hjust = 0.5))
grid.arrange(p1, p2, ncol = 2)
# 注意以上过程与方差分解不同
# rdacca.hpf变差分解
# library(rdacca.hp)
# rdacca.hp(otu.helli, env[env.fac], method = "RDA", type = "adjR2", var.part = F)
# $Method_Type
# [1] "RDA" "adjR2"
# $Total_explained_variation
# [1] 0.313
# $Hier.part
# Unique Average.share Individual I.perc(%)
# AN 0.0150 0.0867 0.1017 32.49
# AP 0.0275 0.0188 0.0463 14.79
# TP 0.0306 0.0130 0.0436 13.93
# AK 0.0413 -0.0010 0.0403 12.88
# TN 0.0251 0.0563 0.0814 26.01
# attr(,"class")
# [1] "rdaccahp"
# vegan变差分解
vp_1 <- varpart(otu.helli, env[env.fac[1]], env[env.fac[-1]])
plot(vp_1, digits = 2, Xnames = c(env.fac[1], 'others'), bg = c('blue', 'red'), main="variation partitioning")
title(main = "variation partitioning")
### 6. 基本绘图
plot(fit, type="n")
text(fit, dis="cn")
points(fit, pch=21, col="blue", bg="blue", cex=2)
text(fit, "sites", col='red', cex=0.5, adj = 0.5, pos=2) 
### 7. 使用ggplot2绘图
# 获取样点坐标
sites <- total$sites[,1:2] %>% data.frame() %>% merge(grp[,1:2],by = 'row.names')
# 获取指示变量坐标
biplot <- total$biplot[,1:2] %>% data.frame()
biplot$env <- rownames(biplot)
# 定义样点组合在图例中的出场顺序
sites$treatment <- factor(sites$treatment, levels = c('S','SH','H'))
# 重命名，须注意与上面的对应关系
levels(sites$treatment) = c('水田','水旱轮作','旱田')
rda.contri <- total$cont$importance[2,1:2]*100 #提取贡献率
labs.x <- sprintf("RDA1 (%.2f%%)", rda.contri[1])
labs.y <- sprintf("RDA1 (%.2f%%)", rda.contri[2])
ggplot(sites, aes(x = RDA1, y = RDA2, color = treatment)) +
geom_point(size = 3) +
stat_ellipse(show.legend = F) + #添加置信区间，点的包围线
geom_segment(data = biplot,
aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
arrow = arrow(length = unit(1/2, 'picas')), lwd = 1,
colour = "blue") + #添加箭头
guides(color = guide_legend(override.aes = list(size=5)))+ #增加图例中点的大小
geom_hline(yintercept=0, linetype=2,color='grey') + #添加经过坐标原点的横线
geom_vline(xintercept=0, linetype=2,color='grey') + #添加经过坐标原点的纵线
geom_text_repel(data = biplot,
aes(x=RDA1,y=RDA2,label=env),
size= 5, fontface='bold',color='black')+ #添加指示变量文本
labs(x = labs.x, y = labs.y, color = '') +
theme_bw() +
theme(legend.position = "inside",
legend.position.inside = c(0.9,0.2),
legend.background = element_blank(),
legend.text = element_text(face = 'bold',color='black',size=12),
axis.title = element_text(face = 'bold',color='black',size=14),
axis.text = element_text(face = 'bold',color='black',size=12),
panel.grid = element_blank())
© 2021-2025, Lcpmgh, All rights reserved.