I haven't been able to gain access to the Simon and Makuch article mentioned above, but having researched the topic I found:
Steven M Snapinn, Qi Jiang & Boris Iglewicz (2005) Illustrating the Impact of a Time-Varying Covariate With an Extended Kaplan-Meier Estimator, The American Statistician, 59:4, 301-307.
That article proposes a time-dependent Kaplan-Meier plot (KM) by simply updating the cohorts at all event times. It also cites the Simon and Makuch article for proposing a similar idea. The regular KM does not allow this, it only allows a fixed division into groups. The proposed method actually splits survival time according to covariate status - just like one could do when estimating a Cox model with piecewise constant covariates. For the Cox model this is a viable idea, and a standard one. It is, however, more intricate when doing a KM plot. Let me illustrate it with a simulation example.
Let's assume that we have no censoring, but some event (e.g., giving birth) that might or might not occur before time of death. Let's also assume constant hazards for the sake of simplicity. We will also assume that giving birth does not alter the hazard of dying. We will now follow the procedure prescribed in the above article. The article clearly states how this is done in R, simply split your subjects on time of giving birth such that they're constant in your grouping variable. Then use the counting process formulation in the Surv function. In code
library(survival)
library(ggplot2)
n <- 10000
data <- data.frame(id = seq(n),
preg = rexp(n, 1),
death = rexp(n, .5),
enter = 0,
per = NA,
event = 1)
data$exit <- data$death
data0 <- data
data0$exit <- with(data, pmin(preg, death))
data0$per <- 0
data0$event[with(data0, preg < death)] <- 0
data1 <- subset(data, preg < death)
data1$enter <- data1$preg
data1$per <- 1
data <- rbind(data0, data1)
data <- data[order(data$id), ]
Sfit <- survfit(Surv(time = enter, time2 = exit, event = event) ~ per, data = data)
autoplot(Sfit, censSize = 0)$plot
I'm more or less splitting it "by hand". We could use survSplit as well. The procedure actually gives me a very nice estimate.
We get almost identical estimates for the two groups as we should. But actually, my simulation was perhaps a bit unrealistic. Let's say a woman cannot give birth in the first two time units for some reason. This is at least reasonable in your example: there'll be some time between two pregnancies corresponding to the same woman. Making a small addition to the code
data <- data.frame(id = seq(n),
preg = rexp(n, 1) + 2,
death = rexp(n, .5),
enter = 0,
preg = NA,
event = 1)
we get the following plot:
The same thing would be happening with your data. You will not see any 3rd pregnancies for at least some initial period of time, meaning that your estimate will be 1 for that group and that period of time. This is in my opinion a misrepresentation of your data. Consider my simulation. The hazards are identical, but for every time point the per1 estimate is greater than the per0 estimate.
You could consider different remedies to this problem. You propose pasting them together at some point (let the per1-curve start from a certain point on the per0-curve). I like this idea. If I do it on the simulation data, we get:
In our specific case, I think this represents data way better, but I don't know of any published results that support this approach. Heuristically, one can use the argument I presented in another answer:
KM plot with time-varying coefficient
Semi-parametric models such as Cox regression are not easily applied in the presence of interval censoring. In this situation the default choice is to use parametric models.
On the other hand, time-dependent covariates are easily included in a right-censored setting and not with interval censoring. To include time-dependent covariates in an interval censoring scenario few methods have been proposed, and this is probably the best available approach. However, this method (as far as I know) has only been implemented in SAS. The Weibull model is among the possible methods they discuss, so if you are interested in a proportional hazard model you will be fine.
For your analyses I see three options:
- Choose this latter approach based on parametric survival, which would be ideal for your data, and use SAS to fit your models.
- If you prefer to use R, make an assumption on when events are occurring within the interval (normally distributed for example), and do a classical right-censored analysis using Cox regression (in this case you can refer to this R tutorial)
- Seek alternative approaches such as the one you mentioned based on logistic regression.
Best Answer
Survival analysis is performed at event times. What matters at each event time in a Cox proportional hazards regression is the set of current values of covariates for the case that had the event, versus current covariate values for those still at risk who didn't have the event.
A regression coefficient and interaction terms involving a time-varying covariate are thus no different from those for a time-invariant covariate. The regression coefficients and associated hazard ratios are still constant over time. That means that the effect on survival associated with any particular value of the time-varying covariate is constant over time. All that differs is that the values of the time-varying covariates can vary with time.
Once you have a model you certainly can plot predicted survival for members of any particular groups or combinations of covariate values, and you can even allow for time-varying values of covariates in the new data that you provide to the model if you format the data properly. This can be done, for example, with the
survfit.coxph()function in the Rsurvivalpackage. The help page includes the following critical warning, however:Things are more complicated if you are considering time-varying coefficients for covariates. The time-dependent vignette provided with the
survivalpackage is an introduction to issues with respect to time-dependent covariates and coefficients.