# Test for proportional hazards assumption in coxme model - Error and abortion of R session

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.

• 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 4, 2021 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. Jun 7, 2021 at 6:28

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.