Interpret survival curve for multiple-event Cox proportional hazard model I am using the bladder2 dataset, which is setup in "conditional" long format:
data(bladder)

> head(bladder2)
  id rx number size start stop event enum
1  1  1      1    3     0    1     0    1
2  2  1      2    1     0    4     0    1
3  3  1      1    1     0    7     0    1
4  4  1      5    1     0   10     0    1
5  5  1      4    1     0    6     1    1
6  5  1      4    1     6   10     0    2

I then build a simple coxph model with stratification of the recurring event (denoted by enum):
library(survminer)

cox.model.cond = coxph(
  Surv(start, stop, event) ~ rx + number + size + strata(enum),
  data=bladder2
)

Then I plot the resultant survival curves:
ggsurvplot(survfit(cox.model.cond), data=bladder2)


How should I interpret these curves? 
 Specifically:


*

*Even though events are sequential, does the plot indicate the average patient is at risk of 4x events at (almost) any time t?

*If an average patient were currently at time t=10 and they had not experienced any events, what would we say their probabilities of experiencing the events are in the next instant? the next 10 'steps' out?

*Same question as above but what about a patient who has already experienced 1x event?

*The lines for enum=1 and enum=2 cross.  How should this be interpreted, that the chance of 1x event occurring becomes less probable than the chance of the second event occurring?  That seems silly; I'm sure I have it wrong :)

 A: I'm not sure that we are reading the data situation similarly. As I understand it ... all of these patients have had at least one diagnosed episode of bladder cancer. The coxph-model you created with a strata()-constructed term yields this result using its print output:
cox.model.cond
#-----------
Call:
coxph(formula = Surv(start, stop, event) ~ rx + number + size + 
    strata(enum), data = bladder2)

            coef exp(coef)  se(coef)      z      p
rx     -0.333489  0.716420  0.216168 -1.543 0.1229
number  0.119617  1.127065  0.053338  2.243 0.0249
size   -0.008495  0.991541  0.072762 -0.117 0.9071

Likelihood ratio test=6.51  on 3 df, p=0.08928
n= 178, number of events= 112 

Notice that there are no results for the individual strata. That's because the strata function basically separates the analysis into groups with separate baseline curves, but makes the assumption that the other covariates will have the same effects across strata relative to their separate baseline curves. This notion of homogeneity of effects across these strata seems to me  suspect and requiring great care before acceptance.
To the point of comparing different approaches to modeling recurrent events (also with bladder cancer data, but 10 times as much) see the paper, "Prediction of Multiple Recurrent Events: A Comparison of Extended Cox Models in Bladder Cancer", by Hilde Smedinga in American Journal of Epidemiology.
To your separate questions:

*

*Yes, I think that these results indicate an ongoing risk for all comers, However, I think it also indicates increasing risks for persons after a second or third episode.


*(& 3) I'm going to assume that by an "event" you mean a second or third cancer. And an average patient would be one with  ... average size? average treatment? average "number" of sites? Or do you mean simply that we should consider all patients in the stratum? I'm going to assume the second option, since I think the answers to the first set of questions raise too many messy imponderables.
To get the risk estimate for a member of the stratum you would take the decrement in survival divided by the time over which the decrement occurred and divide that by the starting survival fraction.
$$
\frac{ -\frac{d(S(t))}{dt} }{ S(t) }
$$
Generally one would smooth the survival fraction as a function of time ( S(t) ), and take the slope of the smoothed estimate as the numerator and the value of the smoothed estimate as the denominator. There are packages that will do that for you, and the muhaz-package for which I am the maintainer is one such. It would require data in the form of (time, status) rather than (start, stop, status). There is a nice review of several packages for hazard estimation in Hagar & Dukic, Comparison of hazard rate estimation in R.


*The crossing of the lines indicates that they could not be properly analyzed with a proportional hazards model that had a simple enum covariate. It is in point of fact the reason that a strata analysis needed to be done.

