如何在R中使用Copula模型生成相关生存终点（如PFS和OS）
2026-6-15
多生存终点（PFS/OS）模拟与Copula模型实现（R语言）
一、实现逻辑
PFS（无进展生存期）与OS（总生存期）存在天然正相关性，Copula模型可以通过构建联合分布生成符合这种相关性的成对生存数据。这里以Clayton Copula为例（适合刻画正相关的生存事件），核心步骤如下：
- 定义PFS和OS的边缘生存分布（沿用原函数的指数分布）；
- 用Copula生成两个相关的均匀分布随机数；
- 将均匀分布转换为对应的生存时间（基于边缘分布的分位数函数）；
- 结合入组时间、失访时间和研究时长，计算最终的观测时间与删失状态。
需提前安装并加载依赖包：
install.packages(c("copula", "survival")) library(copula) library(survival)
二、扩展后的多终点模拟函数
基于你提供的单终点函数修改，支持生成带相关性的PFS和OS数据：
# 多生存终点（PFS/OS）模拟函数 generate_multiple_survival_data <- function( seed = NULL, mt_pfs_trt = NULL, # 试验组PFS中位时间（向量，对应各试验臂） mt_os_trt = NULL, # 试验组OS中位时间（向量，对应各试验臂） mt_pfs_ctrl = NULL, # 对照组PFS中位时间 mt_os_ctrl = NULL, # 对照组OS中位时间 copula_param = 2, # Copula参数（Clayton Copula，>0表示正相关，值越大相关性越强） n_trt = NULL, # 每个试验组样本量 n_ctrl = NULL, # 对照组样本量 enrl_rate = NULL, # 入组速率 drpot_hr = NULL, # 失访率（指数分布参数） nsim = NULL, # 模拟次数 study_duration = NULL,# 研究总时长 arms = NULL # 试验臂数量 ) { set.seed(seed) # 初始化Clayton Copula clayton_cop <- claytonCopula(param = copula_param, dim = 2) results <- data.frame() for (j in 1:nsim) { # 生成各试验臂数据 for (arm in 1:arms) { for (i in 1:n_trt) { subjid <- paste(j, arm, i, sep = "-") enrl_tm <- runif(1, 0, min(n_trt / enrl_rate, study_duration)) # 用Copula生成相关的均匀随机数 u <- rCopula(1, clayton_cop) # 转换为PFS和OS的生存时间（指数分布分位数函数） pfs_tm <- qexp(u[1], rate = log(2)/mt_pfs_trt[arm]) os_tm <- qexp(u[2], rate = log(2)/mt_os_trt[arm]) # 失访时间 drpot_tm <- rexp(1, rate = drpot_hr) # 计算PFS的观测时间与删失状态 if (pfs_tm > drpot_tm) { obs_pfs_tm <- drpot_tm cal_pfs_tm <- enrl_tm + obs_pfs_tm cnsr_pfs <- ifelse(cal_pfs_tm > study_duration, 1, 1) } else { obs_pfs_tm <- pfs_tm cal_pfs_tm <- enrl_tm + obs_pfs_tm cnsr_pfs <- ifelse(cal_pfs_tm > study_duration, 1, 0) } # 计算OS的观测时间与删失状态（OS不能早于PFS） os_actual_tm <- pfs_tm + os_tm if (os_actual_tm > drpot_tm) { obs_os_tm <- drpot_tm cal_os_tm <- enrl_tm + obs_os_tm cnsr_os <- ifelse(cal_os_tm > study_duration, 1, 1) } else { obs_os_tm <- os_actual_tm cal_os_tm <- enrl_tm + obs_os_tm cnsr_os <- ifelse(cal_os_tm > study_duration, 1, 0) } # 修正：PFS删失时，OS起始时间为观测到的PFS时间 if (cnsr_pfs == 1) { os_actual_tm <- obs_pfs_tm + os_tm if (os_actual_tm > drpot_tm) { obs_os_tm <- drpot_tm cal_os_tm <- enrl_tm + obs_os_tm cnsr_os <- ifelse(cal_os_tm > study_duration, 1, 1) } else { obs_os_tm <- os_actual_tm cal_os_tm <- enrl_tm + obs_os_tm cnsr_os <- ifelse(cal_os_tm > study_duration, 1, 0) } } results <- rbind(results, data.frame( nsim = j, subjid = subjid, arm = as.character(arm), enrl_tm = enrl_tm, pfs_tm = pfs_tm, os_tm = os_actual_tm, drpot_tm = drpot_tm, obs_pfs_tm = obs_pfs_tm, cal_pfs_tm = cal_pfs_tm, cnsr_pfs = cnsr_pfs, obs_os_tm = obs_os_tm, cal_os_tm = cal_os_tm, cnsr_os = cnsr_os )) } } # 生成对照组数据 for (i in 1:n_ctrl) { subjid <- paste(j, "CTRL", i, sep = "-") enrl_tm <- runif(1, 0, min(n_ctrl / enrl_rate, study_duration)) u <- rCopula(1, clayton_cop) pfs_tm <- qexp(u[1], rate = log(2)/mt_pfs_ctrl) os_tm <- qexp(u[2], rate = log(2)/mt_os_ctrl) drpot_tm <- rexp(1, rate = drpot_hr) # PFS观测与删失 if (pfs_tm > drpot_tm) { obs_pfs_tm <- drpot_tm cal_pfs_tm <- enrl_tm + obs_pfs_tm cnsr_pfs <- ifelse(cal_pfs_tm > study_duration, 1, 1) } else { obs_pfs_tm <- pfs_tm cal_pfs_tm <- enrl_tm + obs_pfs_tm cnsr_pfs <- ifelse(cal_pfs_tm > study_duration, 1, 0) } # OS观测与删失 os_actual_tm <- pfs_tm + os_tm if (os_actual_tm > drpot_tm) { obs_os_tm <- drpot_tm cal_os_tm <- enrl_tm + obs_os_tm cnsr_os <- ifelse(cal_os_tm > study_duration, 1, 1) } else { obs_os_tm <- os_actual_tm cal_os_tm <- enrl_tm + obs_os_tm cnsr_os <- ifelse(cal_os_tm > study_duration, 1, 0) } if (cnsr_pfs == 1) { os_actual_tm <- obs_pfs_tm + os_tm if (os_actual_tm > drpot_tm) { obs_os_tm <- drpot_tm cal_os_tm <- enrl_tm + obs_os_tm cnsr_os <- ifelse(cal_os_tm > study_duration, 1, 1) } else { obs_os_tm <- os_actual_tm cal_os_tm <- enrl_tm + obs_os_tm cnsr_os <- ifelse(cal_os_tm > study_duration, 1, 0) } } results <- rbind(results, data.frame( nsim = j, subjid = subjid, arm = "0", enrl_tm = enrl_tm, pfs_tm = pfs_tm, os_tm = os_actual_tm, drpot_tm = drpot_tm, obs_pfs_tm = obs_pfs_tm, cal_pfs_tm = cal_pfs_tm, cnsr_pfs = cnsr_pfs, obs_os_tm = obs_os_tm, cal_os_tm = cal_os_tm, cnsr_os = cnsr_os )) } } # 按模拟次数、PFS观测时间排序 results <- results[order(results$nsim, results$cal_pfs_tm, results$cnsr_pfs), ] return(results) }
三、函数调用示例
# 模拟参数设置 sim_data <- generate_multiple_survival_data( seed = 123, mt_pfs_trt = c(8, 10), # 2个试验臂的PFS中位时间 mt_os_trt = c(18, 22), # 2个试验臂的OS中位时间 mt_pfs_ctrl = 6, # 对照组PFS中位时间 mt_os_ctrl = 15, # 对照组OS中位时间 copula_param = 3, # Copula参数，相关性较强 n_trt = 100, # 每个试验组100人 n_ctrl = 200, # 对照组200人 enrl_rate = 20, # 每月入组20人 drpot_hr = log(2)/30, # 失访中位时间30个月 nsim = 5, # 模拟5次 study_duration = 24, # 研究时长24个月 arms = 2 ) # 查看前5行数据 head(sim_data, 5)
四、后续log-rank检验示例
以第一次模拟的数据为例，分别对PFS和OS做log-rank检验：
# 提取第一次模拟的数据 sim1_data <- subset(sim_data, nsim == 1) # PFS的log-rank检验 pfs_surv <- Surv(time = sim1_data$cal_pfs_tm, event = 1 - sim1_data$cnsr_pfs) pfs_test <- survdiff(pfs_surv ~ arm, data = sim1_data) print(pfs_test) # OS的log-rank检验 os_surv <- Surv(time = sim1_data$cal_os_tm, event = 1 - sim1_data$cnsr_os) os_test <- survdiff(os_surv ~ arm, data = sim1_data) print(os_test)
五、关键说明
- Copula选择：若需不同相关性结构，可替换为
frankCopula或gumbelCopula，仅需修改Copula初始化代码； - 边缘分布：示例用指数分布，可根据需求替换为Weibull、对数正态等分布，修改对应分位数函数即可（如
qweibull）； - OS与PFS关联：示例假设OS是PFS之后的生存时间，可根据实际研究设计调整两者的关联逻辑。
内容的提问来源于stack exchange，提问作者sina wang
相关产品推荐
相关产品推荐
