*************************************************************************
*            CALCULATION of CONFIDENCE INTERVALS FOR NNT                *
*                  VIA CONFIDENCE INTERVALS FOR ARR                     *
* METHOD BASED ON WILSON SCORE INTERVALS WITHOUT CONTINUITY CORRECTION  *
*    REFERENCES: NEWCOMBE (1998), STATISTICS IN MEDICINE 17, 873-890    *
*           BENDER (2001), CONTROL. CLIN. TRIALS 22, 102-110            *
*************************************************************************;

  *------------------------------------------------------------------*
  |                           QUESTIONS ?                            |
  |                                                                  |
  | --->  Dr. Ralf BENDER  Phone:  +49 221 35685-451                 |
  |                        Fax:    +49 221 35685-891                 |
  |                        E-Mail: Ralf@rbsd.de                      |
  *------------------------------------------------------------------*;
  *------------------------------------------------------------------*
  !!                         NOTE:                                  !!
  !!   You have to enter the actual results of the two groups       !!
  !!   in lines 30 to 33.                                           !!
  !!   e1 is the number of adverse events in the high risk group    !!
  !!      (e.g. standard treatment)                                 !!
  !!   n1 is the sample size of the high risk group                 !!
  !!   e2 is the number of adverse events in the low risk group     !!
  !!      (e.g. new treatment)                                      !!
  !!   n2 is the sample size of the low risk group                  !!
  *------------------------------------------------------------------*;

options linesize=95;
proc iml;
    e1 =  10;             /* Number of events in high-risk group A */
    n1 = 200;             /* Sample size of high risk group A      */
    e2 =   1;             /* Number of events in low-risk group B  */
    n2 = 100;             /* Sample size of low risk group B       */

*--------------------------*First Checks*------------------------------*;
  if e1>n1 then do;
    print "ERROR !!!" ,,,,,
          "The number of events (e1) in the first group is larger",
          "than the sample size (n1) !!",,
          "Please re-enter the data and try again.";
    abort;
  end;

  if e2>n2 then do;
    print "ERROR !!!" ,,,,,
          "The number of events (e2) in the second group is larger",
          "than the sample size (n2) !!",,
          "Please re-enter the data and try again.";
    abort;
  end;

  if e1<0 | e2<0 then do;
    print "ERROR !!!" ,,,,,
          "The number of events in the first (e1) or second group (e2)",
          "is negative !!",,
          "Please re-enter the data and try again.";
    abort;
  end;

  if n1<=0 | n2<=0 then do;
    print "ERROR !!!" ,,,,,
          "The sample size in the first (n1) or second group (n2)",
          "is not positive !!",,
          "Please re-enter the data and try again.";
    abort;
  end;

  n = n1+n2;
  if n<10 then do;
    print "SORRY !!!" ,,,,,
          "The total sample size is extremely low (n<10)!!",
          "Please use another software providing exact methods.",,
          "Sorry for any inconvenience.";
    abort;
  end;

  print "CONFIDENCE INTERVALS FOR THE NUMBER NEEDED TO TREAT (NNT)" ,,
        "Based on Wilson score intervals without continuity correction" ,,
        "References: NEWCOMBE (1998): Stat. Med. 17, 873-890" ,
        "BENDER (2001): Control. Clin. Trials 22, 102-110" ,,, ;

*--------------------------*Descriptive*-------------------------------*;
  start Descript;
    r1 = e1/n1;
    r2 = e2/n2;
    ARR = r1-r2;
    if ARR=0  then NNT="infinity";
    else  NNT=1/ARR;
    z90 = probit(1-(0.10/2));
    z95 = probit(1-(0.05/2));
    z99 = probit(1-(0.01/2));
    print "Events and sample sizes of the groups: " e1 n1 e2 n2 ,
          "Risks of the groups and ARR:           " r1 r2 ARR ,
          "Number needed to treat (NNT):          " NNT ,,;
  finish;

*------------------------------*Wald Method*----------------------------*;
  start Wald;
    print "Standard Wald method  (frequently inadequate)";
    al90 = ARR - z90*SQRT((e1*(n1-e1)/n1**3)+(e2*(n2-e2)/n2**3));
    au90 = ARR + z90*SQRT((e1*(n1-e1)/n1**3)+(e2*(n2-e2)/n2**3));
    al95 = ARR - z95*SQRT((e1*(n1-e1)/n1**3)+(e2*(n2-e2)/n2**3));
    au95 = ARR + z95*SQRT((e1*(n1-e1)/n1**3)+(e2*(n2-e2)/n2**3));
    al99 = ARR - z99*SQRT((e1*(n1-e1)/n1**3)+(e2*(n2-e2)/n2**3));
    au99 = ARR + z99*SQRT((e1*(n1-e1)/n1**3)+(e2*(n2-e2)/n2**3));
    nl90 = 1/au90;    nu90 = 1/al90;
    nl95 = 1/au95;    nu95 = 1/al95;
    nl99 = 1/au99;    nu99 = 1/al99;
    print "90% confidence intervals for ARR and NNT: "  al90 au90 "  " nl90 nu90;
    print "95% confidence intervals for ARR and NNT: "  al95 au95 "  " nl95 nu95;
    print "99% confidence intervals for ARR and NNT: "  al99 au99 "  " nl99 nu99 ,,,;
  finish;

*-----------------------*Wilson Score Method*---------------------------*;
  start Wilson;
    print "Confidence intervals based upon Wilson score intervals";
    l190 = (2*e1+z90**2) / (2*(n1+z90**2)) - SQRT(((2*e1+z90**2)/(2*(n1+z90**2)))**2-e1**2/(n1**2+n1*z90**2));
    u190 = (2*e1+z90**2) / (2*(n1+z90**2)) + SQRT(((2*e1+z90**2)/(2*(n1+z90**2)))**2-e1**2/(n1**2+n1*z90**2));
    l290 = (2*e2+z90**2) / (2*(n2+z90**2)) - SQRT(((2*e2+z90**2)/(2*(n2+z90**2)))**2-e2**2/(n2**2+n2*z90**2));
    u290 = (2*e2+z90**2) / (2*(n2+z90**2)) + SQRT(((2*e2+z90**2)/(2*(n2+z90**2)))**2-e2**2/(n2**2+n2*z90**2));

    l195 = (2*e1+z95**2) / (2*(n1+z95**2)) - SQRT(((2*e1+z95**2)/(2*(n1+z95**2)))**2-e1**2/(n1**2+n1*z95**2));
    u195 = (2*e1+z95**2) / (2*(n1+z95**2)) + SQRT(((2*e1+z95**2)/(2*(n1+z95**2)))**2-e1**2/(n1**2+n1*z95**2));
    l295 = (2*e2+z95**2) / (2*(n2+z95**2)) - SQRT(((2*e2+z95**2)/(2*(n2+z95**2)))**2-e2**2/(n2**2+n2*z95**2));
    u295 = (2*e2+z95**2) / (2*(n2+z95**2)) + SQRT(((2*e2+z95**2)/(2*(n2+z95**2)))**2-e2**2/(n2**2+n2*z95**2));

    l199 = (2*e1+z99**2) / (2*(n1+z99**2)) - SQRT(((2*e1+z99**2)/(2*(n1+z99**2)))**2-e1**2/(n1**2+n1*z99**2));
    u199 = (2*e1+z99**2) / (2*(n1+z99**2)) + SQRT(((2*e1+z99**2)/(2*(n1+z99**2)))**2-e1**2/(n1**2+n1*z99**2));
    l299 = (2*e2+z99**2) / (2*(n2+z99**2)) - SQRT(((2*e2+z99**2)/(2*(n2+z99**2)))**2-e2**2/(n2**2+n2*z99**2));
    u299 = (2*e2+z99**2) / (2*(n2+z99**2)) + SQRT(((2*e2+z99**2)/(2*(n2+z99**2)))**2-e2**2/(n2**2+n2*z99**2));

    delta90=SQRT((e1/n1-l190)**2+(u290-e2/n2)**2); epsil90=SQRT((u190-e1/n1)**2+(e2/n2-l290)**2);
    delta95=SQRT((e1/n1-l195)**2+(u295-e2/n2)**2); epsil95=SQRT((u195-e1/n1)**2+(e2/n2-l295)**2);
    delta99=SQRT((e1/n1-l199)**2+(u299-e2/n2)**2); epsil99=SQRT((u199-e1/n1)**2+(e2/n2-l299)**2);

    al90 = ARR - delta90;  au90 = ARR + epsil90;
    al95 = ARR - delta95;  au95 = ARR + epsil95;
    al99 = ARR - delta99;  au99 = ARR + epsil99;

    nl90 = 1/au90;    nu90 = 1/al90;
    nl95 = 1/au95;    nu95 = 1/al95;
    nl99 = 1/au99;    nu99 = 1/al99;

    print "90% confidence intervals for ARR and NNT: "  al90 au90 "  " nl90 nu90;
    print "95% confidence intervals for ARR and NNT: "  al95 au95 "  " nl95 nu95;
    print "99% confidence intervals for ARR and NNT: "  al99 au99 "  " nl99 nu99 ,,;

  finish;

  run Descript;
  run Wald;
  run Wilson;
quit;