Utilise Surv object in ggplot or lattice

Question

Anyone knows how to take advantage of ggplot or lattice in doing survival analysis? It would be nice to do a trellis or facet-like survival graphs.

---

So in the end I played around and sort of found a solution for a Kaplan-Meier plot. I apologize for the messy code in taking the list elements into a dataframe, but I couldnt figure out another way.

Note: It only works with two levels of strata. If anyone know how I can use x<-length(stratum) to do this please let me know (in Stata I could append to a macro-unsure how this works in R).

ggkm<-function(time,event,stratum) {

    m2s<-Surv(time,as.numeric(event))

    fit <- survfit(m2s ~ stratum)

    f$time <- fit$time

    f$surv <- fit$surv

    f$strata <- c(rep(names(fit$strata[1]),fit$strata[1]),
            rep(names(fit$strata[2]),fit$strata[2])) 

    f$upper <- fit$upper

    f$lower <- fit$lower

    r <- ggplot (f, aes(x=time, y=surv, fill=strata, group=strata))
        +geom_line()+geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.3)

    return(r)
}

Answer 1

I have been using the following code in lattice. The first function draws KM-curves for one group and would typically be used as the panel.group function, while the second adds the log-rank test p-value for the entire panel:

 km.panel <- function(x,y,type,mark.time=T,...){
     na.part <- is.na(x)|is.na(y)
     x <- x[!na.part]
     y <- y[!na.part]
     if (length(x)==0) return()
     fit <- survfit(Surv(x,y)~1)
     if (mark.time){
       cens <- which(fit$time %in% x[y==0])
       panel.xyplot(fit$time[cens], fit$surv[cens], type="p",...)
      }
     panel.xyplot(c(0,fit$time), c(1,fit$surv),type="s",...)
}

logrank.panel <- function(x,y,subscripts,groups,...){
    lr <-  survdiff(Surv(x,y)~groups[subscripts])
    otmp <- lr$obs
    etmp <- lr$exp
    df <- (sum(1 * (etmp > 0))) - 1
    p <- 1 - pchisq(lr$chisq, df)
    p.text <- paste("p=", signif(p, 2))
    grid.text(p.text, 0.95, 0.05, just=c("right","bottom"))
    panel.superpose(x=x,y=y,subscripts=subscripts,groups=groups,...)
}

The censoring indicator has to be 0-1 for this code to work. The usage would be along the following lines:

library(survival)
library(lattice)
library(grid)
data(colon)  #built-in example data set
xyplot(status~time, data=colon, groups=rx, panel.groups=km.panel, panel=logrank.panel)

If you just use 'panel=panel.superpose' then you won't get the p-value.

Answer 2

I started out following almost exactly the approach you use in your updated answer. But the thing that's irritating about the survfit is that it only marks the changes, not each tick - e.g., it will give you 0 - 100%, 3 - 88% instead of 0 - 100%, 1 - 100%, 2 - 100%, 3 - 88%. If you feed that into ggplot, your lines will slope from 0 to 3, rather than remaining flat and dropping straight down at 3. That might be fine depending on your application and assumptions, but it's not the classic KM plot. This is how I handled the varying numbers of strata:

groupvec <- c()
for(i in seq_along(x$strata)){
    groupvec <- append(groupvec, rep(x = names(x$strata[i]), times = x$strata[i]))
}
f$strata <- groupvec

For what it's worth, this is how I ended up doing it - but this isn't really a KM plot, either, because I'm not calculating out the KM estimate per se (although I have no censoring, so this is equivalent... I believe).

survcurv <- function(surv.time, group = NA) {
    #Must be able to coerce surv.time and group to vectors
    if(!is.vector(as.vector(surv.time)) | !is.vector(as.vector(group))) {stop("surv.time and group must be coercible to vectors.")}

    #Make sure that the surv.time is numeric
    if(!is.numeric(surv.time)) {stop("Survival times must be numeric.")}

    #Group can be just about anything, but must be the same length as surv.time
    if(length(surv.time) != length(group)) {stop("The vectors passed to the surv.time and group arguments must be of equal length.")}

    #What is the maximum number of ticks recorded?
    max.time <- max(surv.time)  

    #What is the number of groups in the data?
    n.groups <- length(unique(group))

    #Use the number of ticks (plus one for t = 0) times the number of groups to
    #create an empty skeleton of the results.
    curves <- data.frame(tick = rep(0:max.time, n.groups), group = NA, surv.prop = NA)

    #Add the group names - R will reuse the vector so that equal numbers of rows
    #are labeled with each group.
    curves$group <- unique(group)

    #For each row, calculate the number of survivors in group[i] at tick[i]
    for(i in seq_len(nrow(curves))){
      curves$surv.prop[i] <- sum(surv.time[group %in% curves$group[i]] > curves$tick[i]) /
          length(surv.time[group %in% curves$group[i]])
    }

  #Return the results, ordered by group and tick - easier for humans to read.
  return(curves[order(curves$group, curves$tick), ])   

}