R语言ggplot2堆叠柱状图添加误差棒、样本量标注及统计检验实现

以下代码可直接复现所有需求，统计方法适配重复测量流行病学随访数据的特征。 #

依赖包加载与数据预处理 #

先加载分析所需包，统一数据格式：

library(tidyverse)
library(binom)
library(lme4)
library(car)
library(emmeans)

# 输入示例数据
ID = c(001, 001, 001, 001, 001, 002, 002, 002, 003, 003, 004, 004, 004, 004, 005, 005, 006, 007, 007, 007, 008, 008, 008, 008, 008, 009, 009, 009, 009, 009)
Visit = c(00, 01, 02, 03, 04, 00, 01, 02,  01, 02, 00, 01, 02, 03, 00, 02, 00, 01, 02, 04, 00, 01, 02, 03, 04, 00, 01, 02, 03, 04)
CVD = c(0, 0, 1, 1, 1, 1, 1, 0, 1, 1,  1, 0, 0, 1, 1,  0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0)
TRT= c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2)
BIO=c(12.00, 11.9, 15.24, 13.10, 30.01, 45.90, 20.09, 23.45, 14.78, 18.05,
      24.23, 12.34, 80.01, 13.98, 12.50, 36.95, 29.00, 39.87, 19.03, 11.48,
      14.14, 28.06, 12.22, 72.08, 15.00, 11.33, 58.00, 17.71, 52.08, 15.25)
df<-data.frame(ID,Visit, CVD, TRT, BIO) %>%
  mutate(
    Visit = factor(Visit),
    TRT = factor(TRT, levels = c(1,2), labels = c("治疗组1", "治疗组2")),
    CVD = factor(CVD, levels = c(0,1), labels = c("未患病", "患病"))
  )

核心统计量计算 #

同步计算每个访视-治疗组的样本量、CVD患病人数、患病率点估计、95%置信区间（采用Wilson区间法，小样本表现优于正态近似），同时生成堆叠图所需的百分比标签：

# 患病率与置信区间统计量
sum_stats <- df %>%
  group_by(TRT, Visit) %>%
  summarise(N = n(), CVD_case = sum(CVD == "患病"), .groups = "drop") %>%
  bind_cols(
    binom.confint(x = .$CVD_case, n = .$N, conf.level = 0.95, methods = "wilson") %>% 
      select(mean, lower, upper)
  ) %>%
  rename(prevalence = mean)

# 堆叠柱百分比标签
percentData <- df %>% 
  group_by(Visit, TRT) %>% 
  count(CVD) %>%
  mutate(ratio=scales::percent(n/sum(n), accuracy = 0.1)) %>% 
  ungroup()

可视化绘制 #

满足横轴标注样本量、添加患病率点估计与95%CI误差棒的要求：

ggplot(sum_stats, aes(x = Visit)) +
  geom_bar(data = df, aes(fill = CVD), position = "fill", alpha = 0.7) +
  geom_text(
    data = percentData,
    aes(y = after_stat(count)/tapply(after_stat(count), after_stat(x), sum)[after_stat(x)], 
        label = ratio, group = CVD),
    position = position_fill(vjust = 0.5), size = 3.5
  ) +
  # 患病率点与误差棒
  geom_point(aes(y = prevalence), color = "#c0392b", size = 2.5, position = position_dodge(0.5)) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, color = "#c0392b", 
                linewidth = 0.8, position = position_dodge(0.5)) +
  # 横轴下方标注样本量
  geom_text(aes(y = -0.08, label = paste0("N=", N)), size = 3, hjust = 0.5) +
  facet_wrap(~TRT) +
  scale_y_continuous(
    name = "CVD患病率", labels = scales::percent,
    expand = expansion(mult = c(0.1, 0.05))
  ) +
  scale_fill_manual(values = c("#bdc3c7", "#3498db"), name = "CVD状态") +
  theme_bw() +
  theme(strip.background = element_rect(fill = "#f0f0f0"), panel.grid.minor = element_blank())

置信区间结果表导出 #

直接导出所有分组的统计结果为独立表格：

write.csv(
  sum_stats %>% mutate(across(c(prevalence, lower, upper), ~scales::percent(., accuracy = 0.1))),
  "CVD患病率_95%CI结果表.csv", row.names = F, fileEncoding = "UTF-8"
)

统计检验方法选择 #

由于数据为同一研究对象的重复随访测量，存在个体内相关性，不可使用普通卡方检验，推荐采用二项分布广义线性混合效应模型（GLMM） 完成两类检验，将个体ID作为随机截距控制非独立误差：

a) 患病率随时间变化的显著性检验 #

• 分别拟合包含Visit固定效应的全模型、不包含Visit的空模型，通过似然比检验比较两个模型的拟合优度，得到时间因素对患病率影响的整体p值
• 若需比较各随访时间点与基线的患病率差异，采用emmeans进行事后两两比较，用Bonferroni法校正多重检验p值
  示例代码：

model_full <- glmer(CVD == "患病" ~ Visit + TRT + (1|ID), data = df, family = binomial)
model_null <- glmer(CVD == "患病" ~ TRT + (1|ID), data = df, family = binomial)
time_effect_p <- anova(model_full, model_null, test = "LRT")$`Pr(>Chisq)`[2]
time_pairwise <- emmeans(model_full, pairwise ~ Visit, adjust = "bonferroni")$contrasts

b) 治疗组间患病率差异的显著性检验 #

• 拟合包含Visit*TRT交互项的模型，通过III型方差分析得到治疗主效应、时间-治疗交互效应的p值：交互效应显著说明两组患病率随时间的变化趋势存在差异，治疗主效应显著说明两组整体患病率存在差异
• 若需比较每个访视下的组间患病率差异，采用emmeans按访视分层做组间比较，校正多重检验p值
  示例代码：

model_inter <- glmer(CVD == "患病" ~ Visit * TRT + (1|ID), data = df, family = binomial)
model_anova_res <- Anova(model_inter, type = 3)
trt_main_p <- model_anova_res["TRT", "Pr(>Chisq)"]
inter_p <- model_anova_res["Visit:TRT", "Pr(>Chisq)"]
trt_compare_per_visit <- emmeans(model_inter, pairwise ~ TRT | Visit, adjust = "bonferroni")$contrasts

注意：若样本量足够、个体随时间的变化趋势存在明显异质性，可将随机效应升级为随机斜率形式(1 + Visit | ID)，模型拟合效果更优；若存在失访，建议先做多重插补再建模，避免选择偏倚。

内容的提问来源于stack exchange，提问作者R. Simian