How to calculate the heterogeneity P value for a variable in two Cox regression models

In my opinion, the heterogeneity P value is used to determine whether a variable has a differential effect on survival in the first 5 years versus the latter 5 years. But I don’t know how to calculate it.

It is extremely unlikely that something magical happens at exactly 5 years, so the researchers should have tested for a smooth departure from a constant hazard ratio. This is readily done with the smoothed scaled Schoenfeld residual “correlation with time” test as implemented in the R survival package cox.zph function. As for getting the test with the artificial dichotomization represented by the table, if we can assume conditional independence and you can the the standard errors of the log HRs, you can get the test by computing the difference in log hazard ratios divided by its standard error, where the SE is the square root of the sum of squares of the two SEs. This provides a z-test. Note that this should not be labeled ‘heterogeneity’ but instead ‘non-proportional hazards’ or ‘non-constant effect over time’.

Thank you, professor. But I just want to know if a variable has a different effect on survival in the two dichotomized time periods. What can I do?

1 Like

Don’t do anything until you confirm with the smooth scaled Schoenfeld residual plot that the change in effect can be well approximated by two piecewise flat relationships where the point of discontinuity happens to be at 5y. Then you can do a time-stratified analysis, i.e., two separate Cox model analyses with “fake censoring” of the first and a landmark analysis for the second where follow-up starts over at the original 5y point.

What is often better is to choose an accelerated failure time parametric survival model that makes a different assumption than proportional hazards, if its new assumptions are well-enough satisfied. And don’t discretize time.

Dr. Frank Harrell,

Complementing this topic, I have one question regarding interpretation of the proportional hazards assumption with cox.zph:

  • The output includes a chisq, df and p-value for every variable and also for the GLOBAL model. How should I interpret a test that gives me a GLOBAL proportional models p-value > 0.5 but one of the variables presented are in fact < 0.05?

The question may be silly, but I understand that multiple tests for non-proportional hazards can have p-values less than 0.05 by random chance (as would normally happen in multiple testing analysis) and the GLOBAL value reflect an overall model fit for non-proportionally (similar to a correction for multiple testing).

Adding on this: in this situation of a divergent GLOBAL and specific variable p-value for cox.xph, how would you report in an article?

Thanks in advance,

Felippe

The global test has a near-perfect multiplicity adjustment, so I usually would not look at the individual tests if the global test p-value > 0.2 unless an individual \rho value (correlation between Schoenfeld residual and time) is very high.

Thanks Dr. Frank Harrell for the answer. It clarifies my question!

Best,

Felippe