Utilise Surv object in ggplot or lattice

analysis, ggplot2, graph, lattice, r

Solution

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.

Problem

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

Original source