# coxph ran out of iterations and did not converge

Yes, I have checked that previous answers to "Ran out of iterations..." questions do not solve my problem.

I have fault data on Firefox, 899 faults and 1395 (estimated) censored faults. The censoring all happens on one of half a dozen start days and half a dozen end days (the initial/final release of a version).

library(survival)

f_sur=Surv(ff_usage$start, ff_usage$end, event=ff_usage$event) plot(survfit(f_sur ~ 1)) f_cox=coxph(f_sur ~ total_usage+cluster(fault_id), data=ff_usage) The Kaplan-Meier curve looks about right. total_usage is an estimate of the number of Firefox users up until the fault is reported. This is very time dependent and so each fault timeline is broken up into 7 day intervals clustered on fault_id; unsplit original. The dependency on total_usage (or its log) could be close to 1 (I am hoping for one or the other). I have tried setting init and increasing iter.max; also strata(src_id) and subsetting on src_id. Most of the start/end times are estimated and have a regular interval, I have tried adding some randomization, e.g., runif(n, -3, 3). No change. All I ever see is: Warning message: In fitter(X, Y, strats, offset, init, control, weights = weights, : Ran out of iterations and did not converge Suggestions for things to try welcome. • You should have posted the fit. It is a warning rather than an error so there was a result and if may be informative to look at the coefficients. In this case the se.est is zero. – DWin Aug 10 '13 at 16:25 ## 2 Answers This may be a case where, as the coxph() documentation page puts it, "the actual MLE estimate of a coefficient is infinity" so that "the associated coefficient grows at a steady pace and a race condition will exist in the fitting routine." In particular, close interrelations of the start / end times with the total_usage variable may be the problem here. When I have problems with a continuous predictor variable like your total_usage in survival analysis, I examine a split of the continuous variable at the median. Look at survival curves from your data based on a split of total_usage at its median value of$5866.2\$ (the coxph() for this simple analysis also didn't converge):

plot(survfit(f_sur~(total_usage > 5866.2),data=ff_usage))

Looks like almost all censoring times and events for the low total_usage cases are before something like time=700, while almost all events and censoring times for the high total_usage subset are greater than that time. Also, examining:

summary(survfit(f_sur~(total_usage > 5866.2),data=ff_usage))

may provide some insight. My data sets are typically much smaller than this, but I have run into related problems in Cox analysis with "a dichotomous variable where one of the groups has no events," so that hazard ratios are ill-defined.

Hope this helps point you in the right direction.

• +1, nice answer! Welcome to CV, we hope we'll see more of these in the future. – gung - Reinstate Monica Aug 6 '13 at 17:22
• Thanks for the suggestions @gung. Randomising the fault report times over a time span (instead of reguar intervals) does not make the problem go away. A histogram of the censored fault total_usage looks roughly Poisson with a peak around 3,000, while for actual faults the Poisson look is a bit rougher and peaks around 4,000. – Derek Jones Aug 7 '13 at 14:15
• Except the result is not congruent with the notion that a coefficient goes to infinity. – DWin Sep 20 '18 at 4:37

The result:

f_cox
Call:
coxph(formula = f_sur ~ total_usage + cluster(fault_id), data = ff_usage)

coef exp(coef) se(coef) robust se    z p
total_usage -0.00407     0.996        0         0 -Inf 0

This is the warning message. Notice it is not a failure to converge:

Warning message:
In fitter(X, Y, strats, offset, init, control, weights = weights,  :
Loglik converged before variable  1 ; beta may be infinite.

When considering problems with survival analyses where estimates blow up, it's often useful to look at tabular displays. (in this case the "explosion" is to the small side rather than the high side.) Consider looking at event crossed with your clustering variable:

with(ff_usage, table(event, fault_id))

Every cluster (all 2294 of them) has event count either 0 or 1, so the algorithm ends up doing a simple tabulation and a zero estimate for sd.coef. It's pretty clearly fake data with not much effort at inserting randomness at least for the counts. The numbers at risk ascend in integer sequences along with "fault_id".

with(ff_usage, table(event, fault_id))
fault_id
event   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23
0   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  15  16  17  18  19  20  21
1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
fault_id
event  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46
0  22  23  24  25  26  27  28  29  30  31  31  32  33  34  35  36  37  38  39  40  41  42  43
1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
fault_id
event  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69
0  44  45  46  47  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  63  64
1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
fault_id
event  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92
0  65  66  67  68  69  70  71  72  73  74  33  35  36  37  39  40  42  43  45  46  47  49  50
1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
#--- truncated output which otherwise goes on for pages since there are 2294 "clusters".

Running this in 2018 produces a mildly different result, although the cause of this issue remains inadequate counts within clusters:

Call:
coxph(formula = f_sur ~ total_usage + cluster(fault_id), data = ff_usage)

coef exp(coef)  se(coef) robust se     z      p
total_usage -1.67e-03  9.98e-01  3.84e-05  1.05e-04 -15.9 <2e-16

Likelihood ratio test=6641  on 1 df, p=<2e-16
n= 174353, number of events= 899
• I think you have missed the point @DWin, total_usage is time dependent and slicing up the total time into intervals, all censored until the last one, is the solution. Of course fault_id is an ascending integer sequence, it was created to identify the clusters created by the slicing process. – Derek Jones Aug 10 '13 at 18:57
• I may have "missed the point" but I think I identified why the design was producing a degenerate result. – DWin Aug 10 '13 at 19:07
• The data is certainly not fake. So the question of how to get get coxph to converge still stands. – Derek Jones Aug 10 '13 at 19:40
• It did converge. The warning is just issued because it converged "too well". You said you were expecting a result of 0 or 1 and got a result of -0.00407 with exp(estimate) of 0.996, albeit with an estimate for se(estimate) of 0, so I remain unclear where the problem lies. – DWin Aug 10 '13 at 22:41
• "Ran out of iterations and did not converge" does not sound like convergence to me. No idea why this is a warning and not an error. – Derek Jones Aug 11 '13 at 1:38