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.