Germán Rodríguez
Demographic Methods Princeton University

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

```. infile group weeks relapse using ///
>         http://data.princeton.edu/eco572/datasets/gehan.dat

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