%macro smooth (data=_last_, time=, width=, survival=survival);
/*********************************************************************
(updated 29 July 2014)
This update works with SAS 9.3 or 9.4. It uses PROC SGPLOT
rather than PROC GPLOT. For that reason, it requires that ODS graphics
be turned ON, and it will work with the SAS University Edition.
MACRO SMOOTH produces graphs of smoothed hazard functions using output
from either LIFETEST or PHREG. With LIFETEST, it uses the data set
produced by the OUTSURV option on the PROC statement. With PHREG, it uses
the data set produced by the BASELINE statement. If either procedure
includes a STRATA statement, there will be more than one
survival curve in the data set, and SMOOTH will produce multiple smoothed
hazard curves on the same axes.
Note that PROC LIFETEST now has a built-in option for producing smoothed hazard
functions, with many more features than the SMOOTH macro. But PHREG does not
a built-in option, which is why the SMOOTH macro is still usefl.
SMOOTH employs a kernel smoothing method described by H. Ramlau-Hansen (1983),
"Smoothing Counting Process Intensities by Means of Kernel Functions,"
The Annals of Statistics 11, 453-466. It uses a Gaussian kernel rather than
the default Epanichnikov kernel used by the PLOT option in PROC LIFETEST.
There are four parameters:
DATA is the name of the data set containing survivor function estimates.
Default is the most recently created data set.
TIME is name of the variable containing event times.
SURVIVAL is the name of a variable containing survivor function
estimates (default is SURVIVAL, which is the automatic name in
the OUTSURV data set produced by PROC LIFETEST)
WIDTH is bandwidth of smoothing function. Default is 1/6 of the range
of event times.
Example of usage:
%smooth(data=my.data,time=duration,width=8,survival=s)
Author: Paul D. Allison, University of Pennsylvania
allison@soc.upenn.edu
*******************************************************************/
data _inset_;
set &data end=final;
retain _grp_ 0;
t=&time;
survival=&survival;
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 survival in (0,1,.) then delete;
run;
proc iml;
use _inset_;
read all var {t _grp_};
%if &width ne %then %let w2=&width;
%else %let w2=(max(t)-min(t))/5;
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;
%if &nset = yes %then %let c=group=group;
%else %let c=;
proc sgplot data=_plt_;
series x=s y=lambda / &c ;
yaxis label="Hazard" min=0;
xaxis label="Time (bandwidth=&width)";
run;
quit;
%mend smooth;