The Gehan Survival Data
The WWS 509 datasets page has a copy of the Gehan survival data which were analyzed by Cox in his original proportional hazards paper, see http://data.princeton.edu/wws509/datasets/#gehan.
The data shows the length of remission in weeks for two groups of leukemia patients, treated and control. You can read the data directly from the WWS509 website using the command
. infile group weeks relapse using /// > http://data.princeton.edu/wws509/datasets/gehan.raw (42 observations read)
As usual the first task is to stset the data:
. 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
As you can see, 30 of the 42 patients had a relapse, with a total observation time of 541 weeks. The weekly relapse rate is 7.8%. We will also label the groups:
. label define group 1 "Control" 2 "Treated" . label values group group
Kaplan-Meier
The Kaplan-Meier estimators for the two groups are easily plotted
using the command sts graph with the by
option to indicate the groups:
. sts graph, by(group)
failure _d: relapse
analysis time _t: weeks
This shows that after 23 weeks all patients in the control group had relapsed but about half those in the treated group remained in remission. This reproduces Figure 1 in the notes.
You can list the estimates using sts list:
. sts list, by(group)
failure _d: relapse
analysis time _t: weeks
Beg. Net Survivor Std.
Time Total Fail Lost Function Error [95% Conf. Int.]
-------------------------------------------------------------------------------
Control
1 21 2 0 0.9048 0.0641 0.6700 0.9753
2 19 2 0 0.8095 0.0857 0.5689 0.9239
3 17 1 0 0.7619 0.0929 0.5194 0.8933
4 16 2 0 0.6667 0.1029 0.4254 0.8250
5 14 2 0 0.5714 0.1080 0.3380 0.7492
8 12 4 0 0.3810 0.1060 0.1831 0.5778
11 8 2 0 0.2857 0.0986 0.1166 0.4818
12 6 2 0 0.1905 0.0857 0.0595 0.3774
15 4 1 0 0.1429 0.0764 0.0357 0.3212
17 3 1 0 0.0952 0.0641 0.0163 0.2612
22 2 1 0 0.0476 0.0465 0.0033 0.1970
23 1 1 0 0.0000 . . .
Treated
6 21 3 1 0.8571 0.0764 0.6197 0.9516
7 17 1 0 0.8067 0.0869 0.5631 0.9228
9 16 0 1 0.8067 0.0869 0.5631 0.9228
10 15 1 1 0.7529 0.0963 0.5032 0.8894
11 13 0 1 0.7529 0.0963 0.5032 0.8894
13 12 1 0 0.6902 0.1068 0.4316 0.8491
16 11 1 0 0.6275 0.1141 0.3675 0.8049
17 10 0 1 0.6275 0.1141 0.3675 0.8049
19 9 0 1 0.6275 0.1141 0.3675 0.8049
20 8 0 1 0.6275 0.1141 0.3675 0.8049
22 7 1 0 0.5378 0.1282 0.2678 0.7468
23 6 1 0 0.4482 0.1346 0.1881 0.6801
25 5 0 1 0.4482 0.1346 0.1881 0.6801
32 4 0 2 0.4482 0.1346 0.1881 0.6801
34 2 0 1 0.4482 0.1346 0.1881 0.6801
35 1 0 1 0.4482 0.1346 0.1881 0.6801
-------------------------------------------------------------------------------
Note that Stata reports the K-M survival function just after each distinct failure time. In the control group there are no censored observations, and the Kaplan-Meier estimate is simply the proportion alive after each distinct failure time. You can also check that the standard error is the usual binomial estimate. For example just after 8 weeks there are 8 out of 21 alive, the proportion is 8/21=0.381 and the standard error is sqrt(0.381(1-0.381)/21) = 0.106.
In the treated group 12 cases are censored and 9 relapse.
Can you compute the estimate by hand?
Try tab weeks relapse if group==2 to get the
basic counts. Only the distinct times of death
--6, 7, 10, 13, 16, 22 and 23-- are relevant.
The counts of deaths are 3,1,1,1,1,1,1.
When there are ties between event and censoring time it is
customary to assume that the event occurred first; that is,
observations censored at t are assumed to be exposed
at that time, effectively censored just after t.
The counts of censored observations after each death time
(but before the next) are 1,1,2,0,3,0 and 5.
Note that when the last observation
is censored the K-M estimate is greater than zero,
and is usually considered undefined from that point on.
Stata can compute the Nelson-Aalen estimate of the cumulative hazard
if you use the na option of sts list.
The command sts generate can create Stata
variables containing the Kaplan-Meier or Nelson-Aalen estimates.
Here's a comparison of the two methods of estimating a survival function for the overall sample
. sts generate kmS = s . sts generate naH = na . gen naS = exp(-naH) . twoway (scatter kmS _t, c(J) sort) (scatter naS _t, c(J) sort)
Mantel-Haenszel
Stata can compute the Mantel-Haenszel test using the
sts test command with the logrank option:
. sts test group
failure _d: relapse
analysis time _t: weeks
Log-rank test for equality of survivor functions
| Events Events
group | observed expected
--------+-------------------------
Control | 21 10.75
Treated | 9 19.25
--------+-------------------------
Total | 30 30.00
chi2(1) = 16.79
Pr>chi2 = 0.0000
We see that the treated have fewer deaths than would be expected (and thus the controls have more) if the two groups had the same survival distribution: 9 instead of 19.25. The Mantel-Haenszel chi-squared satistics of 16.79 on one d.f. is highly significant.
It is instructive to compute this statistic "by hand".
The Mantel-Haenszel test is just one of several statistics that Stata can compute. These are all based on a comparison of observed and expected counts of deaths at each distinct time of death but differ on the weights they assign to each time as shown below:
| Name | Weight | Option |
|---|---|---|
| Wilcoxon | ni | wilcoxon |
| Tarone-Ware | sqrt(ni) | tware |
| Peto-Peto-Prentice | S(ti) | peto |
| Fleming-Harrington | S(ti)p (1-S(ti))q | fh(p,q) |
The most popular of these is the Wilcoxon test (actually an extension of Wilcoxon's well-known non-parametric test due to Gehan and Breslow), which gives more weight to early failures.