Schoenfeld test (cox.zph) shows no covariate violates PH assumption but global test suggests whole model does (p<0.001). What to do?

Question

Using cox.zph (part of the survival package), I find that no covariate violates the PH assumption in a Cox regression model (p > 0.05, all). However, the global test suggests the model violates the PH assumption (GLOBAL, p < 0.001).

I've attached the output of cox.zph as an image (see below). I have already added interaction terms in order to remove violations of the PH assumption for individual covariates. Now, I'm stuck. What should I do?

• Temper interpretations/conclusions drawn from HRs?
• Include time-variation in the model somehow?

Addendum:

I input the following:

k <-coxph(Surv(int_risk, dead) ~ 
            age +
            race +
            sex + 
            bmi +
            lab1 + 
            lab2 +
            lab3 + 
            lab4 +
            class1 +
            event +
            event2 +
            event3 + 

            age:event +
            age:lab1 +
            age:lab2 +
            age:lab3 +
            age:lab4 +
            age:event2 + 
            age:class1 +
            age:event3 + 

            race:bmi+

            class1:lab3, 
          data = df)
kph <- cox.zph(k)

I'm looking at 15,003 observations with 5,001 events. I omit observations if missing data (~1200 observations).

Answer 1

Note that the p values from cox.zph() shouldn't be thought of in the same way that you use p values in standard tests of null hypotheses. The burden is on you to demonstrate that PH holds; it's not enough to show that you can't rule out the PH null hypothesis because p > 0.05. A p value on the order of 0.1, as you have for some of your cox.zph() results, is certainly consistent with a violation of the PH assumption.

On the other hand, with 5001 events it's quite possible to find a statistically significant non-proportionality that is of little practical significance. As Therneau and Grambsch put it in Section 6.5 of their book:

A "significant" nonproportionality may make no difference to the interpretation of a data set, particularly for large sampIe sizes.

They go on to describe a study with 771 events, in which the cox.zph() test gave significant p values for 2 out of 3 predictors with Global p = 1.95e-09. This result might be taken to require time-dependent estimates of each coefficient, $\hat\beta(t)$, instead of a time-independent estimate $\beta$. Examining plots of scaled Schoenfeld residuals over time, however, led to the following conclusion about the association between age and risk of death:

There is a profound effect of age on the risk of death, and the test shows the effect to be significantly nonproportional. However, the smoothed residual plot shows that the variation in $\hat\beta(t)$ is small relative to $\hat\beta$.

Looking at the results instead of just accepting a p value is always a good idea. Pull back to a model omitting all interaction terms except for those that are of substantive interest to the goals of your study. See whether the "significant" trends in residuals over time are of sufficient magnitude to matter.

Here are a few more items to consider if you still find the lack of PH to be an obstacle to your study.

First, from their names it looks like some of your predictors are dummy-coded. For example, you have lab1, lab 2, lab3, and lab 4. If these represent mutually exclusive binary representations of a categorical predictor called lab and lab has 4 levels, then you might be getting into problems with linearly dependent predictors. I don't know whether that would cause a problem with coxph() but it's possible that cox.zph() can't handle that and leads to a large global $\chi^2$. The coding of event, however, suggests that you have taken that into account (perhaps using using event1 as the reference for event, with similar precautions for lab and other potential categorical predictors) so I note this mostly for others who might read this thread.

If that's not the case, next examine whether the linearity assumptions for your continuous predictors, evidently age and bmi, are valid. Failure to satisfy the linearity assumption can lead to problems with satisfying the proportional hazards (PH) assumption, too. See this page for one example with a plausibility argument for why this can be the case. So attend to the linearity assumption. Modeling such a predictor with restricted cubic splines, for example, might solve the PH problem and make further adjustments unnecessary.

Failing that, consider stratification on predictors that don't meet the PH assumption instead of adding interaction terms.* That is trivial for categorical predictors. For a continuous predictor like age, if age still doesn't meet the PH assumption despite your attempts at transforming it to meet the linearity assumption, and age is included as a control covariate that is not of primary interest to you, you could subset by age into a reasonable number (order of 5) of categories and include strata(age.category) as a predictor instead of age. You won't get a hazard ratio for age but this can solve the PH problem if age is to blame.

If the predictor giving you problems with PH is of primary interest then you can model it with a time-dependent coefficient. Make sure, however, that you do this properly as you can find seemingly reputable references to completely inappropriate ways to proceed. See this page for some discussion and links to correct details.

Finally, for further details on Cox PH modeling and how to deal with problems that arise, get and read carefully a copy of Harrell's Regression Modeling Strategies. You will find that very valuable if you will be doing a lot of regression modeling (Cox or otherwise), as your comment suggests will be the case. It ties in nicely to the modeling tools provided by Harrell's rms package in R.

---

*At first thought, it would seem that if you have a continuous predictor that isn't meeting the PH assumption, and then you add interactions of that predictor with others, you would make the overall PH problem worse rather than better.