Interpretation of interaction between variable and time in cox regression

Question

I was running a model to check the survival rate among different treatments. Basically, the treatment has four groups. I firstly plot a K-M curve and noticed clear interactions among the four treatment groups over the time. So, in the cox regression model, I put the interaction between treatment and time in the model: such as treat1 * log(time) treat2 * log(time) treat3*log(time). And treatment 1 was used as the reference group. I obtained the result like belows:

variable	Estimate	Hazard ratio	P value
treatment 2	0.29	1.36	<.0001
treatment 3	-1.20	0.35	<.0001
treatment 4	-1.41	0.38	<.0001
treat2*time	-0.25	0.84	<.0001
treat3*time	0.37	1.21	<.0001
treat4*time	0.28	1.25	<.0001

How can I interpret these results? I thought I could explain it like this, for example, treatment 2 increases the risk of death by 36% in the early stage, and reduces the risk of death by 16% in the late stage. Only treatment 2 seems consistent with the K-M plot, but for the rest of the two treatments, when compared with treatment 1, they don't seem to reduce the risk of death so much in the early stage, and they don't seem to increase the risk of death so much comparing to treatment 1 in the late stage. But if I explain the result like this: treatment 2 increases the risk of death by 36% before interaction, and also increases the risk of death by 14%(1.36*0.74=1.14), it totally wrong.

So I am wondering if there is an appropriate way to explain the interaction results? Or should I run a model without the interaction, and compare it with the interaction model? Thank you!

And below is the K-M plot, the K-M plot seems different from the Cox-regression estimate. For example, blue line is treatment 2, and in the KM plot, it seems have better survival rate than red. But in the Cox table, it increased the risk of death by 36%, and after interaction, the risk of death also increased 14%(1.36*0.84). Is it normal or something wrong with the analysis?

Answer 1

Unless you took special precautions to extend your data set to include separate rows for each individual at each time the individual was at risk (which your coding doesn't seem to have done), your attempt at creating an interaction between treatment and time isn't correct.

If the time value is the event or censoring time for a subject, then as the R survival time dependence vignette says your interaction term is just a:

time-static value for each subject based on their value for the covariate time...This variable most definitely breaks the rule about not looking into the future, and one would quickly find the circularity: large values of time appear to predict long survival because long survival leads to large values for time.

In R the tt() function allows you to model these types of interactions with time correctly. Otherwise there really is nothing that can be interpreted meaningfully.

In response to comments and updated question

You need to decide what you want to study and report. As the Kaplan-Meier plots clearly cross, the proportional hazard (PH) assumption with respect to treatment (as a single predictor) does not hold. There is no single answer to which group has "better survival." The blue group might have the best survival at late times, but it seems to have the worst survival at early times. The overall survival curves might be enough to report on their own.

The tt() function and a related method for SAS deal with the time-varying coefficient problem by translating it to a time-varying covariate problem. At each event time, a new time-varying covariate is generated for each individual at risk, equal to the product of the individual's covariate value and the (function of) that event time value (not the individual's own event/censoring time, which is what your code did). That becomes unwieldy and fills up a computer when you have 200,000 individuals: at early times you have hundreds of thousands of new data points generated for each of tens to hundreds of thousands of event times.

If there is reason to think the there are different "early" or "middle" or "late" effects of the treatments on survival, you might make the analysis feasible by using step functions instead of continuous functions of time for the coefficients. That's illustrated in Section 4.1 of the R survival time dependence vignette. Instead of having to duplicate all cases at risk at each event time, you only have to do that for each time period that you choose for the step functions.

Instead of using treatment as the sole predictor in your Cox model, you might find better results if you included all potentially outcome-associated covariates in your model. You certainly have enough data to allow for very flexible modeling, and your use of propensity score weighting (as stated in a comment) indicates that you have other covariates available. If your propensity score model isn't well specified, then you might make up for that with a "doubly robust" model that also includes outcome-associated covariates in the Cox model.

Also, instead of trying to fit a particular function form for the coefficients as a function of time, you might just let smoothed scaled Schoenfeld residual plots illustrate how the association of each treatment with outcome changes over time. With this large a data set you almost certainly will find some violation of PH; those plots can let you and your audience know how big those violations are, over what time courses.