*************************************************************************
*            CALCULATION OF THE NUMBER NEEDED TO BE EXPOSED (NNE)       *
*                     AND THE EXPOSURE IMPACT NUMBER (EIN)              *
*                       AND 95% CONFIDENCE INTERVALS                    *
*               WITH ADJUSTMENT FOR CONFOUNDING VARIABLES               *
*                                                                       *
*      REFERENCE:  BENDER ET AL. (2007), STAT. MED. 26: 5586-5595       *
*************************************************************************;

  *------------------------------------------------------------------*
  |                           QUESTIONS ?                            |
  |                                                                  |
  | --->  Prof. Dr. Ralf BENDER  Phone:  +49 221 35685-451           |
  |                              Fax:    +49 221 35685-891           |
  |                              E-Mail: Ralf@rbsd.de                |
  *------------------------------------------------------------------*;
  *------------------------------------------------------------------*
  !!  SYNTAX (last row)                                             !!
  !!  %nne(data= , response= , exposure= , covar= )                 !!
  !!                                                                !!
  !!  where                                                         !!
  !!    data     specifies the data set you are using               !!
  !!    response specicifies the event of interest (1=yes, 0=no)    !!
  !!    exposure specifies the binary exposure (1=yes, 0=no)        !!
  !!    covar    specifies the list of confounders in the model     !!
  *------------------------------------------------------------------*;


options linesize=105 pagesize=60;
proc datasets LIBRARY=WORK KILL;

*--------------------------------------------------------------------*
|                       Simulation of Cohort Data                    |
|    --->  REPLACE THE FOLLOWING DATA STEP WITH YOUR OWN DATA !      |
*--------------------------------------------------------------------*;
title4 'Simulated Cohort Data (n=5000, OR=3, 1 Covariate)';
title5 ' ';
data sim;
  do i=1 to 5000;
    No=i;
    alpha=-10;
    gamma=1.0986123;
    beta1=0.1823216;
    Z=RANBIN(0,1,0.5);
    if Z=0 then X1=7*RANNOR(0)+45;
    if Z=1 then X1=7*RANNOR(0)+40;
    prob=exp(alpha+gamma*Z+beta1*X1)/(1+exp(alpha+gamma*Z+beta1*X1));
    Y=RANBIN(0,1,prob);
    output;
  end;
  drop i alpha gamma beta1;
run;


%macro nne(data=,response=,exposure=,covar=);
title1 'CALCULATION OF ADJUSTED NNE AND ADJUSTED EIN WITH 95% CONFIDENCE INTERVALS';
*--------------------------------------------------------------------*
|                       Creation of Data Set                         |
*--------------------------------------------------------------------*;
data cohort;
   set &data;
   response=&response;
   exposure=&exposure;
   %let i=1;
   %let x=%qscan(&covar,&i,%str( ));
   %do %while (%length(%trim(&x))>0);
     covar_&i=&x;
     %let i=%eval(&i+1);
     %let x=%qscan(&covar,&i,%str( ));
     %let K=%eval(&i-1);
  %end;
run;

*--------------------------------------------------------------------*
|                       Descriptive Anlysis                          |
*--------------------------------------------------------------------*;
title4 'Descriptive Anlysis';
title5 ' ';
title6 ' ';

proc means data=cohort n min max mean std sum maxdec=3 fw=8;
  var RESPONSE EXPOSURE COVAR_1-COVAR_&K;
proc means data=cohort n min max mean std sum maxdec=3 fw=8;
  class EXPOSURE;
  var RESPONSE COVAR_1-COVAR_&K;
run;

*--------------------------------------------------------------------*
|                         Simple 2x2 Table                           |
*--------------------------------------------------------------------*;
title4 'Simple 2x2 Table with Crude Estimation of OR';
title5 ' ';
title6 ' ';
proc freq data=cohort;
  tables EXPOSURE*RESPONSE / nopercent nocol cmh exact;
run;

*--------------------------------------------------------------------*
|                     Multiple Logistic Regression Model             |
*--------------------------------------------------------------------*;
title4 'Estimation of Adjusted Odds Ratio (Multiple Logistic Regression)';
title5 ' ';
title6 ' ';
ods output ConvergenceStatus=convstat;
proc logistic data=cohort descending outest=est covout;
  model RESPONSE = EXPOSURE COVAR_1-COVAR_&K / rl lackfit rsq;
  output out=risks p=risk l=low u=up xbeta=linpred;
run;

*--------------------------------------------------------------------*
|                      Data Sets Used for SAS/IML                    |
*--------------------------------------------------------------------*;
data UNEXP;   set cohort;
  keep COVAR_1-COVAR_&K;
  if EXPOSURE=0;
run;

data EXPON;   set cohort;
  keep COVAR_1-COVAR_&K;
  if EXPOSURE=1;
run;

data beta;   set est;   if _type_='PARMS';
  drop _LINK_ _TYPE_ _NAME_ _LNLIKE_;
run;

data beta2_K;   set beta;
  keep COVAR_1-COVAR_&K;
run;

data COVB;   set est;   if _type_='COV';
  drop _LINK_ _TYPE_ _NAME_ _LNLIKE_;
run;

data covb2_k;   set covb;
  keep COVAR_1-COVAR_&K;
run;

*--------------------------------------------------------------------*
|                           PROC IML                                 |
*--------------------------------------------------------------------*;
title4 'Crude and Adjusted Results Calculated with SAS/IML';
title5 'Method: Average Risk Difference Approach (Bender et al., 2007) ';
title6 ' ';

proc IML;

use convstat; read all var{Status} into STATUS;
if STATUS ^= 0 then do;
  print "The logistic regression model did not converge!",
        "Thus, it is not possible to calculate meaningful adjusted NNEs or EINs based upon this model.",,
        "Please build a new model and try again.";
  abort;
end;

use cohort;   read all var{RESPONSE}                    into RESPONSE;
use cohort;   read all var{RESPONSE} where (EXPOSURE=0) into RESP_UNE;
use cohort;   read all var{RESPONSE} where (EXPOSURE=1) into RESP_EXP;
use cohort;   read all var{EXPOSURE}                    into EXPOSURE;
use unexp;    read all var _all_                        into UNEXP;
use expon;    read all var _all_                        into EXPON;
use beta;     read all var{Intercept}                   into BETA0;
use beta;     read all var{EXPOSURE}                    into BETA1;
use beta2_k;  read all var _all_                        into BETA2_k;
use covb;     read all var{Intercept EXPOSURE}          into COV1;
use covb2_k;  read all var _all_                        into COV2;

  COVB = (COV1 || COV2);

  N     = NROW(RESPONSE);
  EVENT = SUM(RESPONSE);
  Ne    = NROW(RESP_EXP);
  Nu    = NROW(RESP_UNE);
  Ee    = SUM(RESP_EXP);
  Eu    = SUM(RESP_UNE);
  K     = NCOL(BETA2_k);

  UER   = SUM(RESP_UNE)/NU;
  EER   = SUM(RESP_EXP)/NE;
  OR    = (EER*(1-UER))/(UER*(1-EER));
  RD    = EER-UER;
  NNE   = 1/RD;

  print "GENERAL INFORMATION";
  print "Sample size and number of events: "
         N " " EVENT;
  print "Number of unexposed and exposed subjects with events: "
         NU EU " " NE EE ,,,,;

  print "CRUDE (UNADJUSTED) RESULTS";
  if RD>0 then do;
    print "Crude results refer to a harmful effect of exposure";
    print "Crude event rates, OR, risk difference and NNE:    "
           UER[format=6.4] "  " EER[format=6.4] "    "
           OR [format=6.4] "    " RD[format=6.4] "  " NNE[format=8.2] ,,,,;
  end;

  if RD<0 then do;
    RD  = -RD;
    NNE = 1/RD;
    print "Crude results refer to a beneficial effect of exposure";
    print "Crude event rates, OR, risk difference and NNE:    "
           UER[format=6.4] "  " EER[format=6.4] "    "
           OR [format=6.4] "    " RD[format=6.4] "  " NNE[format=8.2] ,,,,;
  end;

  if RD=0 then do;
    NNE="infinity";
    print "Crude event rates, risk difference and NNE:  "
           UER[format=6.4] "  " EER[format=6.4] "  "
           RD [format=6.4] "  " NNE ,,,,;
  end;

  print "ADJUSTED RESULTS IN TERMS OF NNE (AVERAGE RISK DIFFERENCE APPROACH, BASIS: UNEXPOSED PERSONS)";
  print "Note: This method describes the exposure effect for a population with a confounder distribution like the unexposed persons.";

  LINPRED1 = BETA0 + UNEXP*BETA2_k`;
  LINPRED2 = BETA0 + BETA1 + UNEXP*BETA2_k`;
  LOGIST1  = EXP(LINPRED1) / ((1+EXP(LINPRED1))##2);
  LOGIST2  = EXP(LINPRED2) / ((1+EXP(LINPRED2))##2);
  d1 = SUM(LOGIST2-LOGIST1)/Nu;
  d2 = SUM(LOGIST2)/Nu;

  DO J=1 to K;
    LOGCOVJ = UNEXP[,J] # (LOGIST2-LOGIST1);
    LOGCOV  = (LOGCOV || LOGCOVJ);
  END;

  d3 = LOGCOV[+,]/Nu;
  d  = (d1 || d2 || d3);

  RISK1 = EXP(LINPRED1)/(1+EXP(LINPRED1));
  RISK2 = EXP(LINPRED2)/(1+EXP(LINPRED2));
  UER   = SUM(RISK1)/Nu;
  EER   = SUM(RISK2)/Nu;
  ARD   = EER-UER;

  if ARD=0 then do;
    NNE="infinity";
    print "Average risk difference (ARD):" ARD[format=6.4];
    print "Number needed to be exposed (NNE):"
          NNE   ,,,;
    print "Sorry, but because ARD=0 the delta method is not applicable to calculate the standard error of ARD!";
    abort;
  end;

  else NNE = 1/ARD;
  VAR  = d*COVB*d`;
  SE   = SQRT(VAR);

  sn     = 0.05;
  Zsn    = PROBIT(1-SN/2);
  ARD_L  = ARD - Zsn#SE;
  ARD_U  = ARD + Zsn#SE;
  NNE_L = 1/ARD_U;
  NNE_U = 1/ARD_L;

  print "Observed unexposed (UER) and hypothetical exposed event rate (EER): "
         UER[format=6.4] " " EER[format=6.4];

  if ARD>0 then do;
    print "Adjusted results refer to a harmful effect of exposure";
    print "Average risk differencee (ARD) with SE and 95% CI: "
          ARD[format=6.4] " " SE[format=8.6] " "
          ARD_L[format=6.4] ARD_U[format=6.4];
    print "Number needed to be exposed (NNE) with 95% CI: "
          NNE[format=8.2] " " NNE_L[format=8.2] NNE_U[format=8.2];
  end;

  if ARD<0 then do;
    ARD  = -ARD;
    NNE  = 1/ARD;
    ARD_L  = ARD - Zsn#SE;
    ARD_U  = ARD + Zsn#SE;
    NNE_L = 1/ARD_U;
    NNE_U = 1/ARD_L;
    print "Adjusted results refer to a beneficial effect of exposure";
    print "Average risk difference (ARD) with SE and 95% CI: "
          ARD[format=6.4] " " SE[format=8.6] " "
          ARD_L[format=6.4] ARD_U[format=6.4];
    print "Number needed to be exposed with 95% CI: "
         NNE[format=8.2] " " NNE_L[format=8.2] NNE_U[format=8.2];
  end;

  if NNE_U<0 then do;
    print "Note: A negative upper confidence limit for NNE means that the confidence interval contains 2 areas:" ,
          "[lower limit to infinity[ and ]-infinity to upper limit]!";
  end;


  print "ADJUSTED RESULTS IN TERMS OF EIN (AVERAGE RISK DIFFERENCE APPROACH, BASIS: EXPOSED PERSONS)";
  print "Note: This method describes the exposure effect for a population with a confounder distribution like the exposed persons.";

  LINPRED1 = BETA0 + EXPON*BETA2_k`;
  LINPRED2 = BETA0 + BETA1 + EXPON*BETA2_k`;
  LOGIST1  = EXP(LINPRED1) / ((1+EXP(LINPRED1))##2);
  LOGIST2  = EXP(LINPRED2) / ((1+EXP(LINPRED2))##2);
  d1 = SUM(LOGIST2-LOGIST1)/Ne;
  d2 = SUM(LOGIST2)/Ne;

  DO J=1 to K;
    LOGCOVJ = EXPON[,J] # (LOGIST2-LOGIST1);
    LOGCOV2  = (LOGCOV2 || LOGCOVJ);
  END;

  d3 = LOGCOV2[+,]/Ne;
  d  = (d1 || d2 || d3);

  RISK1 = EXP(LINPRED1)/(1+EXP(LINPRED1));
  RISK2 = EXP(LINPRED2)/(1+EXP(LINPRED2));
  UER   = SUM(RISK1)/Ne;
  EER   = SUM(RISK2)/Ne;
  ARD   = EER-UER;

  if ARD=0 then do;
    EIN="infinity";
    print "Average risk difference (ARD):" ARD[format=6.4];
    print "Exposure impact number (EIN):"
          EIN ,,,;
    print "Sorry, but because ARD=0 the delta method is not applicable to calculate the standard error of ARD!";
  abort;
  end;

  else EIN = 1/ARD;
  VAR  = d*COVB*d`;
  SE   = SQRT(VAR);

  sn     = 0.05;
  Zsn    = PROBIT(1-SN/2);
  ARD_L  = ARD - Zsn#SE;
  ARD_U  = ARD + Zsn#SE;
  EIN_L = 1/ARD_U;
  EIN_U = 1/ARD_L;

  print "Hypothetical unexposed (UER) and observed exposed event rate (EER): "
         UER[format=6.4] " " EER[format=6.4];

  if ARD>0 then do;
    print "Adjusted results refer to a harmful effect of exposure";
    print "Average risk difference (ARD) with SE and 95% CI: "
          ARD[format=6.4] " " SE[format=8.6] " "
          ARD_L[format=6.4] ARD_U[format=6.4];
    print "Exposure impact number (EIN) with 95% CI: "
          EIN[format=8.2] " " EIN_L[format=8.2] EIN_U[format=8.2];
  end;

  if ARD<0 then do;
    ARD  = -ARD;
    EIN  = 1/ARD;
    ARD_L  = ARD - Zsn#SE;
    ARD_U  = ARD + Zsn#SE;
    EIN_L = 1/ARD_U;
    EIN_U = 1/ARD_L;
    print "Adjusted results refer to a beneficial effect of exposure";
    print "Average risk difference (ARD) with SE and 95% CI: "
          ARD[format=6.4] " " SE[format=8.6] " "
          ARD_L[format=6.4] ARD_U[format=6.4];
    print "Exposure impact number (EIN) with 95% CI: "
          EIN[format=8.2] " " EIN_L[format=8.2] EIN_U[format=8.2];
  end;

  if EIN_U<0 then do;
    print "Note: A negative upper confidence limit for EIN means that the confidence interval contains 2 areas:" ,
          "[lower limit to infinity[ and ]-infinity to upper limit]!";
  end;

quit;
run;
%mend nne;

%nne(data=sim,response=Y,exposure=Z,covar=X1);