Home | GLMs | Multilevel | Survival | Demography | Stata | R

Cox's Proportional Hazards Model

We continue our analysis of the cancer relapse data used for Kaplan-Meier. This is the dataset used as an example in Cox's original paper: Cox, D.R. (1972) Regression Models and Life tables, (with discussion) Journal of the Royal Statistical Society, 34: 187--220. (This is a slightly revised and shorter version of a handout used in my survival analysis short course.)

The first task is to read and stset the data:

. infile group weeks relapse using ///
>         http://data.princeton.edu/eco572/datasets/gehan.dat
(42 observations read)

. stset weeks, failure(relapse)

     failure event:  relapse != 0 & relapse < .
obs. time interval:  (0, weeks]
 exit on or before:  failure

------------------------------------------------------------------------------
       42  total obs.
        0  exclusions
------------------------------------------------------------------------------
       42  obs. remaining, representing
       30  failures in single record/single failure data
      541  total analysis time at risk, at risk from t =         0
                             earliest observed entry t =         0
                                  last observed exit t =        35

Relative Risk

To fit a Cox model we use the stcox with all the defaults:

. gen treated = group == 2

. stcox treated

         failure _d:  relapse
   analysis time _t:  weeks

Iteration 0:   log likelihood =  -93.98505
Iteration 1:   log likelihood = -86.385606
Iteration 2:   log likelihood = -86.379623
Iteration 3:   log likelihood = -86.379622
Refining estimates:
Iteration 0:   log likelihood = -86.379622

Cox regression -- Breslow method for ties

No. of subjects =           42                     Number of obs   =        42
No. of failures =           30
Time at risk    =          541
                                                   LR chi2(1)      =     15.21
Log likelihood  =   -86.379622                     Prob > chi2     =    0.0001

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     treated |   .2210887   .0905501    -3.68   0.000     .0990706    .4933877
------------------------------------------------------------------------------

Stata reports hazard ratios (exponentiated coefficients) by default. Here we see that treatment reduces the risk of relapse by 78% at any duration. To obtain the coefficients you can use the nohr option, or with just one type

. di _b[treated]
-1.5091914

Treatment of Ties

Stata has several ways of handling ties. Cox's original proposal is called exactp for exact partial likelihood. An alternative is exactm for the exact marginal likelihood. Both are computationally intensive. A good approximation is efron, due to Efron. The default is breslow, due to Breslow and Peto. Let us compare all four

. estimates store breslow 

. quietly stcox treated, efron

. estimates store efron

. quietly stcox treated, exactm

. estimates store exactm

. quietly stcox treated, exactp

. estimates store exactp

. estimates table breslow efron exactm exactp

------------------------------------------------------------------
    Variable |  breslow       efron        exactm       exactp    
-------------+----------------------------------------------------
     treated | -1.5091914   -1.5721251   -1.5981915    -1.628244  
------------------------------------------------------------------

As you can see, Efron's approximation is pretty good. Cox reports -1.65 in his paper, a value he obtained by evaluating the partial likelihood over a grid of points. Stata's more exact calculation yields -16.3, so Cox did pretty well with a grid.

Proportionality of Hazards

There are many ways of testing proportionality of hazard. Here we will let the trestment indicator interact with time linearly. We use tvc(treatment) to indicate that we want treatment as a time-varying covariate, and texp(_t) to specify a linear interaction with time

. stcox treated, tvc(treated) texp(_t)

         failure _d:  relapse
   analysis time _t:  weeks

Iteration 0:   log likelihood =  -93.98505
Iteration 1:   log likelihood = -86.600602
Iteration 2:   log likelihood = -86.370984
Iteration 3:   log likelihood = -86.370763
Refining estimates:
Iteration 0:   log likelihood = -86.370763

Cox regression -- Breslow method for ties

No. of subjects =           42                     Number of obs   =        42
No. of failures =           30
Time at risk    =          541
                                                   LR chi2(2)      =     15.23
Log likelihood  =   -86.370763                     Prob > chi2     =    0.0005

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
rh           |
     treated |   .2026644   .1573446    -2.06   0.040     .0442509    .9281822
-------------+----------------------------------------------------------------
t            |
     treated |   1.008168   .0617825     0.13   0.894     .8940658    1.136832
------------------------------------------------------------------------------

Note: Second equation contains variables that continuously vary with respect to time; variables are interacted with
      current values of _t.

We see no evidence that the treatment effect varies linearly with time. Cox used this same test in his paper, except that he worked with t-10 to achieve approximate orthogonality.

An alternative test is to allow a different treatment effect after 10 weeks:

. stcox treated, tvc(treated) texp(_t>10)

         failure _d:  relapse
   analysis time _t:  weeks

Iteration 0:   log likelihood =  -93.98505
Iteration 1:   log likelihood = -86.336856
Iteration 2:   log likelihood = -86.164454
Iteration 3:   log likelihood = -86.163329
Iteration 4:   log likelihood = -86.163329
Refining estimates:
Iteration 0:   log likelihood = -86.163329

Cox regression -- Breslow method for ties

No. of subjects =           42                     Number of obs   =        42
No. of failures =           30
Time at risk    =          541
                                                   LR chi2(2)      =     15.64
Log likelihood  =   -86.163329                     Prob > chi2     =    0.0004

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
rh           |
     treated |   .2801237   .1479049    -2.41   0.016     .0995215     .788466
-------------+----------------------------------------------------------------
t            |
     treated |   .5841662   .4777255    -0.66   0.511     .1176067    2.901622
------------------------------------------------------------------------------

Note: Second equation contains variables that continuously vary with respect to time; variables are interacted with
      current values of _t>10.

We estimate a 72% reduction in risk in the first ten weeks and an additional 42% reduction later (for a total of 84%). However, the difference is not significant; we have no evidence that the treatment effect is not proportional.

Baseline Survival

Stata can compute the baseline hazard using a Kaplan-Meier-like estimate (using the estimated relative risks as weights) using the basesurv option. (It can also compute an estimate of the baseline hazard due to Nelson and Aalen. This is avaiable via the basech option. This is very similar to negative log of the Kaplan-Meier estimate.)

We use the option to compute the baseline survival, save it in a variable called S0, and then raise it to the relative risk to get survival in the treated group, which we call S1

. stcox treated, basesurv(S0)

         failure _d:  relapse
   analysis time _t:  weeks

Iteration 0:   log likelihood =  -93.98505
Iteration 1:   log likelihood = -86.385606
Iteration 2:   log likelihood = -86.379623
Iteration 3:   log likelihood = -86.379622
Refining estimates:
Iteration 0:   log likelihood = -86.379622

Cox regression -- Breslow method for ties

No. of subjects =           42                     Number of obs   =        42
No. of failures =           30
Time at risk    =          541
                                                   LR chi2(1)      =     15.21
Log likelihood  =   -86.379622                     Prob > chi2     =    0.0001

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     treated |   .2210887   .0905501    -3.68   0.000     .0990706    .4933877
------------------------------------------------------------------------------

. gen S1 = S0^exp(_b[treated])

Before we plot these we will generate two separate Kaplan-Meier estimates

. sts gen KM = s, by (treated)

Now we plot them

. twoway (scatter S0  _t, c(J) sort)      /// baseline
>          (scatter S1 _t , c(J) sort)     /// treated
>          (scatter KM _t if treated, msymbol(circle_hollow)) /// KM treated
>          (scatter KM _t if !treated, msymbol(X)) /// KM base
>          , legend(off) ytitle(Survival) xtitle(Weeks) ///
>          title(Relapse of Leukemia Patients in Treated and Control Groups) ///
>          subtitle(Kaplan-Meir and Proportional Hazard Estimates) 

. graph export coxsurvival.png, replace
(file coxsurvival.png written in PNG format)

This reproduces Figure 1 in Cox's original paper.

Hazard Plot

One can also use graphs to test proportionality of hazards. This is best done by plotting the estimated cumiulative hazards, which should be parallel. (This is because the eye is much better at judging whether lines are parallel than at judging whether they are proportional.) Stata can do this automatically via the stphplot command

. stphplot, by(treated) legend(off) title(Plot of -log(-log(S(t))))

         failure _d:  relapse
   analysis time _t:  weeks

. graph export coxhazard.png, replace
(file coxhazard.png written in PNG format)

The two lines look quite parallel indeed.

Notes

There are several tests for comparing survival in different groups, the best known being Hantel-Maenzsel, which is closely related to Cox regression with a dummy variable.

The partial likelihood estimator illustrated here is appropiate for continuous data, with few if any ties. Many demographic applications use survival models appropriate for discrete data (using a logit model) or for grouped continuous data (using a Poisson model). These approaches are easy to implement except for computing events and exposure, and we will have opportunity to illustrate those calculations using demographic survey data.