# Does fitting Cox-model with strata and strata-covariate interaction differ from fitting two Cox models?

In Regression Modeling Strategies by Harrell (second edition) there is a section (S. 20.1.7) discussing Cox models including an interaction between a covariate whose main effect on survival we want to estimate as well (age in the example below) and a covariate whose main effect we do not want to estimate (gender in the example below).

Concretely: suppose that in a population the (unknown, true) hazard $$h(t)$$ follows the model

$$h(t) = \begin{cases} h_f(t) \exp(\beta_1 \textrm{age}), & \textrm{for female patiens} \\ h_m(t) \exp((\beta_1 + \beta_2) \textrm{age}), & \textrm{for male patiens} \end{cases}$$ where $$h_f$$, $$h_m$$ are unknown, true, not to be estimated baseline hazard functions and $$\beta_1$$, $$\beta_2$$ are unknown, true parameters to be estimated from the data.

(This example is taken almost literally from the book.)

Now Harrell remarks that the above situation can be rewritten as the stratified Cox model model 1:

$$h(t) = h_{\textrm{gender}}(t) \exp(\beta_1 \textrm{age} + \beta_2 X)$$ where the 'interaction term' $$X$$ equals zero for females and age for males. This is convenient because it means we can use the standard technique for estimating $$\beta_1$$ and $$\beta_2$$.

Now for the question. Suppose that two researchers A and B are given the same sample of patients drawn from the population described above. Researcher A fits model 1, obtaining estimates $$\hat{\beta}_1$$, $$\hat{\beta}_2$$ for the true parameters $$\beta_1, \beta_2$$ together with confidence intervals.

Researcher B takes the more naive approach of fitting two ordinary (i.e. unstratisfied) Cox-models: model 2a: $$h(t) = h_f(t)\exp(\gamma_1 \textrm{age})$$ on the female patients in the sample only and model 2b: $$h(t) = h_m(t)\exp(\gamma_2 \textrm{age})$$ on the male patients in the sample only. Thus obtaining estimates $$\hat{\gamma_1}$$, $$\hat{\gamma_2}$$ of the true parameters $$\beta_1, \beta_1 + \beta_2$$ respectively, together with confidence intervals.

Question:

• Are these estimates necessarily the same (in the sense that $$\hat{\beta}_1 = \hat{\gamma}_1$$, $$\hat{\beta}_2 = \hat{\gamma}_2 - \hat{\gamma}_1$$)? (Recall that both researchers look at the same data.)
• Are the confidence intervals necessarily the same?
• Does it make any sense to say that researcher A has a psychological advantage over researcher B in the case that $$\beta_2 = 0$$, because researcher A is then more likely to suspect that and switch to estimating the more parsimonious model $$h(t) = h_{\textrm{gender}}(t)\exp(\beta_1 \textrm{age})$$?

With models where each parameter has to be estimated (like Ordinary Least Squares), it is possible to create a situation where two separate models have the same estimates of a single one with an interaction term. For example, we could have: $Y_M=\alpha_M+\beta_M*age$, $Y_F=\alpha_F+\beta_F*age$ summarized by: $Y=\lambda+\lambda_F*F+\gamma*age+\gamma_F*F*age$, so that you could directly estimate the gender difference both in intercept and in slope. In fact: $\alpha_M=\lambda, \beta_M=\gamma, \alpha_F-\alpha_M=\lambda_F,\beta_F-\beta_M=\gamma_F$. In that case, I agree with you that the unique model would allow to have an immediate idea on the gender difference (given by the interaction parameters, $\lambda_F$, since the slope difference has a clearer interpretation, and your question refers to that). However, with the Cox model things are different. First of all, if we don't include gender in the regression there may be a reason, i.e. that it doesn not fulfill the proportional hazard assumption. Also, if we build a unique model with gender as an interaction term, we are assuming a common baseline hazard function (unless I misunderstood the meaning of $h_{\textrm{gender}}(t)$), while the two-separate-models approach allow for two separate baseline hazard functions, thus different models are implied.

See, for example, the Chapter "Survival Analysis" from Kleinbaum and Klein, 2012, Part of the series Statistics for Biology and Health.

If you write it exactly as you specified the models in the post, then yes, they will see the same, i.e. Researcher A:

 coxph(Surv(time_1, event)~ age_ + age_*gender + strata(gender), data =df)


In my simulation that gave age coefficient of 0.10 and age x gender of 0.11 (vs the modelled parameters of 0.1 and 0.2)

Researcher B:

coxph(Surv(time_1, event)~ age_ , data =df[df$gender==0, ]) coxph(Surv(time_1, event)~ age_ , data =df[df$gender==1, ])


Here, it gives 0.10 for the first and 0.21 for the second

However, if researcher A just writes something like below, without "strata",

coxph(Surv(time_1, event)~ age_ + age_*gender, data =df)


then no - this is what @Federico mentioned above, that the baseline hazard function will be assumed the same for both genders in this case.

Full code with simulated example is here:

#simulating population where gender =0 corresponds to Cox beta of 0.2 for age #variable and gender =1 corresponds to Cox beta of 0.1
N=5000
df = data.frame(age = round(runif(N, 25,75),1),
gender = rbinom(N,1,0.5))
df$$age_ = (df$$age - 50) /10
df$$event_time = 0.01+ rexp(N, exp(0.2*df$$age_*(df$$gender==1) + 0.1* df$$age_*(df$$gender==0))) #event is an exponential variable df$$event = ifelse(df$$event_time>1, 0,1) #all censored at time =1 sum(df$$event)/1000
df$$time_1 = ifelse(df$$event_time>1, 1, df$event_time) #time of event or 1 for censored #this is your first equation, which allows for separate baseline function by gender # AND it allows for a different coefficient for the second gender # (without age_*gender the model will allow different baselines, but not the coefficients coxph(Surv(time_1, event)~ age_ + age_*gender + strata(gender), data =df) #this is what researcher B will do: coxph(Surv(time_1, event)~ age_ , data =df[df$$gender==0, ]) coxph(Surv(time_1, event)~ age_ , data =df[df$$gender==1, ]) #results A: coxph(Surv(time_1, event)~ age_ + age_*gender + strata(gender), data =df) Call: coxph(formula = Surv(time_1, event) ~ age_ + age_ * gender + strata(gender), data = df) coef exp(coef) se(coef) z p age_ 0.10 1.11 0.02 6 0.00000002 gender NA NA 0.00 NA NA age_:gender 0.11 1.11 0.02 4 0.00001901 #results B - gender ==0 Call: coxph(formula = Surv(time_1, event) ~ age_, data = df[df$gender ==
0, ])

coef exp(coef) se(coef) z          p
age_ 0.10      1.11     0.02 6 0.00000002

Likelihood ratio test=32  on 1 df, p=0.00000002
n= 2465, number of events= 1519

#results B, gender ==1
coxph(formula = Surv(time_1, event) ~ age_, data = df[df\$gender ==
1, ])

coef exp(coef) se(coef)  z                   p
age_ 0.21      1.23     0.02 12 <0.0000000000000002

Likelihood ratio test=144  on 1 df, p=<0.0000000000000002