Interpretation of cox.zph Output with Smoothing Splines in R

External question — Cross Validated Stack Exchange Author: Brandon Original post: https://stats.stackexchange.com/questions/400815 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. In Luke Keele's paper, "Covariate Functional Form in Cox Models", Dr. Keele carries out two Grambsch and Therneau non-proportionality tests, that is, one for a model without splines and one for a model with splines, using the cox.zph function in R. Interpretation of the results of cox.zph for the model without splines is straight forward, variables found to be statistically significant have a time varying effects on the hazard, however, for a model with splines, the output from cox.zph may have multiple entries per spline variable. How do you interpret these results? Here is an example to add context to my question. The following is a cox proportional hazards model without splines: coxMod1 <- coxph(Surv(starta, stopa, dispute) ~ nudem + nugrow + allies + contig + nucapab + trade + sumdisp, data = war, robust = FALSE, method = "efron") zph1 <- cox.zph(coxMod1) From the results of cox.zph , we see that there are only two variables, allies and sumdisp (Previous Disputes), which seems to have time varying effects on the hazard. If we add smoothing splines for the continuous variables nugrow (Economic Growth), nucapab (Capability Ratio), and trade we get: coxMod2 <- coxph(Surv(starta, stopa, dispute) ~ nudem + pspline(nugrow, df=4) + allies + contig + pspline(nucapab, df=4) + pspline(trade, df=4) + sumdisp, data = war, robust = TRUE, method = "efron") zph2 <- cox.zph(coxMod2) The interpretation remains the same for the variables which do not have splines, but for the three that do, how do you interpret these results? If I were to take a stab at interpreting the spline variables which have all rho as not statistically significant, like nugrow, I would say that since all rho are not statistically significant, the variable nugrow does not have a time varying effect on the hazard. Likewise, if the spline variable had all rho statistically significant, intuition would lead me to say that the variable has a time varying effect on the hazard. Is this safe to say? What happens when the spline variable has rho which are both statistically significant and not statistically significant, such as trade. Trade has ps(trade)3, 4, 5, & 6 being statistically significant and the rest are not statistically significant. How do you interpret this? In the paper, Dr. Keele calculates the average rho, chi square, and p-value for each of the spline variables and then assesses the time varying effect on the hazard using the statistical significance of each variable's average rho. Does this seem like a valid approach? The paper can be found here: https://pdfs.semanticscholar.org/0610/022cd1176e21820f992d0f181fe6d2a6ce6d.pdf The data for this paper and this example can be found here: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/VJAHRG

External answer — Cross Validated Stack Exchange Author: EdM Original post: https://stats.stackexchange.com/a/654451 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. The paper linked from the question (evidently a working paper) no longer seems to be available, but there is a final published version: L. Keele, Proportionally Difficult: Testing for Nonproportional Hazards in Cox Models , Political Analysis, Volume 18, Issue 2, Spring 2010, pp. 189 - 205. With respect to output from the cox.zph() function for testing proportional hazards, Footonote 11 says: "For the variables with spline fits, the Therneau-Grambsch test calculates a test statistic and p-value for each knot placement. I report the average p-value across the spline fit." As the OP recognizes, this is not a good approach. It's also no longer necessary to deal with separate cox.zph() results for the multiple spline terms associated with a single continuous predictor. In early versions of the survival package, as used in that study, cox.zph() used an approximation based on the correlation coefficient ( rho ) for the scaled Schoenfeld residuals versus the chosen transformation of time (default, a transformation based on the overall Kaplan-Meier curve). As explained in this answer , starting with Version 3.0-10 of that package a score test has been used. Furthermore, cox.zph() in recent versions of the package returns score-test chi-square values values that (by default) combine all the terms associated with a predictor. You can in principle override that default by setting terms=FALSE as an argument to the function, but in my quick tests that didn't work for individual pspline() terms. More generally, when a predictor has multiple associated terms (splines as here, a multi-level categorical predictor, etc.) it's risky to try to interpret the "significance" of any individual term whether for a regression coefficient or a test like cox.zph() . That "significance" will typically depend on details of how the predictor was coded. You need a test that combines all the terms, as is now provided by default for cox.zph() .