# Interpretation of R output for stratified cox-ph model

For the following model:

model <- coxph(formula = Surv(time, status) ~ treatment * sex + strata(sex), data = data)


this is (part of) the model summary in R:

coef exp(coef)
treatment1 -0.15 0.86
sex_f NA NA
treatment1:sex_f -0.35 0.70

If I understand correctly, in a stratified model, the coefficients cannot be interpreted anymore as a hazard ratio, since the baseline is not the same for the two sexes. So based on the table alone, I cannot tell the hazard ratio of a female receiving treatment compared to a male receiving treatment? And how to get the female-specific log hazard ratio for females with treatment vs. females without treatment, is it log(HR) = -0.15 + -0.35, or is it log(HR) = -0.35?

Many thanks!

When you use strata(sex), coxph should estimate a separate hazard for males and females. There is no need to include sex in the model formula if you are going to estimate the hazard within strata defined by sex.

To get the hazard ratio for females receiving treatment over males receiving treatment, you can just omit strata(sex) from the model.

Here is an example using lung from the survival library.

library(survival)
library(marginaleffects)
set.seed(0)
lung$$sex <- factor(lung$$sex,labels = c("Male", "Female"))

lung\$treatment <- rbinom(nrow(lung), 1, 0.5)

model <- coxph(Surv(time, status) ~ treatment*sex, data=lung)
summary(model)
#> Call:
#> coxph(formula = Surv(time, status) ~ treatment * sex, data = lung)
#>
#>   n= 228, number of events= 165
#>
#>                         coef exp(coef) se(coef)      z Pr(>|z|)
#> treatment            0.05069   1.05200  0.18980  0.267   0.7894
#> sexFemale           -0.54869   0.57771  0.24389 -2.250   0.0245 *
#> treatment:sexFemale  0.02753   1.02792  0.33499  0.082   0.9345
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#>                     exp(coef) exp(-coef) lower .95 upper .95
#> treatment              1.0520     0.9506    0.7252    1.5261
#> sexFemale              0.5777     1.7310    0.3582    0.9318
#> treatment:sexFemale    1.0279     0.9728    0.5331    1.9820
#>
#> Concordance= 0.578  (se = 0.023 )
#> Likelihood ratio test= 10.78  on 3 df,   p=0.01
#> Wald test            = 10.22  on 3 df,   p=0.02
#> Score (logrank) test = 10.47  on 3 df,   p=0.01


Created on 2024-06-23 with reprex v2.1.0

From here, you can use the log hazard ratios (coef in the output) to create the whatever contrast you like. The linear predictor for this model looks like

$$\beta_1I(Sex=F) + \beta_2I(Txt=1) + \beta_3 I(Sex=F \cap Txt=1)$$

Here, $$I$$ is an indicator function which is equal to 1 when the argument is true, else 0. The linear predictor for a male recieving treatment would be

$$\beta_2$$

The linear predictor for a remale recieving treatment would be

$$\beta_1 + \beta_2 + \beta_3$$

Hence the log odds ratio for females recieving treatment vs males recieving treatment would be

$$(\beta_1 + \beta_2 + \beta_3) - (\beta_2) = \beta_1 + \beta_3$$

We can use marginaleffects to get an estimate and confidence interval for this contrast in the log hazard ratio scale.

nd <- datagrid(model=model, sex=c("Male","Female"), treatment=1)

avg_comparisons(model,
newdata = nd,
variables = 'sex',
type='lp'
)
#>
#>  Term      Contrast Estimate Std. Error     z Pr(>|z|)   S  2.5 %  97.5 %
#>   sex Female - Male   -0.521       0.23 -2.27   0.0234 5.4 -0.972 -0.0704
#>
#> Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
#> Type:  lp
$$$$
`
• Thank you for your reply! The reason for the way I set up the model is to allow sex-specific baseline hazard AND sex-specific treatment effect. (And then I'm just not sure about what the table represents.) In the model you propose, the baseline hazards are the same, no? Commented Jun 23 at 22:02
• @user412691 ah, if you want a separate hazard within strata defined by sex, and a separate treatment effect in each strata, what you’ve done is fine. You can use the second code chunk I’ve provided to get the contrast you want. Commented Jun 23 at 22:13
• Thanks! And to compare females with treatment vs. females without treatment, is it log(HR) = treatment:sexFemale or log(HR) = treatment + treatment:sexFemale ? Commented Jun 24 at 6:24