# 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?

• The explanation is a general one: one factor is a special case of the interaction between factors. Commented Jan 14, 2017 at 11:12

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.

• In practice, how do we decide between these two approaches? stratify or adjust.
– hehe
Commented Sep 29, 2022 at 15:02
• A reason to stratify is to (hopefully) correct violations of the PH assumption. In my practice, I also implement stratification when I don’t need to explicitly measure the effect of a variable but I know its values have vastly different underlying risks (eg, magnitude or “shape”; another form of PH violation). I could not find a related post with more detailed information. If you need a detailed answer, you may want to post a new question. Commented Sep 30, 2022 at 21:49
• @ToddD, in the scenario where a variable isn't of interest as a covariate but might work as a stratifying variable, what criteria do you use to compare (a) coxph(Surv(futime, fustat) ~ age) with (b) coxph(Surv(futime, fustat) ~ age + strata(rx))? Do you fit both models and compare their AIC? Or does it make more sense to make that call based on context (e.g., use rx as a stratifying variable when we have biomedical reasons to expect that baseline hazard differs by rx)? Commented Dec 12, 2022 at 14:48

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

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), data=ovarian)


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