在运行生存后，如何将数据重组回到有意义的单个生存曲线中并从上述曲线中进行预测？

Question

I有一个COX模型，该模型具有时间依赖性，该模型通过将曲线拆分为3片而解决了。这样做并确认比例恢复后，我想对Cox模型进行预测。生存程序包文档显示了这一点，但仅用于非拆分数据。

如何将其转换为更容易解释的单个情节，可以在其中检查治疗与控制等...

最小示例

# Clean your workspace environment
# rm(list = ls())

# Clear the plots tab (if any plots are open)
# if (!is.null(dev.list())) dev.off()

# Clear the console (equivalent to Ctrl+L shortcut)
# cat("\014")

library("survival") 
library("survminer")

vet2 <-survSplit(Surv(time, status) ~., data= veteran, cut=c(90, 180), episode= "tgroup", id="id")
vfit2 <-coxph(Surv(tstart, time, status) ~ trt + prior + karno:strata(tgroup),data=vet2)
scurve0 <- survfit(vfit2, newdata =  vet2)
#plot(scurve0)
ggsurvplot(scurve0, data=vfit2,legend = "none" ) 

plotdata <- data.frame( 
                          trt = c(0,0,0,0,0,0),
                          prior = c(0,10,0,10,0,10),
                          karno = c(60,99,60,99,60,99),
                          tgroup = c(1,2,3)
)
scurve1 <- survfit(vfit2, newdata =  plotdata)
ggsurvplot(scurve0, data=vfit2,legend = "none" ) 
scurve1 <- survfit(vfit2, newdata =  plotdata)

这两个**绘图调用都会产生以下图。

我可以将其作为单个曲线恢复吗？我需要关闭传奇，因为地层会产生很多情况。当我尝试使用新数据绘制模型时，我会得到完全相同的图，因此我认为它不起作用

Answer 1

评论：

”特里·塞诺（Terry Terry Therneau）说他不知道如何汲取预测的生存 使用与时变系数的模型时的曲线。搜索 R-Help档案备反复讨论。这就是什么 他2年前说：“生存曲线”一段时间依赖 协变量是不容易定义的东西。阅读第10.2.4章 泰诺和葛兰布奇的书进行了讨论（在很大程度上很大程度上 通过我自己犯的许多错误所告知。）“

，我想您可以使用@EDM概述的方法：//stats.stackexchange.com/questions/questions/534516/survival-curve-curve-curve-curve-curve-fromtime-dipperent-dipperent-coefficients-in-the-cox-model

作为潜在的解决方法（基于

https：//cran.r-project.org/web/packages/survival/vignettes/timedep.pdf

）中描述的“日志（时间+20）”方法。 e.g。对于您的示例数据：

library(survival)
library(survminer)

vfit <- coxph(Surv(time, status) ~ ., data = veteran)
prop <- cox.zph(vfit)
prop
#>            chisq df       p
#> trt       0.2644  1 0.60712
#> celltype 15.2274  3 0.00163
#> karno    12.9352  1 0.00032
#> diagtime  0.0129  1 0.90961
#> age       1.8288  1 0.17627
#> prior     2.1656  1 0.14113
#> GLOBAL   34.5525  8 3.2e-05
plot(prop[3])
abline(h = vfit$coefficients[5], lwd = 2, lty = 2, col = "firebrick3")

vet2 <-survSplit(Surv(time, status) ~.,
                 data = veteran, cut = c(30, 60),
                 episode = "tgroup", id = "id")
vfit2 <-coxph(Surv(tstart, time, status) ~ trt + prior + karno:strata(tgroup),
              data = vet2)
vfit2
#> Call:
#> coxph(formula = Surv(tstart, time, status) ~ trt + prior + karno:strata(tgroup), 
#>     data = vet2)
#> 
#>                                   coef exp(coef)  se(coef)      z        p
#> trt                          -0.010158  0.989894  0.190075 -0.053 0.957381
#> prior                        -0.006735  0.993288  0.020252 -0.333 0.739479
#> karno:strata(tgroup)tgroup=1 -0.052164  0.949173  0.008086 -6.451 1.11e-10
#> karno:strata(tgroup)tgroup=2 -0.040127  0.960668  0.011541 -3.477 0.000507
#> karno:strata(tgroup)tgroup=3 -0.008322  0.991712  0.008725 -0.954 0.340150
#> 
#> Likelihood ratio test=57.13  on 5 df, p=4.759e-11
#> n= 305, number of events= 128
prop2 <- cox.zph(vfit2)
prop2
#>                      chisq df     p
#> trt                  1.217  1 0.270
#> prior                3.898  1 0.048
#> karno:strata(tgroup) 0.235  3 0.972
#> GLOBAL               5.791  5 0.327
plot(prop2[3])

plot(survfit(vfit2, newdata = vet2))

## tt() "log(time + 20)" approach
dtimes <- sort(unique(veteran$time[veteran$status==1]))
vet3 <-survSplit(Surv(time, status) ~.,
                 data = veteran, cut = dtimes,
                 episode = "tgroup", id = "id")
vet3$time_var_karno <- vet3$karno * log(vet3$time+20)
vfit4 <- coxph(Surv(tstart,time,status) ~ trt + prior + karno + time_var_karno, data = vet3)

plotdata <- data.frame( 
  tstart = c(0, dtimes[1:96]),
  time = dtimes,
  trt = rep(0, 97),
  prior = rep(c(0, 5, 10), length.out = 97),
  karno = rep(30, 97),
  tgroup = rep(c(1,2,3), length.out = 97),
  status = rep(0, 97)
)

plotdata$time_var_karno <- plotdata$karno * log(plotdata$time + 20)
plotdata2 <- plotdata
plotdata2$karno <- 60
plotdata2$time_var_karno <- plotdata2$karno * log(plotdata2$time + 20)
plotdata3 <- plotdata
plotdata3$karno <- 90
plotdata3$time_var_karno <- plotdata3$karno * log(plotdata3$time + 20)

plotdata$id <- "1"
plotdata2$id <- "2"
plotdata3$id <- "3"

col_pal <- viridisLite::viridis(n = 4)
plot(survfit(vfit4, newdata = rbind(plotdata,
                                    plotdata2,
                                    plotdata3),
             id = id, se.fit = FALSE),
     bty = "n", xlab = "Days",
     ylab = "Survival probability",
     col = col_pal[1:3])
text(160,0.15,"karno = 30", col = col_pal[1])
text(200,0.4,"karno = 60", col = col_pal[2])
text(240,0.6,"karno = 90", col = col_pal[3])

# overlay
plot(survfit(vfit2, newdata = vet2))
lines(survfit(vfit4, newdata = rbind(plotdata,
                                     plotdata2,
                                     plotdata3),
              id = id, se.fit = FALSE),
      bty = "n", xlab = "Days",
      ylab = "Survival probability",
      col = col_pal[1:3],
      lwd = 3)

在2025-02-08创建

reprexv2.1.1

^{NB：老实说，我不确定这个答案是否为“正确”。它看起来不错，即Time_var_karno校正线可提供原始GGSURVPLOT的非常接近的近似值，但是如果您在期刊上发布此内容，我建议您要求统计教授事先批评它。我还认为值得在https：//stats.stackexchange.com上寻求建议，看看是否有任何“重型击球手”愿意进一步讨论。}