如何从R语言的rda对象提取PC解释率并添加至图表？

你可以通过以下步骤从vegan的rda对象中提取PC轴的解释百分比，并将其添加到图表的坐标轴标签中：

1. 提取特征值并计算解释百分比 #

因为你是对标准化后的数据做无约束RDA（本质等同于PCA），相关信息存储在all.rda$CA部分：

方法一：直接从特征值计算 #

# 提取所有主成分的特征值
eig <- all.rda$CA$eig
# 计算每个轴的解释方差百分比，保留两位小数
pc_percent <- round(eig / sum(eig) * 100, 2)
# 生成PC1和PC2的标签文本
pc1_label <- paste0("PC1 (", pc_percent[1], "%)")
pc2_label <- paste0("PC2 (", pc_percent[2], "%)")

方法二：使用summary()函数提取 #

summary()函数会返回RDA分析的详细统计结果，其中包含方差解释率：

# 生成RDA分析的摘要
rda_summary <- summary(all.rda)
# 提取"Proportion of Variance"行（即解释比例），转换为百分比并保留两位小数
pc_percent <- round(rda_summary$cont$importance[2, ] * 100, 2)
# 生成坐标轴标签
pc1_label <- paste0("PC1 (", pc_percent[1], "%)")
pc2_label <- paste0("PC2 (", pc_percent[2], "%)")

2. 将百分比添加到图表中 #

下面提供Base R和ggplot2两种绘图方式的完整示例：

示例1：ggplot2绘图 #

library(vegan)
library(tidyverse)

# 示例数据（你的原始数据）
all.data <- data.frame(
  Al.308.215 = c(1933.110774, 2124.933186, 2697.142811, 2454.099034, 2416.374341, 2730.405337,2846.957744, 2957.419342, 2959.348363, 2853.258817),
  Al.396.152 = c(1955.47859, 2142.707892, 2719.348431, 2475.607238, 2428.867771, 2759.503064,2875.096669, 2990.04952, 2988.810866, 2880.507931),
  Ba.233.527 = c(4.758488077, 4.873488218, 5.496073131, 5.322604493, 5.334528828, 5.519915915,5.715194792, 5.774621594, 7.152332657, 7.157559445),
  Ba.455.404 = c(4.676267981, 4.748689692, 5.424586334, 5.246081426, 5.233062678, 5.41001526,5.639145829, 5.687752522, 7.029890126, 7.058378483),
  Ca.183.801 = c(659.0148918, 655.0216397, 705.7125714, 699.3610886, 685.2562524, 682.351838,667.2459663, 666.5119829, 451.5861085, 430.0641182),
  Ca.317.933 = c(663.2478342, 660.1603433, 706.2024694, 704.5687867, 685.9959824, 687.0625639,670.5308695, 667.5984284, 448.1139829, 428.0179807),
  Fe.259.941 = c(902.8585447, 911.2101809, 1017.666402, 1022.594904, 1048.51115, 1002.091143,970.8258109, 905.2628883, 997.8037721, 1046.38376),
  K.766.491 = c(314.1726307, 309.757608, 389.1136882, 376.3265692, 373.4714432, 413.0336628,420.9373827, 437.5226928, 475.3400126, 477.8208541),
  Mg.285.213 = c(510.5124552, 596.0268573, 829.8716737, 768.2083155, 744.1670151, 784.4661692,796.7464128, 837.234733, 684.5431608, 707.3844484),
  Mn.257.611 = c(15.91777975, 15.88328987, 17.680152, 17.57143193, 18.39642397, 17.96018218,17.69058768, 17.24667962, 19.92596564, 20.23646264),
  row.names = c("2A", "4A", "5A", "6A", "7A", "8A", "9A", "10A", "11A", "12A")
)

# 模拟metadata（替换为你的真实metadata）
icp.metadata <- data.frame(
  sample.id = rownames(all.data),
  group = rep(c("Group1", "Group2"), each = 5)
)

# 数据标准化与RDA分析
all.std <- na.omit(decostand(all.data, method = "normalize"))
all.rda <- rda(all.std)

# 提取百分比标签
eig <- all.rda$CA$eig
pc_percent <- round(eig / sum(eig) * 100, 2)
pc1_label <- paste0("PC1 (", pc_percent[1], "%)")
pc2_label <- paste0("PC2 (", pc_percent[2], "%)")

# 提取绘图数据
uscores.all <- data.frame(all.rda$CA$u) %>%
  rownames_to_column("sample.id")
uscores1.all <- inner_join(icp.metadata, uscores.all, by = "sample.id")
vscores.all <- data.frame(all.rda$CA$v) %>%
  rownames_to_column("variable")

# 绘图
ggplot() +
  geom_point(data = uscores1.all, aes(x = PC1, y = PC2, color = group), size = 3) +
  geom_segment(data = vscores.all, aes(x = 0, y = 0, xend = PC1, yend = PC2),
               arrow = arrow(length = unit(0.2, "cm")), color = "darkred") +
  geom_text(data = vscores.all, aes(x = PC1 * 1.1, y = PC2 * 1.1, label = variable),
            color = "darkred", size = 3.5) +
  labs(x = pc1_label, y = pc2_label) +
  theme_bw()

示例2：Base R绘图 #

# 承接上述数据处理步骤，直接绘图
plot(all.rda$CA$u, main = "PCA 可视化", xlab = pc1_label, ylab = pc2_label,
     pch = 16, col = uscores1.all$group)
# 添加变量载荷向量
arrows(0, 0, all.rda$CA$v[,1]*0.8, all.rda$CA$v[,2]*0.8, col = "darkred", length = 0.1)
# 添加变量名称
text(all.rda$CA$v[,1]*1.1, all.rda$CA$v[,2]*1.1, rownames(all.rda$CA$v), col = "darkred")

内容的提问来源于stack exchange，提问作者Geomicro