* Program cstatmacs.v1.sas; * Created by Nancy Cook (BWH, HMS, HSPH); * Macros for c-statistics for BINARY Data; * DETAIL = 2 for detailed printout, 1 for limited, 0 for none; * EXAMPLE USAGE OF MACROS; * %cstat(probs,1,pdxy,outxy); * %corrc(probsxy,1,predx,predxy,outxy); * %ctest(probsxy,1,predx,predxy,outxy); * %delong(probsxy,1,outxy,predx,predxy); ****************************************************************; %macro CSTAT(dat,detail,prob,statvar); proc freq data=&DAT noprint ; tables &PROB*&STATVAR/ measures; output out=meas smdrc; data meas; set meas; * _smdrc_ = Somers D statistic; cstat=(1+_smdrc_)/2; sec=e_smdrc/2; %if &detail>0 %then %do; proc print data=meas; var cstat sec _smdrc_ e_smdrc ; title2 "Measures of Discrimination For &PROB "; title3 "In Data Set = &DAT "; run; %end; * Usage: * %cstat(probs,1,pdxy,outxy); %mend CSTAT; ****************************************************************; %macro CORRC(DSNAME,DETAIL,PROB1,PROB2,OUTC); * Macro to compute SE of difference in c-statistics; * Based on Rosner & Glynn (draft 7/07); * Unstratified; * Variables: * DSNAME = dataset name; * DETAIL = 2 for detailed printout, 1 for limited, 0 for none; * PROB1 = predicted prob model 1; * PROB2 = predicted prob model 2; * OUTC = outcome variable (1=event); title2 "Comparison of c-statistics for &prob1 and &prob2"; proc sort data=&dsname; by &outc; proc rank data=&dsname out=ranks nplus1 ties=mean; by &outc; var &prob1 &prob2; ranks rank1 rank2; run; data ranks; set ranks; by &outc; probit1=probit(rank1); probit2=probit(rank2); proc corr data=ranks outp=outcorr noprint; by &outc; var probit1 probit2 ; run; %if &detail=2 %then %do; proc print data=outcorr; run; %end; data corrp; set outcorr end=eof; by &outc; retain ncase ncontr corr_case corr_contr; if &outc=0 and _TYPE_="N " then ncontr=probit2; else if &outc=1 and _TYPE_="N " then ncase=probit2; if &outc=0 and _NAME_="probit1" then corr_case=probit2; else if &outc=1 and _NAME_="probit1" then corr_contr=probit2; if eof then output; keep ncase ncontr corr_case corr_contr; run; %if &detail>0 %then %do; proc print data=corrp; run; %end; run; title2; run; * Usage: * %corrc(probsxy,1,predx,predxy,outxy); %mend CORRC; ****************************************************************; %macro CTEST(probsxy,detail,predx,predxy,outxy); * Macro to test difference in c-statistics; * Based on Rosner & Glynn (Biometrics 2009); * Computes c-statistics from proc freq; * Uses macros CSTAT and CORRC; * Unstratified; * Variables: * DSNAME = dataset name; * DETAIL = 2 for detailed printout, 1 for limited, 0 for none; * PROB1 = predicted prob model 1; * PROB2 = predicted prob model 2; * OUTC = outcome variable (1=event); %cstat(&probsxy,&detail,&predx,&outxy); data measdx; set meas; keep cstat sec probitcdx varcdx0 ccildx0 cciudx0; rename cstat=cstatdx sec=secdx0; probitcdx=probit(cstat); varcdx0=sec**2; ccildx0= cstat-1.96*sec; cciudx0= cstat+1.96*sec; run; %cstat(&probsxy,&detail,&predxy,&outxy); data measdxy; set meas; keep cstat sec probitcdxy varcdxy0 ccildxy0 cciudxy0; rename cstat=cstatdxy sec=secdxy0; probitcdxy=probit(cstat); varcdxy0=sec**2; ccildxy0= cstat-1.96*sec; cciudxy0= cstat+1.96*sec; run; %corrc(&probsxy,&detail,&predx,&predxy,&outxy); data diffpdc; merge measdx (keep= cstatdx probitcdx varcdx0 secdx0 ccildx0 cciudx0) measdxy (keep= cstatdxy probitcdxy varcdxy0 secdxy0 ccildxy0 cciudxy0) corrp; diffpdc=cstatdxy-cstatdx; * Recompute variance of c from Bernies formula; * Should be higher than null ver from proc freq; varcdx=(cstatdx*(1-cstatdx)+(ncontr+ncase-2)* (probbnrm(probitcdx,probitcdx,0.5)-cstatdx**2))/(ncontr*ncase) ; secdx=sqrt(varcdx); ccildx=cstatdx-1.96*secdx; cciudx=cstatdx+1.96*secdx; varcdxy=(cstatdxy*(1-cstatdxy)+(ncontr+ncase-2)* (probbnrm(probitcdxy,probitcdxy,0.5)-cstatdxy**2))/(ncontr*ncase) ; secdxy=sqrt(varcdxy); ccildxy=cstatdxy-1.96*secdxy; cciudxy=cstatdxy+1.96*secdxy; covc=(probbnrm(probitcdxy,probitcdx,(corr_case+corr_contr)/2) - cstatdxy*cstatdx + (ncase-1)*(probbnrm(probitcdxy,probitcdx,corr_contr/2) - cstatdxy*cstatdx) + (ncontr-1)*(probbnrm(probitcdxy,probitcdx,corr_case/2) - cstatdxy*cstatdx)) / (ncase*ncontr); vardiffpdc=varcdxy + varcdx - 2*covc; sediffpdc=sqrt(vardiffpdc); zdiffpdc=diffpdc/sediffpdc; pdiffpdc=2*(1-probnorm(abs(zdiffpdc))); diffpdccil=diffpdc-1.96*sediffpdc; diffpdcciu=diffpdc+1.96*sediffpdc; run; %if &detail>0 %then %do; proc print data=diffpdc; title2 "Test of c-statistics for &predx and &predxy"; run; %end; title2; run; * Usage: * %ctest(probsxy,1,predx,predxy,outxy); %mend CTEST; ****************************************************************; %macro DELONG(DSNAME,DETAIL,OUTC,PROB1,PROB2); * Macro to compute DeLong test of difference in c-statistics; * Variables: * DSNAME = dataset name; * DETAIL = 2 for detailed printout, 1 for limited, 0 for none; * OUTC = outcome variable (1=event); * PROB1 = predicted prob model 1; * PROB2 = predicted prob model 2; ods listing exclude all; ods html exclude all; ods output ROCAssociation=rocest; ods output ROCContrastEstimate=roccontrast; proc logistic data=&dsname rocoptions(nodetails); model &outc(event='1') = &prob1 &prob2 / nofit; roc 'Prob1' &prob1; roc 'Prob2' &prob2; roccontrast reference('Prob1') / estimate; run; data ests; set rocest; retain roc1 rocse1 roc1cil roc1ciu; if _n_=1 then do; roc1=Area; rocse1=StdErr; roc1cil=LowerArea; roc1ciu=UpperArea; end; else if _n_=2 then do; roc2=Area; rocse2=StdErr; roc2cil=LowerArea; roc2ciu=UpperArea; output; end; keep roc1 rocse1 roc1cil roc1ciu roc2 rocse2 roc2cil roc2ciu; run; data contrast; set roccontrast; diffroc=Estimate; diffse=StdErr; diffcil=LowerCL; diffciu=UpperCL; diffchisq=ChiSquare; diffrocp=ProbChiSq; keep diffroc diffse diffcil diffciu diffchisq diffrocp; run; data delong; merge ests contrast; run; ods listing exclude none; ods html exclude none; %if &detail>0 %then %do; proc print data=delong; title2 "DeLong Test of c-statistics for &prob1 and &prob2 in Data=&dsname"; run; %end; title2; run; * Usage: * %delong(probsxy,1,outxy,predx,predxy); %mend DELONG; ****************************************************************;