%macro smooth (data=_last_, time=, width=, survival=survival);

/*********************************************************************
MACRO SMOOTH produces graphs of smoothed hazard functions using output
from either PROC LIFETEST or PROC PHREG. With PROC LIFETEST, it uses the
data set produced by the OUTSURV option in the PROC statement. With PROC
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. The 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 (the default is SURVIVAL, which is the automatic name in
         PROC LIFETEST).

WIDTH    is bandwidth of smoothing function. The default is 1/5 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=log(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;





/* SAS code to for Kaplan-Meier estimation
   of surviror functions to times from the
   VA lung cancer trial of 137 male patients
   with inoperable lung cancer */

/* Variables
     Treatment: 1=standard, 2=test (chemotherapy)
     Celltype: 1=squamous, 2=smallcell, 3=adeno, 4=large
     Survival in days
     Status: 1=dead, 0=censored
     Karnofsky score
     Months from Diagnosis
     Age in years
     Prior therapy: 0=no, 10=yes
*/


/*
data va;
    infile 'c:\mydocuments\courses\st565\data\va.dat';
    input rx cellt time status karno months age prior_rx;

    prior_rx = prior_rx/10;

    if (cellt=1) then celltype= 'squamous';
    if (cellt=2) then celltype= 'smallcell';
    if (cellt=3) then celltype= 'adeno';
    if (cellt=4) then celltype= 'large';

proc lifetest method=KM plots=(s,ls, lls) graphics outs=su data=va;
 time time*status(0);
 strata rx;
 symbol1 v=none color=black line=1;
 symbol2 v=none color=black line=2;
run;


proc print data=su; run;

* Delete censored cases ;
data su2;
 set su;
 if _CENSOR_=0;
run;



proc gplot data=su2;
 plot SURVIVAL*time SDF_LCL*time SDF_UCL*time/overlay;
 symbol1 v=none interpol=step color=black line=1 w=4;
 symbol2 v=none interpol=step color=black line=2 w=4;
 symbol3 v=none interpol=step color=black line=2 w=4;
 by rx;
run;



%smooth(data=su2, time=time, width=50)

*/