/*******************************************************************/ /* */ /* NAME: BINOMCI */ /* TITLE: Confidence Interval Estimation for Biomial Proportions */ /* PRODUCT: BASE */ /* KEYS: */ /* PROCS: MACRO */ /* VERSION: 8.1 */ /* */ /*******************************************************************/ /*----------------------------------------------------------- TITLE ----- BINOMCI: A SAS macro for confidence interval estimation for binomial proportions. BACKGROUND ---------- It is known for long time that the standard Wald confidence interval has bad coverage probability even for moderate and large sample sizes. A solution (since version 8.1 also implemented in PROC FREQ, use the BINOMIAL-Option, vgl. http://www.sas.com/service/techsup/faq/stat_proc/freqproc922.html) uses the exact confidence but these show unnecessarily conservative behaviour. Brown, Cai and DasGupta show the erratic behaviour of the different standrad methods and recommend three solutions (Wilson's, Agresti-Coull's and Jeffrey's interval for this situation. These procedures have been implemented in this macro along with three standard methods (Wald CI, Wald Ci with continuity correction and the exact interval). The macro is based on an old macro (exactpci.sas) to compute exact confidence intervals delivered by SAS Institute DESCRIPTION ----------- The BINOMCI-Macro is a simple SAS/BASE-macro that uses only data step programming and a PROC PRINT-Statement for printing of results. DETAILS ------- For observed proportions exactly equal to zero or to one only the exact and the Jeffrey's interval are computed, the other CI are not defined in this case. In every case all CIs (except Exact and Jeffreys) are truncated at zero or at one. Formulas for the intervals can be found in the references below. REFERENCES ---------- Snedecor, G.W. and Cochran, W.G. (1980), Statistical Methods, Iowa State University Press, pg 121. Collett, D. (1991), Modelling Binary Data, London: Chapman & Hall, pp 23-25. Brown, L.D., Cai, T.T., and DasGupta, A. (2001): Interval Estimation for a Binomial Proportion, Statistical Science, 16, 101-133. Brown, L.D., Cai, T.T., and DasGupta, A. (2002): Confidence Intervals for a Binomial Proportion and Asymptotic Expansions, Annals of Statistics, 30, 688-707. SYNTAX ------ %binomci (data=_last_, events=, trials=, level=, printdata=1 ) where data= specifies the data set you are using. events= specifies the number of "successes". This is the numerator of your observed proportion. trials= Specifies the number of binomial trials. This is the denominator of your observed proportion. level= Specifies the desired confidence level of the intervals and is a decimal value between 0 and 1. printdata= Controls printing of output. Printdata=0 suppresses all printing. Printdata=1 is the default. AUTHOR ------ Comments, suggestions and error messages are very much appreciated. Send them to: Oliver Kuss University of Halle-Wittenberg Institute of Medical Epidemiology, Biometry, and Informatics Magdeburger Str. 27 06097 Halle/Saale, Germany Tel.: +49-345-5573582 Fax: +49-345-5573580 e-mail: Oliver.Kuss@medizin.uni-halle.de -------------------------------------------------------------*/ %macro binomci(data=,events=,trials=,level=,printdata=1); %if %length(&data)=0 %then %do; %put '*****************ERROR: NUMBER OF EVENTS NOT SPECIFIED********************';%goto exit; %end; %if %length(&events)=0 %then %do; %put '*****************ERROR: NUMBER OF EVENTS NOT SPECIFIED********************';%goto exit; %end; %if %length(&trials)=0 %then %do; %put '*****************ERROR: NUMBER OF EVENTS NOT SPECIFIED********************';%goto exit; %end; data temp1; set &data; conf=&level; events=&events; trials=&trials; call symput('conf',trim(left(put(conf*100,2.)))); phat=events/trials; if phat in (0,1) then do; put "NOTE: Only exact limits are computed when observed p=0 or 1."; goto exact; end; conf=(1-conf)/2 + conf; crit=probit(conf); halfwid=crit*sqrt(phat*(1-phat)/trials); adjhwid=crit*sqrt(phat*(1-phat)/trials)+1/(2*trials); /* Large sample confidence limits */ lslcl=max(phat-halfwid,0); lsucl=min(phat+halfwid,1); /* Large sample limits with continuity correction */ lsadjlcl=max(phat-adjhwid,0); lsadjucl=min(phat+adjhwid,1); /* Wilson Interval */ wilsonlcl=max((events + crit**2*0.5)/(trials + crit**2) - (crit*sqrt(trials))/(trials + crit**2)* sqrt(phat*(1-phat)+ crit**2/(4*trials)),0); wilsonucl=min((events + crit**2*0.5)/(trials + crit**2) + (crit*sqrt(trials))/(trials + crit**2)* sqrt(phat*(1-phat)+ crit**2/(4*trials)),1); /* The Agresti-Coull Interval */ agcevents=events + crit**2*0.5; agctrials=trials + crit**2; agcp=agcevents/agctrials; agrcoulllcl=max(agcp - crit*sqrt(agcp*(1-agcp)/agctrials),0); agrcoullucl=min(agcp + crit*sqrt(agcp*(1-agcp)/agctrials),1); exact: /* Exact und Jeffrey's confidence limits */ if phat>0 then do; exlcl=events / (events+(trials-events+1)*finv( conf,2*(trials-events+1),2*events ) ); jefflcl=betainv(1-conf,events+0.5,trials-events+0.5); testlcl=betainv(1-conf,events,trials-events+1); end; else do; exlcl=0; jefflcl=0; testlcl=0; end; if phat<1 then do; exucl =(events+1) /(events+1+(trials-events)/finv( conf,2*(events+1),2*(trials-events) ) ); jeffucl=betainv(conf,events+0.5,trials-events+0.5); testucl=betainv(conf,events+1,trials-events); end; else do; exucl=1; jeffucl=1; testucl=0; end; run; %if &printdata=1 %then %do; proc print noobs split='/'; var phat exlcl exucl testlcl testucl jefflcl jeffucl wilsonlcl wilsonucl agrcoulllcl agrcoullucl lsadjlcl lsadjucl lslcl lsucl; label phat = "Observed/proportion" exlcl = "Exact/lower/limit" exucl = "Exact/upper/limit" wilsonlcl = "Wilson/lower/limit" wilsonucl = "Wilson/upper/limit" agrcoulllcl = "Agresti-Coull/lower/limit" agrcoullucl = "Agresti-Coull/upper/limit" jefflcl = "Jeffreys/lower/limit" jeffucl = "Jeffreys/upper/limit" lsadjlcl = "Adjusted/approx./lower/limit" lsadjucl = "Adjusted/approx./upper/limit" lslcl = "Approx./lower/limit" lsucl = "Approx./upper/limit"; title "&conf% Confidence limits for Binomial probability"; title2 "Adjusted limits use 0.5 continuity correction"; run; %end; %exit: %mend binomci; data test; input COUNT SUMINZ; cards; 2 50 4 100 10 250 20 500 40 1000 60 1500 80 2000 100 2500 120 3000 140 3500 160 4000 ;run; %binomci(data=test,events=count,trials=suminz,level=0.95);