Does violating the proportional hazards assumption always matter? I've seen a couple of posts saying that violating the proportional hazards assumption (PHA) doesn't always matter, but several more suggesting it does and that the model needs to be modified accordingly.
I'd like to know whether I need to worry about it in my case. I'm testing whether the latency of a species to eat depends on which treatment group they are in (a binary category). I'm using coxme in R, as there are repeated measurements for each individual. I have included three other variables - two categorical and one continuous. Running cox.ph() to check the PHA shows that my variable of interest (the treatment group) violates it (p = 0.04). Another variable also does to a greater extent (tiny p-value), so I've stratified by that. Running cox.ph() again shows treatment is now at p = 0.02.
I'm wondering whether this matters. I'm a novice to survival analyses, but my inclination is "no". There is no overall difference between the treatment groups (as expected), the plots of the survival curves show they are so close over the whole time frame that their overlapping isn't weird, and I was only testing to see if there was a difference, not to estimate the hazard. I would have used a linear regression had there not been censored data (the subjects were given 600 seconds to eat, and most ate within a few seconds).
Does it sound like it is OK to proceed as is?
 A: It seems strange that you would find a "significant" violation of proportional hazards (PH) if there weren't some association of treatment with time-to-eat. The key here might be in your statement:

the subjects were given 600 seconds to eat, and most ate within a few seconds.

This suggests that you might have separate states of rapid-eating versus slow-or-never-eating, and that your treatment might have different effects on the two states. The danger in just accepting your initial finding of no treatment difference is that you might miss an interesting true difference between early and late times. If the PH assumption doesn't hold, then your coefficient estimates are event-averaged values that might not represent the data well.
A smoothed plot of scaled Schoenfeld residuals (provided by the R cox.zph() function that tests for PH) could be very informative, as it shows directly how the estimated hazard changes over time. Something as simple as a step-function in time for the treatment coefficient might then show something useful. Section 4.1 of the R time dependence vignette demonstrates how to implement that.
Other thoughts:
Your data might not meet PH but instead follow an accelerated failure time (AFT) model, like a lognormal or log-logistic model.
Depending on the distribution of event times, you might consider a type of hurdle model that models both the probability of being in the rapid-eating group and the time-to-event within that group.
There are other ways to deal with simple censoring like this, in particular tobit models. This page shows how to include a random effect. But that wouldn't deal with the possible early/late difference directly; you would see that difference only in the patterns of residuals. A proportional-odds model also might be considered with this type of censoring, but again a true early/late difference wouldn't be found except perhaps as a violation of the proportional-odds assumption.
