欢迎关注R语言数据分析指南
❝本节继续来介绍之前的统计分析系列，此次介绍rda分析的过程及绘图，个人观点仅供参考。本数据+详细代码注释稍后会上传到交流群内，购买过小编绘图文档的朋友可在所加的交流群内获取下载，有需要的朋友可关注文末介绍购买小编的R绘图文档。购买前请咨询，零基础不要买。
结果图
加载R包
library(tidyverse)
library(vegan)
library(broom)
导入数据
spe <- read_tsv("sp.txt") %>% column_to_rownames(var="label")
env <- read_tsv("env.txt") %>% column_to_rownames(var="label")
group <- read_tsv("group.txt")
数据标准化
spe.hel <- decostand(spe, method = "hellinger")
env.s <- env %>% mutate(across(where(is.numeric),~
decostand(.,method = "standardize")))
❝标准化的情形
•例如：Temperature（℃）、pH（无单位）、Nutrient（mg/L） 具有不同的度量尺度。
•不标准化的话，大数值变量（如温度）可能主导分析结果，影响环境变量的贡献计算。
•若某些变量有很大的数值范围差异（如 0.01 vs. 100），标准化有助于均衡变量的影响。
•例如：想看 Temperature 和 pH 对物种丰度的相对重要性，则标准化后可以直接比较它们的回归系数。
RDA分析
spe.rda <- rda(spe.hel~.,env.s) # 对全部因子进行rda分析
#逐步前进变量选择（forward selection）
fwd.sel <- ordiR2step(rda(spe.hel ~ 1, data = env.s),
scope = formula(spe.rda),
direction = "forward", R2scope = TRUE,
pstep = 1000, trace = FALSE)
fwd.sel$call #查看最终选中的环境变量
rda_tb <- rda(formula = spe.hel ~ Temperature + DPP + TP, data = env.s)
# 计算 调整 R²（Adjusted R²），评估 环境变量对物种丰度变化的解释力。
RsquareAdj(rda_tb)$adj.r.squared
# 检验 RDA 模型是否显著。
anova.cca(rda_tb, step = 1000)
# 检验各个环境变量的贡献
df <- anova.cca(rda_tb, by = "terms", step = 1000)
# 整理结果
df_tidy <- tidy(df) %>% drop_na() %>%
mutate(
sig = case_when(
p.value <= 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*",
TRUE ~ ""))
结果整理
❝scaling = 2：指定 RDA 中 “双重缩放”（Scaling 2），这种缩放方式特别适用于物种数据。它的目的是同时考虑物种和环境因子的变异，并且使 物种点（site scores） 和 环境变量（environmental variables） 都尽可能地平衡。
rda_tb.scaling2 <- summary(rda_tb,scaling=2)
rda_tb.site <- data.frame(rda_tb.scaling2$sites)[1:2] %>%
rownames_to_column(var="sample") %>%
left_join(.,group,by="sample")
rda_tb.env <- data.frame(rda_tb.scaling2$biplot)[1:2] %>%
rownames_to_column(var="sample") %>%
left_join(.,df_tidy %>% dplyr::select(1,sig),by=c("sample"="term"))
rda_tb.species <- data.frame(rda_tb.scaling2$species)[1:2] %>%
rownames_to_column(var="sample")
计算贡献度
rda_explained <- eigenvals(rda_tb) / sum(eigenvals(rda_tb))
rda1_explained <- round(rda_explained[1] * 100, 2) # 转换为百分比
rda2_explained <- round(rda_explained[2] * 100, 2)
数据可视化
ggplot(rda_tb.site, aes(RDA1, RDA2,color=group)) +
geom_vline(xintercept = 0, color = 'black', size = 0.5) +
geom_hline(yintercept = 0, color = 'black', size = 0.5) +
geom_point(size=2) +
geom_segment(data=rda_tb.env,
aes(x = 0, y = 0, xend = RDA1*1.5,yend = RDA2*1.5),
arrow = arrow(length = unit(0.2, 'cm')),
size = 0.4, color = '#4876FF')+
geom_text(data = rda_tb.env %>% filter(sample !="Temperature"),
aes(RDA1 * 1.55, RDA2 * 1.55,label = sample),
color = '#4876FF',size = 4)+
geom_text(data = rda_tb.env %>% filter(sample =="Temperature"),
aes(RDA1 * 1.55, RDA2 * 1.55,label = sample),
color = '#4876FF',size = 4,vjust=1,hjust=1) +
geom_text(data = rda_tb.env,
aes(RDA1 * 1.55, RDA2 * 1.55,label = sig),
color = 'red',size = 4,vjust=-0.1)+
geom_segment(data=rda_tb.species,
aes(x = 0, y = 0, xend = RDA1,yend = RDA2),
arrow = arrow(length = unit(0.2, 'cm')),
size = 0.4, color = '#FDB462') +
geom_text(data = rda_tb.species,
aes(RDA1 * 1.1, RDA2 * 1.1,
label = sample), color = '#FDB462', size = 4)+
scale_color_manual(values=c('#E31A1C','#228B22','#1F78B4', '#FDB462',
'#8B658B','#4876FF')) +
labs(x = paste0("RDA1 (", rda1_explained, "%)"),
y = paste0("RDA2 (", rda2_explained, "%)")) +
theme_test()
关注下方公众号下回更新不迷路
购买介绍
❝本节介绍到此结束，有需要学习R数据可视化的朋友欢迎到淘宝店铺：R语言数据分析指南，购买小编的R语言可视化文档，2025年购买将获取2025年更新的内容，同时将赠送2024年的绘图文档内容。
更新的绘图内容包含数据+代码+注释文档+文档清单，小编只分享案例文档，不额外回答问题，无答疑服务，更新截止2025年12月31日结束，零基础不推荐买。
案例特点
❝所选案例图绝大部份属于个性化分析图表，数据案例多来自已经发表的高分论文，并会汇总整理分享一些论文中公开的分析代码。
2025年起提供更加专业的html文档，更加的直观易学。文档累计上千人次购买拥有良好的社群交流体验，R代码结构清晰易懂.
目录大纲展示
群友精彩评论
淘宝店铺
2025年更新案例图展示
2024年案例图展示
