将Weibull模型拟合并绘制成生存数据

问题描述

我想实现与这个问题完全相同的事情: How to plot the survival curve generated by survreg (package survival of R)?

除了我不希望数据按变量分层(在上面的问题中按性别分层)之外。

我只希望整个治疗患者的无进展生存期。

因此,当我从另一个问题复制代码时,这就是我遇到的问题:

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)  # in my case here I would replace as.factor(sex) by 1

pred.sex1 = predict(sWei,newdata=list(sex=1),type="quantile",p=seq(.01,.99,by=.01)) #Since I don't want to stratify,what do I do with these 2 lines of code?
pred.sex2 = predict(sWei,newdata=list(sex=2),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))

我尝试将as.factor(sex)替换为1,然后剩下的代码就没有意义了,有人可以帮我吗?

非常感谢!

解决方法

如果您只想绘制整体经验生存曲线,则可以执行以下操作:

library(survival)
library(survminer)
library(tidyr)

s   <- with(lung,Surv(time,status))
fKM <- survfit(s ~ 1,data = survival::lung)
ggsurvplot(fKM,ggtheme = theme_bw())

但是,如果您想在没有预测变量的情况下拟合Weibull模型,那么您的公式就可以了。

sWei  <- survreg(s ~ 1,dist = 'weibull',data = lung) 
probs <- seq(0.01,1,by = 0.01)
time  <- predict(sWei,type = "quantile",se = TRUE,p = probs)

唯一的问题是time现在是两个矩阵的命名列表:fitse.fit。两者都具有与lung相同的行数,但是所有行都是相同的,因此我们只从每一行中取一个并计算数据帧中的置信区间,然后可以使用该置信区间来创建ggplot

ggplot(data = data.frame(p     = 1 - probs,time  = time$fit[1,],upper = time$fit[1,] + 1.96 * time$se.fit[1,lower = time$fit[1,] - 1.96 * time$se.fit[1,])) + 
  geom_step(aes(p,time,colour = "All"),size = 1) +
  geom_ribbon(aes(p,ymin = lower,ymax = upper,fill = "All"),alpha = 0.2) +
  coord_flip(ylim = c(0,1000)) +
  scale_fill_discrete(name = "Strata") +
  scale_color_discrete(name = "Strata") +
  theme_bw() +
  theme(legend.position = "top")

enter image description here

我们可以看到看起来非常合适。

如果您希望两者都在同一个情节中,则可以执行以下操作:

df <- data.frame(p     = 1 - probs,])

ggsurvplot(fKM,ggtheme = theme_bw())$plot +
  geom_line(data = df,aes(time,p),linetype = 2,size = 1) +
  geom_line(data = df,aes(upper,aes(lower,size = 1)

enter image description here

reprex package(v0.3.0)于2020-08-18创建