0

I am having trouble understanding how to interpret a Schoenfeld plot and how to correct for the violation of PH assumption in Cox PH model when using continuous numeric covariates.

Here is the code I am using to call the model, the variable I am investigating is average spring/summer precipitation (AS_avg_total_p_mm) over the study period (10 years):

The image and plot below is the output of the cox.zph() test, and the test graphed using ggcoxzph().

While the statistics indicate the model violates the PH assumption, I am unsure how to correct for this. I have read about introducing a time interaction variable into the model (either by var:time, or tt(var) and adding a tt function), but I am unsure how to do this properly. My knowledge of statistics and modeling is very basic.

When I add the var:time variable to my model, the covariate and the interaction are significant, but even further fail the cox.zph() significance test (image below). I have also read that this might be an inappropriate approach.

Introducing a tt(var) hinders me from using cox.zph(), as the function does not support the use of this.

Therefore:

  • How can I implement a time interaction and still test that my model is following the rules?

  • How do I interpret a Shoenfeld plot that looks like mine? Most that I see example wise are scattered or follow some specific pattern. What should a "good" plot look like?

  • Is the failure of this assumption drastic enough to warrant changes to my model? Is this up to my interpretation?

  • Does my sample size play a role in how these tests are behaving? It is quite large.

Please let me know if any additional information is needed here and forgive me if my question is redundant. I have tried many examples to solve my problem. Thank you for your time!

Tyler Brown
  • 3
  • 1
    Welcome to Cross Validated! Try fitting your continuous predictor flexibly, like with a regression spline. It’s unlikely that there’s a strictly linear association between that predictor and the log-hazard of an event, and a poorly specified form for a continuous predictor can show up like a PH violation. Also, you might be using an old version of ggcoxzph, with a serious coding error that makes its plots nearly useless. Use the standard survival function tools. – EdM Dec 18 '23 at 22:47
  • See section 4.2 page 21 of https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf for why your time:covariate interaction term is problematic. – EdM Dec 18 '23 at 23:03
  • Thank you for your comments, I actually just saw the post where you mentioned issues with ggcoxzph and have been plotting in accordance to the survival package. Additionally, looking at the pdf, specifically the spline viginette, I am having a bit of trouble understanding the effective df. From my understanding, pspline creates 3 df for each variable? How do I choose df for the model? I did a test with the smoothHR package indicating my df = 3.06, however when I use 3 4 or 5 df my model fails the cox.zph test, passing when df = 6. Should I be concerned the model is not passing still? – Tyler Brown Dec 19 '23 at 01:49
  • and to clarify this is without a time interaction, only the non-linear covariate by itself. Should I try pspline and a time interaction? – Tyler Brown Dec 19 '23 at 01:50

1 Answers1

0

Section 6.6 of Therneau and Grambsch discusses causes of non-proportionality. A couple of issues there should be dealt with before you jump to assuming that there your predictor requires a time-varying coefficient.

First, your model allows for only a linear association between your continuous predictor AS_avg_total_p_mm and the log-hazard of an event. Strictly linear associations are seldom correct. Mis-specifying the functional form of a continuous predictor's association with outcome can lead to an apparent violation of proportional hazards (PH), while there might be no violation if that functional form were specified properly. See this page, for example. If you perform a flexible fit in the Cox model,* for example with a spline, your problem might go away.

Second, your model only includes that one predictor. Omitting any outcome-associated predictor from a Cox model can lead to bias in the coefficient estimates for included predictors and to apparent violations of PH, violations that would also disappear if you didn't omit that predictor. If you have information on other variables that might be associated with outcome, based your understanding of the subject matter, include them.

With 259 events you have a lot of potential flexibility in dealing with these issues. Frank Harrell recommends, as a general strategy in regression modeling, to decide how many degrees of freedom (df) you can afford to spend (coefficients to estimate, adjusted for penalization), to decide where to spend them, and then to spend them. With 259 events, the rule of thumb of about 15 events per df in survival models means that you could probably afford to spend 15 df. That allows you to include several additional predictors while including a very flexible fit of AS_avg_total_p_mm.

There are several different approaches to smooth fitting of continuous predictors, outlined on this page. The ps() function provided in the survival package is a particular type of penalized smoothing spline. The penalization means that the effective number of degrees of freedom used up is smaller than the number of coefficients. An example in Section 5.8.3 of Therneau and Grambsch uses up only 4 degrees of freedom for 13 coefficients. You also could consider an unpenalized regression spline; see Section 2.4 of Harrell's Regression Modeling Strategies.

Even if PH is violated, there might be no need to find the function of time that fits the predictor's coefficient. First, as AdamO frequently points out (for example here), the coefficient of a Cox model is then a type of time-averaged value. If you use robust standard errors (cluster() term), you can still do valid inference. Second, the plot of the smoothed Schoenfeld residuals over time is the estimated form of the coefficient's value over time. That is probably more useful to the reader than the results of a numerical fit. And, as you've learned, there's no good way to document that your numerical fit solved the PH problem except for showing that the shape of the numerical fit matches the shape of the Schoenfeld residual plot.

*I mean the Cox model of survival here without any time-varying covariate, not the model that tries to accommodate a general form of a time-varying coefficient as in Section 4.2 of the time-dependence vignette. With a continuous predictor in a Cox model, it's critical to distinguish the functional form of its association with log-hazard from the functional form in time of its potentially time-varying coefficient. You appreciate that, as indicated in a comment, but this needs to be written for those that might later come upon this page.

EdM
  • 92,183
  • 10
  • 92
  • 267
  • thank you for the explanation, I am extremely grateful as it has helped me tremendously. I have been implementing more variables into my model, fitting the non-linear variables with pspline and assigning df appropriate to # of events. The issue I am running into is validating the model with cox.zph() and schoenfeld plots. Some of my variables look very flat and in my opinion acceptable, however the test indicates some variables and the model are still significantly non linear. Should I take this test with a grain of salt and go off of the plot? How does this effect the power/validity? – Tyler Brown Dec 20 '23 at 19:49
  • Additionally, some variables fit with pspline indicate that the linear form is non-significant in the model while the non-linear form is. Does this show that the spline is ineffective in correcting the proportionality issue? Should I remove variables that are unable to be corrected due to a lack of power? I do not really want to do this as they are important to my analysis. Should I provide model and test outputs for clarification on my question? – Tyler Brown Dec 20 '23 at 19:52
  • and just to add, I have been plotting the schoenfeld residuals with plot() rather than the previous method where issues were identified with the code. thanks again! – Tyler Brown Dec 20 '23 at 19:55
  • @TylerBrown unless your model is using up more than about 15 degrees of freedom with your 259 events, don't omit predictors or spline df. If the fit of a continuous predictor is improved by a spline over a linear fit in in the Cox model, that's good. That's not guaranteed to fix PH, however. – EdM Dec 20 '23 at 21:17
  • @TylerBrown If PH is violated, see this answer or several others by AdamO on proportional hazards. You get a sort of time-averaged coefficient estimate, with perhaps some loss of power. The Schoenfeld residual plot shows directly how the estimate changes over time. Show that plot and use a cluster term in the model to get robust standard errors that allow for valid inference if PH is violated. – EdM Dec 20 '23 at 21:21