嘿，针对你的问题，我来拆解一下生存分析的适用性，再给出具体的R代码方案～

首先明确：累计发病率（IR）本质是固定随访期内新发病例占暴露人群的比例，而生存分析（尤其是Kaplan-Meier法）计算的是累计发病风险——当随访时间一致、无失访时，两者几乎等价；但如果存在失访、个体随访时间差异大的情况，生存分析能更准确地估计发病风险，因为它会对失访个体进行删失处理。

所以：

• 如果你的数据是个体水平的随访数据（包含每个研究对象的入组时间、随访终止时间、是否发病、失访状态等），生存分析完全适用，且是更严谨的选择。
• 如果你的数据是汇总统计数据（已经整理好各组的发病数、人年数/暴露人数），直接用流行病学方法计算IR和CI会更高效，不需要用到生存分析。

对于D2的基线现患情况：生存分析可以通过标记基线已发病的个体，但如果只是需要图表不从0开始，更简单的方式是在计算中加入基线现患率，直接调整可视化的y轴起点即可。

我会分两种数据类型（汇总数据/个体数据）给出方案，同时覆盖D1和D2的需求。

所需包 #

先安装并加载必要的包：

install.packages(c("epiR", "ggplot2", "dplyr", "survival", "survminer"))
library(epiR)
library(ggplot2)
library(dplyr)
library(survival)
library(survminer)

方案1：汇总数据（D1示例） #

假设你的D1是长格式汇总数据，包含分组、事件类型、发病数、人年数：

步骤1：模拟/准备数据 #

set.seed(123) # 固定随机种子，方便复现
D1 <- expand.grid(group = paste0("G", 1:5), event = paste0("E", 1:8)) %>%
  mutate(
    cases = sample(5:20, nrow(.), replace = TRUE), # 发病数
    person_years = sample(800:1500, nrow(.), replace = TRUE) # 随访人年数
  )

步骤2：计算IR与95%CI #

用epiR包的epi.incidence函数，基于泊松分布计算置信区间：

D1_ir <- D1 %>%
  rowwise() %>%
  mutate(
    # 计算每1000人年的发病率，方便解读
    incidence = epi.incidence(cases, person_years, method = "poisson")$est * 1000,
    ir_lower = epi.incidence(cases, person_years, method = "poisson")$lower * 1000,
    ir_upper = epi.incidence(cases, person_years, method = "poisson")$upper * 1000
  ) %>%
  ungroup()

步骤3：分组绘制IR图表 #

为每个人群组单独生成图表，展示8类事件的IR及CI：

# 批量生成图表
plots_D1 <- D1_ir %>%
  group_by(group) %>%
  group_map(function(data, group_info) {
    ggplot(data, aes(x = event, y = incidence)) +
      geom_col(fill = "#2c3e50", width = 0.7) +
      geom_errorbar(aes(ymin = ir_lower, ymax = ir_upper), width = 0.2) +
      labs(
        title = paste("累计发病率（每1000人年） - 组", group_info$group),
        x = "事件类型",
        y = "累计发病率（95%CI）"
      ) +
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))
  })

# 查看第一组的图表
plots_D1[[1]]

方案2：汇总数据（D2示例，含基线现患） #

假设D2在D1的基础上，增加了基线现患病例数和基线暴露人数：

步骤1：模拟/准备数据 #

set.seed(456)
D2 <- expand.grid(group = paste0("G", 1:5), event = paste0("E", 1:8)) %>%
  mutate(
    cases = sample(5:20, nrow(.), replace = TRUE),
    person_years = sample(800:1500, nrow(.), replace = TRUE),
    prev_cases = sample(2:10, nrow(.), replace = TRUE), # 基线现患数
    prev_n = person_years # 假设基线暴露人数与随访人数一致，可根据实际调整
  )

步骤2：计算调整后的病例率 #

我们计算基线患病率 + 随访累计发病率，并保留随访IR的置信区间：

D2_ir <- D2 %>%
  rowwise() %>%
  mutate(
    prev_rate = (prev_cases / prev_n) * 1000, # 每1000人的基线患病率
    incidence = epi.incidence(cases, person_years, method = "poisson")$est * 1000,
    ir_lower = epi.incidence(cases, person_years, method = "poisson")$lower * 1000,
    ir_upper = epi.incidence(cases, person_years, method = "poisson")$upper * 1000,
    total_rate = prev_rate + incidence # 总病例率（现患+新增）
  ) %>%
  ungroup()

步骤3：绘制非零起点的图表 #

图表y轴从基线患病率附近开始，同时标注基线现患的参考线：

plots_D2 <- D2_ir %>%
  group_by(group) %>%
  group_map(function(data, group_info) {
    ggplot(data, aes(x = event, y = total_rate)) +
      geom_col(fill = "#e74c3c", width = 0.7) +
      geom_errorbar(aes(ymin = prev_rate + ir_lower, ymax = prev_rate + ir_upper), width = 0.2) +
      geom_hline(aes(yintercept = prev_rate), linetype = "dashed", color = "#95a5a6", size = 0.8) +
      labs(
        title = paste("累计病例率（含基线现患） - 组", group_info$group),
        x = "事件类型",
        y = "累计病例率（每1000人，95%CI）"
      ) +
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
      # 手动设置y轴范围，确保基线现患线可见
      ylim(min(data$prev_rate)*0.8, max(data$total_rate + data$ir_upper)*1.1)
  })

# 查看第一组的图表
plots_D2[[1]]

方案3：个体水平数据（生存分析实现） #

如果你的数据是个体水平（每行一个研究对象），用survival包计算累计发病风险：

示例数据 #

# 模拟个体数据：group, event（是否发生E1）, time（随访时间）, status（1=发病，0=失访/未发病）
set.seed(789)
D1_individual <- expand.grid(group = paste0("G", 1:5), id = 1:200) %>%
  mutate(
    time = sample(1:10, nrow(.), replace = TRUE),
    status = rbinom(nrow(.), 1, 0.1) # 10%的发病概率
  )

计算累计发病风险（KM法） #

# 按分组和事件（这里以E1为例，多事件可循环处理）
km_fit <- survfit(Surv(time, status) ~ group, data = D1_individual)

# 绘制KM曲线（即累计发病风险曲线）
ggsurvplot(
  km_fit,
  data = D1_individual,
  conf.int = TRUE, # 显示95%CI
  risk.table = TRUE, # 显示风险人数表
  title = "累计发病风险（E1） - 各人群组",
  xlab = "随访时间",
  ylab = "累计发病风险"
)

对于多事件（E1-E8），可以通过循环为每个事件单独计算并绘图。

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