# Cox regression with multiple factors [R]

I'm running a Cox regression with three different groups...

> data <- source("https://pastebin.com/raw/padE6C4W")\$value
> coxph(Surv(time, event) ~ group, data = data)
# Call:
# coxph(formula = Surv(time, event) ~ group, data = data)
#
#          coef exp(coef) se(coef)     z       p
# group2 0.2936    1.3412   0.2249 1.306 0.19170
# group3 0.8417    2.3203   0.3099 2.716 0.00661
#
# Likelihood ratio test=6.89  on 2 df, p=0.03193
# n= 329, number of events= 139
#    (46 observations deleted due to missingness)


... and, it's expected, my group1 is used as a reference. I'd like to obtain a figure like this, considering all groups (something similar to their group ISS):

I don't know how they consider the reference group, so I tried this approach to calculate HR for each group:

> data$$group_1 <- factor(ifelse(data$$group == 1, 1, 0))
> data$$group_2 <- factor(ifelse(data$$group == 2, 1, 0))
> data$$group_3 <- factor(ifelse(data$$group == 3, 1, 0))

> coxph(Surv(time, event) ~ group_3 + group_2 + group_1, data = data)
# Call:
# coxph(formula = Surv(time, event) ~ group_3 + group_2 + group_1,
#     data = data)
#
#            coef exp(coef) se(coef)     z       p
# group_31 0.8417    2.3203   0.3099 2.716 0.00661
# group_21 0.2936    1.3412   0.2249 1.306 0.19170
# group_11     NA        NA   0.0000    NA      NA
#
# Likelihood ratio test=6.89  on 2 df, p=0.03193
# n= 329, number of events= 139
#    (46 observations deleted due to missingness)

> coxph(Surv(time, event) ~ group_1 + group_2 + group_3, data = data)
# Call:
# coxph(formula = Surv(time, event) ~ group_1 + group_2 + group_3,
#     data = data)
#
#             coef exp(coef) se(coef)      z       p
# group_11 -0.8417    0.4310   0.3099 -2.716 0.00661
# group_21 -0.5481    0.5780   0.2572 -2.131 0.03307
# group_31      NA        NA   0.0000     NA      NA
#
# Likelihood ratio test=6.89  on 2 df, p=0.03193
# n= 329, number of events= 139
#    (46 observations deleted due to missingness)


But I've got some questions:

• Is this the best approach taking account that I've split my variable and performed a multivariable instead a univariable?
• Why, depending of groups order (1:3 or 3:1), the HR changes? I'd expect to keep it constant...
• Why there's always a group that never gets any statistics?

Sorry if these are a very naïve questions. Thanks in advance.

• For Q2 try adding the coefficients for group 2 from your two models (having made both positive) and compare with the ones for group 1 and group 3. For Q3 it is the same reason as for your first model. For Q1 it depends on your scientific question. Jan 28 '21 at 15:06

The coefficients of regression models are typically expressed in terms of some baseline condition. In ordinary least squares, the estimated value of the outcome variable when all predictors are at their baseline values (e.g., the reference level of a categorical predictor, or 0 for a continuous predictor) is the intercept. With a typical "treatment" coding of a regression model, the regression coefficients represent differences from the intercept value associated with each of the predictors. So the number of coefficients returned by a regression model for a categorical predictor is one less than the number of levels of the predictor. The situation with the reference level of the predictor is included in the intercept, which gets to your third question. You do get "statistics" for that reference level, combined with all the other baseline predictor values, but that's in terms of the point estimate and variance of the model intercept and the covariance of that intercept estimate with the estimates of the other coefficients.

The same is true of Cox regressions, except that instead of an intercept you have a baseline hazard function, a function that represents survival as a function of time at reference values of predictors. The Cox regression coefficients represent differences in log-hazard from that baseline, and the hazard ratios are those coefficients exponentiated. Continuing from the argument above, if you change the reference level of a categorical predictor then you are changing the baseline hazard against which the effects of the other levels of the predictor are estimated. The same thing happens in linear regression, too. Thus the hazard ratios relative to a baseline hazard change when the baseline hazard changes. That gets to your second question.

In terms of your first question, you seem to be misunderstanding the plot that you show. The hazard ratios labeled "Sex" aren't those for male versus female. They are hazard ratios for "MRD undetectable" versus "MRD persistent," broken down into male and female subgroups. As another answer points out and as discussed above here, you can't specify separate coefficients for all levels of a categorical predictor in any regression; its reference level is subsumed in the intercept in ordinary least squares or in the baseline hazard in a Cox model.

It's still possible, however, to get predictions (e.g., 3-year survival) for all levels of a categorical predictor, based on the baseline hazard and specified values of other covariates in the model. That might accomplish what you want, with the prediction functions available for survival models. The emmeans package in R provides a convenient way to get estimated outcomes associated with all values of categorical predictors for a wide range of regression models, which might simplify things somewhat for you.

You have too many parameters in the model. If you have 3 group categories, you can only make 2 indicator variables. If you want the reference group to be group 1, then create two indictor variables (group ==2) and (group==3). The baseline hazard function will apply to group 1, the proportional hazards model will then estimate how groups 2 and groups 3 change from the baseline hazard group 1.

Thanks for your comments, they help me to clarify a bit my ideas and see that my questions was wrongly asked... because, as @John L says, to see the proportional hazard ratio I need a baseline (control) group.

That author of this image (I guess) is doing (here is the source from I took the graph, sorry) is determine how the presence of a feature is affecting to another clinical ones, that is, see if, for example, undetectable MRD is favouring or protecting the event.