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.


1 Answer 1


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
 # 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)

 # 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)

 # 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))


 #      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


 # 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.

Additional notes

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.

At one level this is a terminology issue: the "stratification" achieved in Models 2/3 is akin to what we loosely call stratified estimates in epidemiology (multiple estimates of effect, specific to each stratum), whereas the stratification in Model 1 is perhaps better considered as stratified estimation of a single effect (that is, the estimation proceeds within strata and is then combined).

For more detail on these issues you might like to see Kleinbaum and Klein's Survival analysis: a self learning text, which is reasonably accessible to a non-mathematician and covers most of these points.

  • $\begingroup$ thank you for your answer, however it would be useful (to the purpose of this question which will may become a standard for other similar in the future) to add, if you want, a brief interpretation of the coefficients in Model 2 - in other words, the unsignificant p in the interaction means that the effect do not differ across sexes; but how to interpret, in this scenario, the coefficient of age alone? $\endgroup$ Nov 13, 2021 at 10:21
  • $\begingroup$ I will add a brief note when time permits. The interpretation of coefficients in interaction terms should be covered elsewhere in Stack Exchange, and so this answer concentrated on the differences between these three scenarios. Also this particular analysis example was intended to be reproducible and exemplary of the method, rather than chosen as a prime example of the interaction phenomenon, so it's perhaps not helpful as a canonical example of that issue. $\endgroup$ Nov 14, 2021 at 21:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.