In addition to the good answer from Todd B, here's a process you can follow, and some more of the background math.
A reasonable process to follow
Fit a Cox model (fit.unstrat) with all of the covariates and no stratifying variables.
Use cox.zph(fit.unstrat) to check for violations of the proportional hazards assumption. If you see e.g. rx has p<0.05, then rx violates the proportional hazards assumption and should not be included as a covariate in the Cox model. A reasonable follow-up is to try changing it to a stratifying variable, e.g., fit.strat <- update(fit.unstrat, . ~ . - rx + strata(rx)).
Next, it's appropriate to evaluate whether multiple models improve over the stratified model. The stratified Cox model allows the two rx groups to have different baseline hazards / baseline survivals, but it still only calculates a single effect $\beta_\text{age}$ of age on survival time. You can think of this as being the mean effect of age averaged across both rx groups. This raises a question: do you get a statistically significantly better fit when you allow the two rx groups to have different covariate effects ($\beta_\text{age, rx=2}$ and $\beta_\text{age, rx=1}$)? To answer this question, you:
- fit the stratified model,
- fit two separate Cox models for datasets subsetted by
rx, and
- compare them in a likelihood ratio test, where the p-value is $P\left(\chi^2_{(k-1)(q)}\right)>\ell_\text{stratified model} - \sum \ell_\text{unstratified models}$, with
$k$ being the number of stratifying groups (just two here--
rx=1 and rx=2), $q$ being the number of covariates (just 1 here--$z_\text{age}$), and $\ell$ is a model's log-likelihood.
Here's some sample code for the likelihood ratio tests. In this example (pretending rx violates the proportional hazards assumption), we wouldn't choose the multiple models, because while they improve the log likelihood, it's not a statistically significant improvement in the fit (p=0.1575):
library(survival)
fit.strat <- coxph(Surv(futime, fustat) ~ age + strata(rx), data=ovarian)
fit.grp1 <- coxph(Surv(futime, fustat) ~ age, subset=(rx==1), data=ovarian)
fit.grp2 <- coxph(Surv(futime, fustat) ~ age, subset=(rx==2), data=ovarian)
LL.strat <- fit.strat$loglik[2]
LL.unstrat <- fit.grp1$loglik[2] + fit.grp2$loglik[2]
X2 <- -2*(LL.strat - LL.unstrat)
n_groups <- length( unique(ovarian$rx) )
n_params <- length( fit.strat$coef )
p_val <- 1 - pchisq(X2, df=(n_groups-1)*n_params)
Background math
Standard Cox proportional hazards model:
$$h(t|Z) = h_0 \exp(\beta_\text{age}Z_\text{age} + \beta_\text{rx}Z_\text{rx})$$
fit.unstrat <- coxph(Surv(futime, fustat) ~ age + rx, data=ovarian)
The quantity $h_0$ is called the baseline hazard. It corresponds with the baseline survival, which is the survival of the reference group whose covariates are all equal to 0.
For a Cox model, the baseline hazard is estimated from the data non-parametrically (without assuming any distribution).
The equation above implies:
(hazard for rx==2 group)/(hazard for rx==1 group) is a constant. In other words, group 2's hazard is always higher/lower than group 1's hazard by a constant factor.
The two rx groups have the same baseline hazard and baseline survival.
Stratified Cox proportional hazards model (you're actually fitting two equations that share a single $\beta_\text{age}$):
$$h_\text{rx=2}(t|Z) = h_{0,\text{rx=2}} \exp(\beta_\text{age}Z_\text{age}) \\
h_\text{rx=1}(t|Z) = h_{0,\text{rx=1}} \exp(\beta_\text{age}Z_\text{age})$$
fit.strat <- coxph(Surv(futime, fustat) ~ age + strata(rx), data=ovarian)
These equations imply:
- The two rx groups have different baseline hazards / baseline survivals.
- The two rx groups share a single factor, $\exp(\beta_\text{age})$, by which the hazard is always higher/lower for a 1-unit increase in age. In other words, we're specifying that both
rx groups have the same $\beta_\text{age}$.
Multiple models (this is just a term I'm making up):
$$h_\text{rx=2}(t|Z) = h_{0,\text{rx=2}} \exp(\beta_\text{age, rx=2}Z_\text{age}) \\
h_\text{rx=1}(t|Z) = h_{0,\text{rx=1}} \exp(\beta_\text{age, rx=1}Z_\text{age})$$
fit.grp1 <- coxph(Surv(futime, fustat) ~ age, subset=(rx==1), data=ovarian)
fit.grp2 <- coxph(Surv(futime, fustat) ~ age, subset=(rx==2), data=ovarian)
Short note about time-dependent covariates
In addition to fitting a stratified model, another typical way to remediate a violation of the proportional hazards assumption is to (1) convert the violating variable into a time-dependent variable, and then (2) fit a new Cox PH model and check if the time-dependent covariate satisfies the PH assumption. Time-dependent covariates are an entirely different topic, but they can get invoked for the same reason as a stratified Cox model. Picking between these two often depends on context. If you think the two rx groups truly have different baseline survivals, then a stratified model makes sense. If there's a time component to rx (e.g., the data also contain a variable indicating how much time passed before they started treatment, or we think rx might have a cumulative effect such that the treatment only provides a benefit after some amount of time has passed), then a time-dependent variable might make more sense.
