我正在使用survival
包中的survplot功能。具有置信区间的生存图很好地生成,但现在我遇到了将图转换为累积发生率曲线的问题。曲线本身正确生成,但在使用conf = "bars"
函数时,置信区间仍保留在生存设置中。然而,"bands"
和"diffbands"
工作正常。
我将为您带来一个简单的可重复示例:
library(survival)
library(rms)
Data <- data.frame("time" = sample(1:500), "death" = sample(c(TRUE, FALSE), 500, replace = TRUE))
Data$SurvObj <- with(Data, Surv(Data$time, Data$death == 1))
km.as.one <- npsurv(SurvObj ~ 1, data = Data, conf.type = "log-log")
这是问题所在:
survplot(km.as.one, fun=function(y) 1 - y, conf = "bars")
但是,这些工作正常:
survplot(km.as.one, conf = "bars")
survplot(km.as.one, fun=function(y) 1 - y, conf = "bands")
这个问题有什么可能的解决方案吗?我想ggplot2
包可以正确地做到这一点,但是我已经使用survival
包生成了相当多的生存图,所以现在更改包会导致很多额外的工作。
如果你厌倦了等待rms的mod到CRAN,你可以只考虑所考虑的值:
km.as.one$surv <- 1-km.as.one$surv
km.as.one$lower <- 1-km.as.one$lower
km.as.one$upper <- 1-km.as.one$upper
survplot(km.as.one, fun=function(y) y, conf = "bars")
嘿,刚注意到(0,1)处的“点”。
如果不使用对象中的三个向量的黑客,而是使用fun
参数:
survplot(km.as.one, fun=function(y) 1-y, conf = "bars")
...一个人发现线路被转换但点和错误栏没有。如果修改代码:
getAnywhere(survplot.npsurv)
....并复制到代码编辑器,然后将fun
的值应用于该相当长的函数的主体末端附近的点和误差栏:
surv.plot.npsurv <- <- function (fit, xlim, ylim, xlab, ylab, time.inc, state = NULL,
# lines 2-289 of original code suppressed
ss <- fun(v$surv[j]) # lines 290-292 when doing this in Rstudio's code editor.
lower <- fun(v$lower[j])
upper <- fun(v$upper[j])
# rest of original code
}
我知道你说使用ggplot2
包可能不是最好的开关,但它确实产生了一个答案,最小的额外编码:
library(survival)
library(rms)
library(broom)
library(dplyr)
set.seed(123)
Data <- data.frame("time" = sample(1:500), "death" = sample(c(FALSE, TRUE), 500, replace = TRUE))
Data$SurvObj <- with(Data, Surv(time, death == 1))
km.as.one <- npsurv(SurvObj ~ 1, data = Data, conf.type = "log-log")
survplot(km.as.one, fun=function(y) 1 - y, conf = "bars")
tidydf <- tidy(km.as.one) %>%
mutate(estimate = 1- estimate,
#invert estimates
conf.low = 1- conf.low,
conf.high = 1- conf.high,
#get points and CIs at specific timepoints
pointest = ifelse(row_number()%%50 != 0, NA,estimate),
confestlow = ifelse(row_number()%%50 != 0, NA,conf.high),
confesthigh = ifelse(row_number()%%50 != 0, NA,conf.low))
#plot
ggplot(tidydf)+
geom_line(aes(x=time, y = estimate,group = 1))+
geom_point(aes(x=time, y = pointest))+
geom_errorbar(aes(x=time, ymin = confestlow, ymax = confesthigh))