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

1. Fit a Cox model (fit.unstrat) with all of the covariates and no stratifying variables.

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

3. 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 test:

library(survival)
ovarian.surv <- Surv(ovarian$futime, ovarian$fustat)
fit.strat <- coxph(ovarian.surv ~ age + strata(rx), data=ovarian)
fit.grp1 <- coxph(ovarian.surv ~ age, subset=(rx==1), data=ovarian)
fit.grp2 <- coxph(ovarian.surv ~ 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( cox_fit.stratified$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)

• 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:

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

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

These equations imply:

1. The two rx groups have different baseline hazards / baseline survivals.
2. 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))
fit.grp2 <- coxph(Surv(futime, fustat) ~ age, subset=(rx==2))

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.