*************************************************************************
*             CALCULATION OF THE NUMBER NEEDED TO TREAT (NNT)           *
*                     AND 95% CONFIDENCE INTERVALS                      *
*             WITH ADJUSTMENT FOR BALANCED COVARIATES IN RCTS           *
*                                                                       *
*      REFERENCES:  BENDER ET AL. (2007), STAT. MED. 26: 5586-5595      *
*                   BENDER & VERVÖLGYI (2010), CCT 31: in press         *
*************************************************************************;

  *------------------------------------------------------------------*
  |                           QUESTIONS ?                            |
  |                                                                  |
  | --->  Prof. Dr. Ralf BENDER  Phone:  +49 221 35685-451           |
  |                              Fax:    +49 221 35685-891           |
  |                              E-Mail: Ralf@rbsd.de                |
  *------------------------------------------------------------------*;
  *------------------------------------------------------------------*
  !!  SYNTAX (last row)                                             !!
  !!  %nnt(data= , response= , treatment= , covar= )                !!
  !!                                                                !!
  !!  where                                                         !!
  !!    data      specifies the data set you are using              !!
  !!    response  specicifies the event of interest (1=yes, 0=no)   !!
  !!    treatment specifies the binary treatment (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 RCT Data                       |
|    --->  REPLACE THE FOLLOWING DATA STEP WITH YOUR OWN DATA !      |
*--------------------------------------------------------------------*;
title4 'Simulated Data';
title5 ' ';
data sim;
  do i=1 to 1000;
    No=i;
    alpha=-10;
    gamma=log(0.50);
    beta1=log(1.2);
    X1=5*RANNOR(0)+50;
	if No<=500 then Z=0;
    if No> 500 then Z=1;
    prob = exp(alpha+gamma*Z+beta1*X1)/(1+exp(alpha+gamma*Z+beta1*X1));
	p_1  = exp(alpha+gamma+beta1*X1)/(1+exp(alpha+gamma+beta1*X1));
	p_0  = exp(alpha+beta1*X1)/(1+exp(alpha+beta1*X1));
	RD   = p_0-p_1;
    Y=RANBIN(0,1,prob);
    output;
  end;
run;

proc gplot data=sim;
  plot prob*X1;
  by Z;
run;

*--------------------------------------------------------------------*
|                               True NNT                             |
*--------------------------------------------------------------------*;
title4 'Hypothetical CER, TER, True OR, True RD, and True NNT';
title5 'Basis: Population of All Patients';
title6 ' ';

proc means data=sim noprint;
  var p_0 p_1 RD;
  output out=SimM mean=p_1 p_0 RD;
run;

data SimM;  set SimM;
  OR  = 0.50;
  NNT = 1/RD;
proc print data=SimM;
  var OR p_1 p_0 RD NNT;
  format NNT_True 16.8;
run;


%macro nnt(data=,response=,treatment=,covar=);
title1 'CALCULATION OF ADJUSTED NNT WITH 95% CONFIDENCE INTERVALS';
*--------------------------------------------------------------------*
|                       Creation of Data Set                         |
*--------------------------------------------------------------------*;
data RCT;
   set &data;
   response  = &response;
   treatment = &treatment;
   %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=RCT n min max mean std sum maxdec=3 fw=8;
  var RESPONSE TREATMENT COVAR_1-COVAR_&K;
proc means data=RCT n min max mean std sum maxdec=3 fw=8;
  class TREATMENT;
  var RESPONSE COVAR_1-COVAR_&K;
run;

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

*--------------------------------------------------------------------*
|                     Simple Logistic Regression Model             |
*--------------------------------------------------------------------*;
title4 'Estimation of Crude Odds Ratio (Simple Logistic Regression)';
title5 ' ';
title6 ' ';
proc logistic data=RCT descending;
  model RESPONSE = TREATMENT / rl rsq;
run;

*--------------------------------------------------------------------*
|                     Multiple Logistic Regression Model             |
*--------------------------------------------------------------------*;
title4 'Estimation of Adjusted Odds Ratio (Multiple Logistic Regression)';
title5 ' ';
title6 ' ';
ods output ConvergenceStatus=convstat;
proc logistic data=RCT descending outest=est covout;
  model RESPONSE = TREATMENT 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 ALL;   set RCT;
  keep COVAR_1-COVAR_&K;
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 NNTs based upon this model.",,
        "Please build a new model and try again.";
  abort;
end;

use RCT;      read all var{RESPONSE}                     into RESPONSE;
use RCT;      read all var{RESPONSE} where (TREATMENT=0) into RESP_CTR;
use RCT;      read all var{RESPONSE} where (TREATMENT=1) into RESP_TRT;
use RCT;      read all var{TREATMENT}                    into TREATMENT;
use ALL;      read all var _all_                         into ALL;
use beta;     read all var{Intercept}                    into BETA0;
use beta;     read all var{TREATMENT}                    into BETA1;
use beta2_k;  read all var _all_                         into BETA2_k;
use covb;     read all var{Intercept TREATMENT}          into COV1;
use covb2_k;  read all var _all_                         into COV2;

  COVB = (COV1 || COV2);

  N     = NROW(RESPONSE);
  EVENT = SUM(RESPONSE);
  Nt    = NROW(RESP_TRT);
  Nc    = NROW(RESP_CTR);
  Et    = SUM(RESP_TRT);
  Ec    = SUM(RESP_CTR);
  K     = NCOL(BETA2_k);

  CER   = SUM(RESP_CTR)/Nc;
  TER   = SUM(RESP_TRT)/Nt;
  OR    = (TER*(1-CER))/(CER*(1-TER));
  RD    = TER-CER;
  NNT   = 1/RD;

  sn     = 0.05;
  Zsn    = PROBIT(1-SN/2);
  SE     = SQRT((Et*(Nt-Et)/Nt**3)+(Ec*(Nc-Ec)/Nc**3));

  SE_LOR = SQRT((1/Et)+(1/(Nt-Et))+(1/Ec)+(1/(Nc-Ec)));
  OR_L   = exp(log(OR)-Zsn#SE_LOR); 
  OR_U   = exp(log(OR)+Zsn#SE_LOR); 

  print "GENERAL INFORMATION";
  print "Sample size and number of events: "
         N " " EVENT;
  print "Number of unexposed and exposed subjects with events: "
         Nc Ec " " Nt Et ,,,,;

  print "CRUDE (UNADJUSTED) RESULTS";
  if RD>0 then do;
    RD_L  = RD - Zsn#SE;
    RD_U  = RD + Zsn#SE;
    NNT_L = 1/RD_U;
    NNT_U = 1/RD_L;
    print "Crude results refer to a harmful effect of treatment";
    print "Crude event rates: "
           CER[format=6.4] "  " TER[format=6.4],,;
    print "Crude Odds Ratio with SE(log OR) and 95% CI: "
           OR [format=6.4] " " SE_LOR[format=8.6] " " OR_L[format=6.4] OR_U[format=6.4],,;
    print "Crude risk difference (RD) with SE and 95% CI:"
           RD[format=6.4] " " SE[format=8.6] " " RD_L[format=6.4] RD_U[format=6.4],,;
    print "Crude Number needed to treat (NNT) with 95% CI: "
           NNT[format=8.2] " " NNT_L[format=8.2] NNT_U[format=8.2],,,,;
  end;

  if RD<0 then do;
    RD  = -RD;
    NNT = 1/RD;
    RD_L  = RD - Zsn#SE;
    RD_U  = RD + Zsn#SE;
    NNT_L = 1/RD_U;
    NNT_U = 1/RD_L;
    print "Crude results refer to a beneficial effect of treatment";
    print "Crude event rates: "
           CER[format=6.4] "  " TER[format=6.4],,;
    print "Crude Odds Ratio with SE(log OR) and 95% CI: "
           OR [format=6.4] " " SE_LOR[format=8.6] " " OR_L[format=6.4] OR_U[format=6.4],,;
    print "Crude risk difference (RD) with SE and 95% CI:"
           RD[format=6.4] " " SE[format=8.6] " " RD_L[format=6.4] RD_U[format=6.4],,;
    print "Crude Number needed to treat (NNT) with 95% CI: "
           NNT[format=8.2] " " NNT_L[format=8.2] NNT_U[format=8.2],,,,;
  end;

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

  if RD=0 then do;
    NNT="infinity";
    print "Crude event rates, risk difference and NNT:  "
           CER[format=6.4] "  " TER[format=6.4] "  "
           RD [format=6.4] "  " NNT ,,,,;
  end;


  print "ADJUSTED RESULTS FOR NNT (AVERAGE RISK DIFFERENCE APPROACH, BASIS: ALL PATIENTS)";

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

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

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

  RISK1 = EXP(LINPRED1)/(1+EXP(LINPRED1));
  RISK2 = EXP(LINPRED2)/(1+EXP(LINPRED2));
  CER   = SUM(RISK1)/N;
  TER   = SUM(RISK2)/N;
  ARD   = TER-CER;

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

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

  print "Hypothetical control (CER) and treatment event rate (TER): "
         CER[format=6.4] " " TER[format=6.4];

  if ARD>0 then do;
    ARD_L  = ARD - Zsn#SE;
    ARD_U  = ARD + Zsn#SE;
    NNT_L = 1/ARD_U;
    NNT_U = 1/ARD_L;
    print "Adjusted results refer to a harmful effect of treatment";
    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 treat (NNT) with 95% CI: "
          NNT[format=8.2] " " NNT_L[format=8.2] NNT_U[format=8.2];
  end;

  if ARD<0 then do;
    ARD  = -ARD;
    NNT  = 1/ARD;
    ARD_L  = ARD - Zsn#SE;
    ARD_U  = ARD + Zsn#SE;
    NNT_L = 1/ARD_U;
    NNT_U = 1/ARD_L;
    print "Adjusted results refer to a beneficial effect of treatment";
    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 treat (NNT) with 95% CI: "
         NNT[format=8.2] " " NNT_L[format=8.2] NNT_U[format=8.2];
  end;

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

quit;
run;
%mend nnt;

%nnt(data=sim,response=Y,treatment=Z,covar=X1);