#  Plotting Survival Curves with Uncertainty Estimates 

 





May 07, 2008

 

 

One of the pesky things I've found in my (limited) experience with survival analysis is that it's almost impossible to plot several survival curves in the same space *and* include measures of uncertainty without the entire plot becoming incomprehensible. So, to build on the great R discussions [Ellie](http://www.iq.harvard.edu/blog/sss/archives/2008/05/tuesday_tips_tr.shtml) and [Andy](http://www.iq.harvard.edu/blog/sss/archives/2008/04/google_charts_f.shtml#more) have provided in recent blog posts, I'd like to offer an extension of my own. I've created a fairly flexible function that allows one to plot several survival curves along with estimation uncertainty from [Zelig's](http://gking.harvard.edu/zelig/) Cox proportional hazards output (which was developed by Patrick Lam). Here are two examples of what my surv.plot() function can provide:

[](/file_url/165)

Hopefully this will be of some interest to a few readers. More details and example code below.

Here is the syntax for the command:  
s.out: Simulated output from Zelig for each curve organized as a list()  
duration: Surival time  
censor: Censoring indicator  
type: Display type for confidence bands. The default is "line" but "poly" is also supported (to create the shaded region in the right plot above).  
plotcensor: Creates rug() plot indicating censoring times (Default is TRUE)  
plottimes: plots a point for each survival time in the step function (Default is TRUE)  
int: Desired uncertainty interval (Default is c(0.025,0.975) which corresponds to a 95% interval)

[Here's](/files/sss_blog/files/survplot.R) the plot.surv() source code, and below I've copied the R code I used to create the plots above:

library(Zelig)  
data(coalition)

\# Fit the Cox Model  
z.out1 - zelig(Surv(duration,ciep12)~invest+numst2+crisis,  
 robust=TRUE,cluster="polar",model="coxph",data=coalition)

\# Set Low and High Quantities of Interest   
 low - setx(z.out1,numst2=0)  
 high - setx(z.out1,numst2=1)

\# Simulate for Each   
 s.out1 - sim(z.out1,x=low)  
 s.out2 - sim(z.out1,x=high)

\# Create list output that contains both simulations  
 out - list(s.out1,s.out2)

  
\# Plot the results  
 par(mfrow=c(1,2))  
 surv.plot(s.out=out,duration=coalition$duration,censor=coalition$ciep12,type="line",plottimes=TRUE)  
 surv.plot(s.out=out,duration=coalition$duration,censor=coalition$ciep12,type="poly",plotcensor=TRUE)

Posted by [John Graves](http://www.iq.harvard.edu/blog/sss/archives/author/john-graves/) at May 7, 2008 11:48 AM



 

 

 



 

 

 Share on:- [     Facebook ](#)
- [     Twitter ](#)
- [     Linkedin ](#)