Cox model with factors

Question

I have a simple Cox proportional hazards model, with a small made-up data set. There are three groups, and the data is generated so that the time-to-event is larger in the third group than the other two. Please note that this code is copied from another poster's example (see the following link: KM versus Cox model):

times <- rexp(30)
status <- rep(1, 30)
groups <- rep(c(1,2,3), 10)
dat <- data.frame(groups, times, status)
dat[dat$groups==1, "times"] <-  dat[dat$groups==1, "times"]+6
mod1<- coxph(Surv(times, status) ~ groups, dat)
summary(mod1)

That is fine, but the problem appears when I explicitly make the 'groups' variable a categorical variable:

dat$groups <- factor(dat$groups)
mod1<- coxph(Surv(times, status) ~ groups, dat)
summary(mod1)

I get sensible results when the group variable is numeric, but not when they are converted to factor. It seems like the model, if anything, should be MORE sensible when the groups are categorical.

What is happening? Why when I run the model with categorical 'group' variable do I get the following warning:

In coxph.fit(X, Y, istrat, offset, init, control, weights = weights,
:
Loglik converged before variable  1,2 ; coefficient may be infinite. 

Answer 1

This is similar to perfect separation in logistic regression. When I followed your code after using set.seed(102), the largest value of "times" was 3.35. Adding 6 to those in group1 put their survival times longer than any of the individuals in the other 2 groups.

Treating group as a linear continuous predictor assumes that each extra level of the numeric group has the same extra association with the log-hazard of an event. With 3 groups that poses no problem in solving the score equation to get a regression coefficient. The software starts with an estimated regression coefficient value of 0 (hazard ratio of 1) for the numeric group predictor, then iteratively adjusts that value to best correspond to a single linear association between log-hazard and the 3 numeric values of group. In this case, the result ends up a compromise between the (infinite) log-hazard difference between group1 versus the others and the (lack of) log-hazard difference between the second two groups.

When you converted to a factor and (implicitly) used group1 as the reference level, you created a problem. The regression coefficients for the other 2 groups then represent the log-hazard associated with membership in those two groups, relative to that in group1. Similarly to the above, the software will start with dummy coefficients for the other 2 groups of 0, and then interactively try to find values for each of them that best represent the observed log-hazard relative to that of group1

But all of the events in the other 2 groups happened before any of the events in group1. Those 2 groups have infinitely higher hazards than the reference group, group1. So as the software tries to adjust the initial coefficient estimates to match the data, the coefficients run off toward infinity.

You can get results more like you were expecting if you add a value less than 6 to group1 (adding 3 worked for me) or (maybe) if you set a different group as the reference level of the factor. Alternatively, if you only had 2 groups with non-overlapping survival times and treated them as numeric rather than categorical, you would end up with the same lack of convergence as there is no longer a third group that requires compromise to provide a single regression coefficient.