如何用ggplot绘制agricolae包Tukey检验结果，X轴按年份倒序并显示标准误？

我正在分析2000-2020年样带的植物与土壤数据，用agricolae包完成了单因素ANOVA及Tukey检验。目前用base plot绘制探索性图形，但样式调整困难，需要用ggplot绘制高质量图，要求：

• X轴按2020至2000的时间倒序排列（而非默认的均值降序）
• 误差线显示**标准误（SE）**而非标准差（SD）

数据集（biochemA） #

Year Transect Depth NH4_(ug/L)  
2000 1 1 391,777  
2009 1 1 471,616 
2011 1 1 362,964
2013 1 1 348,208 
2015 1 1 234,341 
2017 1 1 768,844
2019 1 1 599,063 
2020 1 1 912,435
2000 2 1 452,272 
2009 2 1 391,134
2011 2 1 285,286 
2013 2 1 755,01 
2015 2 1 376,022 
2017 2 1 1205,095
2019 2 1 2940,163 
2020 2 1 298,487 
2000 3 1 1409,322 
2009 3 1 847,658 
2011 3 1 332,635 
2013 3 1 487,695 
2015 3 1 337,721 
2017 3 1 1702,21 
2019 3 1 2684,409 
2020 3 1 448,644

现有ANOVA及Tukey检验代码 #

modelNH4 <- aov(biochemA$NH4_(ug/L) ~ Year, data = biochemA)
out_NH4 <- HSD.test(modelNH4, "Year", group=TRUE, console=TRUE)
plot(out_NH4, main="NH4", xlab="Year", ylab="μg/g soil", las=2)

Tukey检验输出结果（中文翻译） #

> out_bioANH4
$statistics
   均方误差 自由度  总体均值   变异系数 最小显著差
  173.8115     38 11.83804 111.3678 26.72761

$parameters
   检验方法 分组变量 组数 学生化极差 显著性水平
  Tukey     Year      8     4.53321      0.05

$means
     NH4_(μg/g土壤)       标准差 重复数       最小值       最大值
2000        7.662795  6.182201     5  3.000779 18.067119
2009        6.314628  3.589682     5  2.905835 12.042757
2011        4.866242  1.180330     5  3.698008  6.192181
2013        6.401148  3.330437     5  3.628023 10.752502
2015        4.286512  1.424005     5  2.777662  6.101243
2017       14.385232  8.981692     5  6.405698 29.624701
2019       42.317093 18.145402     5 12.244628 57.310143
2020        8.470637  4.223832     5  4.573570 15.409000
           下四分位数 中位数 上四分位数
2000  3.547272  5.327341  8.371463
2009  3.728181  6.008561  6.887809
2011  4.097738  4.257121  6.086162
2013  3.684275  4.766714  9.174225
2015  3.089635  4.135425  5.328595
2017  9.569202 12.620927 13.705629
2019 38.758582 50.761796 52.510317
2020  5.761486  7.738720  8.870408

$comparison
NULL

$groups
     NH4_(μg/g土壤) 分组
2019       42.317093    a
2017       14.385232    b
2020        8.470637    b
2000        7.662795    b
2013        6.401148    b
2009        6.314628    b
2011        4.866242    b
2015        4.286512    b

attr("class")
[1] "group"

步骤1：数据预处理 #

首先处理NH4_(ug/L)列的逗号分隔符，转为数值型；同时将Year设为因子类型，方便后续排序。

library(tidyverse)
library(agricolae)

# 数据预处理
biochemA <- biochemA %>%
  mutate(NH4 = as.numeric(str_replace(`NH4_(ug/L)`, ",", ""))) %>%  # 替换逗号并转数值
  mutate(Year = as.factor(Year))

步骤2：计算统计量并合并Tukey分组 #

从out_NH4中提取分组信息，同时按年份计算均值和标准误：

# 提取Tukey分组信息
tukey_groups <- out_NH4$groups %>%
  rownames_to_column("Year") %>%
  rename(NH4_mean = `biochemA1$NH4_.ug.g_soil._`, group = groups) %>%
  mutate(Year = as.factor(Year))

# 计算均值和标准误
summary_stats <- biochemA %>%
  group_by(Year) %>%
  summarise(
    NH4_mean = mean(NH4),
    se = sd(NH4)/sqrt(n())  # 标准误计算公式：标准差/根号样本量
  ) %>%
  left_join(tukey_groups, by = c("Year", "NH4_mean"))  # 合并分组信息

步骤3：ggplot绘图（满足X轴倒序+标准误要求） #

ggplot(summary_stats, aes(x = fct_rev(Year), y = NH4_mean)) +
  geom_col(fill = "#4285F4", width = 0.7) +  # 柱状图
  geom_errorbar(aes(ymin = NH4_mean - se, ymax = NH4_mean + se), 
                width = 0.2, color = "black") +  # 标准误误差线
  geom_text(aes(label = group), vjust = -0.5, size = 4, color = "red") +  # 添加Tukey分组标记
  labs(
    title = "NH4含量随年份变化（Tukey检验）",
    x = "年份",
    y = "NH4含量（μg/g土壤）"
  ) +
  theme_bw() +  # 简洁黑白主题
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )

关键说明 #

• fct_rev(Year)：将年份因子倒序，实现X轴从2020到2000的排列
• se = sd(NH4)/sqrt(n())：手动计算标准误，替代默认的标准差
• geom_text：将Tukey检验的分组标记添加到柱状图顶部，直观展示差异显著性

内容的提问来源于stack exchange，提问作者Christoffer Bugge Harder