使用 Cox 回归绘制 Kaplan-Meier 图表

Question

使用 Cox 回归绘制 Kaplan-Meier 图表

8

我在R中使用以下代码设置了一个Cox比例风险模型，用于预测死亡率。变量A、B和C仅被添加以避免混淆（例如年龄、性别、种族），但我们真正感兴趣的是预测因子X。X是一个连续变量。

cox.model <- coxph(Surv(time, dead) ~ A + B + C + X, data = df)

现在，我在绘制Kaplan-Meier曲线方面遇到了麻烦。我一直在寻找如何创建此图形的方法，但是我没有太多的运气。我不确定是否可以为Cox模型绘制Kaplan-Meier曲线？ Kaplan-Meier是否会调整我的协变量或者它不需要它们？

我尝试过以下方法，但被告知这是不正确的。

plot(survfit(cox.model), xlab = 'Time (years)', ylab = 'Survival Probabilities')

我也尝试绘制一个显示累计死亡风险的图表。我不知道自己是否做得正确，因为我试过几种不同的方法并得到了不同的结果。理想情况下，我想绘制两条线，一条显示X的第75个百分位数的死亡风险，另一条显示第25个百分位数的死亡风险。我该怎么做？

我可以列出我尝试过的所有其他方法，但我不想让任何人感到困惑！

非常感谢。

- Hims

1

并不是“KM曲线调整协变量”，而是可以从模型拟合中构建预测的阶跃函数生存曲线。大多数人会使用KM曲线来指代未经调整的生存曲线。您还需要指定所有变量以进行预测。请参见下面的编码示例。 - IRTFM

3个回答

4

首先需要确定X的25和75百分位数，并指定用作newdata参数的数据框中所有其他协变量的值，以便请求估计的生存曲线。

可以使用来自Fox网站的其他回复者建议的数据，尽管在我的机器上它需要构建一个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

这将是模型拟合和survfit调用的代码：

 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"))

- IRTFM

在红色曲线和蓝色曲线之间进行log-rank检验是否可行？ - jiggunjer

您可以在构建和估计数据的条件下，获取mod.allison模型中prio变量的对数秩检验结果。这符合您的需求吗？ - IRTFM

我认为是这样的。如果我对模型进行summary()，它会返回所有系数的单个logrank结果。但我想要在两个优先级四分位之间进行测试。 - jiggunjer

我认为没有一种方法可以按照您描述的方式进行统计测试。 - IRTFM

重新审视这个问题，我认为足够的做法是将带有 prio 协变量和不带 prio 协变量的两个模型进行比较，然后取出差异的偏差值并与自由度为一的卡方统计量进行比较。这是一个对数似然比检验。 - IRTFM

0

将混淆因素的值设置为固定值，并在给定X值的多个时间点上绘制预测的生存概率（正如@IRTFM在他的答案中建议的那样），会得到一个条件效应估计。这不是标准Kaplan-Meier估计器的用途，我认为这也不是原始帖子的意图。通常我们对平均因果效应感兴趣。换句话说：如果在整个样本中将X设置为某个特定值x，生存概率会是多少？

我们可以使用拟合的cox模型和g-computation来获得这个概率。在g-computation中，我们将X的值设置为整个样本中的x，然后使用cox模型预测每个个体在t时刻的生存概率，过程中使用他们观察到的协变量值。然后我们只需取这些预测的平均值即可得到最终估计值。通过在一系列时间点和可能的X值范围内重复此过程，我们可以获得一个三维生存曲面。然后我们可以使用颜色比例尺来可视化这个曲面。

这可以通过我开发的contsurvplot R包来实现，如在先前的回答中所讨论的：将连续变量的生存分析转换为分类变量或在该软件包的文档中。有关此策略的更多信息可以在我关于此主题的预印本文章中找到：https://arxiv.org/pdf/2208.04644.pdf

- Denzo

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- nathanesau · Accepted Answer

这是一个来自这篇文章的示例。

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。就像在这个文档中提到的那样，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%。