How to handle multiple violations of the Schoenfeld residuals assumption in cox proportional hazard models? below you can see the cox.zph function to test for proportional hazards used on my data. 
Issue number 1: there are multiple violations of the proportional hazards assumption (these are tree_species 3 and 4, Skewness, Area 2 and Cavity height). Do I prioritise stratifying/making an interaction with time for the most significant variable first and then check the assumption again (and again until at least the global assumption is met)? Or do I correct them all at the same time?
Issue number 2: the variable tree_species has 4 levels. Just tree_species 3 and 4 violate the PH assumption. Do I correct this variable despite tree_species 2 not having any issues? Same goes for the variable "area".


 A: Be cautious about assessing model fit statistics with tests. The statistical significance of a test is usually more of a reflection of the sample size. If there are non-constant HRs, the choice boils down to your desired application and the power of the sample.

*

*Use robust standard errors. Within reason, when the hazard ratios are non-constant across time, what the Cox model estimates is the time-averaged hazard ratio which is also useful provided the design is representative of the population of interest. The only issue is that the Cox model does not calculate the correct standard errors using the model based approach. Robust standard errors will generally calculate the same coefficients but different standard errors.
Since they have deprecated the robust=TRUE command, it is done by specifying cluster. For independent data, each observation is its own cluster.


*Increase the models sophistication. Adjust for the interaction between (the log of) time (at risk) and coefficients. This explicitly models the hazard ratio over various levels of time and you can assess and describe how it changes.
