My line of thought goes like this. As per the AG model, each individual is represented by a counting process $N_i(t)$ with intensity $\lambda_i(t)$ which can be written as $\lambda_i(t) = \lambda_0(t) \exp(\beta' x_i(t))$. By this it is implicit that only the "current" value of $x_i$ matters for the intensity. The counting process grows when an individual has an event. I also assume independent censoring (actually stopping of the process). Note that these assumptions are implicit from the way that you fitted the model.
The probability of no events in $(a,b)$ (improperly, a "survival") is then
$$
S_i(t|s) = \exp \left( -\int_s^t \lambda_0(u) \exp(\beta'x_i(u)) du \right)
$$
Note that the probability of no events does not directly relate to the probability of one event, as the complementary event to "no events in a period" is "at least one event during a period".
Then say you are interested in $S(t|s_0)$. Then the idea would be to fix $x_i = x_i(s_0)$ (i.e. assume the covariates don't change after $s_0$), and we would get
$$
S_i(t|s_0) = S_i(t) / S_i(s_0)
$$
Since the part before $s_0$ cancels out, together with my assumption that only the current value of $x_i$ matters for the intensity, means that this is equal to
$$
S_i(t|s_i) = \exp \left( -\exp(\beta'x_i(s_0)) \int_{s_0}^t \lambda_0(u)du \right)
$$
Then in R you can do it like that (I give an example with a data set):
library(frailtypack)
data(readmission)
mod1 <- coxph(Surv(t.start, t.stop, event) ~ sex + charlson + cluster(id), data = readmission)
# set the covariate values at s_0
mycov <- data.frame(sex = "Female", charlson = "1-2")
sf <- survfit(mod1, newdata = mycov)
# with different s_0
par(mfrow=c(2,2))
time <- c(300, 500, 700, 1000)
for(i in time) {
pos <- which.max(sf$time[sf$time <= i])
S_s0 <- sf$surv[pos]
with(sf, plot(time[pos:length(time)], surv[pos:length(surv)] / S_s0), type = "l")
}
Here you get the plots of the survival curves, which correspond to the values
$(t, S(t|s_0))$ for $t\geq s_0$.
The other comments that I would give on this are the following: it is difficult to talk about the distribution of the next event. This is because in the AG formulation the time since previous event does not play any role. In other words, if you would like to take that into account, more complicated stochastic models should be used, where for example you include the "previous number of events" as a time-dependent covariate. This does complicate things a lot and the interpretation of the estimated quantities is most likely very difficult.
The second comment I have is about the nature of the time dependent covariates. Mostly, the AG works nicely with "external" covariates, such as air pollution, or something which is not directly measured on the subject ("external" of the recurrent event process). This is mostly because the first expression I wrote here is the probability of no events during $(a,b)$ relies on the assumption that the number of events in any given interval is Poisson distributed. This is true if the covariates are external. A discussion on this can be found in several textbooks, for example in Cook & Lawless at section 2.5. If your time-dependent covariates do depend on the recurrent event process, then it should be modeled jointly with the recurrent events process.
You need to go beyond the standard output from coxph() and analyze the data in a way that accomplishes two things. For your first question, what you presumably want is a more general idea whether any of your treatments affect developmental times significantly, including both their direct effects and their interactions with shipments and types of events. For your second question, you need flexibility to examine particular combinations (contrasts) of factor levels other than the comparisons of each individually against the reference level (as you correctly take the coxph() output to represent.)
Although it's possible to do this yourself by hand, the rms package in R provides these capabilities directly. There is a bit of learning curve in terms of pre-processing the data (datadist is important) and using the cph() function with appropriate parameter values to allow useful downstream processing, but it's well worth it. Then the anova() function applied to the cph() output provides information about whether the treatments differ among themselves (as main effects or including interaction terms) and the contrast() function allows for great flexibility in terms of examining pre-specified combinations of factor levels.
The rms package also provides ways to calibrate and validate your model. You need to know that your data adequately meet the proportional hazards assumption before you report results that depend on the assumption, and you want to make sure that your results aren't too dependent on the particular data sample you have.
If you are publishing in a journal that provides for supplemental material outside of the main text, showing the anova tables and specific contrast results in the main text and putting more detailed results in the supplement might work. In clinical survival analysis, usual practice is to present hazard ratios (exp(coef) in the output) along with their confidence intervals. As you might be more interested in whether the treatments delay rather than accelerate the developmental steps, you might consider using exp(-coef) instead as a measure of delay and recalculate the confidence intervals from -coef and robust se.
In this particular application, you will of course have to deal with the biological problem of the differences among shipments being the greatest of all influences. That makes it very hard to go from your specific data to generalizations about these types of insects and treatments. A more than 5-fold difference in rates of developmental progression between 2 shipments of insects (ship2 versus ship1) is pretty remarkable. My guess is that, in the background of those shipment differences and their striking interactions with specific developmental stages, any effects of the treatments themselves will be lost.
Best Answer
Therneau and Grambsch discuss recurrent-event models in Chapter 8.
For the Andersen-Gill model: "This method is the simplest to visualize and set up, but makes the strongest assumptions" (page 185); "This model is ideally suited to the situation of mutual independence of the observations within a subject" (page 186).
For the PWP model: "It assumes that a subject cannot be at risk for event 2 until event 1 occurs; in general, a subject is not at risk for the $k^{th}$ event until he/she has experienced event $k - 1$. To accomplish this... each event is assigned to a separate stratum... The use of time-dependent strata means that the underlying intensity function may vary from event to event, unlike the AG model, which assumes that all events are identical" (page 187).
What you found is what you might expect if all events within an individual are not identical, so that the stronger assumption underlying the Andersen-Gill model isn't met.