Date:         Fri, 15 Sep 2000 17:44:37 +0100
Reply-To:     John Whittington <John.W@MEDISCIENCE.CO.UK>
Sender:       "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From:         John Whittington <John.W@MEDISCIENCE.CO.UK>
Subject:      Re: Why offset in confidence interval for proportions?
In-Reply-To:  <E13ZwHh-00060L-00@relay1.netnames.net>
Content-Type: text/plain; charset="us-ascii"

A number of people have asked me to post the macro I use for caluculating
binary confidence intervals.  I include two versions the below - the first
when one supplies parameters of 'number of successes' (nsuc), 'total sample
size' (ntot) and 'confidence level required' (conf) - and the second when
the proportion of successes (p) is supplied instead of the number of
successes.  Both of these macros write the results to the log, but could
obviously be easily modified to return the results programatically.  I've
been using these for ages, and honestly can't remember whether I wrote them
from scratch myself, or was inspired by (or lifted the code from!) someone
else - so I apologise for anyone who deserves credit, but is not being
given it by me!  The algorithm is quite subtle, in as much as it relies
upon the relationship between the beta distribution and the f distribution
- but it does seem to give correct answers!

%macro bin_CI (nsuc, ntot, conf);
   data null;
      p = &nsuc / &ntot ;
      clevel = &conf + ((1 - &conf) / 2) ;
      f_low = finv(clevel, 2 * (&ntot - &nsuc + 1), 2 * &nsuc) ;
      f_up  = finv(clevel, 2 * (&nsuc + 1), 2 * (&ntot - &nsuc)) ;
      low_CL = &nsuc / (&nsuc + (&ntot - &nsuc + 1) * f_low) ;
      up_CL  = (&nsuc + 1) *f_up / (&ntot - &nsuc + (&nsuc + 1) * f_up) ;
      put ; put "Two-sided Confidence Intervals" ;
      put "Number of successes= &nsuc    Sample size= &ntot" ;
      put p= low_CL= up_CL= ; put ;
   run;
%mend;

*** test example ;
%bin_CI (50, 100, 0.95)

%macro bin_CIp (p, ntot, conf);
   data null;
      nsuc = &p * &ntot ;
      clevel = &conf + ((1 - &conf) / 2) ;
      f_low = finv(clevel, 2 * (&ntot - nsuc + 1), 2 * nsuc) ;
      f_up  = finv(clevel, 2 * (nsuc + 1), 2 * (&ntot - nsuc)) ;
      low_CL = nsuc / (nsuc + (&ntot - nsuc + 1) * f_low) ;
      up_CL  = (nsuc + 1) *f_up / (&ntot - nsuc + (nsuc + 1) * f_up) ;
      put ; put "Two-sided Confidence Intervals" ;
      put "Number of successes= " @ ; put nsuc @ ; put "Sample size= &ntot" ;
      put "p= &p " low_CL= up_CL= ; put ;
   run;
%mend;

* test example ;
%bin_CIp (0.5, 100, 0.95)

Kind Regards,


John

----------------------------------------------------------------
Dr John Whittington,       Voice:    +44 (0) 1296 730225
Mediscience Services       Fax:      +44 (0) 1296 738893
Twyford Manor, Twyford,    E-mail:   John.W@mediscience.co.uk
Buckingham  MK18 4EL, UK             mediscience@compuserve.com
----------------------------------------------------------------