** pass to the macro, TIME , DELTA (censored coded 0), GROUP (label for strata);
** GROUPNUM (number of the group you are running the macro for);
** DATA, BAND (bandwidth for the kernel hazard estimate);

%MACRO HAZARD(TIME=, DELTA=, GROUP=, GROUPNUM=, DATA=_LAST_, BAND=);

PROC SORT DATA=&DATA OUT=NULL_A;
  BY &GROUP &TIME;
RUN;

** Get the maximum value of time, whether censored or not,;
** and regardless of the group, for use in plotting;
PROC MEANS DATA=NULL_A MAX NOPRINT;
  VAR &TIME;
  OUTPUT OUT=MAXTIME MAX=MAXTIME;
RUN;

* limit data to the group of interest;
DATA NULL_B;
  SET NULL_A;
  IF &GROUP = &GROUPNUM;
  TIME=&TIME;
  DELTA=Δ
  GROUP=&GROUP;
RUN;

PROC IML;

USE NULL_B;
  READ ALL VAR {TIME DELTA GROUP} INTO Y;
CLOSE NULL_B;

NUMROW = NROW(Y);    *count number of unique failure time points;
*concat to y, number at risk initialized to 1, number of events initialized to delta;
Y = Y || J(NUMROW,1,1) || Y[,2];
        VARNAME = {TIME DELTA GROUP NRISK NEVENT};
        MATTRIB Y COLNAME=VARNAME;

DO J = 1 TO NUMROW-1;            * process each observation/row;
   Y[J,4] = (NUMROW-J+1);        * decrease number at risk by 1 for each observation;
END;
DO J = 1 TO NUMROW-1;                                       * process each observation/row;
   IF Y[J,1] = Y[(J+1),1] THEN
      DO;                                                   *While time is the same;
        Y[(J+1),5] = SUM(Y[J,5], Y[(J+1),2]);               *add up the number of events;
        Y[(J+1),4] = Y[J,4];                                *retain number at risk from;
      END;                                                  *first obs at that time;
END;


** Reduce the data to rows for each last unique failure time (like LAST. in datasteps);
** Matrix then has the number at risk at the start of the failure time and the total;
** number of events at the same failure time;
** Initialize the new matrix and its subscript count;
Z=Y;                          *Z will eventually have same colums as Y, but fewer rows;
K=1;
DO J=1 TO NUMROW;
  IF J^=NUMROW THEN DO;
        IF (Y[J,1] ^= Y[(J+1),1]) & Y[J,5] > 0 THEN DO;
                Z[K,] =  Y[J,];
                K=K+1;
        END;
  END;   *keep the an obs if its time value is not equal to the next time;
  ELSE IF J=NUMROW THEN DO;
        IF (Y[J,1] ^= Y[(J-1),1]) & Y[J,5] > 0 THEN DO;
                Z[K,] =  Y[J,];
                K=K+1;
        END;
  END;   *keep the last record if its time value is not equal to the previous;
END;

** Keep just the first K-1 rows of Z (k-1 is the count of unique failure times);
Z=Z[1:(K-1),];

********************************************************;
**  Calculate the cum b-line hazard estimator   **;
********************************************************;

USE MAXTIME;
  READ ALL VAR {MAXTIME} INTO MAXTIME;
CLOSE MAXTIME;


DNRATIO = Z[,5] / Z[,4];
DNSRATIO = DNRATIO/(Z[,4]-Z[,5]);
NELSON = CUSUM(DNRATIO);
VARNELS = CUSUM(DNSRATIO);


** Vertically concatenate a value for time =0;
** and a value for the maximum time in the dataset;
** set nelson to 0 at time 0 and repeat the last value;
** of nelson at the maximum time value;
** This way the graphs extend across the observed times;
TIME_   =   0    // Z[,1]   // MAXTIME;
NELSON_ =   0    // NELSON  // NELSON[NROW(Z),1];
VARNELS_=   0    // VARNELS // VARNELS[NROW(Z),1];
GROUP_  = &GROUPNUM // Z[,3]   // &GROUPNUM;
*** Put confidence intervals about the cum hazard function;
UPPERN=NELSON_;
LOWERN=NELSON_;
UPPERN = NELSON_+(1.96*SQRT(VARNELS_));
LOWERN = NELSON_-(1.96*SQRT(VARNELS_));

** Create a SAS dataset with the results with name of dataset indicating the group;
CUMH =  TIME_ || GROUP_ || NELSON_ || VARNELS_ || LOWERN || UPPERN;
        VARNAME = {TIME GROUP NELSON VARNELS LOWERN UPPERN};
        MATTRIB CUMH COLNAME=VARNAME;
CREATE CUMH&GROUPNUM FROM CUMH [COLNAME=VARNAME] ;
*APPEND FROM CUMH [COLNAME=VARNAME];
APPEND FROM CUMH;


**************************************************;
**  get the hazard estimates                    **;
**************************************************;
** Create lots of points across the range of your data.;
** The hazard will be calculated at each of these points;

X = J(200, 1, (MAXTIME/200));
X = .0001 // CUSUM(X);

LAMBDA = J(201,1,0);
BAND=&BAND;
DO J = 1 TO 201;
  DO i = 1 TO NROW(Z);
     T = (X[J]-Z[i,1]) / BAND;
     IF T <= 1 & T >= -1 THEN K = 0.75*(1-T**2);
                         ELSE K = 0;
     LAMBDA[J,1] =  LAMBDA[J,1] + (1/BAND)*K*(DNRATIO[i,1]);
  END;
END;

** Create a SAS dataset with the results;
GROUP=J(201,1,&GROUPNUM) ;
BAND = J(201,1,BAND);
HAZ =  LAMBDA || X || GROUP || BAND;
        VARNAME = {LAMBDA X GROUP BAND};
        MATTRIB HAZ COLNAME=VARNAME;
CREATE HAZ&GROUPNUM FROM HAZ [COLNAME=VARNAME] ;
APPEND FROM HAZ;
*APPEND FROM HAZ [COLNAME=VARNAME];


QUIT;

%MEND HAZARD;