### SAS Textbook Examples Applied Survival Analysis by Hosmer and Lemeshow Chapter 2: Descriptive Methods for Survival Data

In this chapter we will be using the hmohiv and the uis data sets.

Fig. 2.1, p. 28.
The first five observations of the hmohiv data set.
data mini;
set hmohiv;
if id <= 5 ;
run;
proc print data=mini noobs;
var id time censor;
run;
ID    time    censor

1       5        1
2       6        0
3       8        1
4       3        1
5      22        1

Table 2.2 and fig. 2.1, p. 32.

The estimated survivorship function computed for the mini data set in table 2.1, p. 28.
Note: The graph was generated by using the options plot=(s) where s=survival curve.
proc lifetest data=mini plots=(s);
time time*censor(0);
run;

<output omitted>

Product-Limit Survival Estimates
Survival
Standard     Number      Number
time     Survival    Failure      Error       Failed       Left
0.0000       1.0000           0           0       0           5
3.0000       0.8000      0.2000      0.1789       1           4
5.0000       0.6000      0.4000      0.2191       2           3
6.0000*           .           .           .       2           2
8.0000       0.3000      0.7000      0.2387       3           1
22.0000            0      1.0000           0       4           0

NOTE: The marked survival times are censored observations.

Fig. 2.2, p. 34.
The figure shows the Kaplan-Meier estimate of the survivorship function.
Note: We use the ods output system to create a data set est containing the Kaplan-Meier estimates and the components used in calculating the estimates.
ods listing close;
ods output  ProductLimitEstimates=est;
proc lifetest data=hmohiv plots=(s);
time time*censor(0);
run;
ods listing;


Table 2.3, p. 35.
Note: The output from proc lifetest has one line of output for each subject. In other words there would be 15 lines of output before the line containing the first estimate. We eliminate this extra output by including a where statement in the proc print stipulating that we do not wish to see the output where there is no Kaplan-Meier estimate.
proc print data=est noobs;
where survival ne . ;
var time left failed survival;
run;

<output omitted>

time    Left    Failed    Survival
0.0000    100        0        1.0000
1.0000     85       15        0.8500
2.0000     78       20        0.7988
3.0000     63       30        0.6894
4.0000     57       34        0.6442
5.0000     49       41        0.5636

...

57.0000      3       79        0.0584
58.0000      2       80        0.0389

In order to get other variables in the table, we will have to use a couple of data steps.

 ods listing close;
ods output  ProductLimitEstimates=t2;
proc lifetest data=hmohiv outs=t  ;
time time*censor(0);
run;
ods listing;
data t2a;
set t2;
retain c 0;
retain s 1;
by time;
if first.time & _n_ >1 then c = censor;
else if _n_ > 1 then c = c + censor;
if survival ~=. then s = survival;
keep time failed left c s ;
if last.time;
run;
data t2b;
set t2a;
flag = lag(failed);
d = failed-flag;
llag = lag(left);
s1 = (llag -d)/llag;
drop flag llag;
run;
proc print data = t2b noobs;
var time left d c s1 s;
run;
    time    Left     d    c       s1         s

0.0000    100      .    0     .         1.00000
1.0000     83     15    2    0.85000    0.85000
2.0000     73      5    5    0.93976    0.79880
3.0000     61     10    2    0.86301    0.68937
4.0000     56      4    1    0.93443    0.64417
5.0000     49      7    0    0.87500    0.56365
6.0000     46      2    1    0.95918    0.54064
7.0000     39      6    1    0.86957    0.47012
8.0000     35      4    0    0.89744    0.42190
9.0000     32      3    0    0.91429    0.38574
10.0000     28      3    1    0.90625    0.34958
11.0000     25      3    0    0.89286    0.31212
12.0000     21      2    2    0.92000    0.28715
13.0000     20      1    0    0.95238    0.27348
14.0000     19      1    0    0.95000    0.25981
15.0000     17      2    0    0.89474    0.23246
19.0000     16      0    1    1.00000    0.23246
22.0000     15      1    0    0.93750    0.21793
24.0000     14      0    1    1.00000    0.21793
30.0000     13      1    0    0.92857    0.20236
31.0000     12      1    0    0.92308    0.18680
32.0000     11      1    0    0.91667    0.17123
34.0000     10      1    0    0.90909    0.15566
35.0000      9      1    0    0.90000    0.14010
36.0000      8      1    0    0.88889    0.12453
43.0000      7      1    0    0.87500    0.10896
53.0000      6      1    0    0.85714    0.09340
54.0000      5      1    0    0.83333    0.07783
56.0000      4      0    1    1.00000    0.07783
57.0000      3      1    0    0.75000    0.05837
58.0000      2      1    0    0.66667    0.03892
60.0000      0      0    2    1.00000    0.03892

Table 2.4, Figure 2.3 and Figure 2.4 on page 38 and page 39.
 ods output  LifetableEstimates=est2;
proc lifetest data=hmohiv  method=lt intervals=(0 6 to 72 by 6);
time time*censor(0);
run;
proc print data = est2 noobs;
var lowertime uppertime failed censored effsize survival;
run;
   Lower       Upper                            Eff
Time        Time    Failed    Censored     Size    Survival

0           6      41         10        95.0      1.0000
6          12      21          3        47.5      0.5684
12          18       6          2        24.0      0.3171
18          24       1          1        16.5      0.2378
24          30       0          1        14.5      0.2234
30          36       5          0        14.0      0.2234
36          42       1          0         9.0      0.1436
42          48       1          0         8.0      0.1277
48          54       1          0         7.0      0.1117
54          60       3          1         5.5      0.0958
60          66       0          2         1.0      0.0435
Fig. 2.3, p. 38.
 symbol i = steplj c=red;
axis1 label=( a= 90 "Proportion Surviving") minor=none order =(0 to 1 by .5);
proc gplot data = est2;
format survival 2.1;
plot survival*lowertime / vaxis=axis1;
run;
quit;
 symbol i = join c=red value=dot;
axis1 label=( a= 90 "Proportion Surviving") minor=none order =(0 to 1 by .5);
proc gplot data = est2;
format survival 2.1;
plot survival*lowertime / vaxis=axis1;
run;
quit;


Fig. 2.5 on page 46. The code below works in SAS 9.1.3.
proc lifetest data=hmohiv noprint;
time time *censor(0);
survival out=Out1 confband=hw;
run;

symbol i=join v=circle r=5 ;
axis1 order = (0 to 1 by .2) label=(a=90 "Survival Probability") minor=none;
axis2 order = (0 to 60 by 10) label = ("Survival Time (Months)") minor=none;
legend label=none value=(h=2 font=swiss 'Kaplan-Meyer' 'HW lower' 'HW upper'
'Pointwise lower' 'Pointwise upper')
position=(top right inside) mode=share cborder=black;

proc gplot data = out1;
plot (survival hw_lcl hw_ucl sdf_lcl sdf_ucl)*time
/legend=legend1 vaxis=axis1 haxis=axis2 overlay;
run;
quit;

Table 2.5, p. 50.

proc lifetest data=hmohiv;
time time*censor(0);
run;

Quartile Estimates

Point     95% Confidence Interval
Percent    Estimate      [Lower      Upper)
75     15.0000     10.0000     34.0000
50      7.0000      5.0000      9.0000
25      3.0000      2.0000      4.0000

Table 2.6, p. 52.
proc lifetest data=hmohiv;
time time*censor(0);
survival out=out2;
run;
proc print data = out2 noobs;
where 4<=time<=9 & conftype~="";
run;
time    _CENSOR_    SURVIVAL    CONFTYPE    SDF_LCL    SDF_UCL

4         0        0.64417     LOGLOG     0.53868    0.73150
5         0        0.56365     LOGLOG     0.45636    0.65769
6         0        0.54064     LOGLOG     0.43343    0.63609
7         0        0.47012     LOGLOG     0.36440    0.56876
8         0        0.42190     LOGLOG     0.31832    0.52174
9         0        0.38574     LOGLOG     0.28454    0.48579

Fig. 2.7, p. 58.
Estimated survivorship functions for each group of drug use.
goptions reset = all;
proc lifetest data=hmohiv plots=(s) noprint;
time time*censor(0);
strata drug;
test age;
run;

Table 2.8, p. 63.
Inputting the mini1 dataset.
data mini1;
input time drug event;
cards;
3 0 1
4 0 0
5 0 1
22 0 1
34 0 1
2 1 1
3 1 1
4 1 1
7 1 0
11 1 1
;
run;
proc print data=mini1;
run;

time    drug    event
3       0       1
4       0       0
5       0       1
22       0       1
34       0       1
2       1       1
3       1       1
4       1       1
7       1       0
11       1       1

Table 2.9 on page 64. We will skip the data step for creating the d's and n's. Instead, we will focus on the calculation of other quantities in the table.
data table2_9;
retain survival 1;
input time d1 n1 d n;
lagn = lag(n);
lagd = lag(d);
if _n_> 1 then
survival = survival*(lagn+1-lagd)/(lagn+1);
e = n1*d/n;
v = n1*(n-n1)*d*(n-d)/(n**2*(n-1));
tw = sqrt(n);
pp = survival*n/(n+1);
lr = 1;
wl = n;
datalines;
2 0 5 1 10
3 1 5 2 9
4 0 4 1 7
5 1 3 1 5
11 0 2 1 3
22 1 2 1 2
34 1 1 1 1
;
run;
proc print data = table2_9 noobs;
var time d1 n1 d n e v lr wl tw pp;
run;
time    d1    n1    d     n       e          v       lr    wl       tw         pp

2      0     5    1    10    0.50000    0.25000     1    10    3.16228    0.90909
3      1     5    2     9    1.11111    0.43210     1     9    3.00000    0.81818
4      0     4    1     7    0.57143    0.24490     1     7    2.64575    0.63636
5      1     3    1     5    0.60000    0.24000     1     5    2.23607    0.53030
11      0     2    1     3    0.66667    0.22222     1     3    1.73205    0.39773
22      1     2    1     2    1.00000    0.00000     1     2    1.41421    0.26515
34      1     1    1     1    1.00000     .          1     1    1.00000    0.13258

Table 2.10, p. 64.
Log-rank and Wilcoxon tests of equality of the survivorship functions across the two drug strata.

proc lifetest data=mini1;
time time*event(0);
strata drug /test=(logrank wilcoxon tarone peto modpeto);
run;
        Test of Equality over Strata

Pr >
Test         Chi-Square      DF    Chi-Square

Log-Rank         1.5118       1      0.2189
Wilcoxon         1.2500       1      0.2636
Tarone           1.3632       1      0.2430
Peto             1.4229       1      0.2329
Modified Peto    1.3695       1      0.2419

Table 2.11, p. 65.
Testing for equality of survivorship functions across drug strata for the hmohiv data.
proc lifetest data=hmohiv;
time time*censor(0);
strata drug /test=(logrank wilcoxon tarone peto);
run;
       Test of Equality over Strata

Pr >
Test      Chi-Square      DF    Chi-Square

Log-Rank     11.8556       1      0.0006
Wilcoxon     10.9104       1      0.0010
Tarone       12.3359       1      0.0004
Peto         11.4974       1      0.0007

Table 2.12 , Figure 2.8, Table 2.14 and Table 2.15 on page 65, page 69 and page 70.

Creating the categorical variable for age.

data agecat;
set hmohiv;
if age le 29 then agecat = 1;
else if age le 34 then agecat = 2;
else if age le 39 then agecat = 3;
else agecat = 4;
run;
goptions reset=all;

proc lifetest data=agecat plots=(s);
time time*censor(0);
strata agecat /trend test=(logrank wilcoxon tarone peto);
symbol1 c=blue l=12 h=.5;
symbol2 c=red l=1 h=.5;
symbol3 c=green l=8 h=.5;
symbol4 c=purple l=2 h=.5;
run;

Summary of the Number of Censored and Uncensored Values

Percent
Stratum          agecat       Total  Failed    Censored    Censored
1               1          12       8           4       33.33
2               2          34      29           5       14.71
3               3          25      20           5       20.00
4               4          29      23           6       20.69
-------------------------------------------------------------------
Total                         100      80          20       20.0

Quartile Estimates

Stratum 1: agecat = 1

Point     95% Confidence Interval
Percent    Estimate      [Lower      Upper)
75     58.0000     43.0000       .
50     43.0000      8.0000       .
25      8.0000      5.0000     54.0000

Stratum 2: agecat = 2

Point     95% Confidence Interval
Percent    Estimate      [Lower      Upper)
75     15.0000     10.0000     32.0000
50      9.0000      7.0000     12.0000
25      5.0000      2.0000      8.0000

Stratum 3: agecat = 3

Point     95% Confidence Interval
Percent    Estimate      [Lower      Upper)
75     11.0000      7.0000     34.0000
50      7.0000      3.0000      9.0000
25      3.0000      2.0000      5.0000

Stratum 4: agecat = 4

Point     95% Confidence Interval
Percent    Estimate      [Lower      Upper)
75      5.0000      4.0000     13.0000
50      4.0000      3.0000      5.0000
25      2.0000      1.0000      3.0000
       Test of Equality over Strata

Pr >
Test      Chi-Square      DF    Chi-Square

Log-Rank     19.9058       3      0.0002
Wilcoxon     14.1433       3      0.0027
Tarone       16.9556       3      0.0007
Peto         15.6650       3      0.0013
                        Trend Tests

Test      Standard
Test       Statistic         Error       z-Score    Pr > |z|

Log-Rank     34.3443        8.0563        4.2630      <.0001
Wilcoxon   1919.0000      517.5932        3.7075      0.0002
Tarone      243.3701       60.4948        4.0230      <.0001
Peto         18.8199        4.8370        3.8908      <.0001

Fig. 2.8, p. 69.

Table 2.17, p. 76.
proc phreg data = hmohiv noprint;
model time*censor(0)=;
baseline  out=na cumhaz = h survival=s;
run;
data na;
set na;
shat = exp(-h);
run;
proc print data = na noobs;
run;
time       s          h         shat

0     1.00000    0.00000    1.00000
1     0.85000    0.15000    0.86071
2     0.79880    0.21024    0.81039
3     0.68937    0.34723    0.70664
4     0.64417    0.41280    0.66179
5     0.56365    0.53780    0.58403
6     0.54064    0.57862    0.56067
7     0.47012    0.70905    0.49211
8     0.42190    0.81162    0.44414
9     0.38574    0.89733    0.40766
10     0.34958    0.99108    0.37118
11     0.31212    1.09822    0.33346
12     0.28715    1.17822    0.30783
13     0.27348    1.22584    0.29351
14     0.25981    1.27584    0.27920
15     0.23246    1.38111    0.25130
22     0.21793    1.44361    0.23608
30     0.20236    1.51503    0.21980
31     0.18680    1.59196    0.20353
32     0.17123    1.67529    0.18725
34     0.15566    1.76620    0.17098
35     0.14010    1.86620    0.15471
36     0.12453    1.97731    0.13844
43     0.10896    2.10231    0.12217
53     0.09340    2.24517    0.10591
54     0.07783    2.41183    0.08965
57     0.05837    2.66183    0.06982
58     0.03892    2.99517    0.05003

Fig. 2.10, p. 77.
Graphs of the Nelson-Aalen and Kaplan-Meier estimators of the survivorship function from the hmohiv data.

symbol1 value = dot h=.8  i = stepjl color = black;
symbol2 value = triangle h=.8  i = stepjl color = red;
legend1 label=none value=(height=1 font=swiss 'Kaplan-Meier' 'Nelson-Aalen' )
position=(top right inside) mode=share cborder=black;
axis1 label=(a=90 'Survival Probability') order=(0 to 1.0 by .25);

proc gplot data=na;
plot (s shat )*time /overlay legend=legend1 vaxis=axis1;
run;
quit;



The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.