********************************************************************
*          ESTIMATION OF BENCHMARK VALUES (VARL, VARG)             *
*            AND CALCULATION OF CONFIDENCE INTERVALS               *
*        BASED UPON THE SIMPLE LOGISTIC REGRESSION MODEL           *
*              (BENDER 1999, BIOMETRICAL JOURNAL 41)               *
********************************************************************;
*------------------------------------------------------------------*
|                           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.          !!
!!   THE DATA SET 'DATA' MUST CONTAIN AT LEAST TWO VARIABLES:     !!
!!               - Y = BINARY RESPONSE VARIABLE                   !!
!!               - X = CONTINUOUS COVARIATE                       !!
*------------------------------------------------------------------*;

options linesize=90;

*------------------------------------------------------------------*
|                         Data: Simulation                         |
*------------------------------------------------------------------*;
data DATA;
  do i=0 to 10 by 0.01;
    X=i;
    a=-5;  b=1;
    prob=exp(a+b*x)/(1+exp(a+b*x));
    Y=RANBIN(0,1,prob);
    output;
  end;
  label Y='Response';
  label X='Covariate';
  keep Y X;
run;


*------------------------------------------------------------------*
|                         Logistic Regression                      |
*------------------------------------------------------------------*;
title1 'SIMPLE LOGISTIC REGRESSION';
proc logistic data=DATA descending simple outest=est covout;
  model Y = X  / lackfit risklimits;
  output out=out(keep=X l p u) l=l p=p u=u / alpha=0.05;
  quit;
run;

proc sort data=out out=out;
  by X;
run;

data theta;    set est;    if _type_='PARMS';
  alpha=intercep;   beta=X;
  keep alpha beta;
run;

data sigma;    set est;    if _type_='COV';
  alpha=intercep;   beta=X;
  if _name_='INTERCPT' then par='ALPHA';
  if _name_='X'        then par='BETA ';
  keep par alpha beta;
run;

*------------------------------------------------------------------*
|                    Plot of Logistic Response Curve               |
*------------------------------------------------------------------*;
goptions device=win target=winprtg rotate=landscape
         horigin=0.5cm vorigin=0.5cm hsize=26.5cm vsize=17.5cm
         ftext=swissl htext=2.0;
symbol1 v=none i=join l=1 w=6;
symbol2 v=none i=join l=2 w=6 r=2;
axis1 offset=(0,0) label=(a=90 f=zapf h=2.2 'Probability')
      order=0 to 1 by 0.1  minor=(N=1 w=4) major=(w=4) width=4;
axis2 offset=(0,0) label=(f=zapf h=2.2 'X')
      major=(w=4) width=4;
title f=centx  h=3.0 'Logistic Regression Response Curve';
proc gplot data=out;
  plot (p l u)*X / overlay vaxis=axis1 haxis=axis2;
run;

*------------------------------------------------------------------*
|                 IML: Benchmarks VARL and VARG                    |
*------------------------------------------------------------------*;
title1 'CALCULATION OF BENCHMARK VALUES BASED UPON THE LOGISTIC CURVE';
title3 '(BENDER 1999, Biometrical Journal 41)';
title4 ' ';

proc IML;

print
"*----------------------------------------------------------------------*",
"! Note: 1) The following results are appropriate only if the logistic  !",
"!          model describes the data adequately.  I strongly recommend  !",
"!          to consider at first the goodness-of-fit criteria given by  !",
"!          PROC LOGISTIC (AIC, -2LogL, Hosmer-Lemeshow-Test).          !",
"!       2) It is assumed that the association between the response and !",
"!          the explanatory variable is positive.  If this is not the   !",
"!          case in your data, please reverse the sign of X (X -> -X)   !",
"!          or the definition of the event (Y -> 1-Y).                  !",
"!       3) This concept of benchmark values requires the specification !",
"!          of acceptable limits for p and d.                           !",
"!          These quantities are defined by                             !",
"!          p = Probability (Y=1|X=x)                                   !",
"!          d = Increase of p for a unit increase of X (gradient).      !",
"!       4) For an acceptable risk level p the corresponding X          !",
"!          value of an acceptable risk level (VARL) is calculated.     !",
"!       5) For an acceptable risk gradient d the corresponding X       !",
"!          value of an acceptable risk gradient (VARG) is calculated.  !",
"*----------------------------------------------------------------------*",,;

  use data;   read all var{X}          into X;
  use data;   read all var{Y}          into Y;
  use theta;  read all var{alpha}      into ALP;
  use theta;  read all var{beta}       into BET;
  use sigma;  read all var{alpha beta} into SIG;

  NY = NROW(Y);
  do k=1 to NY;
    if Y[k]^=0 then do;
      if y[k]^=1 then do;
      print "ERROR-MESSAGE:",,
            "Sorry, your response variable Y is not binary!  For Y  ",
            "only the values 0 (no event) and 1 (event) are allowed.",
            "Please check data and model and try again.             ",;
      abort;
      end;
    end;
  end;

  if NROW(SIG)=0 then do;
    print "ERROR-MESSAGE:",,
          "Sorry, no computations are possible because ",
          "the results of PROC LOGISTIC are inadequate!",
          "Please check data and model and try again.  ",;
    abort;
  end;

  if BET<0 then do;
    print "ERROR-MESSAGE:",,
          "The association between Y and X is negative in your data!   ",
          "Please reverse the sign of X or the definition of the event.",
          "Check your data and try again.                              ",;
    abort;
  end;

  N   = NROW(X);      EVENTS = SUM(Y);
  MIN = MIN(X);          MAX = MAX(X);
  SDA = SQRT(SIG[1,1]);   SDB = SQRT(SIG[2,2]);

  P0V  = {0.01, 0.05, 0.10, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99};

  D0M = BET/4;

  IF 0.05<=D0M THEN DO;  D0=0.05;  D0V=D0V//D0; END;
  IF 0.10<=D0M THEN DO;  D0=0.10;  D0V=D0V//D0; END;
  IF 0.20<=D0M THEN DO;  D0=0.20;  D0V=D0V//D0; END;

  IF 0.05>D0M THEN DO;
    DK  = D0M/10-0.00001;
    D0  = D0M/10-0.00001;
    D0V = D0;
    DO K=2 TO 10 by 2;
      D0  = K#DK;
      D0V = D0V//D0;
    END;
  END;

  DEV = SQRT(BET##2 - 4#D0V*BET);

  VARLV = 1/BET * (log(P0V/(1-P0V))-ALP);
  VARGV = 1/BET * (log((BET-(2#D0V)-DEV)/(2#D0V))-ALP);

  SD1V = 1/BET # SQRT(SIG[1,1] + 2#VARLV#SIG[1,2] + (VARLV##2)#SIG[2,2]);
  SD2V = 1/BET # SQRT(SIG[1,1] + 2#(VARGV+1/DEV)#SIG[1,2] + ((VARGV+1/DEV)##2)#SIG[2,2]);

  niv  = 0.05;
  Zniv = PROBIT(1-NIV/2);
  LO1V = VARLV - Zniv#SD1V;     UP1V = VARLV + Zniv#SD1V;
  LO2V = VARGV - Zniv#SD2V;     UP2V = VARGV + Zniv#SD2V;

  print "Sample size, events: "  N EVENTS;
  print "Range of X: "           MIN MAX;
  print "Logistic coefficients (alpha, SD, beta, SD): " ALP SDA BET SDB;
  print "Benchmark VARL (p0, VARL, SD, 95% CI): "  P0V[format=4.2] VARLV SD1V LO1V UP1V;
  print "Maximal Value for d0: "  D0M[format=5.3];
  print "Benchmark VARG (d0, VARG, SD, 95% CI): "  D0V[format=5.3] VARGV SD2V LO2V UP2V;

quit;
run;


*********** All Rights by R. BENDER, Duesseldorf, 1998 *****************;