This is a follow-up question to the one posted here.

I have run a proportional hazards model in R using the coxph function. The model includes 3 covariates: treatment (5 levels), sex (2 levels), diet (2 levels), and their interactions:

mdl <- coxph(Surv(Day, Death) ~ treatment*sex + sex*diet + treatment*diet)

My data set contains 700 individuals, with 443 events.

I would like to make pairwise comparisons between survival probabilities of individuals from different treatments for a given sex and diet combination (e.g. compare treatments 2 and 3 for diet = 1, sex = male). I understand that this can be done using the hazard ratio for the two conditions.

The output of the coxph model in R lists hazard ratios (the exp(coef) column) for each condition versus some reference level. However, as my treatment covariate has 5 levels, I'd like to know how to make a direct comparison between any two of these, whilst appropriately adjusting p-values for multiple comparisons. It's not clear how this can be done in R.

I have explored the possibility of using the pairwise_survdiff() function in the survminer package, but this does not allow for interaction terms in the model.

  • $\begingroup$ I assume you've cycled through the levels of the dummy coded reference variable? if not, the code is simply df$variable<-relevel(df$variable, ref="whatever level") $\endgroup$
    – Huy Pham
    Nov 27, 2018 at 10:44
  • $\begingroup$ @HuyPham The problem here (I think) is that I want to compare different treatments (5 levels) for specific values of the diet and sex variables (2 levels). I'm not entirely sure what the output of the coxph model is telling me, specifically in the case of interaction terms. What is the reference case for these? How can I directly compare one combination of variable options to another? $\endgroup$
    – allhands
    Nov 27, 2018 at 11:08
  • $\begingroup$ So the reference category is whatever is NOT in the output. So if you have diet2, diet3,..., diet5 then diet1 is your reference category, that is what all the other categories are being compared to. Same for sex. So for an interaction you might have the following diet2:female, diet3:female,...,diet5:female. This is saying the hazard ratio in/decreases by X% for being in diet2 compared to diet1 IF they are female. The diet reference category's effect is the effect of female because that's when there is no comparison of groups anymore. $\endgroup$
    – Huy Pham
    Nov 27, 2018 at 11:13
  • $\begingroup$ this website should be pretty useful. It's in stata but everything's basically the same. $\endgroup$
    – Huy Pham
    Nov 27, 2018 at 11:19
  • $\begingroup$ @HuyPham thanks for the link. Just on your previous comment - that makes sense. But what if I want to compare diet2:female, diet3:female when treatment = 1 (as opposed to treatment = 2)? There is no 3-way interaction term in the model. $\endgroup$
    – allhands
    Nov 27, 2018 at 11:24


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.