[R] Survival (KM/Cox) plots with lattice?

[Dieter Menne]

Deepayan Sarkar <deepayan.sarkar <at> gmail.com> writes:

> > I am looking for a lattice-panel for survival (KM/Cox) plots. I know it's
> > not standard, but maybe someone has already tried?
> 
> There are some half-formed ideas in
> 
> http://dsarkar.fhcrc.org/talks/extendingLattice.pdf
> 
> but nothing packaged (mostly because I'm not all that familiar with
> survival analysis). If you have some well-defined use cases, I would
> be happy to collaborate on something generally useful.
> 

Deepayan,

Here my attempts. The code is not very generic, cox does not work yet, and I
don't believe the syntax is good. And, there is an environment marked !!! below,
which I could not get around and replaced it by a hardcoded grouping.

But nevertheless, I think it already looks better than standard graphics

Dieter

#------------------------------------------------------------------------------
library(survival)
library(lattice)

ovarian$resfak = factor(ovarian$resid.ds)
levels(ovarian$resfak) =c("resid.no","resid.yes")
ovarian$rxfak = factor(ovarian$rx)
levels(ovarian$rxfak) =c("rx.no","rx.yes")


xyplot.survfit = function(x,form,groups=NULL,...) {
  st = names(x$strata)
  # treat strata like other grouping variable
  st = sub("strata\\(.*)=","",st)
  faks = strsplit(st,"[=,]")[]
  nc = length(faks[[1]])
  # bug (?) in survfit; sometimes has spaces
  faks = gsub(" ","", unlist(faks))
  nam = faks[seq(1,nc,by=2)]
  gr = data.frame(
    matrix(faks[seq(2,length(faks),by=2)],ncol=nc %/%2,byrow=TRUE)) 
  names(gr)= nam          
  dfr = with(x,data.frame(time,n.risk,n.event,surv,std.err,upper,lower,
    gr[rep(1:nrow(gr),x$strata),]))
  # !!! Some env-problem: groups = groups does not work here
  if (!is.null(groups)) {
    xyplot(form, data=dfr, type="s", groups=rxfak, #groups=groups,
      cens = dfr$n.event==0,
      panel = function(x,y,subscripts,groups,cens,...){
        panel.superpose(x,y,subscripts,groups,...)
        ce = cens[subscripts]
        panel.superpose(x[ce],y[ce],ce,groups,type="p")
      })
  }  else {
    xyplot(form,data=dfr,type="s",cens = dfr$n.event==0, subscripts=TRUE,
      panel = function(x,y,cens,subscripts,...){
        panel.xyplot(x,y,...)
        ce = cens[subscripts]
        # dirty. Should pass ... too to set pwd etc.  
        panel.xyplot(x[ce],y[ce],type="p")
      }
      )
  }
}

svfit = survfit( Surv(futime,fustat)~resfak+rxfak,data=ovarian)
cph = coxph( Surv(futime,fustat)~resfak+rxfak,data=ovarian)
#coxfit = survfit(cph,newdata=... )# Not tried yet

# The formula in xyplot is redundant and error prone, since the first term should
# always be surv~time.
# The correct formula could be inferred from svfit$call, but I would be 
# important to distinguish between groups and panels in plotting 
# form = ~|resfak, groups=rxfak # looks minimal, but not really nice
# panels=resfak, groups=rxfak # collides with panel=

xyplot(svfit,surv~time|resfak, groups="rxfak") # looks ok, bad syntax see !!!
xyplot(svfit,surv~time|resfak+rxfak) # looks ok
#  would look ok if the !!! problem above were solved
#xyplot(svfit,surv~time|rxfak, groups="resfak") 

#xyplot(svfit,surv~time)  # plot is not valid, should not be possible

# Does not work yet
#xyplot(coxfit,surv~time|resfak, groups="rxfak")
#xyplot(coxfit,surv~time|resfak+rxfak)