基于RDA、CCA的层次分割分析及绘图

学术   2024-12-24 15:31   云南  

也写在前面

又是思路干涸的一天,写写推文换换脑子,今天就不说废话了,把废话都留给年终总结。我们之前陆续更新了(点击链接就可以跳转):

1. 基于线性回归的层次分割分析及绘图

2. 基于双因素的层次分割分析及绘图

3. 基于混合效应模型的层次分割分析及绘图

这一期,我们来更新一下基于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即可

#调用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")
#简单的RDAmite.rda <- rda(sampledata~., env, scale = FALSE)summary(mite.rda)plot(mite.rda)

运行完毕后,我们就会得到RDA的关键数据和一幅不好看的图,图的美化还是那句话,看之前的推文就可以

Biplot scores for constraining variables
RDA1 RDA2 RDA3 RDA4 RDA5 RDA6TC 0.342161 -0.263115 -0.388125 0.137307 0.251537 -0.02199TN 0.054603 -0.549167 -0.264299 0.343325 -0.316158 0.27765CN 0.253506 0.302074 -0.085245 -0.206632 0.538819 -0.26552TP -0.613183 -0.301599 0.005675 -0.005963 -0.169216 0.06092NO3 -0.548305 -0.475723 0.151527 0.271580 -0.204094 -0.15554NH4 -0.712661 -0.288958 0.349700 0.360944 -0.189774 -0.05820AP -0.666277 0.008007 -0.275832 0.407262 -0.241981 0.06987AK 0.579230 0.331089 -0.503094 -0.252480 0.234978 0.12582pH 0.069097 0.243470 0.111079 0.119347 0.071269 0.14927Moisture -0.257933 -0.271649 0.716006 -0.150313 -0.146433 0.01676BD -0.625180 -0.094909 -0.264975 0.428801 -0.295267 0.09396shoot 0.009318 -0.294962 -0.097992 -0.326072 0.249450 -0.39597root -0.425110 0.474072 -0.088597 0.005317 -0.151714 -0.16808ANPP 0.176580 0.484541 0.165789 0.270270 -0.407461 -0.02564plant.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[1] "RDA" "R2" 
$Total_explained_variation[1] 0.869
$Hier.part Unique Average.share Individual I.perc(%)TC 0.0220 0.0362 0.0582 6.70TN 0.0234 0.0458 0.0692 7.96CN 0.0231 0.0169 0.0400 4.60TP 0.0338 0.0320 0.0658 7.57NO3 0.0301 0.0308 0.0609 7.01NH4 0.0321 0.0438 0.0759 8.73AP 0.0106 0.0748 0.0854 9.83AK 0.0314 0.0341 0.0655 7.54pH 0.0264 0.0059 0.0323 3.72Moisture 0.0272 0.0253 0.0525 6.04BD 0.0139 0.0606 0.0745 8.57shoot 0.0099 0.0199 0.0298 3.43root 0.0245 0.0271 0.0516 5.94ANPP 0.0661 0.0123 0.0784 9.02plant.diversity 0.0204 0.0087 0.0291 3.35
attr(,"class")[1] "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-valueTC              0.17276122   0.208TN              0.29606161   0.039CN              0.13342479   0.315TP              0.38043888   0.011NO3             0.44076671   0.014NH4             0.48038258   0.006AP              0.36832716   0.021AK              0.36418834   0.029pH              0.05936992   0.607Moisture        0.11913897   0.323BD              0.32645575   0.028shoot           0.08417330   0.491root            0.38656397   0.024ANPP            0.24340040   0.074plant.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-valueKK_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-valueKK_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")
#简单的RDAmite.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-valueKK_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-valueKK_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”即可下载完整代码与示例数据




走天涯徐小洋地理数据科学
一个爱生活的地理土博,分享GIS、遥感、空间分析、R语言、景观生态等地理数据科学实操教程、经典文献、数据资源
 最新文章