I've used pspline to add restricted cubic splines to certain continuous covariates that violated the PH assumption upon first analysis. I identified these covarites using cox.zph and also used Schoenfeld plots to plot the offending covariates (p<0.05).
After using pspline and re-running cox.zph I see that my offending variables are no longer significant. For completeness, however, I'd also like to show a new plot that shows the violation has been "fixed" because of the spline (i.e., that time-dependency has been accounted for).
However, when I make Schoenfeld plots for the new model (that also incorporates pspline) the plot looks hardly any different. I guess this is because the spline terms are all associated into one coefficient (fourth comment on the answer) but this then begs the question: how instead can I "show" graphically that the PH violation has been dealt with in a similar way that I "show" PH is not satisfied using Schoenfeld plots. I should be consistent with the reporting of my results.
Edit: Number of events in dataset = 1600, N=100k
pspline()function doesn't provide restricted cubic splines. Instead it provides a penalized smoothing spline whose coefficients can't be readily interpreted on their own. If you want restricted cubic splines, try thercs()function in Frank Harrell'srmspackage. Then specifyterms=FALSEin the call tocox.zph()and you can get plots for each (unpenalized) spline coefficient separately. – EdM Aug 16 '22 at 15:07