也写在前面
又是思路干涸的一天,写写推文换换脑子,今天就不说废话了,把废话都留给年终总结。我们之前陆续更新了(点击链接就可以跳转):
这一期,我们来更新一下基于RDA、CCA的层次分割分析及绘图。我在之前的推文中讲过,当我们的目标变量是独立变量时,我们就可以基于线性回归的层次分割完成分析(基于glmm.hp完成);当我们的变量是群落矩阵的时候,我们就可以基于RDA、CCA的层次分割分析完成分析(基于rdaccaa.hp完成),详细的解释请点击链接跳转查看对应推文。基于RDA、CCA的层次分割分析是配合RDA或CCA分析进行的,怎么做RDA以及绘图,请查看我之前古老的推文,RDA、CCA的绘图结果以及完成层次分割分析后的绘图分别如下。
RDA和CCA的图
RDA+层次分割的图
那么,问题来了——为什么要做个RDA再做个层次分割呢?
答:其实我们进行RDA或CCA分析的时候,主要目的是为了探讨哪些环境因子影响或驱动了群落的变化。但是我们大多数论文其实只展示了一个RDA的图,我们仅通过指标载荷的大小来判断驱动群落变化的主要因子是哪个,显然这样的方式很是粗放,或者说只是定量的描述了一下,并不能通过这个图定性的说就是某某指标驱动了群落的变化。所以,我们见到的大多数论文其实只做了一半的RDA,即只做了定量分析而没有做定性分析。鉴于此,大佬们在RDA中又加入了置换检验结果来定性描述结果(这个推文后边的代码会提到),但是由于数据分析的不平衡设计,置换检验的结果是会随机变化的,即如果将用于分析的因子的顺序改变,那么靠前的因子更容易出现显著的情况,这就让人觉得很魔幻。不过,只要不改变因子的分析顺序,每次的结果还是一致的。为了解决这个问题,赖江山老师开发了rdaccaa.hp包,这个包可以分解每个因子驱动群落变化的百分比,而且这个百分比不会因为因子顺序的变化而变化,再通过置换检验就可以定量且定性的描述驱动群落变化的因子有哪些。
还有一个问题,层次分割分析为什么不用考虑共线性问题?
答:层次分割分析(Hierarchical Partitioning, HP) 提出了一种比较好解决共线性的方案,它基于所有变量子集的全局分解,并通过计算每个变量在不同子集组合中的独立贡献和共享贡献来量化其重要性。所以,大家要是还是不理解的话,直接引用这几篇论文就可以:
Lai Jiangshan,Zou Yi, Zhang Jinlong, Peres-Neto Pedro (2022) Generalizing hierarchical andvariation partitioning in multiple regression and canonical analyses using therdacca.hp R package. Methods in Ecology and Evolution, 13: 782-788.
Jiangshan Lai, Weijie Zhu, Dongfang, Cui, Lingfeng Mao. 2023. Extension of the glmm.hp Package to Zero-Inflated Generalized Linear Mixed Models and Multiple Regression Journal of Plant Ecology, 16(6):rtad038
Jiangshan Lai, Yi Zou, Shuang Zhang, Xiaoguang Zhang,Lingfeng Mao(2022). glmm.hp: an R package for computing individual effect of predictors in generalized linear mixed models. Journal of Plant Ecology,15(6):1302-1307.
准备数据集
数据准备和RDA分析一样的数据就可以,分别是群落数据和环境因子。小编已上传,文末下载。
RDA CCA分析
关于RDA CCA的分析和绘图,这里就不再赘述,以及什么时候该选用RDA或CCA等,请查看之前的推文
这是推文的链接,点击即可跳转。
层次分割分析
我们直接进入层次分割环节,不过在这之前,我们先读取数据做一个简单的RDA。
PS:想做CCA把16行的rda三个字母换成cca即可
library(vegan)
library(ggrepel)
library(ggplot2)
library(ggpubr)
library(rdacca.hp)
sampledata <- read.csv(file.choose(), head = TRUE, row.names=1)
env <- read.csv(file.choose(), header=TRUE, row.names=1)
sampledata <- t(sampledata)
sampledata <- decostand(sampledata,method = "hellinger")
mite.rda <- rda(sampledata~., env, scale = FALSE)
summary(mite.rda)
plot(mite.rda)
运行完毕后,我们就会得到RDA的关键数据和一幅不好看的图,图的美化还是那句话,看之前的推文就可以
Biplot scores for constraining variables
RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
TC 0.342161 -0.263115 -0.388125 0.137307 0.251537 -0.02199
TN 0.054603 -0.549167 -0.264299 0.343325 -0.316158 0.27765
CN 0.253506 0.302074 -0.085245 -0.206632 0.538819 -0.26552
TP -0.613183 -0.301599 0.005675 -0.005963 -0.169216 0.06092
NO3 -0.548305 -0.475723 0.151527 0.271580 -0.204094 -0.15554
NH4 -0.712661 -0.288958 0.349700 0.360944 -0.189774 -0.05820
AP -0.666277 0.008007 -0.275832 0.407262 -0.241981 0.06987
AK 0.579230 0.331089 -0.503094 -0.252480 0.234978 0.12582
pH 0.069097 0.243470 0.111079 0.119347 0.071269 0.14927
Moisture -0.257933 -0.271649 0.716006 -0.150313 -0.146433 0.01676
BD -0.625180 -0.094909 -0.264975 0.428801 -0.295267 0.09396
shoot 0.009318 -0.294962 -0.097992 -0.326072 0.249450 -0.39597
root -0.425110 0.474072 -0.088597 0.005317 -0.151714 -0.16808
ANPP 0.176580 0.484541 0.165789 0.270270 -0.407461 -0.02564
plant.diversity -0.044953 -0.447974 -0.212302 -0.174326 0.009537 0.14304
然后,就是正式的层次分割了,这一步的计算会很慢,大家数据量越大计算就会越慢,耐心等待即可。
PS:这里如果想对CCA进行层次分割,将method = 'RDA'中的RDA改为CCA即可。此外,type = 'R2'可以选择R2(是未校正R方)和adjR2(校正后R方)。到底选择哪个,看大家的结果哪个和预期符合就用哪个,除非审稿人非要让用校正后的
#层次分割(这一步计算很慢)
mite.rda.hp <- rdacca.hp(sampledata, env, method = 'RDA', type = 'R2', scale = FALSE)
mite.rda.hp
我们的计算结果如下,我们可以看见整体的解释率是86.9%,已经很高了。然后就是每个因子的具体解释率。看到这里大家就有疑问了,那有了具体的解释率只是定量分析,定性的显著性去哪里了呢?关于显著性的我们继续往后看
mite.rda.hp
$Method_Type
"RDA" "R2"
$Total_explained_variation
0.869
$Hier.part
Unique Average.share Individual I.perc(%)
TC 0.0220 0.0362 0.0582 6.70
TN 0.0234 0.0458 0.0692 7.96
CN 0.0231 0.0169 0.0400 4.60
TP 0.0338 0.0320 0.0658 7.57
NO3 0.0301 0.0308 0.0609 7.01
NH4 0.0321 0.0438 0.0759 8.73
AP 0.0106 0.0748 0.0854 9.83
AK 0.0314 0.0341 0.0655 7.54
pH 0.0264 0.0059 0.0323 3.72
Moisture 0.0272 0.0253 0.0525 6.04
BD 0.0139 0.0606 0.0745 8.57
shoot 0.0099 0.0199 0.0298 3.43
root 0.0245 0.0271 0.0516 5.94
ANPP 0.0661 0.0123 0.0784 9.02
0.0204 0.0087 0.0291 3.35
attr(,"class")
"rdaccahp"
层次分割分析的显著性检验
关于层次分割分析的显著性检验,rdacca.hp包自带了一个检验功能,但是该检验功能计算速度很慢,这个就属于耐心等一等也没法解决的范畴,所以,我一般不使用这个检验结果,我跑了两小时也没出来。这一步运行不出来没关系,继续往下运行就行。
#显著性检验(计算量很大,一般算出来要很久)
set.seed(123)
permu_hp <- permu.hp(dv = sampledata, iv = env, method = 'RDA', type = 'adjR2', permutations = 999)
permu_hp
因此,我使用的是另一个置换检验的结果。我们来看看我常用的置换检验,这个置换检验就是我说的由于不平衡设计导致改变因子顺序就会改变显著性的那个检验,那我为什么还要用?因为他计算起来快,此外,一般情况下置换检验的结果和我们计算的解释率是成正比的。所以,当我们给期刊提交数据和代码的时候,只要不改变因子的顺序保持和论文里一致的顺序,那么结果也是一致的,这也是被认可的。
#置换检验
envfit <- envfit(mite.rda, env, permutations = 999)
r <- as.matrix(envfit$vectors$r)
p <- as.matrix(envfit$vectors$pvals)
env.p <- cbind(r,p)
colnames(env.p) <- c("r2","p-value")
KK <- as.data.frame(env.p)
KK
置换检验的结果如下。这样我们不但得到了每个因子的具体解释率(定量),也得到了每个因子的显著性(定性),这个RDA分析才算是真正的圆满。
> KK
r2 p-value
TC 0.17276122 0.208
TN 0.29606161 0.039
CN 0.13342479 0.315
TP 0.38043888 0.011
NO3 0.44076671 0.014
NH4 0.48038258 0.006
AP 0.36832716 0.021
AK 0.36418834 0.029
pH 0.05936992 0.607
Moisture 0.11913897 0.323
BD 0.32645575 0.028
shoot 0.08417330 0.491
root 0.38656397 0.024
ANPP 0.24340040 0.074
plant.diversity 0.19299186 0.161
层次分割分析的绘图
分析完了我们就开始绘图,绘图和之前的一样,可以选择玫瑰图或者柱形图,我们先绘制玫瑰图,图的细节大家自己微调
#====================================================玫瑰图绘制==========================================================
# 加载必要的包
library(ggplot2)
library(dplyr)
# 1. 提取 mite.rda.hp 的 I.perc(%) 和指标名称
mite_rda_hp_data <- data.frame(
variable = rownames(mite.rda.hp$Hier.part),
percentage = mite.rda.hp$Hier.part[, "I.perc(%)"]
)
# 2. 提取 KK 中的 r2 和 p-value
KK_data <- data.frame(
variable = rownames(KK),
p_value = KK[, "p-value"]
)
# 3. 合并两个数据框
data_plot <- merge(mite_rda_hp_data, KK_data, by = "variable")
# 4. 设置显著性和颜色条件
data_plot$significance <- ifelse(data_plot$p_value < 0.05, "Significant", "Not Significant")
data_plot$color <- ifelse(data_plot$p_value < 0.05, "#F39B7FFF", "#8491B4FF")
# 5. 为了保持变量的原始顺序
data_plot$variable <- factor(data_plot$variable, levels = unique(data_plot$variable))
# 6. 绘制极坐标柱状图
ggplot(data_plot, aes(x = variable, y = percentage, fill = significance)) +
geom_bar(stat = "identity") +
geom_text(aes(label = round(percentage, 2)), vjust = -0.5, size = 3.5) + # 标注百分比
coord_polar(start = 1.5) + # 转换为极坐标
ylim(-3, max(data_plot$percentage) + 3) + # 自动设置Y轴范围
theme_minimal() +
labs(y = "% Explained Variation", x = "", fill = "Significance") +
scale_fill_manual(values = c("Significant" = "#F39B7FFF", "Not Significant" = "#8491B4FF")) +
theme(
axis.text.x = element_text(size = 10, hjust = 1), # x轴标签设置
panel.background = element_rect(fill = "white")
)
然后,再绘制柱形图看看,还是一样,细节自己微调即可
#=====================================================柱形图绘制==================================================================
# 1. 提取 mite.rda.hp 的 I.perc(%) 和指标名称
mite_rda_hp_data <- data.frame(
variable = rownames(mite.rda.hp$Hier.part),
percentage = mite.rda.hp$Hier.part[, "I.perc(%)"]
)
# 2. 提取 KK 中的 p-value
KK_data <- data.frame(
variable = rownames(KK),
p_value = KK[, "p-value"]
)
# 3. 合并两个数据框
data_plot <- merge(mite_rda_hp_data, KK_data, by = "variable")
# 4. 设置显著性和颜色条件
data_plot$significance <- ifelse(data_plot$p_value < 0.05, "Significant", "Not Significant")
data_plot$color <- ifelse(data_plot$p_value < 0.05, "#F39B7FFF", "#8491B4FF")
# 5. 按 percentage 降序排列
data_plot_sorted <- data_plot[order(-data_plot$percentage), ]
data_plot_sorted$variable <- factor(data_plot_sorted$variable, levels = data_plot_sorted$variable)
# 6. 绘制柱形图
ggplot(data_plot_sorted, aes(x = variable, y = percentage, fill = significance)) +
geom_bar(stat = "identity", width = 0.7) + # 调整柱形宽度
geom_text(aes(label = round(percentage, 2)), vjust = -0.5, size = 3.5) + # 在柱形上方标记数据
theme_minimal() +
ylim(0, max(data_plot_sorted$percentage) + 2) + # 自动设置 y 轴范围,增加上边距
labs(y = "% Explained Variation", x = "", fill = "Significance") + # 添加图例标题
scale_fill_manual(values = c("Significant" = "#F39B7FFF", "Not Significant" = "#8491B4FF")) + # 自定义颜色
theme(
axis.text.x = element_text(size = 10, angle = 45, hjust = 1), # x 轴标签旋转
panel.background = element_rect(fill = "white") # 白色背景
)
最后,大家自行选择使用玫瑰图还是柱形图,然后与RDA或CCA的图组合一下即可。相信这样的分析一定会让审稿人眼前一亮。
完整版代码如下
#调用R包
library(vegan)
library(ggrepel)
library(ggplot2)
library(ggpubr)
library(rdacca.hp)
#读取数据,依次为otu数据、环境因子数据、分组信息
sampledata <- read.csv(file.choose(), head = TRUE, row.names=1)
env <- read.csv(file.choose(), header=TRUE, row.names=1)
sampledata <- t(sampledata)
#对OTU数据进行hellinger转化
sampledata <- decostand(sampledata,method = "hellinger")
#简单的RDA
mite.rda <- rda(sampledata~., env, scale = FALSE)
summary(mite.rda)
plot(mite.rda)
#层次分割(这一步计算也很慢)
mite.rda.hp <- rdacca.hp(sampledata, env, method = 'RDA', type = 'R2', scale = FALSE)
mite.rda.hp
#显著性检验(计算量很大,一般算出来要很久)
set.seed(123)
permu_hp <- permu.hp(dv = sampledata, iv = env, method = 'RDA', type = 'adjR2', permutations = 999)
permu_hp
#置换检验
envfit <- envfit(mite.rda, env, permutations = 999)
r <- as.matrix(envfit$vectors$r)
p <- as.matrix(envfit$vectors$pvals)
env.p <- cbind(r,p)
colnames(env.p) <- c("r2","p-value")
KK <- as.data.frame(env.p)
KK
#====================================================玫瑰图绘制==========================================================
# 加载必要的包
library(ggplot2)
library(dplyr)
# 1. 提取 mite.rda.hp 的 I.perc(%) 和指标名称
mite_rda_hp_data <- data.frame(
variable = rownames(mite.rda.hp$Hier.part),
percentage = mite.rda.hp$Hier.part[, "I.perc(%)"]
)
# 2. 提取 KK 中的 r2 和 p-value
KK_data <- data.frame(
variable = rownames(KK),
p_value = KK[, "p-value"]
)
# 3. 合并两个数据框
data_plot <- merge(mite_rda_hp_data, KK_data, by = "variable")
# 4. 设置显著性和颜色条件
data_plot$significance <- ifelse(data_plot$p_value < 0.05, "Significant", "Not Significant")
data_plot$color <- ifelse(data_plot$p_value < 0.05, "#F39B7FFF", "#8491B4FF")
# 5. 为了保持变量的原始顺序
data_plot$variable <- factor(data_plot$variable, levels = unique(data_plot$variable))
# 6. 绘制极坐标柱状图
ggplot(data_plot, aes(x = variable, y = percentage, fill = significance)) +
geom_bar(stat = "identity") +
geom_text(aes(label = round(percentage, 2)), vjust = -0.5, size = 3.5) + # 标注百分比
coord_polar(start = 1.5) + # 转换为极坐标
ylim(-3, max(data_plot$percentage) + 3) + # 自动设置Y轴范围
theme_minimal() +
labs(y = "% Explained Variation", x = "", fill = "Significance") +
scale_fill_manual(values = c("Significant" = "#F39B7FFF", "Not Significant" = "#8491B4FF")) +
theme(
axis.text.x = element_text(size = 10, hjust = 1), # x轴标签设置
panel.background = element_rect(fill = "white")
)
#=====================================================柱形图绘制==================================================================
# 1. 提取 mite.rda.hp 的 I.perc(%) 和指标名称
mite_rda_hp_data <- data.frame(
variable = rownames(mite.rda.hp$Hier.part),
percentage = mite.rda.hp$Hier.part[, "I.perc(%)"]
)
# 2. 提取 KK 中的 p-value
KK_data <- data.frame(
variable = rownames(KK),
p_value = KK[, "p-value"]
)
# 3. 合并两个数据框
data_plot <- merge(mite_rda_hp_data, KK_data, by = "variable")
# 4. 设置显著性和颜色条件
data_plot$significance <- ifelse(data_plot$p_value < 0.05, "Significant", "Not Significant")
data_plot$color <- ifelse(data_plot$p_value < 0.05, "#F39B7FFF", "#8491B4FF")
# 5. 按 percentage 降序排列
data_plot_sorted <- data_plot[order(-data_plot$percentage), ]
data_plot_sorted$variable <- factor(data_plot_sorted$variable, levels = data_plot_sorted$variable)
# 6. 绘制柱形图
ggplot(data_plot_sorted, aes(x = variable, y = percentage, fill = significance)) +
geom_bar(stat = "identity", width = 0.7) + # 调整柱形宽度
geom_text(aes(label = round(percentage, 2)), vjust = -0.5, size = 3.5) + # 在柱形上方标记数据
theme_minimal() +
ylim(0, max(data_plot_sorted$percentage) + 2) + # 自动设置 y 轴范围,增加上边距
labs(y = "% Explained Variation", x = "", fill = "Significance") + # 添加图例标题
scale_fill_manual(values = c("Significant" = "#F39B7FFF", "Not Significant" = "#8491B4FF")) + # 自定义颜色
theme(
axis.text.x = element_text(size = 10, angle = 45, hjust = 1), # x 轴标签旋转
panel.background = element_rect(fill = "white") # 白色背景
)
脚本与示例数据下载
下载的时候顺便点个赞呗~
关注公众号“生态R学社”回复“层次分割4”即可下载完整代码与示例数据