如何在ggplot2中为嵌套分面图添加误差棒？

我尝试为以下绘图添加误差棒，调整了width和position_dodge参数但都没效果，不确定是不是和mean_cl_boot有关。代码如下：

library(ggplot2)
library(dplyr)
library(tidyr)
library(ggh4x)

df <- expand.grid(  
  group = c("Children", "Teenagers"),  
  task = c("A", "B", "C"),  
  hemisphere = c("Left", "Right"),  
  region = c("Frontal", "Motor"),  
  frequency = c("Alpha", "Beta", "Theta")  
)  

df$power <- runif(nrow(df))  # 生成演示用的随机power值  

avg_power <- df %>%  
  group_by(task, group, frequency) %>%  
  summarize(mean_power_left = mean(power[hemisphere == "Left"]),   
            mean_power_right = mean(power[hemisphere == "Right"]),   
            mean_power_frontal = mean(power[region == "Frontal"]),   
            mean_power_motor = mean(power[region == "Motor"]),   
            .groups = 'drop')  

avg_power_long <- avg_power %>%  
  pivot_longer(cols = starts_with("mean_power"),  
               names_to = "area",   
               names_pattern = "mean_power_(.*)") %>%  
  mutate(area = factor(area, levels = c("left", "right", "frontal", "motor"), 
                       labels = c("Left", "Right", "Frontal", "Motor"))) %>%
  arrange(task, group, frequency, area) %>%
  select(group, task, frequency, area, value) %>%
  rename(power = value)

avg_power_long <- avg_power_long %>%
  mutate(region = case_when(area == c("Left", "Right") ~ "Hemisphere",
                            area == c("Frontal", "Motor") ~ "Area"),
         .before = power)

interaction_plot <- avg_power_long %>% 
  mutate(frequency = factor(frequency, levels = c("Theta", "Alpha", "Beta"))) %>%
  ggplot(aes(x = task, y = power, color = group)) +
  stat_summary(fun = mean, geom = "line", aes(group = group), linewidth = 1.5) + 
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.3, 
               position = position_dodge(0.4)) +
  facet_nested(frequency ~ region + area) +
  labs(title = "",
       x = "Task", y = "Mean EEG Band Power",
       color = "Group") +
  theme_minimal(base_size = 15) +
  theme(strip.text.x = element_text(size = 14),
        strip.placement = "outside")

print(interaction_plot)

绘图效果：

你提前把数据汇总成了单个均值，mean_cl_boot需要基于重复的原始样本计算置信区间，但avg_power_long里每个task-group-frequency-area组合只有1个power值，没有重复数据，所以生成不了误差棒。另外，stat_summary的position_dodge需要和x轴分组的偏移匹配，但这不是核心问题。

提供两种可行方案：

方案1：基于原始数据直接统计（推荐） #

跳过提前汇总的步骤，让stat_summary直接对原始数据计算均值和置信区间，误差棒会自动生成：

library(ggplot2)
library(dplyr)
library(tidyr)
library(ggh4x)

df <- expand.grid(  
  group = c("Children", "Teenagers"),  
  task = c("A", "B", "C"),  
  hemisphere = c("Left", "Right"),  
  region = c("Frontal", "Motor"),  
  frequency = c("Alpha", "Beta", "Theta")  
)  

df$power <- runif(nrow(df))  

# 整理原始数据，拆分出area和region字段
df_long <- df %>%
  pivot_longer(cols = c(hemisphere, region), names_to = "region_type", values_to = "area") %>%
  mutate(region = case_when(
    region_type == "hemisphere" ~ "Hemisphere",
    region_type == "region" ~ "Area"
  )) %>%
  select(-region_type)

# 基于原始数据绘图
interaction_plot <- df_long %>% 
  mutate(frequency = factor(frequency, levels = c("Theta", "Alpha", "Beta"))) %>%
  ggplot(aes(x = task, y = power, color = group)) +
  # 均值线：自动计算分组均值，position_dodge和误差棒保持一致
  stat_summary(fun = mean, geom = "line", aes(group = group), linewidth = 1.5, 
               position = position_dodge(0.4)) + 
  # 误差棒：用mean_cl_boot计算置信区间
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.3, 
               position = position_dodge(0.4)) +
  facet_nested(frequency ~ region + area) +
  labs(title = "",
       x = "Task", y = "Mean EEG Band Power",
       color = "Group") +
  theme_minimal(base_size = 15) +
  theme(strip.text.x = element_text(size = 14),
        strip.placement = "outside")

print(interaction_plot)

方案2：汇总时手动计算置信区间 #

如果必须用提前汇总的数据，需要在汇总阶段就计算好每个分组的置信区间，再用geom_errorbar绘制：

library(ggplot2)
library(dplyr)
library(tidyr)
library(ggh4x)
library(boot) # 用于bootstrap置信区间计算

df <- expand.grid(  
  group = c("Children", "Teenagers"),  
  task = c("A", "B", "C"),  
  hemisphere = c("Left", "Right"),  
  region = c("Frontal", "Motor"),  
  frequency = c("Alpha", "Beta", "Theta")  
)  

df$power <- runif(nrow(df))  

# 自定义函数：计算均值的bootstrap置信区间
boot_ci <- function(x) {
  boot_mean <- function(data, i) mean(data[i])
  boot_result <- boot(data = x, statistic = boot_mean, R = 1000)
  ci <- boot.ci(boot_result, type = "perc")$percent[4:5]
  tibble(mean = mean(x), lower = ci[1], upper = ci[2])
}

# 按分组计算均值和置信区间
avg_power_with_ci <- df %>%
  pivot_longer(cols = c(hemisphere, region), names_to = "region_type", values_to = "area") %>%
  mutate(region = case_when(
    region_type == "hemisphere" ~ "Hemisphere",
    region_type == "region" ~ "Area"
  )) %>%
  select(-region_type) %>%
  group_by(task, group, frequency, region, area) %>%
  summarise(boot_ci(power), .groups = "drop")

# 绘图
interaction_plot <- avg_power_with_ci %>% 
  mutate(frequency = factor(frequency, levels = c("Theta", "Alpha", "Beta"))) %>%
  ggplot(aes(x = task, y = mean, color = group)) +
  geom_line(aes(group = group), linewidth = 1.5, position = position_dodge(0.4)) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.3, position = position_dodge(0.4)) +
  facet_nested(frequency ~ region + area) +
  labs(title = "",
       x = "Task", y = "Mean EEG Band Power",
       color = "Group") +
  theme_minimal(base_size = 15) +
  theme(strip.text.x = element_text(size = 14),
        strip.placement = "outside")

print(interaction_plot)

内容的提问来源于stack exchange，提问作者always.learning