I have tried to visualise the proportional hazard (ph) assumption for my cox-ph model by plotting the smooth scaled schoenfeld residual plot. I'm working in R.
The residual plot is computed in R like this
test_ph <- cox.zph(my_cox_model, transform="identity")
plot(test_ph, resid=FALSE)
The output I get from the above code, is the y-axis representing the Beta(t) for the covariate over time (on the x-axis). However, I'm interested in having the y-axis showing the HR. I have tried to exp() the y values like in the following codes, and then plotting, but it did not work
hr <- exp(test_ph$y)
time <- test_ph$time #X-axis is names time in the test_ph
plot(time, hr)
How can I convert the values on the y-axis to HR for the residual plot? Is it even possible to do that, or should I code a plot for plotting the HR over time instead converting the y-axis on the residual plot? If so, I still need help to code the latter in R as well.
Thank you in advance!
Does your model have multiple predictors? Are you getting this error?
Error in xy.coords(x, y, xlabel, ylabel, log) :
'x' and 'y' lengths differ
If so, the problem is that test_ph$y has multiple columns, one for each predictor, and plot can only do one at a time. Try:
hr <- exp(test_ph$y)
time <- test_ph$time #X-axis is names time in the test_ph
plot(time, hr[,1])
To see the plot for the first predictor.