使用CoxPH预测和绘制生存曲线

Question

我正在尝试在R中预测并绘制一个（估计的）新观察值的生存曲线。使用“survival”库和“lung”数据集，我首先对数据进行了Cox比例风险模型拟合。然后，我试图预测并绘制一个假想的新观察值的生存曲线（我在“list”命令中输入了此假想新观察值的详细信息）。但是，这不起作用。

以下是我的代码：

#load library
library(survival)
data(lung)

#create survival object
s <- with(lung,Surv(time,status))

#create model
modelA  <- coxph(s ~ as.factor(sex)+age+ph.ecog+wt.loss+ph.karno,data=lung, model=TRUE)
summary(modelA)

#plot
plot(survfit(modelA), ylab="Probability of Survival",
     xlab="Time", col=c("red", "black", "black"))

#predict for a hypothetical NEW observation (here is where the error is)
lines(predict(modelA, newdata=list(sex=1, 
                                  age = 56, 
                                  ph.ecog = 1, 
                                  ph.karno = 50,
                                  wt.loss = 11),
               type="quantile",
               p=seq(.01,.99,by=.01)),
       seq(.99,.01,by=-.01),
       col="blue")
## Error in match.arg(type) : 
##  'arg' should be one of “lp”, “risk”, “expected”, “terms”, “survival”

有人知道我做错了什么吗？谢谢

- stats_noob

Answer 1

这就是survfit函数的作用。在您的例子中，您绘制了模型的survfit，但是您可以将newdata参数输入此函数，并生成这些数据的估计生存率。

如果我们重现您的例子：

library(survival)

s <- with(lung, Surv(time, status))

modelA  <- coxph(s ~ as.factor(sex) + age + ph.ecog + wt.loss + ph.karno,
                 data = lung, model = TRUE)

plot(survfit(modelA), ylab = "Probability of Survival",
     xlab = "Time", col = c("red", "black", "black"))

然后，我们可以根据您指定的协变量创建生存曲线，如下所示：

est <- survfit(modelA, newdata = data.frame(sex      = 1,
                                            age      = 56,
                                            ph.ecog  = 1, 
                                            ph.karno = 50,
                                            wt.loss  = 11))

现在，est是一个包含time和survival成员的S3对象，因此我们可以绘制一条蓝色线来跟踪具有给定协变量的个体的估计生存率，如下所示：

lines(est$time, est$surv, col = 'blue', type = 's')

或者单独绘制，并附带95%置信区间：

plot(est, ylab = "Probability of Survival",
     xlab = "Time", col = c("red", "black", "black"))

^{2022年5月26日由reprex包（v2.0.1）创建}

- Allan Cameron

Answer 2

请参考predict()函数的描述（您可以通过运行?predict.coxph在R帮助中打开它，或者在这里查看示例）： type - 预测值的类型。可选项为线性预测器("lp")、风险得分exp(lp)("risk")、给定协变量和随访时间的预期事件数("expected")以及线性预测器的术语("terms")。一个受试者的生存概率等于exp(-expected)。

您可以看到，您的type="quantile"与期望输入不符。如果您调用predict()而没有type参数，在您的情况下，它将默认使用lp (线性预测器)。
当您为对象modelA调用predict()函数时，它确定它是coxph类，因此应用predict.coxph()函数。像type="quantile"和p=seq(.01,.99,by=.01)这样的参数对于predict.coxph()是不可接受的（p被忽略，type引发错误）。它们在另一个函数predict.survreg()中使用 - 为了调用它，您的modelA对象必须是survreg类，即它应该是使用survreg()调用而不是coxph()调用创建的。

- Vasily A

---

谢谢您的回复！我想知道这个参数是否需要指定，还是可以留空让计算机自动赋值。我会尝试为这个参数设置另一个值。您是否熟悉统计学中的“生存分析”领域？我的整体方法看起来合理吗？谢谢！ - stats_noob

如果您将其留空，它将被视为“lp”。我对生存分析有基本的概念，不确定您的方法 - 可能详细解释一下您想要得到什么？ - Vasily A

当然，我明天早上会解释。但基本上我想在一些数据上制作一个生存Cox回归模型，然后使用该模型预测新观察值（例如医疗患者）的估计生存函数（Breslow方法）。根据新患者的估计生存曲线的外观，我将使用这些信息来决定谁先接受治疗（例如，如果患者A的估计生存曲线比患者B“陡峭”，那么患者A将首先接受治疗）。 - stats_noob

谢谢您的回复！我认为survreg()只能进行加速失效时间（AFT）模型（即参数回归）。是否可以使用survreg()进行Cox比例风险回归？ - stats_noob

据我理解，survreg()函数不运行比例风险模型，在这种情况下应该使用coxph()函数。 - Vasily A