CODEWITHSUNDEEP
rsurvival-analysisconfidence-interval
Joonas Lehto
Joonas Lehto

Reputation: 23

R survplot function with confidence bars not working correctly

I'm using survplot function from the survival package. The survival plots with confidence intervals are producing nicely, but now I've faced a problem with transforming the plots to cumulative incidence curves. The curve itself is producing correctly, but when using the conf = "bars" function the confidence intervals remain in the survival setting. The "bands" and "diffbands", however, are working correctly.

I'll bring you a simple reproducible example:

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

Here's the problem:

survplot(km.as.one, fun=function(y) 1 - y, conf = "bars")

However, these are working correctly:

survplot(km.as.one, conf = "bars")
survplot(km.as.one, fun=function(y) 1 - y, conf = "bands")

Is there any possible solutions for this problem? I guess the ggplot2 package would do this correctly, but I've produced quite a number of survival plots already with the survival package, so changing the package now would cause a lot of extra work.

Upvotes: 2

Views: 847

Answers (2)

IRTFM
IRTFM

Reputation: 263342

If you are tired of waiting for teh mods to rms to make it to CRAN, you can just take the converse of the values under consideration:

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

enter image description here

Heh, just noticed the "dot" at (0,1).

If instead of using the hacking of the three vectors in the object one instead uses the fun parameter:

survplot(km.as.one, fun=function(y) 1-y, conf = "bars")

... one discovers that the line gets transformed but the points and errorbars do not. If you modify the code:

 getAnywhere(survplot.npsurv)

.... and copy to a code editor, then apply the value of fun to the points and error bars near the end of the body of that rather long function:

    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
     }

enter image description here

Upvotes: 1

Mike
Mike

Reputation: 4370

I know you said using the ggplot2 package might not be the best switch but it does produce an answer with minimal extra coding:

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

Is this the plot you are looking for? enter image description here

Upvotes: 1

Related Questions