/**********************************************************************************************************************
	Used in output for calculating proportion CIs (called by HCW_VaccFluSumCMS.sas).  prop_l and prop_u are the output parameters.

										Mid-P exact method


Developed by Minn M. Soe, MBBS,MPH ( Epidemiologist/Biostatistican, DHQP, CDC )

Reference 
1. Dean AG, Sullivan KM, Soe MM. OpenEpi: Open Source Epidemiologic Statistics for Public Health, Version 2.3.1. www.OpenEpi.com
2. Rothman KJ,  Boice JD  Jr:  Epidemiologic analysis with a programmable calculator.   NIH Pub No. 79-1649.  Bethesda, MD:  National Institutes of Health, 1979;31-32.
3. Fleiss JL.  Statistical Methods for Rates and Proportions, 2nd Ed. John Wiley & Sons, New York, 1981.
4. Newcombe RG.  Two-sided confidence intervals for the single proportion: comparison of seven methods.  Statistics in Medicine 1988;17:857-872.
Input values:
x=numerator,  n=denominator.

Output values: 
Calculate p-value (one-sided and two-sided) for testing the null hypothesis of observed propportion similar to a known proportion
Calculate 95%CI of single proportion.

**********************************************************************************************************************/

%macro proportionCI(vX,vN);

%macro BinP(x1,x2);
q=p2/(1-p2);  k=0;  v = 1;  s=0;  tot=0;
    do while(k<=&vN) ;
        tot=tot+v;
        if(k>=&x1 and k<=&x2) then s=s+v ;
        if(tot>10**30) then do;
			s=s/10**30; tot=tot/10**30; v=v/10**30;
		end;
        k=k+1; v=v*q*(&vN+1-k)/k;
     end;
    binP=s/tot;
%mend;
%macro BinP2(x1,x2);
q=p2/(1-p2);  k=0;  v = 1;  s=0;  tot=0;
    do while(k<=&vN) ;
        tot=tot+v;
        if(k>=&x1 and k<=&x2) then s=s+v ;
        if(tot>10**30) then do;
			s=s/10**30; tot=tot/10**30; v=v/10**30;
		end;
        k=k+1; v=v*q*(&vN+1-k)/k;
     end;
    binP2=(s/tot)*0.5;
%mend;


	format prop_L prop_U 8.3;

vP=&vX/&vN;
vX=&vX;
vN=&vN; 

if (&vX=0) then prop_L = .; 
else do;
        p2=vP/2; vsL=0; vsH=vP;  p=2.5/100;
        do while((vsH-vsL)>10**-5);
			%BinP(x1=&vX+1, x2=&vN);
			%BinP2(x1=&vX, x2=&vX);
			 if (BinP+BinP2) >p  then do; vsH=p2; p2=(vsL+p2)/2 ;end;
			else do; vsL=p2; p2=(vsH+p2)/2 ;end;
		end;
        prop_L = p2;
end;       
   

if (&vX=&vN) then prop_U = 1; 
else do;
        p2=(1+vP)/2; vsL=vP; vsH=1; p=2.5/100;
        do while((vsH-vsL)>10**-5);
			%BinP(x1=0, x2=&vX-1);
			%BinP2(x1=&vX, x2=&vX);
			if (BinP+BinP2) <p then do; vsH=p2; p2=(vsL+p2)/2;end;
			else do; vsL=p2; p2=(vsH+p2)/2 ;end;
		end;
        prop_U = p2;
end; 

/*handing missing values/ illogical values*/
if &vX=. or &vN=. then do;
vp=.;prop_L=.;prop_U=.;
end;
 
*run;
	drop vx vn vp p2 vsl vsh p q k v s tot binp binp2; 

%mend;

****************IMPORT AN EXTERNAL DATASET AND RUN THE MACRO*********************;
/*SAS code Example if you download and save macro on your computer*/
options mprint;

data ClabsiExample;
input Numerator Denominator;
cards;
3 50
;
run;

data ClabsiExample_;set ClabsiExample;
%proportionCI(vX=Numerator,vN=Denominator);
run;

proc print;run;