* Program bootitmacro.v2.sas;
* Created by Nancy Cook (BWH, HMS, HSPH);
* Macros to compute reclassification statistics with Bootstrap sampling ;
* For correct standard errors;
* RC stat, NRI, contin NRI, and IDI for BINARY outcomes;
* Includes macro BOOTIT - takes overall bootstrap with no weights;

* Read in reclassification macros - Use correct path;
%include '/home/swift/nancyc/macros/shared/binary/reclass.macros.v4.sas' ;

/********
* Example usage;
%bootit(data,1,100,outxy,pdx,pdxy,4,0.05,0.10,0.20,seed=923);
run;

* Example code to include at end to run and save results in file;
* Change datafile names as necessary;
* To save in files;
data out.rcbootb; set rcbootb; by strat;
data out.nribootb; set nribootb; by strat;
data out.contbootb; set contbootb; by strat;
data out.idibootb; set idibootb; by strat;
run;

* delete to stop appending if you want to rerun the program;
proc datasets;
delete rcbootb nribootb contbootb idibootb;
run;

********/

****************************************************;

%macro BOOTIT(DSNAME,DETAIL,NBOOT,STATVAR,PROB1,PROB2,
       NCAT,C1,C2,C3,C4,C5,C6,C7,C8,C9,SEED=0);
* Macro to compute RC statistic and NRI in bootstrap samples;
* And save results in file;
* Bootstrap of overall sample ;

*  Variables:
*  DSNAME = dataset name;
*  DETAIL = 2 for detailed printout, 1 for limited, 0 for none;
*           Applies to print out - 0 or >0 only (1 or 2);
*           Individual results for reclass not printed;
*  NBOOT = number of bootstrap samples;
*  STATVAR = outcome variable (coded 0,1);
*  PROB1 = probability for model 1;
*  PROB2 = probability for model 2;
*  NCAT = number of categories in classification;
*  C1-C9 = category cutpoints (should have ncat-1 cutpoints);

%put SEED= &seed ;

%let nstrat=1;
* note: leaving nstrat and istrat in although not used here (placeholder);

data _one_; set &DSNAME ;
run;

* Run programs in original data;
%reclass(_one_,&DETAIL,&STATVAR,&PROB1,&PROB2,&NCAT,&C1,&C2,&C3,&C4,&C5,&C6,&C7,&C8,&C9);
data reclassdat; set calstatsb;
strat=1;
run;
%nricat(_one_,&DETAIL,&STATVAR,&PROB1,&PROB2,&NCAT,&C1,&C2,&C3,&C4,&C5,&C6,&C7,&C8,&C9);
data nridat; set _all_;
strat=1;
run;
%contnri(_one_,&DETAIL,&STATVAR,&PROB1,&PROB2);
data contnridat; set _all_;
strat=1;
run;
%idimacro(_one_,&DETAIL,&STATVAR,&PROB1,&PROB2);
data _ididat_; set _idi_;
strat=1;
run;

* options nonotes;

%do istrat= 1 %to &nstrat;
%do iiter= 1 %to &nboot;
%let seed=%eval(&seed+1);

proc surveyselect method=urs data=_one_ out=boot rate=1 seed=&seed outhits noprint;

* Reclass calibration with NCAT categories;
%reclass(boot,0,&STATVAR,&PROB1,&PROB2,&NCAT,&C1,&C2,&C3,&C4,&C5,&C6,&C7,&C8,&C9);
   data calstatsb; set calstatsb;
   strat=&istrat; iter=&iiter;
   if .<pcalib1<0.05 then rcsig1=1;
      else if pcalib1 ne . then rcsig1=0;
   if .<pcalib2<0.05 then rcsig2=1;
      else if pcalib2 ne . then rcsig2=0;
   if .<pcal_adj1<0.05 then rcsig_adj1=1;
      else if pcal_adj1 ne . then rcsig_adj1=0;
   if .<pcal_adj2<0.05 then rcsig_adj2=1;
      else if pcal_adj2 ne . then rcsig_adj2=0;
   * Diff in chisq;
   diffchi=chisq1-chisq2;
   diffadj=chisq_adj1-chisq_adj2;
   keep strat iter ncat maxcat dfmax nccat df chisq1 pcalib1
       chisq_adj1 pcal_adj1 chisq2 pcalib2 chisq_adj2 pcal_adj2 
       rcsig1 rcsig_adj1 rcsig2 rcsig_adj2 diffchi diffadj;
   proc append base=rcbootb data=calstatsb force;
   run;

* NRI with NCAT categories;
%nricat(boot,0,&STATVAR,&PROB1,&PROB2,&NCAT,&C1,&C2,&C3,&C4,&C5,&C6,&C7,&C8,&C9);
   data _all_; set _all_;
   strat=&istrat; iter=&iiter;
   if .<p2nri<0.05 then nrisig=1;
      else if p2nri ne . then nrisig=0;
   keep strat iter ncase ncontr 
   ncat cut1-cut9 up_case down_case up_contr down_contr 
   ri_case ri_contr nri zri_case zri_contr znri p2ri_case p2ri_contr p2nri
   chitrend ptrend chitrend2 ptrend2 nrisig ;
   proc append base=nribootb data=_all_ force;
   run;

* NRI with Infinite Categories;
%contnri(boot,0,&STATVAR,&PROB1,&PROB2);
   data _all_; set _all_;
   strat=&istrat; iter=&iiter;
   if .<p2nri<0.05 then nrisig=1;
      else if p2nri ne . then nrisig=0;
   keep strat iter ncase ncontr 
       up_case down_case up_contr down_contr 
       ri_case ri_contr nri vri_case vri_contr vnri senri nricil nriciu
       zri_case zri_contr znri p2ri_case p2ri_contr p2nri
       chitrend ptrend nrisig;
   proc append base=contbootb data=_all_ force;
   run;

* IDI;
%idimacro(boot,0,&STATVAR,&PROB1,&PROB2);
  data _idi_; set _idi_;
   strat=&istrat; iter=&iiter;
   if .<idip<0.05 then idisig=1;
      else if idip ne . then idisig=0;
   proc append base=idibootb data=_idi_ force;
   run;

%end; * end iiter iteration;
%end; * end istrat iteration;
run;

options notes;

%if &DETAIL>0 %then %do;

/* Output parameter and sample stats */
/* proc print data=rcbootb; by strat; */
title2 "Bootstrap Reclassification Calibration Statistics for &NCAT Categories";
proc means n mean stddev median min max data=rcbootb; by strat;
* format strat stratf.;
run;

/* proc print data=nribootb; by strat; */
title2 "Bootstrap NRI Statistics for &NCAT Categories";
proc means n mean stddev median min max data=nribootb; by strat;
* format strat stratf.;
run;

* Compute CI and p for NRI;
proc univariate data=nribootb noprint; by strat;
var nri ri_case ri_contr;
output out=nridist mean=avenri averi_case averi_contr std=sdnri sdri_case sdri_contr
       pctlpts=2.5 5 10 50 90 95 97.5 pctlpre=nri_q ricase_q ricontr_q;
run;

data nristats; merge nridat nridist; by strat;
znri=nri/sdnri;
pnri=2*(1-probnorm(abs(znri)));
zri_case=ri_case/sdri_case;
pri_case=2*(1-probnorm(abs(zri_case)));
zri_contr=ri_contr/sdri_contr;
pri_contr=2*(1-probnorm(abs(zri_contr)));
run;

proc print data=nristats;
title2 "NRI Bootstrap Statistics and Confidence Interval";
title3 "Note: P-value assumes normality, CI based on &NBOOT bootstrap samples";
var nri sdnri pnri nri_q2_5 nri_q97_5
 ri_case sdri_case pri_case ricase_q2_5 ricase_q97_5
 ri_contr sdri_contr pri_contr ricontr_q2_5 ricontr_q97_5 ;
run;

/* proc print data=contbootb; by strat; */
title2 "Bootstrap NRI Statistics for Infinite Categories (Continuous NRI)";
proc means n mean stddev median min max data=contbootb; by strat;
* format strat stratf.;
run;

proc univariate data=contbootb noprint; by strat;
var nri ri_case ri_contr;
output out=contnridist mean=avenri averi_case averi_contr std=sdnri sdri_case sdri_contr
       pctlpts=2.5 5 10 50 90 95 97.5 pctlpre=nriq ri_caseq ri_contrq;
run;

data contnristats; merge contnridat contnridist; by strat;
znri=nri/sdnri;
pnri=2*(1-probnorm(abs(znri)));
zri_case=ri_case/sdri_case;
pri_case=2*(1-probnorm(abs(zri_case)));
zri_contr=ri_contr/sdri_contr;
pri_contr=2*(1-probnorm(abs(zri_contr)));
run;

proc print data=contnristats;
title2 "Continuous NRI Bootstrap Statistics and Confidence Interval";
title3 "Note: P-value assumes normality, CI based on &NBOOT bootstrap samples";
var nri sdnri pnri nriq2_5 nriq97_5
 ri_case sdri_case pri_case ri_caseq2_5 ri_caseq97_5
 ri_contr sdri_contr pri_contr ri_contrq2_5 ri_contrq97_5 ;
run;

/* proc print data=idibootb; by strat; */
title2 "Bootstrap IDI Statistics";
proc means n mean stddev median min max data=idibootb; by strat;
* format strat stratf.;
run;

proc univariate data=idibootb noprint; by strat;
var idi rel_idi;
output out=ididist mean=aveidi averel_idi std=sdidi sdrel_idi
       pctlpts=2.5 5 10 50 90 95 97.5 pctlpre=idiq rel_idiq;
run;

data idistats; merge _ididat_ ididist; by strat;
zidi=idi/sdidi;
pidi=2*(1-probnorm(abs(zidi)));
zrel_idi=rel_idi/sdrel_idi;
prel_idi=2*(1-probnorm(abs(zrel_idi)));
run;

proc print data=idistats;
title2 "IDI Bootstrap Statistics and Confidence Interval";
title3 "Note: P-value assumes normality, CI based on &NBOOT bootstrap samples";
var idi sdidi pidi idiq2_5 idiq97_5
 rel_idi sdrel_idi prel_idi rel_idiq2_5 rel_idiq97_5 ;
run;

%end;

%mend bootit;

****************************************************;