C2.3 Cox Proportional Hazards Regression

Unadjusted Cox regression

The ‘stcox’ command is used to run a Cox proportional hazards regression model. You must use the ‘stset’ command before using stcox, and due to this you do not need to specify an outcome variable or a time variable again. You only need to specify the explanatory variables that you are examining in your model.

stcox [varlist] [if] [in] [, options]

For example, for frailty, the command would be:

stcox i.frailty

The output looks like this:

Cox regression with Breslow method for ties

No. of subjects =       4,327                           Number of obs =  4,327
No. of failures =       1,526
Time at risk    = 29,476.9825

                                                                            LR chi2(4)    = 577.22
Log likelihood = -12128.151                                                 Prob > chi2   = 0.0000

 ——————————————————————————
          _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
————-+—————————————————————-
     frailty |
          2  |    1.45484   .1522779     3.58   0.000     1.185005    1.786119
          3  |   1.944776    .195232     6.63   0.000     1.597421    2.367664
          4  |   2.884121   .2643057    11.56   0.000     2.409949    3.451588
          5  |   5.740893   .5005441    20.04   0.000     4.839091    6.810752
——————————————————————————

Looking at the output, you can see that the hazard of dying is nearly 6 times greater for participants in frailty level 5 compared with frailty level 1 (HR 5.7, 95% CI:4.8-6.8). This is statistically significant.

• Question C2.3a.i: What is association of vitamin D (using the variable coded above and below the median) with the risk (hazard) of dying?

Multivariable Cox Regression

We can adjust for confounders in a Cox regression as we do in any other regression model. For example, if we want to examine the independent association of vitamin D on the hazard of mortality, adjusted for age group, we type:

stcox vitd_med i.age_grp

The output should look like this:

        Failure _d: death
  Analysis time _t: fu_years
       ID variable: whl1_id

Iteration 0:   log likelihood = -12416.762
Iteration 1:   log likelihood = -12186.532
Iteration 2:   log likelihood = -12137.041
Iteration 3:   log likelihood = -12136.832
Iteration 4:   log likelihood = -12136.832
Refining estimates:
Iteration 0:   log likelihood = -12136.832

Cox regression with Breslow method for ties

No. of subjects =       4,327                           Number of obs =  4,327
No. of failures =       1,526
Time at risk    = 29,476.9825
                                                        LR chi2(4)    = 559.86
Log likelihood = -12136.832                             Prob > chi2   = 0.0000

——————————————————————————
          _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
————-+—————————————————————-
    vitd_med |   .6728851    .035984    -7.41   0.000      .605928    .7472412
             |
     age_grp |
      71-75  |   1.860847   .1792291     6.45   0.000     1.540729    2.247475
      76-80  |   3.128741   .3046026    11.72   0.000     2.585233    3.786513
      81-95  |   5.800825   .5781615    17.64   0.000     4.771463    7.052256
——————————————————————————

We can see that once vitamin D has been adjusted for age group, the hazard ratio weakens from 0.57 (unadjusted HR, from section C2.3a.i.) to 0.67. There is still a statistically significant association (p<0.001), but it is weaker (i.e. closer to the null value of HR=1) since we controlled for the important confounder of age.

• Question C2.3a.ii: What is the impact of frailty on the hazard of dying, once adjusted for age_grp?

Test the proportional hazards assumption

There are many similarities between Cox regression and linear and logistic regression. We have seen that you can adjust for additional variables in the same way. You can also fit interactions using the same commands. However, the one key difference between a Cox regression and other regression models is the proportional hazards assumption.

The proportional hazards assumption states that hazard ratio (comparing levels of an exposure variable) is constant over time. This assumption needs to be tested. The lecture explained 3 ways to do this: fitting an interaction with time, examining the residuals with a log-log plot and testing the residuals for a correlation with time. In practice, generally just choosing one method to use is sufficient.

Fitting an interaction with time

If the PH assumption is violated, it means that the hazard ratio is not constant over time, i.e. there is an interaction between a covariate and time. Thus, we can fit a Cox model with a time interaction in it and see if this is significant. We can also use a likelihood ratio (LR) test to see if a model with a time-varying covariate fits better than one without. The ‘tvc’ and ‘texp’ options of the stcox are used for this.

For example, to test age the command would be:

stcox i.age_grp frailty, tvc(age_grp)

——————————————————————————
          _t | Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]
————-+—————————————————————-
main         |
     age_grp |
      71-75  |   1.534684   .1700356     3.87   0.000     1.235122      1.9069
      76-80  |   2.028394   .2977641     4.82   0.000      1.52124    2.704624
      81-95  |   3.009572   .5683724     5.83   0.000     2.078505     4.35771
             |
     frailty |   1.434477   .0290646    17.81   0.000     1.378627    1.492588
————-+—————————————————————-
tvc          |
     age_grp |   1.026909    .012475     2.19   0.029     1.002747    1.051653
——————————————————————————
Note: Variables in tvc equation interacted with _t.

When we look at the p-value for the coefficient in the ‘tvc’ row, we see that it is significant (p=0.03). This suggests there is a significant interaction with time, which indicates the HR for age group is not constant over time.

Inspecting log-log plots

To compute a log-log plot, the command is ‘stphplot’, which can be adjusted for confounders.

A log-log plot for age group is computed using the command below:

stphplot, by(age_grp) adjust(frailty)

When we inspect a log-log plot, we want to see that the curves in the main portion of the plot (i.e. not the tails of the plot, where there are few events) are generally parallel over time. If the curves touch or cross, this indicates there is a violation of the PH assumption.

Looking at the curves for each level of age group above, they appear to be generally parallel with similar slopes over time. This looks acceptable and does not indicate a meaningful violation of the PH assumption. This may contradict the results of the ‘tvc’ command above as statistical tests will often be significant for small deviations if you have a lot of power (e.g. if you have a large dataset). For this reason, we often prefer to visually inspect the data for obvious violations of the PH assumption.

• Question C2.3b:  Interpret a log-log plot for frailty (adjusted for age group) and assess if it violates the PH assumption. 

Fixing violations in the PH assumption

If the PH assumption has been violated for a covariate, the lecture mentioned 2 ways to deal with the situation. One way is to specify the model using the ‘tvc’ option for the offending variable, which models an interaction with time. The other way, which we are doing here, is to fit a stratified Cox model where we stratify by the time-varying covariate.

We do this by using the strata option of the stcox command.

stcox i.age_grp, strata(frailty)

You do not get HRs or p-values for each stratum. Instead, Stata assumes there is a separate baseline hazard for each stratum in order to estimate the proportional effects of the other variables across strata (which provides only one coefficient since the proportional effects are assumed to be constant, which is the no-interaction assumption [ie no interaction between other explanatory variables and the time-varying coefficient]). Since you do not get HR for levels of a variable you stratify on, do not use this method if you have non-proportionality in your main exposure of interest.

If you run the command you above, you will get this output:

Stratified Cox regression with Breslow    method    for ties
Strata variable: frailty

No. of subjects =       4,327            Number of obs =  4,327
No. of failures =       1,526
Time at risk    = 29,476.9825
            LR chi2(3)    = 280.76
Log likelihood = -9708.127            Prob > chi2   = 0.0000

           
_t  Haz. ratio   Std. err.    z    P>z    [95% conf. interval]
            
age_grp 
71-75     1.756097   .1693938    5.84    0.000    1.453588    2.121562
76-80     2.642535   .2592631    9.90    0.000    2.180259    3.202827
81-95      4.28578   .4339019    14.37    0.000    3.514409    5.226457

Notice that you do not get any output regarding the variable ‘frailty’ that you stratified on, although Stata indicates at the top of the output that you have used this variable in the strata. You would still interpret the HR of age_grp as having been adjusted for frailty, though. So you would say that the hazard of dying was 4.3 times greater in the 81-95 age group compared to the 60-70 age group, once it was adjusted for frailty (as stratifying on a variable is another way to adjust for its association).