R packages (or SAS code) to produce two simultaneous Kaplan-Meier curves? There's a way to do survival analysis of two (or more I suppose) mutually exclusive competing risks as a mixture of two different survival curves. Something like what you see in A.C. Ghani et al. Methods for Estimating the Case Fatality Ratio for a Novel, Emerging Infectious Disease. American Journal of Epidemiology (2005) Vol. 162, No. 5
What I'm looking for is a package that would help produce something like this figure:

Where the survival curve of one outcome, and 1-the survival curve of the other outcome will eventually meet at a particular point that is the mixture of the two outcomes.
 A: In R, a survfit.object---returned by survfit()---stores a fitted survival curve. In particular, this object contains the time points at which the curve has a step and the ordinates at those points. You can therefore construct the survival function, $t\mapsto \hat{S}(t)$, by constant interpolation. Here is the way I would do this:
km <- summary(survfit(Surv(time, event) ~ 1, data=data))
S <- approxfun(km$time, km$surv,
               method="constant", f=0, yleft=1, rule=2)

Now, S can be used as any user-defined function in R: in particular, you can evaluate S(t) at any time t, you can make plots using plot(), and you can superimpose two K-M curves on the same graph using lines(), ...
Hope this helps!
A: What you are asking for is a simultaneous plot of the survival function for one process and the cumulative incidence function (= 1- S(t)) for the competing process. The 'cmprsk' R package should be able to do the plots, but since the usual mode is to display both process as the cumulative incidence, you will need to do some work to transform the data so that one is S(t) and the other is H(t).
A: Wouldn't it be good enough if you could plot two curves using  par(new=T)?
plot(survfit(KMfit1 ~ 1),main="Kaplan-Meier estimate with 95% confidence bounds",xlab="time", ylab="survival function",col="red",xlim=c(0,70))

par(new=T)

plot(survfit(KMfit2 ~ 1),col="green",xlim=c(0,70))


