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

command 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 Mantel-Haenzsel, 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.