/*------------------------------------------------------------------*
   | MACRO NAME  : senspe
   | SHORT DESC  : Calculate sensitivity and specificity
   *------------------------------------------------------------------*
   | CREATED BY  : Cha, Stephen                  (08/18/2005  9:51)
   *------------------------------------------------------------------*
   | PURPOSE
   |
   | Get a) Freq table with missing options
   | b) calculate sensitivity and specificity with 95% CI
   | c) calculate PPV and PNV with 95% CI
   *------------------------------------------------------------------*
   | OPERATING SYSTEM COMPATIBILITY
   |
   | UNIX SAS v8   :
   | UNIX SAS v9   :   YES
   | MVS SAS v8    :
   | MVS SAS v9    :
   | PC SAS v8     :
   | PC SAS v9     :
   *------------------------------------------------------------------*
   | MACRO CALL
   |
   | %senspe  (
   |            data= ,
   |            testvar= ,
   |            truthvar=
   |          );
   *------------------------------------------------------------------*
   | REQUIRED PARAMETERS
   |
   | Name      : data
   | Default   :
   | Type      : Text
   | Purpose   : SAS dataset
   |
   | Name      : testvar
   | Default   :
   | Type      : Text
   | Purpose   : a test variable (Y=1 and N=1)
   |
   | Name      : truthvar
   | Default   :
   | Type      : Text
   | Purpose   : the gold standard variable (Y=1 and N=0)
   |
   *------------------------------------------------------------------*
   | Copyright 2005 Mayo Clinic College of Medicine.
   |
   | This program is free software; you can redistribute it and/or
   | modify it under the terms of the GNU General Public License as
   | published by the Free Software Foundation; either version 2 of
   | the License, or (at your option) any later version.
   |
   | This program is distributed in the hope that it will be useful,
   | but WITHOUT ANY WARRANTY; without even the implied warranty of
   | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
   | General Public License for more details.
   *------------------------------------------------------------------*/
 
/** Macro to compute the
    a) specificity and sensitivity, and their 95% confidence interval
    b) PPV and PNV and their 95% confidence Interval
 
    *****************************************************************;
    *** SAS MACRO spesen                                          ***;
    *** DATE: 9/02/2003                                           ***;
    *** AUTHOR: CHA, STEPHEN & STEPHANIE BAGNIEWSKI               ***;
    *****************************************************************;
 
    THIS MACRO(SPESEN) COMPUTES THE Specifity,Sensitivity, PPV, and PNV
    and their 95% CONFIDENCE INTERVAL. FOR A 2*2 TABLE WHERE THE ROWS
    AND COLUMNS CONTAIN 2 LEVELS (0=NO, 1=YES).
 
    REFERENCE: BISHOP, FIENBERG AND HOLLAND. DISCRETE MULTIVARIATE
               ANALYSIS: THEORY AND PRACTICE. MIT PRESS(1975).
 
    INPUT:  SAS DATA SET WITH N OBSERVATIONS CONTAINING THE ROW VARIABLE
      AND THE COLUMN VARIABLE ON THE SAME OBSERVATION.
 
    MACRO CALL:  %spesen(DATA,testVAR=,TruthVAR=);
 
    OUTPUT: A)FREQUENCY TABLE (WITH MISSING OPTIONS)
            B)Sensitivity and Specivity and 95% CI
            C)PPV and PNV and 95% CI
 
    example: %spesen(data=temp,testvar=xvar,truthvar=yvar);
**/
 
%macro spesen(data=,testvar=,truthvar=);
 
  **************************************;
  * Take out missing and dump freq     *;
  **************************************;
 
  proc format;
    value ansf  0="No "
                1="Yes";
 
  data temp1; set &data;
  format &testvar &truthvar ansf.;
  proc freq; tables &testvar*&truthvar/nopct missing;
 
  data temp2; set temp1;
  if &testvar=. or &truthvar=. then delete;
  proc freq; tables &testvar*&truthvar/noprint sparse out=dd1;
 
  **************************************;
  * get specificity and Sensitivity    *;
  **************************************;
 
    data tmp1; set dd1;
    id=1;
    proc sort; by id;
 
    data tmp2; set tmp1; by id;
    keep a b c d;
    if first.id then do;
      a=0; b=0; c=0; d=0;
      retain a b c d;
    end;
    if &testvar=1 & &truthvar=1 then A=count;
    if &testvar=1 & &truthvar=0 then B=count;
    if &testvar=0 & &truthvar=1 then C=count;
    if &testvar=0 & &truthvar=0 then D=count;
    if last.id then output;
 
    data tmp3; set tmp2;
    cat="Sensitivity";
    rate=A/(A+C);
    LPCT=(1.0-95.0/100)/2;
    UPCT=1-LPCT;
    Phat=rate;
 
   **NORMAL APPROXIMATION;
    L_normal=PHAT+PROBIT(LPCT)*SQRT(PHAT*(1-PHAT)/(A+C));  *LOWER LIMIT;
    U_normal=PHAT+PROBIT(UPCT)*SQRT(PHAT*(1-PHAT)/(A+C));  *UPPER LIMIT;
 
   **POISSON APPROXIMATION;
    DF_LO=2*A;
    DF_UP=2*(A+1);
    IF A GT 0 THEN L_poissn=GAMINV(LPCT,DF_LO/2)/(A+C);  *LOWER LIMIT;
     ELSE L_POI=0;
    U_poissn=GAMINV(UPCT,DF_UP/2)/(A+C);                 *UPPER LIMIT;
 
   *** EXACT BINOMIAL CONFIDENCE LIMITS ***equn 10.7 of Feller*;
     * Upper limit;
    if A<(A+C) then u_exact=1-betainv(lpct,C,A+1);
     else u_exact=1;
     * Lower limit;
    if A>0 then l_exact=1-betainv(upct,C+1,A);
     else l_exact=0;
    keep cat rate L_normal u_normal l_poissn u_poissn l_exact u_exact;
 
    data tmp4; set tmp2;
    cat="Specificity";
    rate=D/(B+D);
    LPCT=(1.0-95.0/100)/2;
    UPCT=1-LPCT;
    Phat=rate;
 
   **NORMAL APPROXIMATION;
    L_normal=PHAT+PROBIT(LPCT)*SQRT(PHAT*(1-PHAT)/(B+D));  *LOWER LIMIT;
    U_normal=PHAT+PROBIT(UPCT)*SQRT(PHAT*(1-PHAT)/(B+D));  *UPPER LIMIT;
 
   **POISSON APPROXIMATION;
    DF_LO=2*D;
    DF_UP=2*(D+1);
    IF D GT 0 THEN L_poissn=GAMINV(LPCT,DF_LO/2)/(B+D);  *LOWER LIMIT;
     ELSE L_POI=0;
    U_poissn=GAMINV(UPCT,DF_UP/2)/(B+D);                 *UPPER LIMIT;
 
   *** EXACT BINOMIAL CONFIDENCE LIMITS ***equn 10.7 of Feller*;
     * Upper limit;
    if D<(B+D) then u_exact=1-betainv(lpct,B,D+1);
     else u_exact=1;
     * Lower limit;
    if D>0 then l_exact=1-betainv(upct,B+1,D);
     else l_exact=0;
    keep cat rate L_normal u_normal l_poissn u_poissn l_exact u_exact;
 
    data tmp5; set tmp2;
    cat="PPV";
    rate=A/(A+B);
    LPCT=(1.0-95.0/100)/2;
    UPCT=1-LPCT;
    Phat=rate;
 
   **NORMAL APPROXIMATION;
    L_normal=PHAT+PROBIT(LPCT)*SQRT(PHAT*(1-PHAT)/(A+B));  *LOWER LIMIT;
    U_normal=PHAT+PROBIT(UPCT)*SQRT(PHAT*(1-PHAT)/(A+B));  *UPPER LIMIT;
 
   **POISSON APPROXIMATION;
    DF_LO=2*A;
    DF_UP=2*(A+1);
    IF A GT 0 THEN L_poissn=GAMINV(LPCT,DF_LO/2)/(A+B);  *LOWER LIMIT;
     ELSE L_POI=0;
    U_poissn=GAMINV(UPCT,DF_UP/2)/(A+B);                 *UPPER LIMIT;
 
   *** EXACT BINOMIAL CONFIDENCE LIMITS ***equn 10.7 of Feller*;
     * Upper limit;
    if A<(A+B) then u_exact=1-betainv(lpct,B,A+1);
     else u_exact=1;
     * Lower limit;
    if A>0 then l_exact=1-betainv(upct,B+1,A);
     else l_exact=0;
    keep cat rate L_normal u_normal l_poissn u_poissn l_exact u_exact;
 
    data tmp6; set tmp2;
    cat="PNV";
    rate=D/(C+D);
    LPCT=(1.0-95.0/100)/2;
    UPCT=1-LPCT;
    Phat=rate;
 
   **NORMAL APPROXIMATION;
    L_normal=PHAT+PROBIT(LPCT)*SQRT(PHAT*(1-PHAT)/(C+D));  *LOWER LIMIT;
    U_normal=PHAT+PROBIT(UPCT)*SQRT(PHAT*(1-PHAT)/(C+D));  *UPPER LIMIT;
 
   **POISSON APPROXIMATION;
    DF_LO=2*D;
    DF_UP=2*(D+1);
    IF D GT 0 THEN L_poissn=GAMINV(LPCT,DF_LO/2)/(C+D);  *LOWER LIMIT;
     ELSE L_POI=0;
    U_poissn=GAMINV(UPCT,DF_UP/2)/(C+D);                 *UPPER LIMIT;
 
   *** EXACT BINOMIAL CONFIDENCE LIMITS ***equn 10.7 of Feller*;
     * Upper limit;
    if D<(C+D) then u_exact=1-betainv(lpct,C,D+1);
     else u_exact=1;
     * Lower limit;
    if D>0 then l_exact=1-betainv(upct,C+1,D);
     else l_exact=0;
    keep cat rate L_normal u_normal l_poissn u_poissn l_exact u_exact;
 
    data tmp7; set tmp2;
    cat="ACCURACY";
    rate=(A+D)/(A+B+C+D);
    LPCT=(1.0-95.0/100)/2;
    UPCT=1-LPCT;
    Phat=rate;
 
   **NORMAL APPROXIMATION;
    L_normal=PHAT+PROBIT(LPCT)*SQRT(PHAT*(1-PHAT)/(C+D));  *LOWER LIMIT;
    U_normal=PHAT+PROBIT(UPCT)*SQRT(PHAT*(1-PHAT)/(C+D));  *UPPER LIMIT;
 
   **POISSON APPROXIMATION;
    DF_LO=2*D;
    DF_UP=2*(D+1);
    IF D GT 0 THEN L_poissn=GAMINV(LPCT,DF_LO/2)/(C+D);  *LOWER LIMIT;
     ELSE L_POI=0;
    U_poissn=GAMINV(UPCT,DF_UP/2)/(C+D);                 *UPPER LIMIT;
 
   *** EXACT BINOMIAL CONFIDENCE LIMITS ***equn 10.7 of Feller*;
     * Upper limit;
    if D<(C+D) then u_exact=1-betainv(lpct,C,D+1);
     else u_exact=1;
     * Lower limit;
    if D>0 then l_exact=1-betainv(upct,C+1,D);
     else l_exact=0;
    keep cat rate L_normal u_normal l_poissn u_poissn l_exact u_exact;
 
    data output; set tmp3 tmp4 tmp5 tmp6 tmp7;
    rate=rate*100;
    L_normal=L_Normal*100;
    u_normal=u_normal*100;
    l_poissn=l_poissn*100;
    u_poissn=u_poissn*100;
    l_exact=l_exact*100;
    u_exact=u_exact*100;
    label L_normal="Normal Approximation (L)"
          U_normal="95% CI (Upper)"
          L_poissn="Poisson Approximation (L)"
          U_poissn="95% CI (Upper)"
          L_exact="Exact Binomial (L)"
          U_exact="95% CI (Upper)";
    proc print label n;
       var cat rate L_normal u_normal l_poissn u_poissn l_exact u_exact;
    run;
 %mend;