R语言循环遍历多列生成带标准误差棒的图表实现问题

需求 #

对数据集指定列分别按年份、季节分组计算均值与标准误，生成带误差棒的散点图，需保留列名作为图表标题、列名作为文件命名前缀的原有规则。

测试数据集 #

data <- structure(list(year = c(2019L, 2019L, 2019L, 2019L, 2019L, 2019L, 
                        2019L, 2019L, 2019L, 2019L, 2020L, 2020L, 2020L, 2020L, 2020L, 
                        2020L, 2020L, 2020L, 2020L, 2020L), season = structure(c(1L, 
                                                                                 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
                                                                                 2L, 2L, 2L), .Label = c("dry", "wet"), class = "factor"), site = c(1L, 
                                                                                                                                                    2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 
                                                                                                                                                    3L, 4L, 5L), temp = c(26.9, 27.8, 27.3, 26.4, 26.8, 29.4, 29.9, 
                                                                                                                                                                          29.9, 29.9, 30, 24, 23.5, 23.9, 23.7, 24.1, 30.2, 30.8, 30.8, 
                                                                                                                                                                          30.3, 30.1), do = c(7.5, 7.2, 7.8, 8.1, 7.4, 3.7, 3.2, 6.4, 4.9, 
                                                                                                                                                                                              5, 6.5, 5.2, 5.9, 6.7, 6.3, 3.63, 1.81, 1.85, 4.25, 0.69), salinity = c(32.29, 
                                                                                                                                                                                                                                                                      30.1, 31.35, 31.93, 30.77, 27.35, 27.34, 28.42, 28.37, 28.24, 
                                                                                                                                                                                                                                                                      32.69, 29.72, 29.15, 28.9, 25.29, 24.37, 23.47, 25.1, 24.79, 
                                                                                                                                                                                                                                                                      23.62), pH = c(8.24, 8.28, 8.37, 8.32, 8.39, 7.85, 7.84, 8.13, 
                                                                                                                                                                                                                                                                                     8.04, 8.06, 8.26, 8.17, 8.18, 8.24, 8.13, 7.8, 7.61, 7.61, 7.95, 
                                                                                                                                                                                                                                                                                     7.41), water_depth = c(95L, 95L, 62L, 63L, 55L, 100L, 107L, 110L, 
                                                                                                                                                                                                                                                                                                            140L, 95L, 85L, 80L, 60L, 53L, 55L, 125L, 135L, 145L, 125L, 100L
                                                                                                                                                                                                                                                                                     ), sed_depth = c(56L, 40L, 22L, 1L, 20L, 60L, 47L, 68L, 40L, 
                                                                                                                                                                                                                                                                                                      20L, 55L, 35L, 20L, 1L, 25L, 45L, 30L, 35L, 1L, 20L), SAV = c(25.5, 
                                                                                                                                                                                                                                                                                                                                                                    41.5, 50, 47.5, 60.1, 53.5, 46.5, 80.5, 20, 32.5, 26.1, 39.5, 
                                                                                                                                                                                                                                                                                                                                                                    29.1, 48.5, 39.6, 63, 50.5, 70, 70, 56)), row.names = c(1L, 155L, 
                                                                                                                                                                                                                                                                                                                                                                                                                            309L, 463L, 617L, 771L, 925L, 1079L, 1233L, 1387L, 1541L, 1695L, 
                                                                                                                                                                                                                                                                                                                                                                                                                            1849L, 2003L, 2157L, 2311L, 2465L, 2619L, 2773L, 2927L), class = "data.frame")

原有问题代码 #

col_names <- colnames(data[,-c(1:3)])
for (i in col_names) {       # for-loop over columns
  cdata2 <- plyr::ddply(i, c("year", "season"), summarise,
                        N    = length(i),
                        n_mean = mean(i),
                        n_median = median(i),
                        sd   = sd(num),
                        se   = sd / sqrt(N))
  
  ggplot(cdata2, aes(x = year, y = n_mean, color = season)) +
    geom_errorbar(aes(ymin=n_mean-se, ymax=n_mean+se), 
                  width=.2, 
                  color = "black") +
    geom_point(color = "black", # Make both seasons have black borders
               shape = 21, 
               size = 3,
               aes(fill = season)) +
    scale_fill_manual(values=c("white", "#C0C0C0")) + 
    scale_x_continuous(breaks=c(2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2018,2019,2020)) +
    labs(x= NULL, y = "Mean count") +
    ggtitle(i)
  
  setwd('D:/.../Trend plots')
  
  ggsave(paste0(i, "- ENV_trend_plot.png"), 
         height = 5, width=7, units = "in")
}

修正后可运行代码 #

library(plyr)
library(ggplot2)

col_names <- colnames(data[,-c(1:3)])
for (i in col_names) {
  # 按年份、季节分组计算统计量
  cdata2 <- plyr::ddply(data, c("year", "season"), summarise,
                        N    = length(.data[[i]]),
                        n_mean = mean(.data[[i]]),
                        n_median = median(.data[[i]]),
                        sd   = sd(.data[[i]]),
                        se   = sd / sqrt(N))
  
  # 生成绘图对象
  p <- ggplot(cdata2, aes(x = year, y = n_mean, color = season)) +
    geom_errorbar(aes(ymin=n_mean-se, ymax=n_mean+se), 
                  width=.2, 
                  color = "black") +
    geom_point(color = "black",
               shape = 21, 
               size = 3,
               aes(fill = season)) +
    scale_fill_manual(values=c("white", "#C0C0C0")) + 
    # 如需保留原全量年份刻度，可替换为你原有scale_x_continuous配置
    scale_x_continuous(breaks=c(2019,2020)) +
    labs(x= NULL, y = "Mean count") +
    # 保留原标题规则
    ggtitle(i)
  
  # 输出文件，路径替换为本地实际路径即可
  # setwd('D:/.../Trend plots')
  ggsave(paste0(i, "- ENV_trend_plot.png"), plot = p,
         height = 5, width=7, units = "in")
}

核心修改点 #

• 修正ddply入参：第一个参数传入完整数据集，通过.data[[i]]读取当前循环的目标列数值，解决原代码传入列名字符串无法计算的问题
• 删除原代码中不存在的变量num引用，改用目标列数值计算标准差
• 显式生成ggplot绘图对象并传入ggsave，避免输出空白文件
• 完全保留原有的ggtitle(i)标题设置和文件命名规则

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