0

I'm estimating a cause-specific hazard function using the coxph() function in R. This the first time I've run such an estimation. In performing model diagnostics, I used cox.zph() and found that the GLOBAL was significant at p == 0.0010 and many of my individual covariates (categorical and continuous) were also significant at p < 0.05. I inspected the residuals using the plot() function on the object created by cox.zph() to visualize the residuals. My main question is: OK - what now?

A little bit about the model. I'm estimating in a competing risks framework using health care data where there are several forms of discharge, plus in-hospital death. I've separately run a crr() regression to estimate the impact of the covariates on relative incidence of the different discharge types. In the current step, I'm hoping to learn the effect of the covariates on the cause-specific hazard (instantaneous rate of discharge) for each discharge type.

Dan
  • 63

1 Answers1

1

I've answered my own question, it turns out. First, the title of the post was a misnomer; I'm actually wondering what to do about time-varying coefficients. This is a much different approach than dealing with time-varying covariates. For dealing with the coefficients, I took the following steps:

  1. Examine the Schoenfeld residual plots associated with coefficients that are time varying, which was determined from running cox.zph() on the coxph() model results. A simple set of R commands for plotting is:
plot(coxzphtest[n])
abline(h = 0, col = 2)

where coxzphtest is the name of the object resulting from the cox.zph() function and n is the index number of the covariate of interest in the list of covariates from the output of cox.zph.

  1. Plots of time-varying coefficients where the red line is within the margin of error on the plot might not need any adaptation. This is a judgement call.

  2. Otherwise, if the coefficients need to be re-estimated with a time interaction, there are multiple ways of doing this. Therneau, Crowson, and Atkinson (2022) have written an excellent, freely available tutorial for how to do this. Start on page 16. Based on the plots and what makes sense for the data, I chose a single cutoff of 7 days. But there are other types of adaptations.

Another resource, Prof. Marin's Stats Lectures, was also helpful for coding the plots and interpretation.

Dan
  • 63
  • I'm not certain that you are interpreting the plots of time-varying coefficients properly. I don't think that the standard survival plots for cox.zph objects print any red lines, and the ggplotcoxzph function from the survminer package (which does show individual values in red) has a major bug, at least as of earlier this year, that makes its plots pretty useless. What's important is the flatness of the curve. Yourcox.zph p-values do suggest that there is something important that you shouldn't ignore. How many cases/events do you have? – EdM May 09 '22 at 20:54
  • I have added the red line myself with the abline(h = 0, col = 2) command below the plot command in my posting. This was a technique demonstrated in the Stats Lecture series that I linked to in this post. I am familiar with the bug in ggplotcoxzph from earlier posts in this forum. I have 783 total observations in my data, 254 cases for this particular event. – Dan May 10 '22 at 13:59