0

I'm trying to do a survival analysis using a mixed-effects cox proportional hazards model in R with coxme. My data is - actually - fairly simple. I followed the survival of insect that were kept in groups and exposed to 2 different chemical treatments (A,B) + their combination (A+B), as well as a control (10 groups per treatment). I include the groups as random term and would like to investigate the effects of my chemicals on survival.

I fitted the following model:

fit <- coxme(surv(days.survived, status) ~ A * B + (1|group), data = mydata)

I then tried to check whether the proportional hazards assumption is met, and wanted to use the function cox.zph(), which should actually work. See e.g. this question: coxme proportional hazard assumption

ph.test <- cox.zph(fit)

Unfortunately, this does not work, since R just produces a fatal error and the session is aborted. I have updated all packages just recently and it doesn't help. Out of curiosity, I removed the random term "group" and run the cox.zph again. The results showed, that for one chemical, the proportional hazards assumption is not met. I guess, however, that this result is not really reliable without the random term.

Another thing I tried is using a frailty model with coxph instead (as suggested also in the above answer):

fit2 <- coxph(Surv(days.survived,status) ~ A * B + frailty(group), data = mydata)

and then test.ph2 <- cox.zph(fit2)

The frailty model itself seems to work, but gives me a warning message saying In coxpenal.fit(X, Y, istrat, offset, init = init, control, weights = weights, : Inner loop failed to coverge for iterations 2 4

The cox.zph(fit2) does not work. It gives an error imatr[kk, kk] : subscript out of bounds

Does anyone have an idea what could be the problem here? I'm really desperate and I think it can't be that difficult, since my data is not really complex at all... It just contains the number of days that each individual insect survived, the event status (1 = censored, 2 = dead), the treatment and the group. No covariates, no NAs.

Thank you very much in advance for your help.

AdamO
  • 52,330
  • 5
  • 104
  • 209
Asuka
  • 21
  • 2
  • Sounds like the mixed-effect models (whether with `coxme` or frailty term) isn't fitting properly. Two questions. First, do you have 40 separate groups and 40 separate values of `group`? Second, do you care about the variance among groups, or are you just trying to account for potential within-group correlations in event times? – EdM Jun 04 '21 at 15:20
  • Dear Ed, thanks for your time. My model contains 40 separate groups (10 groups per treatment) and each group has between 30-50 individual values (individual insects with their survival time). At the end of the experiment, almost all individuals were dead, only very few (~ < 10%) were still alive. I do not per se care about the variance, but have to account for it in the model, since there are apparently some position effects. – Asuka Jun 07 '21 at 06:28

1 Answers1

1

One possibility is that the random-effect variance is close to 0. That can lead to problems with convergence during the maximum-likelihood fitting.

If the random-effect variance isn't itself of interest and all you need to do is to account for within-group correlations of event probabilities, you could use a cluster() term for groups instead of a frailty() term or a mixed-effect model. A cluster() term adjusts the covariance matrix of coefficients in the fixed-effect model for those correlations similarly to generalized estimating equations, while frailty() and mixed modeling use maximum (partial) likelihood. This page is an introduction to these different approaches. I suspect that going to a cluster() term will solve your problem.

Two more things.

First, if you actually have "position effects," as you say in a comment, neither a random-effect model nor a cluster term might be adequate. Those account for random differences among groups, with correlations among group members. You might need to model the "position" directly in some way if the position/treatment assignment combinations aren't adequately random. You'll have to apply your knowledge of the subject matter and your experimental setup to evaluate that.

Finally, with a large study it's possible to find a statistically significant violation of proportional hazards (PH) that isn't large enough to have practical significance. That could be the case with the 1500 or so events that you evidently have in your study. Checking for PH is important, but evaluate the magnitude of any violation thoughtfully.

EdM
  • 57,766
  • 7
  • 66
  • 187
  • I will look into your suggestion about the cluster(), thank you! In my case, the position/treatment combinations are randomized, so a mixed-model should be appropriate. Also, I also found a way in which I can check for proportional hazards via cox.zph() without my R crashing: I coded my random term as numeric instead of as factor and suddenly it works without a problem. The model output is otherwise the same, and no violation of prop. hazards is shown. I have, honestly, no idea why this works, and have never encountered advice to code random terms as numeric in cox models, but am happy now! – Asuka Jun 09 '21 at 07:27