我明白你现在用survival包做竞争风险分析，已经拿到生存曲线，但卡在了获取各竞争状态的单独风险率，还要基于这个做模拟。结合你提到的compete vignette示例，我给你一步步拆解解决方案：

一、获取各竞争状态的单独风险率 #

在竞争风险里，我们需要的是原因别风险率（Cause-Specific Hazard, CSH）——也就是针对每个事件类型，把其他事件当作删失后拟合的风险率。你可以用两种方式实现：

1. 拟合分原因的Cox模型 #

这是最直接的方法，对每个竞争事件单独拟合Cox模型，把其他事件视为删失：

library(survival)
# 用pbc数据集（compete vignette里的示例数据）举例，status=1=死亡，2=移植，0=删失
data(pbc)
pbc$status <- ifelse(pbc$status == 2, 2, pbc$status) # 统一状态编码

# 拟合"死亡"事件的Cox模型（其他事件作为删失）
cox_death <- coxph(Surv(time, status == 1) ~ age + bili + albumin, data = pbc)
# 拟合"移植"事件的Cox模型（其他事件作为删失）
cox_transplant <- coxph(Surv(time, status == 2) ~ age + bili + albumin, data = pbc)

之后你可以提取每个模型的基线累积风险和个体的风险率：

# 获取基线累积风险（时间点+累积风险值）
bh_death <- basehaz(cox_death, centered = FALSE)
bh_transplant <- basehaz(cox_transplant, centered = FALSE)

# 对新个体计算线性预测器，得到个体的累积风险
new_patient <- data.frame(age = 50, bili = 1.5, albumin = 3.5)
lp_death <- predict(cox_death, newdata = new_patient)
lp_transplant <- predict(cox_transplant, newdata = new_patient)

# 个体的累积风险 = 基线累积风险 * exp(线性预测器)
H_death <- bh_death$hazard * exp(lp_death)
H_transplant <- bh_transplant$hazard * exp(lp_transplant)

2. 用多状态Cox模型统一提取 #

如果你想更系统地处理多状态转移，可以用survival包的多状态模型功能，直接提取各转移的风险率：

# 构建多状态生存对象
pbc$state <- factor(pbc$status, levels = c(0,1,2), labels = c("alive", "death", "transplant"))
ms_surv <- Surv(pbc$time, pbc$state, type = "mstate")

# 拟合多状态Cox模型
cox_ms <- coxph(ms_surv ~ age + bili + albumin, data = pbc)

# 针对新个体生成多状态风险拟合结果
ms_fit <- msfit(cox_ms, newdata = new_patient)

# 直接提取各状态的累积风险
H_death_ms <- ms_fit$Haz[, "alive->death"]
H_transplant_ms <- ms_fit$Haz[, "alive->transplant"]

这种方式更适合复杂的多状态场景，结果也更规整。

二、实现你计划的模拟操作 #

拿到各状态的风险率后，就可以按照你的思路做模拟了，步骤如下：

# ----------------------
# 步骤1：生成随机数x，反演生存曲线得到事件时间
# ----------------------
set.seed(123) # 设种子保证可重复
x <- runif(1) # 从[0,1]抽均匀随机数
H_target <- -log(x) # 总累积风险目标值（因为S(t)=exp(-H_total(t))=x → H_total(t)=-log(x)）

# 合并时间点并对齐累积风险（以多状态模型的结果为例）
all_times <- ms_fit$time
H_total <- H_death_ms + H_transplant_ms

# 找到第一个总累积风险≥H_target的时间点，就是事件发生时间
event_time_idx <- which.min(H_total < H_target)
event_time <- all_times[event_time_idx]

# ----------------------
# 步骤2：生成随机数y，确定结局状态
# ----------------------
# 用差分法近似瞬时风险率（h(t)≈ΔH/Δt）
if (event_time_idx == 1) {
  # 第一个时间点直接用累积风险/时间
  h_death <- H_death_ms[event_time_idx] / event_time
  h_transplant <- H_transplant_ms[event_time_idx] / event_time
} else {
  delta_t <- all_times[event_time_idx] - all_times[event_time_idx-1]
  h_death <- (H_death_ms[event_time_idx] - H_death_ms[event_time_idx-1]) / delta_t
  h_transplant <- (H_transplant_ms[event_time_idx] - H_transplant_ms[event_time_idx-1]) / delta_t
}

# 生成y并判断结局
y <- runif(1, 0, h_death + h_transplant)
event_status <- ifelse(y <= h_death, 1, 2) # 1=死亡，2=移植

# 输出结果
cat("模拟事件时间：", event_time, "\n模拟结局状态：", event_status, "\n")

小提示 #

• 如果你的生存曲线是累积发生率函数（CIF），而不是整体生存概率，那反演逻辑需要调整，但你的描述里是用生存曲线（S(t)=P(T>t)），所以上面的代码是适配的。
• 如果你需要批量模拟多个个体，把上面的代码写成循环或者用purrr包的函数批量处理即可。

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