Plotting Survival Curves with Uncertainty Estimates 

 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  and  Andy  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  Cox proportional hazards output (which was developed by Patrick Lam).  Here are two examples of what my surv.plot() function can provide: 

     

 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  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  
			robust=TRUE,cluster="polar",model="coxph",data=coalition) 
			 

 # Set Low and High Quantities of Interest		 
	low  
	high  

 # Simulate for Each  
	s.out1  
	s.out2  

 # Create list output that contains both simulations 
	out  

  
# 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)