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

Posted 2023-02-19

【中文标题】如何绘制survreg生成的生存曲线（R的包生存）？

【英文标题】：How to plot the survival curve generated by survreg (package survival of R)?

【发布时间】：2012-02-27 10:12:02

【问题描述】：

我正在尝试根据生存数据拟合和绘制 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')

【问题讨论】：

如果你发布了一个完整的例子，我会试着为你弄清楚。我们需要 subSetCdm 对象。试试 dput(subSetCdm)

?predict.survreg中有例子。

【参考方案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，来自 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")

这里的技巧是颠倒绘图与预测的分位数顺序。可能有更好的方法来做到这一点，但它在这里有效。祝你好运！

【讨论】：

蒂姆，快速提问。如果您想重新创建上述但不是按性别划分的子集...例如以-- sWei

嗨，克里斯，我很难过，对不起，但也许其他回答者之一知道。如果不是，那么可能是一个新问题。

考虑到你的把戏，颠倒分位数的顺序（我在这里称之为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。

【参考方案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 参数添加置信区间。

【参考方案3】：

这只是一个说明Tim Riffe's answer的注释，它使用以下代码：

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 的逆 CDF 的值(t) - 当生存曲线绘制 1-(f 的 CDF) 与 t 时。换句话说，如果您绘制 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")

【参考方案4】：

如果有人想将 Weibull 分布添加到 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))

【参考方案5】：

如果您想使用生存函数本身 S(t)（而不是此处其他答案中使用的逆生存函数 S^-1(p)），我已经编写了一个函数来实现 Weibull 分布的情况（遵循与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')