Cox Regression
We continue our analysis of the Gehan data fitting a proportional hazards model. This is, of course, 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.
The first task is to read and stset the data:
. infile group weeks relapse using ///
> http://data.princeton.edu/wws509/datasets/gehan.raw, clear
(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
. gen treated = group==2
Here's a run with all the defaults:
. 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, unless you specify the
option nohr.
We see that the treatment reduces the risk of relapse by 78%
at any duration.
To see the unexponentiated coefficient you can type
stcox, nohr; but with only one predictor it is
just as easy to display the coefficient:
. di _b[treated] -1.5091914
Handling Ties
Stata has several options for handling ties, Cox's original
proposal is called exactp for "exact partial",
and relies on the discrete partial likelihood. An alternative
is exactm for exact marginal likelihood.
Both are computationally intensive.
A good approximation is efron, due to Efron.
The default is breslow, which is the Breslow-Peto
approximation. Here are results using 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, t
------------------------------------------------------------------
Variable | breslow efron exactm exactp
-------------+----------------------------------------------------
treated | -1.5091914 -1.5721251 -1.5981915 -1.628244
| -3.68 -3.81 -3.79 -3.76
------------------------------------------------------------------
legend: b/t
As you can see, Efron's approximation to the exact marginal likelihood is pretty good. Cox reports -1.65 in his paper, which he obtained by evaluating the likelihood using a grid of points. A more exact calculation yields -1.63, so he did pretty well by evaluating the likelihood on a grid of points (see page 197).
Proportionality of Hazards
Stata provides several way of testing proportionality of hazards. One general approach is to introduce an interaction between a variable, in this case treatment, and time. This can be done using the options
tvc(varlist), to specify the variable(s) we want to interact with time, andtexp(expression)to specify a function of time_t, typically just time,texp(_t), or log-time,texp(log(_t).
Stata will then create a variable equal to the product of
the variable specified in tvc() by the time expression
specified in texp().
For example to let the treatment effect vary linearly with time we can specify
. 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]
-------------+----------------------------------------------------------------
main |
treated | .2026644 .1573446 -2.06 0.040 .0442509 .9281822
-------------+----------------------------------------------------------------
tvc |
treated | 1.008168 .0617825 0.13 0.894 .8940658 1.136832
------------------------------------------------------------------------------
Note: variables in tvc equation interacted with _t
We see that there's no evidence that the treatment effect
changes linearly with time. In this case we didn't really have to
specify texp(_t) because that's the default.
This is exactly what Cox did in his original paper, except that he worked with t-10 to achieve approximate orthogonality of estimation.
As an alternative we could allow different effects before and after 10 weeks, using an indicator variable in the time expression:
. 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]
-------------+----------------------------------------------------------------
main |
treated | .2801237 .1479049 -2.41 0.016 .0995215 .788466
-------------+----------------------------------------------------------------
tvc |
treated | .5841662 .4777255 -0.66 0.511 .1176067 2.901622
------------------------------------------------------------------------------
Note: variables in tvc equation interacted with _t>10
The point estimates indicate a 72% reduction in risk in the first ten weeks and an additional 42% reduction after ten weeks, for a total of 84% in the later period. However, the difference in treatment effects between the two periods is not significant.
Another way to obtain this result is to stsplit
the dataset. This has the advantage that one could use more than two
categories for time. We illustrate with two, however, to show that
we get the same results. (This method has more general applicability
to handle time-varying covariates.)
Before we split we need an id variable:
. gen id = _n
. streset, id(id)
-> stset weeks, id(id) failure(relapse)
id: id
failure event: relapse != 0 & relapse < .
obs. time interval: (weeks[_n-1], weeks]
exit on or before: failure
------------------------------------------------------------------------------
42 total obs.
0 exclusions
------------------------------------------------------------------------------
42 obs. remaining, representing
42 subjects
30 failures in single failure-per-subject data
541 total analysis time at risk, at risk from t = 0
earliest observed entry t = 0
last observed exit t = 35
. stsplit dur, at(10)
(21 observations (episodes) created)
. gen later = treated * (dur == 10)
. stcox treated later
failure _d: relapse
analysis time _t: weeks
id: id
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 = 63
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]
-------------+----------------------------------------------------------------
treated | .2801237 .1479049 -2.41 0.016 .0995215 .788466
later | .5841662 .4777255 -0.66 0.511 .1176067 2.901622
------------------------------------------------------------------------------
As you can see, we get exactly the same results. Although we could
continue working with the split datatset we will "rejoin" it,
illustrating Stata's stjoin command, which requires
first dropping the variables that vary by episode
. drop dur later . stjoin (option censored(0) assumed) (21 obs. eliminated)
Another way to test for proportionality of hazards is to use Schoenfeld residuals (and their scaled counterparts). You can obtain an overall test using the Schoenfeld residual, or a variable-by-variable test based on the scaled variant. In this case with just one predictor there is only one test.
Note:
in older versions of Stata you had to specify the residuals and
other quantities to be used in post-estimation when the model was fit.
Starting with Stata 11 the commands stcurve and
estat phtest automatically generate the statistics they require.
This log has been modified to use the new commands, although the older
options still work.
. estat phtest
Test of proportional-hazards assumption
Time: Time
----------------------------------------------------------------
| chi2 df Prob>chi2
------------+---------------------------------------------------
global test | 0.00 1 0.9886
----------------------------------------------------------------
The test shows that there is no evidence against the proportional
hazards assumption. If there had been, we could get a hint of the
nature of the time dependence by plotting the (scaled) residuals
against time and using a smoother to glean the trend, if any,
using the plot() option:
. estat phtest, plot(treated) . graph export phplot.png, width(500) replace (file phplot.png written in PNG format)
The residuals show no time trend at all, showing that the treatment hazard ratio is fairly constant over time. (We will confirm this result below with a plot of cumulative hazards that provides more direct evidence.)
I'll skip the Cox-Snell residuals because I am persuaded by Therneau and Grambsch's argument that they are not very informative.
Baseline Survival and Hazard
Stata can compute the baseline survival using a Kaplan-Meier type
estimate, or the baseline cumulative hazard using a Nelson-Aalen
like estimate. You can generate these functions setting the
covariates to zero using sts gen, or plot them
at selected values of the covariates using stcurve
with the at() to set some covariates at selected
values and the rest at their means.
Here are the estimated survival functions for the control and
treated groups done the easy way:
. stcurve , surv at1(treated=0) at2(treated=1) . graph export coxstcurve.png, width 400, replace
c(J) option, which is short for
"stairstep", to plot a step function. In this case the
baseline survival is the survival functin for the control group.
. sts gen S0 = basesurv // control . gen S1 = S0^exp(_b[treated]) // treated . sts gen KM = s, by(treated) // two Kaplan-Meiers . twoway (scatter S0 _t, c(J) ms(none) sort) /// baseline > (scatter S1 _t , c(J) ms(none) sort) /// treated > (scatter KM _t if treated, msymbol(circle_hollow)) /// KM treated > (scatter KM _t if !treated, msymbol(X)) /// KM base > , legend(off) /// > title(Relapse of Leukemia Patients in Treated and Control Groups) /// > subtitle(Kaplan-Meier and Proportional Hazards Estimates) . graph export CoxKM.png, width(500) replace (file CoxKM.png written in PNG format)
The figure looks just like Figure 1 in Cox's paper.
You can get pretty much the same thing using Stata's
stcoxkm, by(treated), although you will have to
tweak the default options.
If the purpose of the graph is to test the proportional hazards assumption
you should try stphplot, by(treated),
which plots -log(-log(S(t)) for each group.
Under the proportional hazards assumption the resulting curves should be
parallel. The eye is much better at judging whether curves are parallel
than whether they are proportional. For this dataset:
. stphplot, by(treated) legend(off) title(Plot of -log(-log(S(t))))
failure _d: relapse
analysis time _t: weeks
id: id
. graph export stphplot.png, width(500) replace
(file stphplot.png written in PNG format)
The two lines look quite parallel indeed.