Combining two separate Cox PH models into one model? I have two Cox proportional hazards models (in R), using same outcomes and predictors, one for $n_m$ males and one for $n_f$ females. Is it possible to combine them into one equivalent model over all $n_m+n_f$ individuals, using for example an interaction between the predictors and the sex?
For example, this doesn't work in the sense of having same coefficients as each of the sex-specific models:
library(survival)
n <- 200
d <- data.frame(time=rexp(n), event=sample(0:1, n, replace=TRUE),
   x=rnorm(n), sex=sample(c("m", "f"), n, replace=TRUE))
l1 <- coxph(Surv(time, event) ~ x, data=subset(d, sex == "m"))
l2 <- coxph(Surv(time, event) ~ x, data=subset(d, sex == "f"))
l3 <- coxph(Surv(time, event) ~ x * sex, data=d)

> coef(l1)
     x
-0.1841861
> coef(l2)
     x
0.01391554
> coef(l3)
       x         sexm       x:sexm
-0.009145576  0.060432017 -0.179186702

 A: The two separate Cox PH models assume
$$
\lambda_{i,j}(t)=\lambda_{0i}(t)\exp(\beta_ix_j),\quad t>0.
$$
for $i=\mathrm{f},\mathrm{m}$ and $j=1,\ldots,n_i$. That is, the hazard ratio of two individuals of the same sex is constant since
$$
\frac{\lambda_{i,j_1}(t)}{\lambda_{i,j_2}(t)}=\exp(\beta_i(x_{j_1}-x_{j_2})),\quad t>0,
$$
for $i=\mathrm{f},\mathrm{m}$. The hazard ratio of two individuals of the opposite sex is not necessarily constant since
$$
\frac{\lambda_{\mathrm{f},j_1}(t)}{\lambda_{\mathrm{m},j_2}(t)}=\frac{\lambda_{0\mathrm{f}}(t)}{\lambda_{0\mathrm{m}}(t)}\exp(\beta_{\mathrm f} x_{j_1}-\beta_{\mathrm{m}}x_{j_2}),\quad t>0,
$$
which varies as $t$ varies unless the ratio of baseline hazards $\lambda_{0\mathrm{f}}(t)/ \lambda_{0\mathrm{m}}(t)$ is constant.
If you do not believe that baseline hazards of males and females are proportional, then you might want to specify a stratified model (now dropping the subscript $j$ for individuals)
$$
\lambda_{i}(t)=\lambda_{0i}(t)\exp(\beta_i x),\quad t>0,\quad i=\mathrm{f},\mathrm{m}.
$$
If you do believe that the baseline hazards of males and females are proportional, then you might want to specify the model (which is the model you have specified in R if I'm not mistaken)
$$
\lambda_{i}(t)=\lambda_{0}(t)\exp(\beta_i x),\quad t>0,\quad i=\mathrm{f},\mathrm{m}.
$$
