在R语言中如何将两个不同模型的Group 0生存曲线绘于同一图

问题背景 #

我使用两个不同的二分类模型，二者测试样本数量不同且存在重叠。以预测得分的中位数为阈值，将每个模型的样本划分为Group 0和Group 1。目前能分别绘制两个模型Group 0的生存曲线，但无法合并到同一张图中展示。

现有代码 #

# `survivalCurve`是自定义函数，可绘制Group 0和Group 1的生存曲线，同时返回所用数据框
tme = survivalCurve(mymodel = model1, model_num = 8, use.TME = TRUE)
nontme = survivalCurve(mymodel = model2, model_num = 29, use.TME = FALSE)

# 提取两个模型的Group 0
tme_0 = tme %>% dplyr::filter(group == 0)
nontme_0 = nontme %>% dplyr::filter(group == 0)

# 尝试合并样本但仅能单独绘制曲线
tme_0$rowname = rownames(tme_0)
nontme_0$rowname = rownames(nontme_0)
zeros = merge(tme_0, nontme_0, by = "rowname", all = TRUE) 

zeros_x = zeros[!is.na(zeros$survival_time.x) & !is.na(zeros$survival_status.x), ]
zeros_y = zeros[!is.na(zeros$survival_time.y) & !is.na(zeros$survival_status.y), ]

# 创建生存对象
survival_object_x = Surv(time = zeros_x$survival_time.x, event = zeros_x$survival_status.x)
survival_object_y = Surv(time = zeros_y$survival_time.y, event = zeros_y$survival_status.y)

# 拟合生存曲线
fit_x = survfit(survival_object_x ~ 1)
fit_y = survfit(survival_object_y ~ 1)

# 分别绘制曲线
ggsurv_x = ggsurvplot(fit_x, data = zeros_x, risk.table = TRUE, pval = TRUE, conf.int = TRUE, legend.title = "Group X", legend.labs = c("Group X"), palette = c("blue"))
ggsurv_y = ggsurvplot(fit_y, data = zeros_y, risk.table = TRUE, pval = TRUE, conf.int = TRUE, legend.title = "Group Y", legend.labs = c("Group Y"), palette = c("red"))

数据结构示例 #

> dput(tme[1:20,])
structure(list(survival_time = c(9, 7, 8, 32, 27, 21, 5, 7, 25, 
15, 14, 30, 28, 6, 25, 27, 37, 63, 63, 22), survival_status = c(1, 
1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1), prediction_scores = c(0.452489739997856, 
0.440583099511596, 0.509308651793677, 0.457296723067907, 0.455967522911745, 
0.456091537587146, 0.52058697349597, 0.508446710238974, 0.468852537302421, 
0.401894923842953, 0.459980058832451, 0.488860355868262, 0.457574356551677, 
0.49155383636979, 0.470455235473811, 0.43061685052414, 0.428729855269241, 
0.442110813635852, 0.413582392521575, 0.407929513543207), group = structure(c(1L, 
1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 
1L, 1L, 1L), levels = c("0", "1"), class = "factor")), row.names = c("Nivolumab_2017_p002_ar_8815", 
"Nivolumab_2017_p005_ar_8883", "Nivolumab_2017_p010_ar_8887", 
"Nivolumab_2017_p031_ar_8913", "Nivolumab_2017_p034_ar_8929", 
"Nivolumab_2017_p037_ar_8900", "Nivolumab_2017_p039_ar_8819", 
"Nivolumab_2017_p046_ar_8904", "Nivolumab_2017_p072_ar_8861", 
"Nivolumab_2017_p077_ar_8846", "Nivolumab_2017_p082_ar_8822", 
"Nivolumab_2017_p085_ar_8829", "Nivolumab_2017_p089_ar_8831", 
"Nivolumab_2017_p090_ar_8866", "Nivolumab_2017_p098_ar_8853", 
"Nivolumab_2017_p101_ar_8834", "EA595720", "EA632149", "EA632174", 
"EA639099"), class = "data.frame")

> dput(nontme[1:20,])
structure(list(survival_time = c(5, 31, 24, 36, 16, 15, 3, 5, 
30, 11, 27, 63, 17, 59, 42, 55, 14, 53, 16, 2), survival_status = c(0, 
0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1), prediction_scores = c(0.399196720335709, 
0.399196712585732, 0.399196718778835, 0.399196729515291, 0.399196719003673, 
0.39919672164811, 0.399196722134395, 0.399196714987811, 0.399196722354385, 
0.399196713792278, 0.39919672386451, 0.395944775789647, 0.395944769442692, 
0.395944782920291, 0.395944768860647, 0.395944773326361, 0.395944782577356, 
0.395944769654135, 0.39594477072273, 0.395944782418024), rank = c(144, 
112, 139, 182, 141, 150, 151, 125, 154, 117, 159, 52, 35, 77, 
32, 46, 75, 38, 40, 74), group = structure(c(1L, 1L, 1L, 2L, 
1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L
), levels = c("0", "1"), class = "factor")), row.names = c("Nivolumab_2017_p024_ar_8906", 
"Nivolumab_2017_p026_ar_8920", "Nivolumab_2017_p028_ar_8922", 
"Nivolumab_2017_p036_ar_8898", "Nivolumab_2017_p062_ar_8856", 
"Nivolumab_2017_p077_ar_8846", "Nivolumab_2017_p078_ar_8864", 
"Nivolumab_2017_p084_ar_8850", "Nivolumab_2017_p085_ar_8829", 
"Nivolumab_2017_p092_ar_8867", "Nivolumab_2017_p101_ar_8834", 
"EA632174", "EA632234", "EA632688", "EA632802", "EA639069", "EA639120", 
"EA639131", "G109543_RCCBMS_00114_T_v1_RNA_OnPrem", "G109543_RCCBMS_00147_T_v1_RNA_OnPrem"
), class = "data.frame")

解决方案 #

核心思路是将两个数据集整理为长格式，添加分组标识区分两个模型的Group 0，再统一拟合生存曲线并绘制。

步骤1：整理数据集 #

给每个模型的Group 0数据添加模型标识，合并为单一数据框：

library(dplyr)
library(survival)
library(survminer)

# 提取并整理模型1的Group 0数据
tme_0_clean = tme_0 %>%
  mutate(model_group = "Model1_Group0") %>%
  select(model_group, survival_time, survival_status)

# 提取并整理模型2的Group 0数据
nontme_0_clean = nontme_0 %>%
  mutate(model_group = "Model2_Group0") %>%
  select(model_group, survival_time, survival_status)

# 合并两个数据集
combined_data = bind_rows(tme_0_clean, nontme_0_clean)

步骤2：拟合生存曲线 #

基于合并后的数据集，按模型分组拟合生存曲线：

# 创建统一的生存对象
surv_obj = Surv(time = combined_data$survival_time, event = combined_data$survival_status)

# 按模型分组拟合生存曲线
fit_combined = survfit(surv_obj ~ model_group, data = combined_data)

步骤3：绘制合并后的生存曲线 #

用ggsurvplot绘制包含两条曲线的图，设置图例、颜色等参数：

ggsurv_combined = ggsurvplot(
  fit_combined,
  data = combined_data,
  risk.table = TRUE,
  pval = TRUE,  # 自动计算两组的log-rank检验p值
  conf.int = TRUE,
  legend.title = "Groups",
  legend.labs = c("Model1 Group0", "Model2 Group0"),
  palette = c("blue", "red"),
  risk.table.y.text.col = TRUE,
  risk.table.y.text = FALSE
)

# 显示图形
print(ggsurv_combined)

重叠样本处理（可选） #

如果同一个样本在两个模型的Group 0中都存在，可根据行名去重避免重复统计：

# 添加样本ID列（行名）
tme_0_clean = tme_0 %>%
  mutate(sample_id = rownames(tme_0), model_group = "Model1_Group0") %>%
  select(sample_id, model_group, survival_time, survival_status)

nontme_0_clean = nontme_0 %>%
  mutate(sample_id = rownames(nontme_0), model_group = "Model2_Group0") %>%
  select(sample_id, model_group, survival_time, survival_status)

# 合并后去重（保留每个样本的一条记录）
combined_data_unique = bind_rows(tme_0_clean, nontme_0_clean) %>%
  distinct(sample_id, .keep_all = TRUE)

# 用去重后的数据拟合和绘图
surv_obj_unique = Surv(time = combined_data_unique$survival_time, event = combined_data_unique$survival_status)
fit_unique = survfit(surv_obj_unique ~ model_group, data = combined_data_unique)

ggsurv_unique = ggsurvplot(
  fit_unique,
  data = combined_data_unique,
  risk.table = TRUE,
  pval = TRUE,
  conf.int = TRUE,
  legend.title = "Groups",
  legend.labs = c("Model1 Group0", "Model2 Group0"),
  palette = c("blue", "red")
)
print(ggsurv_unique)

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