下面是一个取自this paper的示例。

url <- "http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt"
Rossi <- read.table(url, header=TRUE)
Rossi[1:5, 1:10]

#   week arrest fin age  race wexp         mar paro prio educ
# 1   20      1  no  27 black   no not married  yes    3    3
# 2   17      1  no  18 black   no not married  yes    8    4
# 3   25      1  no  19 other  yes not married  yes   13    3
# 4   52      0 yes  23 black  yes     married  yes    1    5
# 5   52      0  no  19 other  yes not married  yes    3    3

mod.allison <- coxph(Surv(week, arrest) ~ 
                        fin + age + race + wexp + mar + paro + prio,
                        data=Rossi)
mod.allison

# Call:
# coxph(formula = Surv(week, arrest) ~ fin + age + race + wexp + 
#    mar + paro + prio, data = Rossi)
#
#
#                   coef exp(coef) se(coef)      z      p
# finyes         -0.3794     0.684   0.1914 -1.983 0.0470
# age            -0.0574     0.944   0.0220 -2.611 0.0090
# raceother      -0.3139     0.731   0.3080 -1.019 0.3100 
# wexpyes        -0.1498     0.861   0.2122 -0.706 0.4800
# marnot married  0.4337     1.543   0.3819  1.136 0.2600
# paroyes        -0.0849     0.919   0.1958 -0.434 0.6600
# prio            0.0915     1.096   0.0286  3.194 0.0014
#
# Likelihood ratio test=33.3  on 7 df, p=2.36e-05  n= 432, number of events= 114

请注意，该模型使用fin, age, race, wexp, mar, paro, prio来预测arrest。如this document中所述，survfit()函数使用Kaplan-Meier估计值来估计生存率。

plot(survfit(mod.allison), ylim=c(0.7, 1), xlab="Weeks",
     ylab="Proportion Not Rearrested")

我们得到了生存率的曲线图（置信区间为95%）。对于累积风险率，你可以

# plot(survfit(mod.allison)$cumhaz)

但这并没有给予置信区间。不过，不用担心！我们知道H（t）= -ln（S（t）），我们有S（t）的置信区间。我们需要做的是

sfit <- survfit(mod.allison)
cumhaz.upper <- -log(sfit$upper)
cumhaz.lower <- -log(sfit$lower)
cumhaz <- sfit$cumhaz # same as -log(sfit$surv)

那就把这些画出来

plot(cumhaz, xlab="weeks ahead", ylab="cumulative hazard",
     ylim=c(min(cumhaz.lower), max(cumhaz.upper)))
lines(cumhaz.lower)
lines(cumhaz.upper)

您可能希望使用survfit(..., conf.int=0.50)获得75%和25%的波段，而不是97.5%和2.5%。

在X的第25和75百分位数处估计生存曲线的请求首先需要确定这些百分位数，并指定数据框中所有其他协变量的值，以用作survfit的新数据参数。：
可以使用来自Fox网站的其他resondent建议的数据，尽管在我的机器上需要构建一个url对象：

url <- url("http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt")
 Rossi <- read.table(url, header=TRUE)

这可能不是这个问题的最佳例子，但它确实有一个数字变量，我们可以计算四分位数：

> summary(Rossi$prio)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.000   2.000   2.984   4.000  18.000

这就是模型拟合和曲面拟合调用：

mod.allison <- coxph(Surv(week, arrest) ~ 
                         fin + age + race + prio ,
                         data=Rossi)
 prio.fit <- survfit(mod.allison, 
                     newdata= data.frame(fin="yes", age=30, race="black", prio=c(1,4) ))
 plot(prio.fit, col=c("red","blue"))

将混杂因素的值设置为固定值，并针对给定的X值绘制多个时间点的预测生存概率（正如@IRTFM在回答中所建议的），结果是一个"条件"效应估计。这不是标准的Kaplan-Meier估计量的用途，我也不认为这是最初的发帖者想要的。通常我们感兴趣的是"平均因果效应"。换句话说：如果在整个样本中X被设置为某个特定值x，那么生存概率是多少？
我们可以使用拟合的cox模型加上g计算来获得这个概率。在g计算中，我们将整个样本中的X值设置为x，然后使用cox模型来预测每个个体在t的生存概率。在过程中使用他们观察到的协变量值。然后我们简单地取那些预测的平均值来获得我们的最终估计。通过对一系列的时间点和X的可能值重复这个过程，我们得到了一个三维的生存曲面，然后我们可以用颜色标度来可视化这个曲面。
这可以使用我开发的contsurvplot R包来完成，如前面的答案所述：Converting survival analysis by a continuous variable to categorical或软件包的文档中。关于这个策略的更多信息可以在我关于这个主题的文章的预印本中找到：https://arxiv.org/pdf/2208.04644.pdf