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)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
