首先明确：完全可以为多变量Cox回归生成"number at risk"表，你遇到的报错是因为用法不对——ggsurvplot()需要的是survfit()生成的生存曲线对象，而不是直接传入coxph()的模型结果。单变量时你用了survfit()，所以能正常运行，但多变量时直接传coxph的fit对象，函数找不到需要的ggsurv结构，就报错了。

下面是具体的解决步骤和代码，适配你的数据集和模型：

步骤1：基于多变量Cox模型生成survfit对象 #

首先，你需要构造一个新数据集，用来指定你想展示的自变量分组（比如你想重点看chr1q21_status的不同水平），同时把其他协变量固定在典型值（比如中位数、均值或参考类别）。然后用这个新数据集生成生存曲线对象：

# 你的多变量Cox模型（保持不变）
fit <- coxph(Surv(time, event) ~ chr1q21_status + CCND1 + CRTM1 + IRF4, data = myeloma)

# 构造新数据集：指定chr1q21_status的所有水平，其他变量取中位数（如果是连续变量）
# 如果某个变量是分类变量，替换成unique(myeloma$变量名)
newdata <- expand.grid(
  chr1q21_status = unique(myeloma$chr1q21_status),
  CCND1 = median(myeloma$CCND1, na.rm = TRUE),
  CRTM1 = median(myeloma$CRTM1, na.rm = TRUE),
  IRF4 = median(myeloma$IRF4, na.rm = TRUE)
)

# 基于Cox模型和新数据集生成survfit对象
surv_fit <- survfit(fit, newdata = newdata)

步骤2：用ggsurvplot绘制带risk table的生存曲线 #

现在把生成的surv_fit对象传入ggsurvplot()，就能正常生成带"number at risk"表的图了：

library(survminer)

ggsurvplot(
  surv_fit, 
  data = myeloma, 
  risk.table = TRUE, 
  break.time.by = 365, 
  xlim = c(0, 4000), 
  risk.table.y.text = FALSE, 
  legend.labs = c("2", "3", "4+"),
  xlab = "Time (days)", # 可选：添加x轴标签
  risk.table.title = "Number at Risk" # 可选：给risk table加标题
)

额外注意事项 #

• 如果你的协变量是分类变量（比如二分类的CCND1），记得把median()替换成unique(myeloma$CCND1)，或者指定你想展示的参考水平。
• 确保你的survminer包是最新版本，旧版本可能存在兼容性问题，可通过update.packages("survminer")更新。
• 如果想展示所有协变量组合的生存曲线（不推荐，变量多的话曲线会过于杂乱），可以直接用surv_fit <- survfit(fit)，但需要调整图例和分组标签。

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