How do I interpret the HRs for each time interval computed by survSplit in R?

I’m testing the association between an exposure and outcome. I have performed a Cox proportional hazard model in R and estimated a HR. Both the exposure and the outcome are binary variables, where 1 = present and 0 = absent. I adjust for age (continuous covariate) and stratify by sex (binary variable, F or M). I calculate standard errors by including robust=T. My Cox model in R looks like this:

cox_model <- coxph(Surv(start, end, outcome)~exposure+age+strata(sex),
robust =T, data=my_data)


By examining the smoothed scaled Schoenfeld residual plot like this:

ph_assump_test <- cox.zph(cox_model)
ph_assump_test
plot(ph_assump_test, resid=FALSE)


I clearly see a violation of the proportional hazard assumption on the residual plot (x-axis shows event time in years and y-axis shows Beta(t) coefficients for the exposure). The cox.zph() test also indicates the violation significantly. But the association between the exposure and outcome is significant p<0.05. So having the exposure, gives you a higher risk of getting the outcome (HR above 1).

So I choose to partition the time on the x-axis using this example from Therneau et al. (section 4.1 Step functions) which used survSplit see the example here: survival vignettes. In my situation, I see the violation occurs at year 18 from the residual plot. Therefore, I make the cut at 18. This is how I compute it in R:

split_time <- survSplit(Surv(start, end, outcome) ~ ., data= my_data, cut=c(18),
episode= "tgroup", id="id")

cox_model_split <- coxph(Surv(start, end, outcome) ~ exposure:strata(tgroup)+
age+strata(sex), robust=T, data=split_time)


Then I compute a proportional hazard test and plot a residual plot:

ph_assump_test <- cox.zph(cox_model, terms=FALSE)
ph_assump_test
plot(ph_assump_test[2], resid=FALSE)    # the residual plot before year 18
plot(ph_assump_test[3], resid=FALSE)    # the residual plot after year 18


After visualizing the proportional hazard test by the residual plot, and observing the p-value, the HRs for both time intervals do not violate the proportional hazard assumption. And I still see a significant association (p<0.05) between the exposure and outcome in each time intervals (HR above 1 in each time interval).

My questions are as follows:

1. How do I interpret the HR for the two time intervals? I know how to interpret the time-averaged HR, which considers the entire follow-up period from year 0 to 50. But how should I interpret the two HRs below 18 and above 18?
2. How does the cut in survSplit work theoretically? I was told that all observations in my sample would still contribute to the risk calculation in both time intervals (maybe this is wrong?). I have a hard time understanding this because individuals in my sample who experienced the event in the first time interval (0-18) are no longer at risk in the second time interval (18-50). So how can they still contribute to the risk in the second time interval?
3. Do I still need to include robust=T in the cox_model_split? As the time splitting cox model did not violate the proportional hazard assumption.

• Could you please help me with these questions, if you can? @EdM ? Commented Apr 29 at 18:17
• For future reference, the callout with @ only works for a set of comments in which the individual is already engaged.
– EdM
Commented Apr 29 at 19:47
• Thanks for letting me know. I had hoped that I could tag you when I needed help. You're the reason I understand what Cox-ph regression is all about Commented Apr 29 at 20:13

Question 1. If you have two time intervals separated out by survSplit() at time = 18, then you effectively have two separate models, one for each of the time intervals. The HR for each time interval only holds within that particular time interval. In your example that means a discontinuity in the HR at time = 18, which some find to be disturbing. The simplicity in describing the results might nevertheless be worth it.

Question 2. Your sense is correct; you might want to get clarification from your source who seemed to suggest otherwise. In a Cox model, hazards are modeled from data at event times, based on the current covariate values of those still at risk. If there can be no more than 1 event per individual, an individual experiencing an event ceases to contribute information after the event time. That's true whether or not the data have been re-formatted by survSplit(). In the vignette that you cite, look at the patient with id = 1 who died within the first time interval of that example; that patient's information was not included in the analyses of subsequent intervals.

Question 3. A robust variance estimator can help with problems beyond violations of proportional hazards. See this page.

• For Q1 @EdM - I was considering reporting the failed averaged time hazard ratio for the whole follow-up, as well as the two HRs for each time interval. From there, I was thinking of being very cautious about how to interpret the HRs. I would like to hear your personal opinion on this. Would you suggest sticking only with the averaged HR in general? Because it shows very clear and interesting results by splitting the data as well. Also, what are your thoughts on using survSplit() for this purpose? Do you have better suggestions? Commented Apr 29 at 20:09
• @Aphi11 you could show the plot of smoothed, scaled Schoenfeld residuals over time for that variable. Say something like "the pattern of residuals could be approximated by 2 separate constant hazard ratios for times before and after 18." Overlay the plot with the corresponding horizontal lines at the HR values for the 2 time periods. If there's some subject-matter knowledge that supports such a 2-stage hazard, explain that too.
– EdM
Commented Apr 30 at 13:57
• Hi @Edm do you have time to look at this post? I hope you can help! stats.stackexchange.com/questions/647967/… Commented May 25 at 15:39
• Hi again @EdM do you know why time 18 occurs twice? When visualising the residual plot. I see the time from 0 to 18 and then 18 to 50. what is the purpose of including time 18 twice in the estimates for each interval? Commented May 30 at 14:53
• @DeviSita it's hard to say without seeing the plot. You do fundamentally have two different models, one from 0 to 18 and another from 18 to 50, so I suspect that's why 18 is included twice: once as the end point and once as the start point.
– EdM
Commented May 31 at 2:20