I am trying to conduct a cox proportional hazard, after I checked the assumptions, it turns out that interaction exists. So, I have to put an interaction into the model. But the time varying covariate is polytomous variable with four categories.

The treatment is time varying covariate with four categories. And I get the results like this:

Parameter      DF      Estimate       Error    Chi-Square    Pr > ChiSq       Ratio    

 treatment 2     1      -0.33139       0.02126     1183.5006        <.0001   0.386    
 treatment 3     1      -0.17445       0.02588       45.4387        <.0001   0.745    
 treatment 4     1      -0.69876       0.03591      278.0054        <.0001   0.646    
 treat_ti        1       0.01096     0.0001625      145.9549         <.0001   1.023

I am not sure if I put the interaction term wrong because the treatment is 4-category variable, but there is only one interaction term which makes the interpretation really hard. When I saw the K-M plot, only two treatments had the interaction. So my question is, did I do it correctly? If I did it correctly, how to interpret these results? Could anyone help me with it? Thank you! I will put the code below, I actually used the SAS to code it, but I think the code also can clearly demonstrate the model.

proc phreg data=one;
class treatment(ref='1');
model time*survival(0) = treatment treat_ti;
weight ps_weight;

BTW, I used the propensity score weight to balance the group.

  • 1
    $\begingroup$ Please edit the question to show the model that gave the problem and the model that produced the results that you show. There's a chance that you fell into a common trap with time-varying effects in survival models, as this question originally did, but it's hard to tell without the code. When you edit, please use the {} tool on the toolbar to display code and associated results properly (or paste in text with each line having 4 leading spaces). $\endgroup$
    – EdM
    May 23, 2022 at 19:14
  • $\begingroup$ Thank you! I edited it with the code provided. I used the SAS, but I think it also can clearly illustrate the model that I used. $\endgroup$
    – NewRUser
    May 23, 2022 at 20:23
  • $\begingroup$ This depends on how SAS interprets your code, and I don't use SAS. Perhaps your treat_ti is just reporting a single "time-static" value for each individual based on the product of its group number and event/censoring time? Study Section 4.2 of the R time dependence vignette, the "dichotomy" between time-static and time-varying covariates in such modeling, and the warning about SAS coding on page 23. Seems like you have produced a "time-static" covariate. Revise the question when you're sure the coding is correct. $\endgroup$
    – EdM
    May 23, 2022 at 21:10
  • $\begingroup$ Thank you for the explanation, I think the code that I used is correct because I learned the coding from link. But when I check the article you sent, I noticed that they put log before time, and in another article linkwhich use SAS and R as example, they also use log, should I do it too? And how do you think I should revise my question? Change time-varying to time-static? Thank you! $\endgroup$
    – NewRUser
    May 23, 2022 at 22:35
  • $\begingroup$ My answer deals with your question. It's possible, however, that you don't need to worry so much about the apparent violation of proportional hazards (PH) or that there are better ways of dealing with it. See this page, this page, and this page, for example. Omitted or improperly modeled predictors can lead to an apparent PH violation. Fix those first. A large study can have a "significant" PH violation that doesn't matter in practice. $\endgroup$
    – EdM
    May 24, 2022 at 15:22

1 Answer 1


It is possible to use the time-transform mechanism to model time-varying coefficients for polytomous predictors in survival models, but it takes some extra effort.

If you try to do something similar in R to what you did with SAS, you get an error message instead. Using the veteran data set with its 4-level celltype predictor, and following the procedure indicated in Section 4.2 of the time dependence vignette:

fitFail <- coxph(Surv(time,status)~celltype + tt(celltype),data=veteran,tt=function(x,t,...) x * log(20+t))
# Error in coxph(Surv(time, status) ~ celltype + tt(celltype), data = veteran,  : 
#   data contains an infinite predictor
# In addition: Warning message:
# In Ops.factor(x, log(20 + t)) : ‘*’ not meaningful for factors

The software doesn't know how to parse the indicated time-transform (tt) function, which tries to multiply a categorical predictor by a function of time.

You can make this work if you recognize that celltype is represented in the model matrix by 3 separate columns beyond the intercept (which corresponds to the reference level of the predictor). If you extract those columns into a matrix and then multiply by the same function of time:

fitPolyt <- coxph(Surv(time,status)~celltype + tt(celltype),data=veteran,tt=function(x,t,...) model.matrix(~x)[,2:4] * log(20+t))
# Call:
# coxph(formula = Surv(time, status) ~ celltype + tt(celltype), 
#    data = veteran, tt = function(x, t, ...) model.matrix(~x)[, 
#         2:4] * log(20 + t))
#                            coef exp(coef) se(coef)      z      p
# celltypesmallcell       0.28513   1.32993  1.31895  0.216 0.8288
# celltypeadeno          -2.47072   0.08452  1.83996 -1.343 0.1793
# celltypelarge          -3.70383   0.02463  1.66983 -2.218 0.0265
# tt(celltype)xsmallcell  0.12856   1.13718  0.29261  0.439 0.6604
# tt(celltype)xadeno      0.82688   2.28617  0.42261  1.957 0.0504
# tt(celltype)xlarge      0.82076   2.27224  0.33824  2.427 0.0152
# Likelihood ratio test=34.08  on 6 df, p=6.484e-06
# n= 137, number of events= 128 

you get tt coefficients for each of the non-reference levels of the predictor. Any similarly constructed function of the polytomous covariate and time could be used, in principle.

Note, however, that the cox.zph() function does not work on models with tt() terms. The choice of tt() function is typically based on visual inspection of the shape of the plot of smoothed scaled Schoenfeld residuals over time.

I suspect that you can do something similar in SAS, but I don't use it and software-specific questions are off-topic on this site.

  • $\begingroup$ Thank you so much! I saw that I might need to log the time since the predictor is categorical, but what if some of my time values are 0? That will make some values of log(time) impossible. In this case, do you suggest that I should still log the time or I can just use time instead? $\endgroup$
    – NewRUser
    May 24, 2022 at 20:40
  • $\begingroup$ @NewRUser you don't need to log the time. The danger is that without the log transform you might end up putting more emphasis than warranted on late times with fewer cases at risk, while log-time tends to weight later times less. Any "0" times are (perhaps silently) ignored by the Cox fitting software. You could use some arbitrary small number like 0.5 days instead of 0, unless the event actually happened at the exact time of study enrollment. Cox models don't work with the actual times, just the ordering of events in time, so any value less than the smallest non-0 time will work. $\endgroup$
    – EdM
    May 24, 2022 at 22:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.