# Difference in stratification, interaction and subsetting Cox Survival Regression

I'm struggling to understand the actual differences in these three approaches that came to my mind.

The problem is: I want to analyze the impact of age on mortality in my cohort. However, I want to explore whether age has different effect on mortality across sexes.

The approaches that came to my mind are:

1. Performing a stratified Cox-Regression Analysis, such as this:

model1 <- coxph(Surv(days,event) ~ age + strata(sex), data=data):

1. Introducing an interaction term into the Cox-PH analysis, such as this:

model2 <- coxph(Surv(days,event) ~ age*sex, data=data)

1. Fit two different Cox model by sex:

model3a <- coxph(Surv(days,event) ~ age, data=data, subset=(sex=="male"))

model3b <- coxph(Surv(days,event) ~ age, data=data, subset=(sex=="female"))

What is the best approach to do that? Note that I am interested in evaluating if and how the effect of age changes across sexes, not the effect of sex itself on mortality.

Any help (And explaination on why one approach is better than the other) will be appreciated.

To re-iterate, your key question is "evaluating if and how the effect of age changes across sexes, not the effect of sex itself on mortality".

I've included some R code and output below that illustrates this using the lung dataset.

Model 1 is thus not of interest here as it is a way to control for confounding. The stratification variable allows for adjusting for covariates (confounders) that do not meet the proportional hazards assumption (while allowing for inclusion of other covariates that do meet the PH assumption). These stratification covariates are effectively treated as nuisance information, with no estimation of hazard ratios.

Model 2 is the better answer to your question: it allows estimation of hazard ratios for age, and allows for these impacts to differ according to the level of sex. Thus we get a direct answer to your question.

Model 3 in a way aims to estimate the same parameters as model 2, but from two separate models. (That is, if you solve Model 2 to estimate the hazard ratios for age for each level of sex, you get approx. the same answers as from Models 3a and 3b -- I suspect it's not exactly the same due to how baseline hazard is estimated in each instance).

However, with Model 3 you do not get a direct statistical estimate (the interaction term from Model 2) of whether the effect size for age differs by level of sex.

Finally, note that the distinction between Model 2 and Model 3 gets more complex when we start including other covariates (without interaction terms): for example, Model 2 will estimate the impact of these confounders over both levels of sex (that is, assumed independent of sex); whereas Models 3a/3b will estimate the impact of these confounders within each level of sex, which is more akin to having a version of Model 2 which interacts additional covariate(s) with sex.

R code

 # this uses the lung cancer data from the survival package
# in this example, there is little to no evidence for interaction
library(survival)
data(lung)
# convert sex to a factor (important for interaction estimation in Model 2)
lung$$sex_fct <- factor(lung$$sex)

model1 <- coxph(Surv(time, status) ~ age + strata(sex_fct), data=lung)
model1

# Call:
# coxph(formula = Surv(time, status) ~ age + strata(sex_fct), data = lung)

#         coef exp(coef) se(coef)     z      p
# age 0.016215  1.016347 0.009187 1.765 0.0776

# Likelihood ratio test=3.18  on 1 df, p=0.07444
# n= 228, number of events= 165

model2 <- coxph(Surv(time, status) ~ age * sex_fct, data=lung)
model2

# Call:
# coxph(formula = Surv(time, status) ~ age * sex_fct, data = lung)

#                 coef exp(coef)  se(coef)      z      p
# age         0.020343  1.020551  0.011395  1.785 0.0742
# sex_fct2    0.095322  1.100013  1.228718  0.078 0.9382
# age:sex_f2 -0.009689  0.990358  0.019414 -0.499 0.6177

# Likelihood ratio test=14.37  on 3 df, p=0.002441
# n= 228, number of events= 165

model3a <- coxph(Surv(time, status) ~ age, data=subset(lung, sex_fct==1))
model3b <- coxph(Surv(time, status) ~ age, data=subset(lung, sex_fct==2))

model3a

#      Call:
# coxph(formula = Surv(time, status) ~ age, data = subset(lung, sex_fct == 1))

#        coef exp(coef) se(coef)     z      p
# age 0.01906   1.01925  0.01135 1.679 0.0931

# Likelihood ratio test=2.89  on 1 df, p=0.0889

model3b

# Call:
# coxph(formula = Surv(time, status) ~ age, data = subset(lung, sex_fct == 2))

#        coef exp(coef) se(coef)     z     p
# age 0.01071   1.01077  0.01565 0.684 0.494


Notes on model interpretation

Note that this example dataset was chosen to illustrate the modelling approaches, rather than as an illustration of an interaction phenomenon.

A brief note of interpretation, focusing on model 2: the hazard ratio (HR) on the line labelled as age is the HR at the reference level of the other variable (i.e. when sex equals 1, which happens to be males).

This suggests an HR of 1.02 for a one-year difference in patient age among males (confidence intervals can be extracted using standard functions like conf.int() or the broom package's tidy function)

The interaction term for age:sex (HR = 0.99) is the modifier for the age HR when sex==2 (i.e. the differential impact of age for females relative to males). Here this is pretty close to 1, so there is limited evidence in this data for survival prospects by age being differential between females and males. (of course, any real differential would be hard to detect in a sample of this size, so I would treat this example as an indeterminate result).

Ideally one would solve out the product of the age and age:sex parameters: this would be exp(0.020343 + -0.009689) = 1.01 as the direct hazard ratio for a one-unit change in age among females (sex == 2).

If you want to get the confidence intervals for that estimate, you have to do a bit of extra work. I won't cover this here, but as one example you can use the rms:: package and use cph() instead of coxph() to fit the model; and then use contrast() to calculate the direct contrast of interest.

Note that these models still have several assumptions that may or may not be tenable/reasonable: for example, that the impact of age as a continuous variable is approximately linear (on the log-scale): in other words, that a fixed difference in age is associated with a fixed hazard ratio.