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