# Survival analyses: how to relevel contrasts for multiple group comparisons?

I have trouble understanding how do I compare multiple groups in a single survival analysis. When I only have two groups of patients, I simply perform the Kaplan-Meier estimator. However, if there are more than 2 exposure levels, the inference is not so straightforward.

If you look at the colon data (available in R, survival plot below), you will find 3 different cohorts (rx column): Obs (observation), Lev (Levamisole treatment) and Lev+5FU (Levamisole + 5-FU treatment).

The plot demonstrates clear evidence that:

• Obs and Lev cohorts have similar survival levels
• Lev+5FU cohort has higher survival level than Obs and Lev ones

However, since a single P is given (a single value is usually given in multiple online tutorials), I guess that this P represents the higher survival in Lev+5FU compared to others, just like the P given with an ANOVA test states that data is, somehow, different and you need to do additional tests to figure out what is different. Is my thought correct? I actually would like a better measurement of this. What test can I perform to measure survival differences between:

• Obs versus Lev
• Obs versus Lev+FU
• Lev versus Lev+FU

Should I just perform multiple Kaplan-Meier estimators and then adjust P for multiple tests (Bonferroni, BH or whatever) or is there any alternate solution?

I also experienced that if I perform a Cox PH, I do not get full information.

coxph(Surv(time,status) ~ rx, data = colon)


OUT:

Call:
coxph(formula = Surv(time, status) ~ rx, data = colon)

coef exp(coef) se(coef)     z       p
rxLev     -0.0209    0.9793   0.0768 -0.27    0.79
rxLev+5FU -0.4410    0.6434   0.0839 -5.26 1.5e-07

Likelihood ratio test=35.2  on 2 df, p=2.23e-08
n= 1858, number of events= 920


I guess that the two given P compare Obs versus Lev and Obs versus Lev+5FU right? I my case, no "control" cohort can be used as a survival reference.

Can someone help me with that?

• With three groups you only have two independent comparisons so you do have "full information". – mdewey Mar 4 '17 at 13:48
• You have to either use some group as a control or reference group or else compare them all to the grand mean. – Peter Flom Mar 4 '17 at 13:58
• @mdewey : can you explain why there only is two independant comparisons with 3 groups? I want to compare all group combinations. – TomL Mar 7 '17 at 10:22
• @peter: the colon data was an exemple and as I mentionned: I have no "control" cohort. Suppose I have a tumor cohort and I want to see survival differences between tumor locations (extremities, internal trunk, head and neck). How do I deal with that? – TomL Mar 7 '17 at 10:29
• Suppose you have three people: Amy is 5cm taller than Barbara, Barbara is 3cm shorter than Chloe. The fact that Amy is 2cm taller than Chloe is not an independent fact as it is already contained in the previous two facts. – mdewey Mar 7 '17 at 11:57

## 3 Answers

You can compare the three Kaplan-Meier survival curves. One of your curves is clearly statistically superior to the other two. Those other two are not distinguishable. Another statistical approach could be to use parametric survival distributions such as Weibull or Negative Exponential et cetera. The alternative approach would be to produce smooth curves but not lead to a statistically significant difference between the lower two curves. The fit would be by maximum likelihood, taking account of any right censoring

When categorical covariates are modeled as exposures, R codes the first category as a referent group. The remaining categories are modeled as contrasts whose coefficients are log hazard ratios. To relevel categorical covariates, factor can be called again to permute the order of factor levels. An example of a clinical trial where the control group is coded as Obs is given here:

coxph(Surv(time,status) ~ rx, data = colon)
levels(colon$rx) # c("Obs","Lev","Lev+5FU") coxph(Surv(time,status) ~ factor(rx,levels=c("Obs","Lev","Lev+5FU")), data = colon)  In the second run of coxph, the factor level is not changed and identical estimates are produced each time. These compare the relative survival of participants in the Lev and Lev+5FU arms to control (Obs). If I wanted contrasts with Lev+5FU as the referent group, and compare these survival outcomes with Lev and Obs : coxph(Surv(time,status) ~ factor(rx,levels=c("Lev+5FU","Lev","Obs")), data = colon)  The log-hazard ratios, 95% CIs, and p-values are calculated with Lev+5FU as the referent group. You may be interested in learning about contrasts. In R, the default is contr.treatment which says: contr.treatment contrasts each level with the baseline level (specified by base): the baseline level is omitted. Note that this does not produce ‘contrasts’ as defined in the standard theory for linear models as they are not orthogonal to the intercept. Which seems to accurately define the source of your confusion here. Contrast (no pun intended) this with the help file for contr.helmert: contr.helmert returns Helmert contrasts, which contrast the second level with the first, the third with the average of the first two, and so on. contr.poly returns contrasts based on orthogonal polynomials. contr.sum uses ‘sum to zero contrasts’. This is very close to what you describe later as what you want. We only conduct two inferences here. First compare Lev and Obs alone. If you don't observe any difference between Lev and Obs, pool them and compare against Lev+5FU. See here: contrasts(colon$rx) <- contr.helmert(3)
coxph(Surv(time, status) ~ rx, data=colon)


gives:

Call:
coxph(formula = Surv(time, status) ~ rx, data = colon)

coef exp(coef) se(coef)     z       p
rx1 -0.0104    0.9896   0.0384 -0.27    0.79
rx2 -0.1435    0.8663   0.0250 -5.75 9.2e-09

Likelihood ratio test=35.2  on 2 df, p=2.23e-08
n= 1858, number of events= 920


Which quantitatively describes exactly what you observed: Lev shows no difference from Obs (p=0.79). But the combined Lev, Obs cohort show significant differences in survival HR=exp(-0.44), p < 0.001. Because of the coding the coefficients are off by a factor of 2 from the contr.treatment. Note the identical z and p values.

You can reproduce very similar output by running:

coxph(Surv(time, status) ~ I(rx == "Lev"), data=colon, subset= rx != "Lev+5FU")
coxph(Surv(time, status) ~ I(rx == "Lev+5FU"), data=colon)


The differences arise because modeling subsets gives rise to a different baseline hazard function, and because there are ties in these data. The helmert contrast is probably a better bet.