经IPTW调整的药物心血管死亡竞争风险分析及Kaplan Meier曲线绘制

• 包名调用错误：crr是cmprsk包下的函数，你写的library(crr)会报错，无名为crr的R包。
• 研究终点设置错误：你要分析心血管死亡（编码为1）的效应，但crr的failcode参数设为了2，实际分析的是非心血管死亡的关联，完全不符合研究目标。
• 权重使用逻辑错误：你将IPTW权重放入了协变量矩阵cov1中，crr函数不会自动识别这部分为权重，正确做法是通过crr的weights参数传入加权值。
• 缺少权重稳健性处理：未对极端IPTW权重做截断，极端值会大幅拉低效应估计的稳定性。
• 绘图逻辑误区：存在竞争事件时，传统Kaplan-Meier曲线会高估目标结局的累积发生风险，针对竞争风险场景应绘制加权累积发生率（CIF）曲线。

# 加载正确的依赖包
library(tableone)
library(cmprsk)
library(ipw)
library(sandwich)
library(survey)
library(ggplot2)
library(survminer)

# 1. 计算IPTW权重，添加截断处理
# 倾向得分模型，可根据研究需求调整协变量
psmodel <- glm(treatment ~ age + sex + BMI, data = df, family = binomial(link = "logit"))
ps <- predict(psmodel, type = "response")

# 计算IPTW并做1%、99%分位数截断，避免极端权重影响
weight <- ifelse(df$treatment == 1, 1/ps, 1/(1-ps))
weight_cut <- pmax(pmin(weight, quantile(weight, 0.99)), quantile(weight, 0.01))

# 2. 协变量平衡检验（可选但推荐，用于验证IPTW效果）
vars <- c("age", "sex", "BMI")
# 加权前平衡
tab_unweighted <- CreateTableOne(vars = vars, strata = "treatment", data = df)
# 加权后平衡
svydes <- svydesign(ids = ~1, weights = ~weight_cut, data = df)
tab_weighted <- svyCreateTableOne(vars = vars, strata = "treatment", design = svydes)
print(tab_unweighted, smd = TRUE)
print(tab_weighted, smd = TRUE)

# 3. 加权竞争风险模型，目标结局为心血管死亡（编码1）
# 协变量矩阵仅放入核心暴露变量treatment
cov_mat <- as.matrix(df$treatment)
competingrisk <- crr(
  ftime = df$Survival, 
  fstatus = df$Outcome, 
  cov1 = cov_mat, 
  failcode = 1, # 正确指定目标结局编码
  cencode = 0, # 指定删失编码为存活（0）
  weights = weight_cut # 传入截断后的IPTW权重
)
# 输出模型结果，得到治疗对心血管死亡的效应值（HR及95%CI）
summary(competingrisk)

# 4. 绘制加权累积发生率（CIF）曲线，为竞争风险场景下的标准生存曲线
# 分别计算治疗组和对照组的CIF
cif_trt <- predict(crr(
  ftime = df$Survival[df$treatment == 1], 
  fstatus = df$Outcome[df$treatment == 1], 
  cov1 = as.matrix(rep(1, sum(df$treatment == 1))), 
  failcode = 1, 
  cencode = 0, 
  weights = weight_cut[df$treatment == 1]
), times = sort(unique(df$Survival)))

cif_ctrl <- predict(crr(
  ftime = df$Survival[df$treatment == 0], 
  fstatus = df$Outcome[df$treatment == 0], 
  cov1 = as.matrix(rep(0, sum(df$treatment == 0))), 
  failcode = 1, 
  cencode = 0, 
  weights = weight_cut[df$treatment == 0]
), times = sort(unique(df$Survival)))

# 整理绘图数据
plot_df <- rbind(
  data.frame(time = as.numeric(names(cif_trt[,1])), cif = cif_trt[,1], group = "治疗组"),
  data.frame(time = as.numeric(names(cif_ctrl[,1])), cif = cif_ctrl[,1], group = "对照组")
)

# 绘制CIF曲线
ggplot(plot_df, aes(x = time, y = cif, color = group)) +
  geom_step(linewidth = 1) +
  labs(x = "生存时间", y = "心血管死亡累积发生率", color = "分组") +
  theme_bw()

# 若确实需要绘制传统加权Kaplan-Meier曲线（将竞争事件当作删失处理），可运行以下代码
df$status_km <- ifelse(df$Outcome == 1, 1, 0) # 仅心血管死亡为事件，其余（竞争事件、存活）为删失
km_fit <- svykm(Surv(Survival, status_km) ~ treatment, design = svydes)
ggsurvplot(km_fit, data = df, risk.table = TRUE, xlab = "生存时间", ylab = "心血管死亡无事件生存率", legend.labs = c("对照组","治疗组"))

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