/***************************************************************************************
# SAS MACRO TO CALCULATE GREENWOOD-NAM-D'AGOSTINO CALIBRATION TEST FOR SURVIVAL MODEL
# The most recent version of this code is available at http://ncook.bwh.harvard.edu/sas-macros.html
# FOR MORE DETAILS SEE Demler, Paynter, Cook "Tests of Calibration and Goodness of Fit 
# in the Survival Setting" DOI: 10.1002/sim.6428
# TO RUN:
# GND.calib(datain, groupvar, timevar, eventvar, predpvar)
# PARAMETERS:
# datain   - dataset with the following variables:
# groupvar - GROUPING ASSIGNMENT FOR EACH OBSERVATION
# timevar  - time of the event or censoring time
# eventvar - EVENT INDICATOR 1=EVENT 0=NONEVENT
# predpvar - PREDICTED PROBABILITIES OF AN EVENT CALCULATED FOR TIME=adm_cens 
# adm_cens - TIME AT WHICH PREDICTED PROBABILIRIES ARE CALCULATED
# REQUIRES AT LEAST 2 EVENTS PER GROUP, AT LEAST 5 EVENTS PER GROUP IS RECOMMENDED
# IF <2 EVENTS PER GROUP THEN QUITS
***************************************************************************************/
%macro gnd_calib(datain, groupvar, timevar, eventvar, predpvar, adm_cens);

   * KM estimates should be calculated for a fixed time=adm.cens;
   data &datain; set &datain;
      if &timevar>&adm_cens then do;
	     &timevar=&adm_cens;
		 &eventvar=0;
      end;
   run;

   %let error2=0;
   %let error5=0;

   title; 
   proc sort data=&datain; by &groupvar ; run;

   proc means data=&datain noprint; by &groupvar ; output out=check sum(&eventvar)=nevents; run;

   proc print data=check(drop=_TYPE_ _FREQ_); run;

   data check; set check; if nevents < 2 then call symput('error2','1'); else if nevents < 5 then call symput('error5','1'); 
   run;

   %if &error2=1 %then %do; %put MACRO GND_CALIB WAS STOPPED: at least one of the groups contains <2 events. Consider collapsing some groups.;
      %return;
   %end;
   %else %if &error5=1 %then %put WARNING: at least one of the groups contains < 5 events. GND can become unstable.
      (see Demler, Paynter, Cook 'Tests of Calibration and Goodness of Fit in the Survival Setting')
      Consider collapsing some groups to avoid this problem;


   proc lifetest data=&datain outsurv=datain_surv stderr noprint; time &timevar*&eventvar(0); by &groupvar; run;

   data datain_surv_noncens; set datain_surv; where _censor_=0 and survival ne .; pr_failure=1-survival; run;

   data datain_surv_noncens_; set datain_surv_noncens; by &groupvar;  if last.&groupvar eq 1; run;

   proc means data=&datain noprint; output out=m mean (&predpvar)=m_&predpvar; by &groupvar;  var &predpvar; run;

   data all; merge datain_surv_noncens_ m(keep=&groupvar m_&predpvar); by &groupvar; run;

   data all; set all; gnd_component=(((1-survival)-m_&predpvar)/sdf_stderr)**2; one=1; run;

   proc print data=all; var &groupvar. pr_failure sdf_stderr m_&predpvar gnd_component; run;

   proc means data=all noprint; output out=result sum (gnd_component one)=chi_square ngroups; run; 

   %global _gnd_chisq;
   %global _gnd_df;
   %global _gnd_pvalue;
   data result; set result; df=ngroups-1; p_value=1-probchi(chi_square, df); 
      call symput('_gnd_chisq',chi_square);
      call symput('_gnd_df',df);
      call symput('_gnd_pvalue',p_value);
   run;

   title "GND test: chi_square=%trim(&_gnd_chisq.) df=%trim(&_gnd_df.) p_value=%trim(&_gnd_pvalue.)";
   proc print data=result(drop=_TYPE_ _FREQ_); run;

%mend gnd_calib;