%macro smooth (data=_last_, time=, width=, survival=survival);
/*********************************************************************
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. 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. If there is more than one survival curve in the input data set, 
SMOOTH will produce multiple smoothed hazard curves on the same axes.

There are four parameters:
 
DATA is the name of the data set containing survivor function estimates.
	Default is the mosts 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 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@ssc.upenn.edu
	    
*******************************************************************/

data _inset_;
 set &data end=final;
 retain _grp_ _censor_ 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 _censor_ = 1 then delete;
 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;
  %else %let c=;

proc gplot data=_plt_;
  plot lambda*s &c / 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;