3
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ Please edit the question to say more about the nature of the censoring. For example, did you just stop recording time-to-eat at some fixed point in time, so that all censored values are simply something like 'greater than 2 minutes'? If so, there might be an alternative to a Cox model. Please provide that information by editing the question, as comments are easy to overlook and can be deleted. $\endgroup$
    – EdM
    Commented Oct 5, 2022 at 16:52
  • $\begingroup$ Thanks - that information is now included. $\endgroup$
    – Picapica
    Commented Oct 5, 2022 at 17:05

1 Answer 1

3
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Thanks for your reply. I arbitrarily divided early and late responders and there is indeed a difference early on. I tried a tobit model but it wouldn't run because it was singular, but, as you say, it wouldn't model the early/late difference anyway. I'll keep looking into the other models. Is a hurdle model not just for count data? $\endgroup$
    – Picapica
    Commented Oct 10, 2022 at 18:53
  • $\begingroup$ @Picapica "cure models" in survival analysis are a conceptual equivalent to hurdle models: a cure model combines the probability of never having the event along with the time-to-event analysis for those who have it. In your case, I suspect that a simple step-function in time for the treatment coefficient, as described in Section 4.1 of the R time dependence vignette, will do what you need. $\endgroup$
    – EdM
    Commented Oct 10, 2022 at 20:12
  • $\begingroup$ This looked really promising! I followed their example, but unfortunately I get a warning message when using my own data: "Error in coxme.fit(X, Y, strats, offset, init, control, weights = weights, : No starting estimate was successful". I'll probably have to make a new question, though I'm not confident I'll be able to sort it. $\endgroup$
    – Picapica
    Commented Oct 11, 2022 at 16:07
  • $\begingroup$ @Picapica playing with the starting estimates in coxme would probably fix this, but it might be simpler to use the standard coxph() function with a cluster() term for the individuals. That provides robust estimates of the coefficient standard errors to account for the repeated measures. It sounds like you have a reasonably complete data set so that missing data (which mixed models like coxme can handle better) isn't a problem. See discussion on GEE vs mixed models. $\endgroup$
    – EdM
    Commented Oct 11, 2022 at 17:00
  • $\begingroup$ I wasn't able to find any helpful starting estimates for the coxme model. Running coxph() works though doesn't show me the pattern I was expecting based on running the separate models (it says there is no difference between early and late - but maybe I am wrong to give weight to that) and the PHA is still violated (this time for TrialType:strata(tgroup) ). Another thing is that I get NAs for both TrialTypeTreatment:tgroup estimates and I'm not sure why - is it supposed to be two baselines? Also I was hoping to get a random effects estimate out. I really appreciate all the help by the way! $\endgroup$
    – Picapica
    Commented Oct 11, 2022 at 18:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.