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.