data visualizationrsurvival
I used the veteran database of the package "survival" to do survival analysis. I tried to check the proportional hazard assumption, with with schoenfied residual plot. I always thought if we had a flat line this mean the covariate effect are constant over time, and therfore the proportional hazard assumption were verified. But I see that in these graphs below all lines are flat, although the p.value show the opposite. Finally I wonder how the shoenfield residuals are visually interpreted.
There is a separate Schoenfeld residual for each covariate for every uncensored individual.
In the case of a categorical covariate with $k$ levels, $k-1$ dummy variates appear in the model.
If the global p-value is significant then the PH assumption does not hold for the model. Schoenfeld is like a Shapiro-Wilk test of normality, if $p<0.05$ then the feature is not normally distributed. If Schoenfeld $p<0.05$, then the model or feature does not meet the PH assumption. You would need the global p-value to not be less than 0.05 before you start dipping into the individual feature p-values.
Also for Cox PH regression and the PH assumption, it's better to have categorical features for which you can plot the survivorship function and the follow-up days (months) on the log scale to see if the lines cross. (although the PH assumption is based on the hazard functions crossing, you can still see the crossings within the survivorship functions). The attached plots show these S(t). vs days (log scale) for a 3-level categorical factor and a 4-level factor. The PH assumption holds for the 3-level factor, but not for the 4-level factor, since one of the lines cross. The global p-value for the model was about 0.4, and the p-value for the 4-level factor was <0.05.
Best Answer
This question illustrates both some weaknesses and a strength of using open-source statistical software. These plots were evidently generated by the
ggcoxzph()function in the Rsurvminerpackage.The main weakness: a coding problem. Although that function is supposed to be a wrapper around the standard R
survival:::plot.cox.zph()function to generateggplot()output, the y-axis limits are much broader than those seen in standard plots; see corresponding plots for similar models on this data set in Section 3.5.2 of the R survival vignette. Those broad limits are needed to include what seem to be displayed as "confidence bands," but they make it impossible to see details of the smoothed line for scaled Schoenfeld residuals versus time.The strength: you can find and fix the coding problem. With open-source R, the code for
ggcoxzph()and forsurvival:::plot.cox.zph()is open for inspection. The wide "confidence bands" focus attention on how their display might differ between the 2 functions. In the recent survminer_0.4.9 version ofggcoxzph(), a critical line of code for calculating those bands differs fromplot.cox.zph()by a multiplicative factor. The offending line:where
dis earlier set to the number of events andd *is not found in thesurvival:::plot.cox.zph()code. Removing that addition to the code allows plots with reasonable y-axis limits and correct confidence bands.The next weakness: Even with that fix, it's hard to see details of the shapes of the smoothed curves, as all y axes are compressed to fit plots for all predictors into a single display. That's evidently the default for that function--it's important to question such default settings even if there aren't errors in code. The original
plot.cox.zph()produces plots for one predictor at a time, which allows for better visualization. Here's the example forcelltype:With the physically larger y axis, details of the curve are much more apparent. Remember that those visual curve details depend on another parameter setting, the transformation of the
timeaxis chosen in thecox.zph()function. The plots shown here are for the default "KM" (Kaplan-Meier) transformation, but there are other options and any function of one argument can be specified.Code: