logo标签列表

如何绘制由R中survival包的survreg函数生成的生存曲线？

rplotsurvival-analysisweibull
37

我正在尝试将Weibull模型拟合和绘制到生存数据中。数据只有一个协变量，即队列，从2006年到2010年。那么，您有什么想法可以添加到以下两行代码中，以绘制2010年队列的生存曲线？

library(survival)
s <- Surv(subSetCdm$dur,subSetCdm$event)
sWei <- survreg(s ~ cohort,dist='weibull',data=subSetCdm)

使用Cox PH模型实现相同的任务相当简单，只需使用以下代码即可。问题在于survfit()不接受survreg类型的对象。

sCox <- coxph(s ~ cohort,data=subSetCdm)
cohort <- factor(c(2010),levels=2006:2010)
sfCox <- survfit(sCox,newdata=data.frame(cohort))
plot(sfCox,col='green')

使用survival包中的lung数据，这是我想要完成的目标。

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

#plot kaplan-meier estimate, per sex
fKM <- survfit(s ~ sex,data=lung)
plot(fKM)

#plot Cox PH survival curves, per sex
sCox <- coxph(s ~ as.factor(sex),data=lung)
lines(survfit(sCox,newdata=data.frame(sex=1)),col='green')
lines(survfit(sCox,newdata=data.frame(sex=2)),col='green')

#plot weibull survival curves, per sex, DOES NOT RUN
sWei <- survreg(s ~ as.factor(sex),dist='weibull',data=lung)
lines(survfit(sWei,newdata=data.frame(sex=1)),col='red')
lines(survfit(sWei,newdata=data.frame(sex=2)),col='red')
2
如果您发布一个完整的示例，我会尝试为您解决问题。我们需要 subSetCdm 对象。请尝试使用 dput(subSetCdm)。 - tim riffe
3
?predict.survreg 中有示例。 - Vincent Zoonekynd
5个回答
26

希望这可以帮到你，我没有犯什么误导性的错误：

复制自上面：

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

    #plot kaplan-meier estimate, per sex
    fKM <- survfit(s ~ sex,data=lung)
    plot(fKM)

    #plot Cox PH survival curves, per sex
    sCox <- coxph(s ~ as.factor(sex),data=lung)
    lines(survfit(sCox,newdata=data.frame(sex=1)),col='green')
    lines(survfit(sCox,newdata=data.frame(sex=2)),col='green')

对于Weibull分布，使用predict函数，参考Vincent的评论：

    #plot weibull survival curves, per sex,
    sWei <- survreg(s ~ as.factor(sex),dist='weibull',data=lung)

    lines(predict(sWei, newdata=list(sex=1),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
    lines(predict(sWei, newdata=list(sex=2),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")

图表输出

这里的技巧是在绘图和预测时反转分位数顺序。可能有更好的方法，但它在这里有效。祝你好运！

Tim，有个问题。如果你想重新创建上面的内容，但不按性别子集...例如从- sWei <- survreg(s ~ 1,dist='weibull',data=lung) 开始。你会如何更改你的预测函数中newdata部分的规定？我试图理解你在上面是如何指定的... - Chris
嗨，Chris，我被难住了，抱歉，也许其他回答者中有人知道。如果没有，那么可能需要提出一个新问题。 - tim riffe
2
考虑到你的技巧，反转分位数序列（我在这里称之为q.seq）：如果我理解正确，predict期望死亡分位数而不是生存分位数。因此，我认为最好提供1-q.seq而不是rev(q.seq)。在你的情况下，这并不重要，因为你的q.seq是对称的，从0.01的死亡概率（= 0.99的生存概率）到0.99（= 0.01的生存概率）。但是，如果你从0.01到0.8的死亡概率，就会有所不同：你的补充序列应该是0.99...0.2，而不是0.8到0.01。 - Igor F.
我在这里扩展了Igor的评论：https://dev59.com/JGox5IYBdhLWcg3whUqm#30215350 - Paul
19

另一种选择是使用包flexsurv。相比于survival包，它提供了一些额外的功能 - 其中参数回归函数flexsurvreg()有一个很好的绘图方法，可以满足您的要求。

使用如上所述的lung数据集；

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

require(flexsurv)
sWei  <- flexsurvreg(s ~ as.factor(sex),dist='weibull',data=lung)
sLno  <- flexsurvreg(s ~ as.factor(sex),dist='lnorm',data=lung)   

plot(sWei)
lines(sLno, col="blue")

plot.flexsurvreg的输出结果

您可以使用type参数在累积风险或风险尺度上绘制图表，并使用ci参数添加置信区间。

10

这只是一份说明，澄清Tim Riffe的回答，该回答使用以下代码：

lines(predict(sWei, newdata=list(sex=1),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
lines(predict(sWei, newdata=list(sex=2),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")

这两个镜像序列 seq(.01,.99,by=.01) 和 seq(.99,.01,by=-.01) 的原因是 predict() 方法提供了事件分布 f(t) 的分位数 - 也就是 f(t) 的反 CDF 值，而生存曲线则是以 1-(f的CDF) 为 y 轴，t 为 x 轴。换句话说，如果你绘制 p 和 predict(p) 的图表，你将得到 CDF，而如果你绘制 1-p 和 predict(p) 的图表，你将得到生存曲线，即 1-CDF。以下代码更加透明且适用于任意的 p 值向量：

pct <- seq(.01,.99,by=.01)
lines(predict(sWei, newdata=list(sex=1),type="quantile",p=pct),1-pct,col="red")
lines(predict(sWei, newdata=list(sex=2),type="quantile",p=pct),1-pct,col="red")
3

如果有人想在 ggplot2 生态系统中将威布尔分布添加到 Kaplan-Meyer 曲线中，可以按照以下步骤进行：

library(survminer)
library(tidyr)

s <- with(lung,Surv(time,status))
fKM <- survfit(s ~ sex,data=lung)
sWei <- survreg(s ~ as.factor(sex),dist='weibull',data=lung)

pred.sex1 = predict(sWei, newdata=list(sex=1),type="quantile",p=seq(.01,.99,by=.01))
pred.sex2 = predict(sWei, newdata=list(sex=2),type="quantile",p=seq(.01,.99,by=.01))

df = data.frame(y=seq(.99,.01,by=-.01), sex1=pred.sex1, sex2=pred.sex2)
df_long = gather(df, key= "sex", value="time", -y)

p = ggsurvplot(fKM, data = lung, risk.table = T)
p$plot = p$plot + geom_line(data=df_long, aes(x=time, y=y, group=sex))

enter image description here

1
如果您想使用生存函数本身S(t)（而不是其他答案中使用的逆生存函数S^{-1}(p)），我编写了一个函数来实现威布尔分布的情况（遵循与pec::predictSurvProb函数系列相同的输入：）
survreg.predictSurvProb <- function(object, newdata, times){
  shape <- 1/object$scale # also equals 1/exp(fit$icoef[2])
  lps <- predict(object, newdata = newdata, type = "lp")
  surv <- t(sapply(lps, function(lp){
    sapply(times, function(t) 1 - pweibull(t, shape = shape, scale = exp(lp)))
  }))
  
  return(surv)
}

您可以执行以下操作：

sWei <- survreg(s ~ as.factor(sex),dist='weibull',data=lung)
times <- seq(min(lung$time), max(lung$time), length.out = 1000)
new_dat <- data.frame(sex = c(1,2))
surv <- survreg.predictSurvProb(sWei, newdata = new_dat, times = times)

lines(times, surv[1, ],col='red')
lines(times, surv[2, ],col='red')

enter image description here

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接