"coefficient may be infinite" even with non-zero events

Question

I'm performing survival analysis for a couple of very simple datasets:

dt.pat1 <- data.frame(
  Marker=c(F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,T,T,T,T,T),
  Event=c(F,T,F,F,F,F,F,F,F,T,T,F,F,F,F,T,F,T,T,T,T,T),
  days=c(30,42,164,168,169,196,197,231,234,249,260,331,370,408,454,486,577,101,183,190,314,328)
)

dt.pat2 <- data.frame(
  Marker=c(F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,T,T,T,T),
  Event=c(F,F,F,F,F,F,F,T,T,T,F,F,F,T,F,F,F,F,F,F,F,F,T,T,F,T),
  days=c(41,111,136,150,163,172,175,209,228,245,246,294,297,298,321,357,358,372,441,447,478,505,106,153,155,189)
)

For the first one, it is very straightforward with surv_fit() and coxph():

dt.pat_this <- dt.pat1
fit_obj  <- surv_fit(Surv(days,Event) ~ Marker, data = dt.pat_this)
ggsurvplot(fit = fit_obj, data = dt.pat_this, risk.table=T)
summary(coxph(fit_obj$call$formula, data=dt.pat_this))
# Call:
#   coxph(formula = fit_obj$call$formula, data = dt.pat_this)
# 
# n= 22, number of events= 9 
# 
# coef exp(coef) se(coef)    z Pr(>|z|)  
# MarkerTRUE 1.8005    6.0528   0.7348 2.45   0.0143 *
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# exp(coef) exp(-coef) lower .95 upper .95
# MarkerTRUE     6.053     0.1652     1.434     25.55
# 
# Concordance= 0.712  (se = 0.09 )
# Likelihood ratio test= 6.14  on 1 df,   p=0.01
# Wald test            = 6  on 1 df,   p=0.01
# Score (logrank) test = 7.76  on 1 df,   p=0.005
#

Reporting: HR=6.1, 95% CI: 1.4-25.6, p=0.0143

For the second set (dt.pat2), the plot looks totally fine but Cox PH summary is not so nice:

# Call:
#   coxph(formula = fit_obj$call$formula, data = dt.pat_this)
# 
# n= 26, number of events= 7 
# 
# coef exp(coef)  se(coef)     z Pr(>|z|)
# MarkerTRUE 2.236e+01 5.119e+09 1.397e+04 0.002    0.999
# 
# exp(coef) exp(-coef) lower .95 upper .95
# MarkerTRUE 5.119e+09  1.953e-10         0       Inf
# 
# Concordance= 0.755  (se = 0.08 )
# Likelihood ratio test= 13.1  on 1 df,   p=3e-04
# Wald test            = 0  on 1 df,   p=1
# Score (logrank) test = 22.01  on 1 df,   p=3e-06
# 
# Warning message:
#   In coxph.fit(X, Y, istrat, offset, init, control, weights = weights,  :
#                  Loglik converged before variable  1 ; coefficient may be infinite. 

I previously thought this happens when there's no events in one of the groups, but apparently even having few events can cause problems? What metrics would you recommend reporting in such case? For now, I resort to using coxphf() instead of coxph() - is there any better approach?

Answer 1

This is discussed in the documentation of the coxph() function.

In certain data cases the actual MLE estimate of a coefficient is infinity, e.g., a dichotomous variable where one of the groups has no events. When this happens the associated coefficient grows at a steady pace and a race condition will exist in the fitting routine: either the log likelihood converges, the information matrix becomes effectively singular, an argument to exp becomes too large for the computer hardware, or the maximum number of interactions is exceeded. (Most often number 1 is the first to occur.) The routine attempts to detect when this has happened, not always successfully. The primary consequence for the user is that the Wald statistic = coefficient/se(coefficient) is not valid in this case and should be ignored; the likelihood ratio and score tests remain valid however.

As suggested in the comments, this is really a statistical issue. In the second dataset all of the deaths in the Marker = TRUE group happen before any deaths in the Marker = FALSE group. (There are individuals in the other group who are censored during this time, but none who experience an event.) Therefore the estimated hazard in the Marker = TRUE group is infinitely higher than the hazard in the other group. The coxph() function successfully detects the problem in this case and issues a warning.

As described in the documentation, the p-value of 0.999 in the summary of the second model is incorrect (this is from the Wald statistic). Instead you can report the likelihood ratio test and conclude that the effect of Marker is statistically significant (p < 0.01). Reporting an estimate of the hazard ratio is nonsensical, but you can conclude the hazard is higher in the Marker = TRUE group.