如何绘制由survreg（R的生存包）生成的生存曲线？

Question

我想对一个生存数据进行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')

使用数据肺（来自生存包），这是我要完成的任务。

#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')

Answer 1

希望能帮到你，我没有犯什么误导性的错误。

复制自上面的内容

    #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，重新从文森特的评论。

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

这里的诀窍是颠倒了绘制与预测的分位数顺序。可能还有更好的方法，但在这里是可行的。祝您好运

Answer 2

另一个选择是使用包中的 flexsurv. 这提供了一些额外的功能，比 survival 包--包括参数回归函数 flexsurvreg() 有一个很好的情节方法，它可以做到你所要求的。

如上图所示，使用肺部。

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

你可以在累积危险或危险等级上绘制，使用的是 type 参数，并添加置信区间与 ci 争论。

Answer 3

这只是一个说明，澄清 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)换句话说，如果你绘制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")

Answer 4

如果有人想在Kaplan -Meyer曲线中加入Weibull分布的话 ggplot2 生态系统，我们可以做到以下几点。

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

Answer 5

如果你想使用生存功能本身的话 S(t) (而不是反生存函数 S^{-1}(p) 我写了一个函数来实现魏布尔分布的情况下的函数（输入与 pec::predictSurvProb 职能系列。

predictSurvProb.survreg <- function(object, newdata, times){
  lp <- predict(object, newdata = newdata, 
                type = "lp")
  wei_shape <- 1/exp(object$icoef[2])

  surv <- sapply(lp, function(lp_i){
    1 - pweibull(times, shape = wei_shape, scale = exp(lp_i))
  })

  return(t(surv))
}