无竞争风险的两组累积发病率分析及绘图技术咨询

修正后的基础数据集代码 #

首先补全数据集定义（原代码缺少ID变量）：

ID <- 1:25
time<-c(1.5989,6.9433, 0.8890, 3.2691, 1.0514, 2.7625, 1.4319, 0.9681, 7.4416, 0.0268, 1.5168, 1.9647, 0.0657, 4.3571, 6.4490, 0.2198, 1.2028, 0.9555, 0.2601, 2.0096, 7.5156, 0.4463, 0.2355, 0.9391, 2.6996)
censor<-c(1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0)
group<-c(1, 2, 1, 1, 2, 2, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 2, 1, 2, 1)
df<-data.frame(ID, time, censor, group)
# 将group转为因子，方便后续分组显示
df$group <- factor(df$group, labels = c("治疗组A", "治疗组B"))

问题1：按治疗组获取随时间变化的累积发病率，并在图下方显示风险人数和事件数 #

使用survival包拟合KM模型，结合survminer包的ggsurvplot函数，可直接生成带风险表的累积发病率曲线：

library(survival)
library(survminer)

# 拟合KM生存模型
km_fit <- survfit(Surv(time, censor) ~ group, data = df)

# 绘制累积发病率曲线+风险表
ggsurvplot(
  km_fit,
  data = df,
  fun = "event", # 直接绘制累积发病率（等价于1-生存概率）
  xlab = "随访时间（月）",
  ylab = "HBV感染累积发病率",
  legend.title = "治疗组",
  risk.table = TRUE, # 启用风险表
  risk.table.title = "风险人数与事件数",
  risk.table.col = "strata", # 风险表按组着色
  risk.table.y.text.col = TRUE,
  risk.table.y.text = FALSE,
  ggtheme = theme_bw()
)

风险表会自动展示每个时间点的风险人数（n.risk）、事件数（n.event）及删失数（n.censor）。

问题2：在累积发病率图上显示log rank检验的p值 #

通过ggsurvplot的pval参数可直接添加log rank检验结果，无需额外计算：

ggsurvplot(
  km_fit,
  data = df,
  fun = "event",
  xlab = "随访时间（月）",
  ylab = "HBV感染累积发病率",
  legend.title = "治疗组",
  risk.table = TRUE,
  risk.table.title = "风险人数与事件数",
  risk.table.col = "strata",
  pval = TRUE, # 添加log rank检验p值
  pval.coord = c(15, 0.8), # 自定义p值显示位置（可选）
  ggtheme = theme_bw()
)

若需手动计算p值，可使用survdiff函数：

lr_test <- survdiff(Surv(time, censor) ~ group, data = df)
p_val <- 1 - pchisq(lr_test$chisq, df = length(lr_test$n)-1)

问题3：获取指定时间点的累积发病率、标准误及95%置信区间 #

使用summary.survfit的times参数指定目标时间点，再通过生存概率推导累积发病率的相关指标：

# 指定需要提取的随访时间点
target_times <- c(0, 6, 12, 18, 24)

# 提取指定时间点的生存分析结果
km_summary <- summary(km_fit, times = target_times, extend = TRUE)

# 整理成结果数据框
ci_table <- data.frame(
  随访时间 = km_summary$time,
  治疗组 = km_summary$strata,
  累积发病率 = 1 - km_summary$surv,
  标准误 = km_summary$std.err, # 累积发病率标准误与生存概率标准误一致（Var(1-S)=Var(S)）
  95%CI下限 = 1 - km_summary$upper, # 反转生存概率置信区间上限
  95%CI上限 = 1 - km_summary$lower  # 反转生存概率置信区间下限
)

# 输出结果
print(ci_table, row.names = FALSE)

• extend=TRUE确保超过最大随访时间的点沿用最后一个时间点的累积发病率结果
• 累积发病率的置信区间通过反转生存概率的置信区间得到，因为两者是互补关系

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