*************************************************************************
*            CALCULATION OF THE NUMBER NEEDED TO BE EXPOSED (NNE)       *
*                       AND 95% CONFIDENCE INTERVALS                    *
*               WITH ADJUSTMENT FOR CONFOUNDING VARIABLES               *
*                                                                       *
* REFERENCE:  BENDER & BLETTNER (2002), J. CLIN. EPIDEMIOL. 55: 525-530 *
*************************************************************************;

  *------------------------------------------------------------------*
  |                           QUESTIONS ?                            |
  |                                                                  |
  | --->  Dr. Ralf BENDER  Phone:  +49 221 35685-451                 |
  |                        Fax:    +49 221 35685-891                 |
  |                        E-Mail: Ralf@rbsd.de                      |
  *------------------------------------------------------------------*;
  *------------------------------------------------------------------*
  !!                            NOTE:                               !!
  !!                                                                !!
  !!   This SAS code is appropriate for the SAS Version 6.          !!
  !!   You have to replace the simulated data set below with your   !!
  !!   own data. The SAS name of your data set should be "cohort".  !!
  !!                                                                !!
  !!   This data set must contain the following variables:          !!
  !!     1) RESPONSE = adverse event of interest (1=yes, 0=no)      !!
  !!     2) EXPOSURE = binary exposure variable (1=yes, 0=no)       !!
  !!     3) COVAR_1  = the first confounding variable               !!
  !!     4) If there are more than 1 confounding variables these    !!
  !!        variables must be named COVAR_2, COVAR_3, and so on.    !!
  !!                                                                !!
  !!   In line 39 you have to enter the number of confounding       !!
  !!   variables (e.g. %LET K=3) !                                  !!
  !!                                                                !!
  !!   It is necessary to investigate the goodness-of-fit of the    !!
  !!   logistic regression model to get useful NNE estimates !      !!
  *------------------------------------------------------------------*;

options linesize=105;

*----->; %LET K=1 ; *** K = Number of covariables (excluding exposure);

*--------------------------------------------------------------------*
|                       Simulation of Cohort Data                    |
*--------------------------------------------------------------------*;
title1 'Calculation of Adjusted NNEs with 95% Confidence Intervals';
title2 'Reference:  BENDER & BLETTNER (2002), J. Clin. Epidemiol. 55: 525-530';
title4 'Simulated Cohort Data (n=1000, OR=2, 1 Covariate)';
title5 ' ';
data sim;
  do i=1 to 1000;
    No=i;
    alpha=-10;
    gamma=0.693147;
    beta1=0.2;
    Z=RANBIN(0,1,0.5);
    if Z=0 then X1=5*RANNOR(0)+45;
    if Z=1 then X1=5*RANNOR(0)+41;
    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;

data UE;  set sim; if Z=0;
  alpha=-10;
  beta1=0.2;
  UER=exp(alpha+beta1*X1)/(1+exp(alpha+beta1*X1));
  keep UER;
run;

proc means data=UE noprint;
  var UER;
  output out=UEmean mean=;
run;

data NNE;  set UEmean;
  OR=exp(0.693147);
  ARI = (OR-1)*UER*(1-UER) / (1+UER*(OR-1));
  EER = ARI + UER;
  NNEH = 1/ARI;
  run;

title4 'True UER, OR, ARI, and NNEH';
title5 ' ';
proc print data=NNE;
  var UER OR ARI NNEH;
run;

data cohort;  set sim;
  RESPONSE = Y;
  EXPOSURE = Z;
  COVAR_1  = X1;
  drop Y Z X1 PROB;
run;

title4 'Descriptive Anlysis of Generated Cohort Data';
title5 ' ';
proc means data=cohort n min max mean std sum maxdec=3;
proc means data=cohort n min max mean std sum maxdec=3;
  class EXPOSURE;
run;

proc corr data=cohort spearman;
  var EXPOSURE COVAR_1;
run;

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

*--------------------------------------------------------------------*
|                     Multiple Logistic Regression Model             |
*--------------------------------------------------------------------*;
title4 'Estimation of Adjusted Exposure Effect (Logistic Regression)';
title5 ' ';
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 COVARIAT;   set cohort;
  keep COVAR_1-COVAR_&K;
  if EXPOSURE=0;
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 ' ';
proc IML;
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 covariat; read all var _all_                        into COVARIAT;
use beta;     read all var{INTERCEP}                    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{INTERCEP EXPOSURE}           into COV1;
use covb2_k;  read all var _all_                        into COV2;

  COVB = (COV1 || COV2);

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

  UER   = SUM(RESP_UNE)/NU;
  EER   = SUM(RESP_EXP)/NE;
  ARI   = EER-UER;
  NNEH  = 1/ARI;

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

  print "CRUDE (UNADJUSTED) RESULTS";
  if ARI>0 then do;
    print "Crude event rates, ARI and NNEH:  "
           EER[format=6.4] "  " UER[format=6.4] "  "
           ARI[format=6.4] "  " NNEH[format=8.2] ,,,,;
  end;

  if ARI<0 then do;
    ARR  = -ARI;
    NNEB = 1/ARR;
    print "Crude event rates, ARR and NNEB:  "
           EER[format=6.4] "  " UER[format=6.4] "  "
           ARR[format=6.4] "  " NNEB[format=8.2] ,,,,;
  end;

  if ARI=0 then do;
    NNEH="infinity";
    print "Crude event rates, ARI and NNEH:  "
           EER[format=6.4] "  " UER[format=6.4] "  "
           ARI[format=6.4] "  " NNEH[format=8.2] ,,,,;
  end;


  D10   = 0;
  D11   = exp(BETA1);
  D12_k = 0#BETA2_k;
  D1    = (D10 || D11 || D12_k);

  LINPRED = BETA0 + COVARIAT*BETA2_k`;
  LOGIST2 = EXP(LINPRED)/((1+EXP(LINPRED))##2);
  DO J=1 to K;
    LOGISTJ = COVARIAT[,J] # EXP(LINPRED)/((1+EXP(LINPRED))##2);
    LOGIST  = (LOGIST || LOGISTJ);
  END;

  D20   = SUM(LOGIST2)/NU;
  D21   = 0;
  D22_k = LOGIST[+,]/NU;
  D2    = (D20 || D21 || D22_k);

  D = (D1 // D2);

  COVF = D*COVB*D`;

  OR   = exp(BETA1);
  RISK = EXP(LINPRED)/(1+EXP(LINPRED));
  UER  = SUM(RISK)/Nu;

  dv1 = UER#(1-UER) / ((1+UER#(OR-1))##2);
  dv2 = (OR-1)#(1-2#UER-(OR-1)#UER##2) / ((1+UER#(OR-1))##2);
  dv  = (dv1 || dv2);

  ARI  = (OR-1)#UER#(1-UER) / (1+UER#(OR-1));
  EER  = ARI+UER;

  print "ADJUSTED RESULTS";

  if ARI=0 then do;
    NNEH="infinity";
    print "Absolute risk increase (ARI):" ARI[format=6.4];
    print "Number needed to be exposed to be harmed (NNEH):"
          NNEH[format=8.2] ,,,;
    print "Sorry, but because ARI=0 the delta method is not applicable
           to calculate the standard error of ARI!";
  abort;
  end;

  else NNEH = 1/ARI;

  VAR  = dv*COVF*dv`;
  SE   = SQRT(VAR);

  sn     = 0.05;
  Zsn    = PROBIT(1-SN/2);
  ARI_L  = ARI - Zsn#SE;
  ARI_U  = ARI + Zsn#SE;
  NNEH_L = 1/ARI_U;
  NNEH_U = 1/ARI_L;

  print "OR, exposed (EER) and unexposed event rate (UER): "
         OR[format=7.3] " " EER[format=6.4] " " UER[format=6.4];

  if ARI>0 then do;
    print "Absolute risk increase (ARI) with SE and 95% CI: "
          ARI[format=6.4] " " SE[format=8.6] " "
          ARI_L[format=6.4] ARI_U[format=6.4];
    print "Number needed to be exposed to harm (NNEH) with 95% CI: "
          NNEH[format=8.2] " " NNEH_L[format=8.2] NNEH_U[format=8.2];
  end;

  if ARI<0 then do;
    ARR  = -ARI;
    NNEB = 1/ARR;
    ARR_L  = ARR - Zsn#SE;
    ARR_U  = ARR + Zsn#SE;
    NNEB_L = 1/ARR_U;
    NNEB_U = 1/ARR_L;
    print "Absolute risk reduction (ARR) with SE and 95% CI: "
          ARR[format=6.4] " " SE[format=8.6] " "
          ARR_L[format=6.4] ARR_U[format=6.4];
    print "Number needed to be exposed to benefit (NNEB) with 95% CI: "
         NNEB[format=8.2] " " NNEB_L[format=8.2] NNEB_U[format=8.2];
  end;

quit;
run;