What are the options in proportional hazard regression model when Schoenfeld residuals are not good?

Question

I am doing a Cox proportional hazards regression in R using coxph, which includes many variables. The Martingale residuals look great, and the Schoenfeld residuals are great for ALMOST all of the variables. There are three variables whose Schoenfeld residuals are not flat, and the nature of the variables is such that it makes sense that they could vary with time.

These are variables I'm not really interested in, so making them strata would be fine. However, all of them are continuous variables, not categorical variables. So I perceive strata to not be a viable route*. I have tried building interactions between the variables and time, as described here, but we get the error:

  In fitter(X, Y, strats, offset, init, control, weights = weights,  :
  Ran out of iterations and did not converge

I'm working with nearly 1000 data points, and am working with half a dozen variables with many factors each, so it feels like we're pushing the limits of how this data can be sliced and diced. Unfortunately, all the simpler models I've tried with fewer included variables are clearly worse (ex. Schoenfeld residuals are crumbier for more variables).

What are my options? Since I don't care about these particular poorly-behaved variables, I'd like to just ignore their output, but I suspect that isn't a valid interpretation!

*One is continuous, one is an integer with a range of over 100, and one is an integer with a range of 6. Perhaps binning?

Answer 1

The most elegant way would be to use a parametric survival model (Gompertz, Weibull, Exponential, ...) if you have some idea what the baseline hazard might look like.

If you want to stay with your Cox model you can take up an extended cox model with time-dependent coefficients. Bear in mind that there are also extended cox models with time depending covariats - these do not solve your problem!

For R see here: http://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf

Answer 2

Couple of ideas -

1) Try the Royston-Parmar modelling approach e.g. http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0047804 and references therein. We've had useful results with it.

2) Centring and standardising continuous variables can be useful numerically.

3) In many models with factors with lots of levels there are a few levels where there are basically no data. Merging levels to remove these, but based on good substantive criteria, can be very helpful.

Good luck!

Answer 3

If using an interaction with the underlying time doesn't work, you can try step functions (for more information see Therneau's 2016 vignette).

Step functions is stratify in specific coefficients at specific intervals. After seeing your plotted Schoenfeld residuals for the problematic covariates (i.e. plot(cox.zph(model.coxph))) you need to visually check where the lines change angle. Try to find one or two points where the beta seems markedly different. Suppose this occurred at time 10 and 20. So we will create data using survSplit() from the survival package which will create a data frame for the specific data model grouping at the aforementioned times:

step.data <- survSplit(Surv(t1, t2, event) ~ 
                      x1 + x2,
                      data = data, cut = c(10, 20), episode = "tgroup")

And then run the cox.ph model with the strata function as interactions with the problematic variables (as with interacting with time, do not add a main effect for time or the strata):

> model.coxph2 <- coxph(Surv(t1, t2, event) ~ 
                          x1 + x2:strata(tgroup), data = step.data)

And that should help.