Table 13.1, page 471Sorting the data to get it in the same order as in the table in the book.
proc sort data=datasets.honking out=sort_honk; by seconds; run; proc print data=sort_honk; id id; var seconds censor; run; ID SECONDS CENSOR 9 1.41 0 54 1.41 1 40 1.51 0 50 1.67 0 60 1.68 0 7 1.86 0 17 2.12 0 31 2.19 0 3 2.36 1 56 2.48 0 5 2.50 0 22 2.53 0 23 2.54 0 45 2.56 0 27 2.62 0 4 2.68 0 49 2.76 1 21 2.78 1 13 2.83 0 ...
Table 13.2, page 477Life table for the horn-honking data with discrete-time and actuarial estimates of the survivor functions. First, we will produce the survivor and hazard functions using actuarial estimates.
The option method allows us to specify that we wish to obtain the actuarial estimates.
ods output LifetableEstimates=act;
proc lifetest data=datasets.honking method=lt intervals=1 to 8 by 1, 18;
time seconds*censor(1);
run;
proc print data=act noobs;
var lowertime uppertime failed censored survival stderr hazard hazstderr;
run;
<output omitted>
Lower Upper HazStd
Time Time Failed Censored Survival StdErr Hazard Err
0 1 0 0 1.0000 0 0 .
1 2 5 1 1.0000 0 0.092593 0.041364
2 3 14 3 0.9115 0.0378 0.329412 0.086837
3 4 9 2 0.6537 0.0643 0.315789 0.103943
4 5 6 4 0.4754 0.0690 0.333333 0.134179
5 6 2 2 0.3396 0.0680 0.181818 0.128033
6 7 2 2 0.2830 0.0674 0.285714 0.199958
7 8 1 0 0.2122 0.0666 0.222222 0.220846
8 18 3 1 0.1698 0.0654 0.15 0.057282
We now produce the estimates for the discrete-time survivor function.
In the data step we are creating a categorical time variable which is used in the following proc lifetest from which we output the product limit estimates into a data set called est. We remove the redundant observations with missing values for the survival function from the data set and then print.
data honking; set datasets.honking; if seconds < 2 then time = 1; else if seconds < 3 then time = 2; else if seconds < 4 then time = 3; else if seconds < 5 then time = 4; else if seconds < 6 then time = 5; else if seconds < 7 then time = 6; else if seconds < 8 then time = 7; else time = 8; run; ods output ProductLimitEstimates=est; ods listing close; proc lifetest data=honking; time time*censor(1); run; ods listing;
data est_dt; set est (where=(survival~=.)); lags = lag(survival); pt =1- survival /lags; lagtime = lag(time); if time = 8 then time = 17; hazard = pt/(time - lagtime); keep time pt survival stderr hazard; run;
proc print data = est_dt noobs; format time 3.0; run;
time Survival StdErr pt hazard
0 1.0000 0 . . 1 0.9123 0.0375 0.08772 0.08772 2 0.6619 0.0632 0.27451 0.27451 3 0.4867 0.0683 0.26471 0.26471 4 0.3597 0.0673 0.26087 0.26087 5 0.3044 0.0674 0.15385 0.15385 6 0.2367 0.0673 0.22222 0.22222 7 0.1894 0.0685 0.20000 0.20000 17 0.0473 0.0444 0.75000 0.07500
Figure 13.1, page 479Left panel based on discrete-time model. The data set created est_dt for the previous example then can be used here for the plots.
symbol i = join v=star;
axis1 order = (0 to 20 by 5) minor = none
label = ('Seconds after light turns green');
axis2 order = (0 to 1 by .25) minor = none ;
axis3 order = (0 to .35 by .05) minor = none;
proc gplot data = est_dt;
format time 4.0 survival 4.2 hazard 4.2;
plot survival*time /haxis=axis1 vaxis = axis2;
plot hazard*time /haxis=axis1 vaxis=axis3;
run;
quit;
Right panel based on the actuarial estimates.
symbol c=blue i=stepjr;
axis1 label=(a=90 'S(tj)') order=(0 to 1.00 by 0.25) minor=none;
axis2 order=(0 to 20 by 5) minor=none label=('Seconds after light turns green');
axis3 order=(0 to .35 by .05) minor = none;
proc gplot data=act;
format survival 4.2;
plot survival*lowertime / vaxis=axis1 haxis=axis2;
plot hazard*lowertime /haxis=axis2 vaxis=axis3;
run;
quit;
Table 13.3, page 484We use the ods output statement to output the Kaplan-Meier estimates into a data set called km. Then we create variables for Table 13.3 and drop the observations where survival is missing.
ods output ProductLimitEstimates=km; ods listing close; proc lifetest data=datasets.honking method=km; time seconds*censor(1); run; ods listing;
data table13_3; set km ; retain myw 0; start = lag(seconds); end = seconds; width = end - start; lagc = lag(censor); if censor = 1 then do; if lagc = 0 then myw = width; else myw = width + myw; end; else if lagc = 1 then myw = width+myw; else myw = 0; lagf = lag(failed); lagl = lag(left); pt = (failed - lagf)/lagl ; lagp = lag(pt); if lagc = 1 and censor = 0 then width = myw; if survival ~=. ; keep survival stderr start end pt width; run;
data table13_3a;
set table13_3;
lagp = lag(pt);
hazard = lagp/width;
surv = lag(survival);
se = lag(stderr);
drop lagp survival stderr;
if _n_>=2;
run;
proc print data = table13_3a noobs;
var start end pt surv se width hazard;
run;
start end pt surv se width hazard
0.00 1.41 0.01754 1.00000 0.000000 1.41 . 1.41 1.51 0.01818 0.98246 0.017389 0.10 0.17544 1.51 1.67 0.01852 0.96459 0.024592 0.16 0.11364 1.67 1.68 0.01887 0.94673 0.029929 0.01 1.85185 1.68 1.86 0.01923 0.92887 0.034283 0.18 0.10482 1.86 2.12 0.01961 0.91100 0.037993 0.26 0.07396 2.12 2.19 0.02000 0.89314 0.041234 0.07 0.28011 2.36 2.48 0.02083 0.87528 0.044109 0.29 0.06897 2.48 2.50 0.02128 0.85704 0.046808 0.02 1.04167 2.50 2.53 0.02174 0.83881 0.049236 0.03 0.70922 2.53 2.54 0.02222 0.82057 0.051432 0.01 2.17391 2.54 2.56 0.02273 0.80234 0.053424 0.02 1.11111 2.56 2.62 0.02326 0.78410 0.055234 0.06 0.37879 2.62 2.68 0.02381 0.76587 0.056880 0.06 0.38760 2.78 2.83 0.02564 0.74763 0.058376 0.15 0.15873 2.83 2.88 0.02632 0.72846 0.059944 0.05 0.51282 2.88 2.89 0.02703 0.70929 0.061355 0.01 2.63158 2.89 2.92 0.02778 0.69012 0.062620 0.03 0.90090 2.92 2.98 0.02857 0.67095 0.063747 0.06 0.46296 3.05 3.14 0.03030 0.65178 0.064744 0.16 0.17857 3.14 3.17 0.03125 0.63203 0.065726 0.03 1.01010 3.17 3.21 0.03226 0.61228 0.066574 0.04 0.78125 3.21 3.22 0.03333 0.59253 0.067292 0.01 3.22581 3.22 3.24 0.03448 0.57278 0.067886 0.02 1.66667 3.46 3.56 0.03704 0.55303 0.068358 0.32 0.10776 3.56 3.57 0.03846 0.53255 0.068826 0.01 3.70370 3.57 3.58 0.04000 0.51206 0.069160 0.01 3.84615 3.58 3.78 0.04167 0.49158 0.069360 0.20 0.20000 4.01 4.10 0.04545 0.47110 0.069429 0.32 0.13021 4.10 4.18 0.04762 0.44968 0.069497 0.08 0.56818 4.30 4.44 0.05263 0.42827 0.069408 0.26 0.18315 4.44 4.51 0.05556 0.40573 0.069318 0.07 0.75188 4.51 4.52 0.05882 0.38319 0.069035 0.01 5.55556 4.71 4.96 0.07143 0.36065 0.068555 0.44 0.13369 5.12 5.39 0.08333 0.33489 0.068327 0.43 0.16611 5.39 5.73 0.09091 0.30698 0.068094 0.34 0.24510 5.88 6.03 0.11111 0.27907 0.067381 0.30 0.30303 6.21 6.30 0.14286 0.24807 0.066648 0.27 0.41152 6.60 7.20 0.20000 0.21263 0.065878 0.90 0.15873 7.20 9.59 0.25000 0.17010 0.064994 2.39 0.08368 9.59 12.29 0.33333 0.12758 0.061094 2.70 0.09259 12.29 13.18 0.50000 0.08505 0.053521 0.89 0.37453
Figure 13.2, page 485All the data sets containing the survival functions (est, act and km) are sorted and merged in order to graph all the survival functions in one graph.
proc sort data=est;
by seconds;
run;
proc sort data=act;
by seconds;
run;
proc sort data=km;
by seconds;
run;
data combo;
merge km act est;
by seconds;
run;
goptions reset=all;
symbol1 c=blue i=stepjll;
symbol2 c=red i=stepjr line=21;
symbol3 c=green i=j line=3;
axis1 label=(a=90 'S(tj)');
axis2 order=(0 to 20 by 5);
legend label=none value=(height=1 font=swiss 'Kaplan-Meier' 'Actuarial' 'Discrete-time' )
position=(top right inside) mode=share cborder=black;
proc gplot data=combo;
format survkm seconds 4.2;
plot survkm*seconds / vaxis=axis1 haxis=axis2;
plot (survkm survact survival)*seconds / overlay vaxis=axis1 haxis=axis2 legend=legend;
run;
quit;
Figure 13.4, page 493The Nelson-Aalen and negative log survivor function estimates of the cumulative hazard function for the horn-honking dataset.
In the output statement of proc phreg the option method is used to specify the Nelson-Aalen estimates. These estimates are outputted into a data set called nelson which is then sorted and all the redundant observations are dropped.
proc phreg data=datasets.honking noprint; model seconds*censor(1) = ; output out=nelson logsurv=logs / method=emp; run; proc sort data=nelson; by seconds; run; data nelson; set nelson; nlogs = -logs; by seconds; if first.seconds; run;
In this proc phreg we use the output statement to output the product-limit estimates of the log survival function to a data set called pl. The pl data set is then sorted by seconds and the redundant observations are dropped. The data sets pl and nelson are then merged and the survival functions are plotted.
proc phreg data=datasets.honking noprint;
model seconds*censor(1) = ;
output out=pl logsurv=logs;
run;
proc sort data=pl;
by seconds;
run;
data pl;
set pl;
nlogs_pl = -logs;
by seconds;
if first.seconds;
run;
data combo;
merge nelson pl;
by seconds;
run;
goptions reset=all;
symbol1 c=blue i=stepjll;
symbol2 c=red i=stepjl line=21;
axis1 label=(a=90 'S(tj)') order=(0 to 3.5 by .5);
axis2 order=(0 to 20 by 5);
legend label=none value=(height=1 font=swiss 'Negative Log survivor' 'Nelson-Aalen' )
position=(top left inside) mode=share cborder=black;
proc gplot data=combo;
plot (nlogs_pl nlogs)*seconds / overlay vaxis=axis1 haxis=axis2 legend=legend;
run;
quit;
Figure 13.5, page 496To generate these Kernel-smoothed hazard functions, we modified a SAS macro supplied by the authors and written by Paul Allison. The full macro can be found here.
proc lifetest data=datasets.honking method=km outsurv = test;
time seconds*censor(1);
run;
data _inset_;
set test end=final;
retain _grp_ _censor_ 0;
t=seconds;
if t=0 and survival=1 then _grp_=_grp_+1;
keep _grp_ t survival;
if final and _grp_ > 1 then call symput('nset','yes');
else if final then call symput('nset','no');
if _censor_ = 1 then delete;
if survival in (0,1) then delete;
run;
%macro smooth(width);
proc iml;
use _inset_;
read all var {t _grp_};
w2=&width;
w=w2;
z=char(w,8,2);
call symput('width',z);
numset=max(_grp_);
create _plt_ var{ lambda s group};
setin _inset_ ;
do m=1 to numset;
read all var {t survival _grp_} where (_grp_=m);
n=nrow(survival);
lo=t[1] + w;
hi=t[n] - w;
npt=50;
inc=(hi-lo)/npt;
s=lo+(1:npt)`*inc;
group=j(npt,1,m);
slag=1//survival[1:n-1];
h=1-survival/slag;
x = (j(npt,1,1)*t` - s*j(1,n,1))/w;
k=.75*(1-x#x)#(abs(x)<=1);
lambda=k*h/w;
append;
end;
quit;
proc gplot data=_plt_;
plot lambda*s / vaxis=axis1 vzero haxis=axis2;
axis1 label=(angle=90 f=titalic 'Hazard Function' ) minor=none ;
axis2 label=(f=titalic "Time (bandwidth=&width)") minor=none;
symbol1 i=join color=black line=1;
symbol2 i=join color=red line=2;
symbol3 i=join color=green line=3;
symbol4 i=join color=blue line=4;
run;
quit;
%mend;
%smooth(1)[http://www.ats.ucla.edu/stat/sas/footer.htm]%smooth(2)
%smooth(3)
