stratification in cox model I am trying to understand difference in output between the following lines of code:
library(survival)

summary(coxph(Surv(futime, fustat) ~ age + strata(rx), data=ovarian))
summary(coxph(Surv(futime, fustat) ~ age, data=ovarian))

What does adding strata(rx) do? 
Explanation ?strata says: 
**Value**
a new factor, whose levels are all possible combinations
of the factors supplied as arguments.

What is the point to take all possible combinations of factor values, when they are mutually exclusive. What am I missing? Why example contains strata(rx) i.e. only one factor, when explanation is talking about interaction of factors?
 A: In a Cox model, stratification allows for as many different hazard functions as there are strata. Beta coefficients (hazard ratios) optimized for all strata are then fitted.
In your example, the model coxph(Surv(futime, fustat) ~ age + strata(rx) will output a hazard ratio for age in the presence of two (or more) hazards intrinsic to the levels of rx. If rx violated the proportional hazards assumption, for example, stratifying may help meet the PH assumption and provide more valid estimates for age. The effect of rx is not explicitly provided as a hazard ratio. Likelihood estimates for the model can be used to assess whether stratification by rx improved the model fit.
coxph(Surv(futime, fustat) ~ age will output a hazard ratio for age only, assuming that the hazard for different levels of rx are the same. In this model, the effect of rx is not explicitly modeled. 
The model coxph(Surv(futime, fustat) ~ age + rx may be useful to consider. This model would provide estimates of the HR for age and rx with both present in the model ("adjusted for one another"). This would be different from coxph(Surv(futime, fustat) ~ age + str(rx) in that the unstratified model provides estimation of effect for both age and rx using a single underlying hazard. 
A: 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.
