Note that the p values from
cox.zph() shouldn't be thought of in the same way that you use p values in standard tests of null hypotheses. The burden is on you to demonstrate that PH holds; it's not enough to show that you can't rule out the PH null hypothesis because p > 0.05. A p value on the order of 0.1, as you have for some of your
cox.zph() results, is certainly consistent with a violation of the PH assumption.
On the other hand, with 5001 events it's quite possible to find a statistically significant non-proportionality that is of little practical significance. As Therneau and Grambsch put it in Section 6.5 of their book:
A "significant" nonproportionality may make no difference to the interpretation of a data set, particularly for large sampIe sizes.
They go on to describe a study with 771 events, in which the
cox.zph() test gave significant p values for 2 out of 3 predictors with Global p = 1.95e-09. This result might be taken to require time-dependent estimates of each coefficient, $\hat\beta(t)$, instead of a time-independent estimate $\beta$. Examining plots of scaled Schoenfeld residuals over time, however, led to the following conclusion about the association between age and risk of death:
There is a profound effect of age on the risk of death, and the test shows the effect to be significantly nonproportional. However, the smoothed residual plot shows that the variation in $\hat\beta(t)$ is small relative to $\hat\beta$.
Looking at the results instead of just accepting a p value is always a good idea. Pull back to a model omitting all interaction terms except for those that are of substantive interest to the goals of your study. See whether the "significant" trends in residuals over time are of sufficient magnitude to matter.
Here are a few more items to consider if you still find the lack of PH to be an obstacle to your study.
First, from their names it looks like some of your predictors are dummy-coded. For example, you have
lab 4. If these represent mutually exclusive binary representations of a categorical predictor called
lab has 4 levels, then you might be getting into problems with linearly dependent predictors. I don't know whether that would cause a problem with
coxph() but it's possible that
cox.zph() can't handle that and leads to a large global $\chi^2$. The coding of
event, however, suggests that you have taken that into account (perhaps using using
event1 as the reference for
event, with similar precautions for
lab and other potential categorical predictors) so I note this mostly for others who might read this thread.
If that's not the case, next examine whether the linearity assumptions for your continuous predictors, evidently
bmi, are valid. Failure to satisfy the linearity assumption can lead to problems with satisfying the proportional hazards (PH) assumption, too. See this page for one example with a plausibility argument for why this can be the case. So attend to the linearity assumption. Modeling such a predictor with restricted cubic splines, for example, might solve the PH problem and make further adjustments unnecessary.
Failing that, consider stratification on predictors that don't meet the PH assumption instead of adding interaction terms.* That is trivial for categorical predictors. For a continuous predictor like
age still doesn't meet the PH assumption despite your attempts at transforming it to meet the linearity assumption, and
age is included as a control covariate that is not of primary interest to you, you could subset by age into a reasonable number (order of 5) of categories and include
strata(age.category) as a predictor instead of
age. You won't get a hazard ratio for
age but this can solve the PH problem if
age is to blame.
If the predictor giving you problems with PH is of primary interest then you can model it with a time-dependent coefficient. Make sure, however, that you do this properly as you can find seemingly reputable references to completely inappropriate ways to proceed. See this page for some discussion and links to correct details.
Finally, for further details on Cox PH modeling and how to deal with problems that arise, get and read carefully a copy of Harrell's Regression Modeling Strategies. You will find that very valuable if you will be doing a lot of regression modeling (Cox or otherwise), as your comment suggests will be the case. It ties in nicely to the modeling tools provided by Harrell's
rms package in R.
*At first thought, it would seem that if you have a continuous predictor that isn't meeting the PH assumption, and then you add interactions of that predictor with others, you would make the overall PH problem worse rather than better.