For the statistical part of the question, if you don't care about modeling the effect of Test_Temp_C
itself, stratification is often a good choice for a binary or factor predictor with a small number of values. You might want to think, however, if there's something important hiding in the observation that this predictor is associated with greater hazard at early times but less hazard at later times. That said, the associate hazards (on the order of +/- 10-20%) might be too small to matter in your study.
For visualization, a cox.zph
object contains a component x
that represents the transformed time axis, a component time
that contains the original event-time values, and a matrix y
of scaled Schoenfeld residuals. You could extract the relevant values from the object into a data frame, pull out the categorical covariate values that apply at each event time from your data, and color points by generating a plot yourself and specifying in the plot command a point color based on the value of your binary covariate, something like col=c("red","blue")[1+Test_Temp_C=="320"]
with basic R graphics.* If there are no tied event times that should be straightforward. If there are ties, you will have to determine which of the cases has which of the residual values by looking directly at the residuals of the coxph
object.
*The col
argument in plot.cox.zph()
is for the curve, not for the points. There are more elegant ways to do such coloring with ggplot()
, and you might consider starting with the survminer
package and its ggcoxdiagnostics()
function to get a plot similar to yours, then do the point coloring with an added geom_point()
command.