Is the proportional hazards assumption violated (interpreting Schoenfeld residuals)? What is my best option if so?

Question

I am using an extension of the Cox model in the counting process format (Andersen-Gill model), with time-dependent covariates. The analysis is investigating the effect of an intervention on hospital admissions, adjusting for known confounders that change over time, through the inclusion of time-dependent covariates in the time-dependent version of the coxph function. In the example below EmergencyAdmission is a binary indicator of whether an admission has occurred at each time interval, and the covariates CareHomeStay and EDAttendance are binary indicators of whether a patient had a stay in a care home or ED attendance at each time interval (all other covariates are binary except for DeprivationIndicator and CharlsonGroupIntegers which have 5 levels and SPARRAScorePeriodPrior which is an integer between 1-50 roughly). These covariates are known confounders, in terms of their prior history. The cluster() term accounts for repeated measures.

results<-coxph(Surv(timeStart,timeStop,EmergencyAdmission,type="counting") ~ InterventionIndicator + DeprivationIndicator + CareHomeStay + CharlsonGroupInteger + DeporEmoConcernInteger + SPARRAScorePeriodPrior + EmergencyAdmissionPeriodPrior + EDAttendancePeriodPrior + CommunityEpisodePeriodPrior + OutpatientAttendancePeriodPrior + DelayedDischargePeriodPrior + cluster(AnonID),data=CohortTimeSeriesMonthly)

The format of the data can be seen in my previous question.

I have three questions:

Question 1 - Do I need to test the proportional hazards assumption in a Cox model with time-dependent covariates?

1. My first question in a clarification of what I believe to be true - I believe am using time-dependent COVARIATES not COEFFICIENTS, hence I still need to test the proportionality assumption. I have done a lot of reading on this and this is asked in over CV questions (which I will reference further on), which give some conflicting answers but I believe it is down to terminology. From my understanding, a time-dependent covariate is simply a variable included in the model whose values change through time, whereas a time-dependent coefficient is an interaction term included in the model for the interaction between time and a given covariate. The two models are given mathematically in Therneau's vignette.

The answers to the other CV questions are conflicting (Question1 suggests PH assumption still needs to be tested with time-dep covariates, whereas Question2 suggests otherwise and Question3 gives two conflicting answers). However, overall they seem to indicate (given the differences are because of terminology) that in a Cox model using time-dependent covariates, the PH assumption still needs to be tested, whereas in a model with time-dependent coefficients, the PH assumption no longer holds. Is this a correct interpretation of these answers? Is it purely down to terminology? What further convinces me is this excerpt from Grant, Chen and May (2014)

A crucial assumption of the PH model is that the effect of a covariate does not change over time (Cox 1972). In other words, β are assumed to be constant for all t. This assumption applies even in the case of time-dependent covariates; though values may change, the effect of the covariate is assumed to be constant.

Question 2 - How do I interpret these results to my test of the PH assumption?

2. If I'm correct in my interpretation to 1) above, I need to test the PH assumption. When I do that I get the following results based on the above model.

> cox.zph(results2)
                                  chisq df       p
InterventionIndicator           6.20e-01  1  0.4311
SIMD16QuintileInteger           4.87e-01  1  0.4853
CareHomeStay                    1.74e-01  1  0.6768
FeedingConcernsInteger          5.43e-01  1  0.4612
CharlsonGroupInteger            3.12e+00  1  0.0773
DeporEmoConcernInteger          3.46e+00  1  0.0628
SPARRAScorePeriodPrior          2.09e+01  1 4.9e-06
EmergencyAdmissionPeriodPrior   1.30e+01  1  0.0003
EDAttendancePeriodPrior         1.58e+00  1  0.2089
CommunityEpisodePeriodPrior     5.94e+00  1  0.0148
OutpatientAttendancePeriodPrior 6.39e-08  1  0.9998
DelayedDischargePeriodPrior     1.24e+00  1  0.2651
GLOBAL                          4.44e+01 12 1.3e-05

>ggcoxzph(cox.zph(results))

The cox.zph documentation states

The plot gives an estimate of the time-dependent coefficient β(t). If the proportional hazards assumption holds then the true β(t) function would be a horizontal line. The table component provides the results of a formal score test for slope=0, a linear fit to the plot would approximate the test.

Although the test indicates significance on several variables and on the global test which I interpret to mean the PH assumption is violated, the plots of the scaled Schoenfield residuals indicate a monotone trend of the residuals over time across all variables. Why is the test statistic significant then? Based on the answer to this question, large sample sizes can cause seemingly strong evidence against the PH assumption, so I am wondering it that's it (over 100,000 rows/observations given 42 monthly observations for about 4000 patients). How do I interpret these test statistics and Schoenfield residuals?

Question 3 - What is my best option for dealing with the violation of the PH assumption?

3. If after all, the PH is violated, how should I proceed? What is my best option? This question ,among others, suggest using a time interaction term (i.e. adding at time-dependent coefficient as I understand it) in the form of a tt() term in the model. The issue with this approach however is that the cox.zph() function cannnot proccess tt() terms so I would be unable to test if the PH assumption is no longer violated (highlighted at the end of this question). Another way to include a time interaction term is by including a variable:time term as suggested in this question, however, I run out of memory when I try this (error message: Cannot allocate vector of size 4.0 GB). I have previously fixed memory problems by getting a bigger machine but not sure it will it this time because previously the errors have only been roughly 200MB big, so this feels like a massive jump which probably won't work. I am also aware of the step function option to stratify on specific intervals as suggested here and given in Therneau's vignette. If this is my best option, how might I do this? Would my model suffer from the stratification? Any suggestions? Even if the PH is not violated in answer to Q2 above, I would appreciate help with this as I will likely run into this again.

Any help appreciated (please!) thank you!

Answer 1

Question 1. If proportional hazards (PH) don't hold, then the concept of a time-independent hazard ratio doesn't hold. I hold to my answer on your linked Question 1, which agrees with your quote from Grant et al. PH are important, and can be checked in a model with time-dependent covariates. As the covariate values used in modeling are instantaneous as of each event time, there is still a meaningful log-hazard difference associated with each change in a time-varying covariate over time, even if there is some question about just what the "baseline hazard" is and thus some potential terminological confusion.

Question 2. Large sample sizes can lead to significant PH violations that aren't of practical importance, just like large sample sizes can lead to significant normality violations that aren't of practical importance. Only you can judge the practical importance, based on your understanding of the subject matter. At first glance your residual plots don't look bad, but there seem to be a lot of multiple cases having the same residual/time combinations so it's hard to tell. Try a loess smoothing of the plots of residuals against (transformed) time (that's done by default in the standard plot.cox.zph() function), and omit the individual points. NOTE: your plots seem to have a reasonably small number of individual time points, suggesting that discrete-time analysis might be more appropriate here. It's not clear why times to EmergencyAdmissions are so discretized; you might need to consider that in your model.

Question 3. There is no one-size-fits all answer. For continuous predictors, proper transformation of the predictor might help. Stratification by a troublesome covariate, time strata, and time-dependent coefficients are standard ways to proceed, explained nicely in the survival package vignette. Just make sure that you handle time-dependent coefficients properly, unlike in this answer to your linked Question 1. If you use time-dependent coefficients you can't thereafter do a cox.zph() test; you simply impose a reasonable estimate of the proper form of the coefficient versus time that matches adequately what you saw with the time-constant-coefficient plot. These lecture notes go into detail about other tests of PH and what to do when they fail.

Time-dependent coefficients get really complicated to think about with recurrent events, however: do you want to have thde time-dependence start from time 0 of entry into the study, or from the time of a previous Admission? Should you be evaluating each subsequent Admission as a different type of event from prior Admissions? Those matters need to be evaluated carefully from combined clinical and statistical perspectives, and might require a change in your model. You don't want to get a seemingly beautiful answer to a meaningless question.