* Program cnrimacro.v1.sas for BINARY data;
* Created by Nancy Cook (BWH, HMS, HSPH);
* Macro CNRICAT to compute Adjusted Conditional NRI (CNRI);
* For use in intermediate categories of risk;
* DETAIL = 2 for detailed output, 1 for limited, 0 for none (eg., to save in a file);

* Example Usage:
* %cnricat(probs,1,outxy,pdx,pdxy,0.05,0.20,4,0.05,0.10,0.20);

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

%macro CNRICAT(DSNAME,DETAIL,STATVAR,PROB1,PROB2,LIM1,LIM2,NCAT,C1,C2,C3,C4,C5);
* Macro to compute Conditional NRI (CNRI) adjusted for expected value;
* Uses up to 6 categories with cutpoints c1-c5;
* Control output of frequency tables and means below;
*  Variables:
*  DSNAME = dataset name;
*  DETAIL = 2 for detailed printout, 1 for limited, 0 for none;
*  STATVAR = outcome variable (coded 0,1);
*  PROB1 = probability for model 1;
*  PROB2 = probability for model 2;
*  LIM1 AND LIM2 = limits to define the intermediate range;
*    (Need to equal category cut points);
*  NCAT = number of categories in classification;
*  C1-C5 = category cutpoints (should have ncat-1 cutpoints);
*  Allows up to 6 categories;

title2 Adjusted CNRI for &NCAT Categories with CutPoints = %do j = 1 %to &ncat; &&c&j %end; and Limits &lim1 to &lim2;

data nri1; set &dsname;
* compute diffs in probs;
stat=&STATVAR;
prob1=&PROB1;
prob2=&PROB2;
diffp=prob2-prob1;
if prob1 ne . and prob2 ne .;

* compute categories - set to above cutpoints;
* define limits in terms of number of categories;
%do i=1 %to &ncat-1;
  cut&i=&&c&i;
%end;
cut&ncat=1.01;
    
array cut {6} cut1-cut6;
array pcat {6} pcat1-pcat6;
if .<prob1<cut1 then pcat1=1;
if .<prob2<cut1 then pcat2=1;
do i=2 to &ncat;
   if cut(i-1)<=prob1<cut(i) then pcat1=i;
   if cut(i-1)<=prob2<cut(i) then pcat2=i;
end;
run;

*** Control output as needed ***;
%if &detail=2 %then %do;
proc means data=nri1 n sum mean median stddev min max;
var prob1 prob2 diffp stat;
run;
proc freq data=nri1;
tables pcat1*pcat2 stat*pcat1*pcat2;
run;
%end;

ods listing exclude all;
ods html exclude all;
ods output CrossTabFreqs=freqs;
run;
proc freq data=nri1;
tables stat*pcat1*pcat2 ;
run;
ods listing exclude none;
ods html exclude none;
* proc print data=freqs;
run;

data freqs; set freqs end=eof; by stat;
retain f11-f16 f21-f26 f31-f36 f41-f46 f51-f56 f61-f66 0;
array freq{6,6} f11-f16 f21-f26 f31-f36 f41-f46 f51-f56 f61-f66;
if pcat1 ne . and pcat2 ne . then freq(pcat1,pcat2)=Frequency;

if last.stat then do;
* Compute expected freqs under H0;
array expf {6,6} e11-e16 e21-e26 e31-e36 e41-e46 e51-e56 e61-e66;
do i=1 to &ncat;
do j=1 to &ncat;
  if i=j then expf(i,j)=freq(i,j);
  else expf(i,j) = (freq(i,j)+freq(j,i))/2;
end; end; 
output;
keep stat f11--f66 e11--e66 ;
end;
%if &detail>0 %then %do;
  proc print data=freqs;
  title3 "Observed and Expected Freqs";
%end;

data cnri; set freqs end=eof;
* compute discordance for sequential NRI;
retain ncase ncontr encase encontr nup_ca ndn_ca nup_co ndn_co 
 enup_ca endn_ca enup_co endn_co cklim1 cklim2 0;
array freq{6,6} f11-f16 f21-f26 f31-f36 f41-f46 f51-f56 f61-f66;
array expf{6,6} e11-e16 e21-e26 e31-e36 e41-e46 e51-e56 e61-e66;

cklim1=1; cklim2=1;
 
if &lim1=0 then do; int1=1; cklim1=0; end;
if &lim2=1 then do; int2=&ncat; cklim2=0; end;
%do i=1 %to &ncat-1;
  if &lim1=&&c&i then do; int1=&i+1; cklim1=0; end;
  if &lim2=&&c&i then do; int2=&i; cklim2=0; end;
%end;

* Check category cutpoints;
if int1>0 and int2>0 then do;

do i=int1 to int2;
do j=1 to &ncat;
  if stat=1 then do;
     ncase=ncase+freq(i,j);
     encase=encase+expf(i,j);
     if j>i then do;
       nup_ca=nup_ca+freq(i,j);
       enup_ca=enup_ca+expf(i,j);
     end;
     else if j<i then do;
       ndn_ca=ndn_ca+freq(i,j);
       endn_ca=endn_ca+expf(i,j);
     end;
  end;
  else if stat=0 then do;
     ncontr=ncontr+freq(i,j);
     encontr=encontr+expf(i,j);
     if j>i then do;
       nup_co=nup_co+freq(i,j);
       enup_co=enup_co+expf(i,j);
     end;
     else if j<i then do;
       ndn_co=ndn_co+freq(i,j);
       endn_co=endn_co+expf(i,j);
     end;
  end;
end;
end;
end;
else put "Check Category Limits";

if eof then do;
  ri_case=(nup_ca-ndn_ca)/ncase;
  ri_contr=(ndn_co-nup_co)/ncontr;
  orig_cnri=ri_case+ri_contr;
  eri_case=(enup_ca-endn_ca)/encase;
  eri_contr=(endn_co-enup_co)/encontr;
  cri_case=ri_case-eri_case;
  cri_contr=ri_contr-eri_contr;
  cnri=cri_case+cri_contr;
  vri_case=(nup_ca+ndn_ca)/ncase**2;
  vri_contr=(nup_co+ndn_co)/ncontr**2;
  vnri=vri_case+vri_contr;
  secnri=sqrt(vnri);
  cnricil=cnri-1.96*secnri;
  cnriciu=cnri+1.96*secnri;
  if vri_case ne 0 then zri_case=cri_case/sqrt(vri_case);
  if vri_contr ne 0 then zri_contr=cri_contr/sqrt(vri_contr);
  if vnri ne 0 then znri=cnri/sqrt(vnri);
  p2ri_case=2*(1-probnorm(abs(zri_case)));
  p2ri_contr=2*(1-probnorm(abs(zri_contr)));
  p2nri=2*(1-probnorm(abs(znri)));
  * Use corrected variance;
  * Base probs on expected under null;
  pup_case=enup_ca/encase;
  pdn_case=endn_ca/encase;
  pup_contr=enup_co/encontr;
  pdn_contr=endn_co/encontr;
  vri2_case=(pup_case*(1-pup_case)+pdn_case*(1-pdn_case)+2*pup_case*pdn_case)/ncase;
  vri2_contr=(pup_contr*(1-pup_contr)+pdn_contr*(1-pdn_contr)+2*pup_contr*pdn_contr)/ncontr;
  vnri2=vri2_case+vri2_contr;
  if vnri2 ne 0 then znri2=cnri/sqrt(vnri2);
  p2nri2=2*(1-probnorm(abs(znri2)));

  *** check for no cases or controls or 0 variance;
  * set cnri to mss if so;
  if ncase=0 or ncontr=0 or vnri=0 then do;
     orig_cnri=.; cnri=.; 
  end;
  ***;

  * compute categories again - set to above cutpoints;
  %do i=1 %to &ncat-1;
    cut&i=&&c&i;
  %end;
  if &ncat<6 then do;
    %do i=&ncat %to 6;
      cut&i=.;
    %end;
  end;
  lim1=&LIM1;
  lim2=&LIM2;
  ncat=&NCAT;

  output;
  keep ncase ncontr encase encontr ncat cut1-cut6 lim1 lim2
       nup_ca ndn_ca nup_co ndn_co ri_case ri_contr orig_cnri 
       enup_ca endn_ca enup_co endn_co cri_case cri_contr cnri 
       vri_case vri_contr vnri secnri cnricil cnriciu
       zri_case zri_contr znri p2ri_case p2ri_contr p2nri
       pup_case pdn_case pup_contr pdn_contr
       vri2_case vri2_contr vnri2 znri2 p2nri2
       cklim1 cklim2;
end;

%if &detail>0 %then %do;  
  proc print data=cnri;
  title3;
  run;
%end;

title2; run;
* Usage:
* %cnricat(probs,1,outxy,pdx,pdxy,0.05,0.20,4,0.05,0.10,0.20);
%exit: %mend CNRICAT;

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