/**********************************************************************************************************************;
February 1, 2014
Used in statistics calculator to compare two IDRs (called by IDRcomp.sas).  IDR_pval is the output parameter.
	Used in output for DU ratios (called by StandardRateTable.sas).  IDR_pval, IDR_pval1, or IDR_pval2 is the output parameter.
*********************COMPARE 2 INCIDENCE RATES by EXACT BINOMIAL TEST and MID-P 95%CI***********************************

Function: compare 2 incidence density rates (based on the fact that the distribution of two poisson variates conditional on their sum is binomial)

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

Reference 
1. Erich L. Lehmann, Joseph P. Romano. Testing Statistical Hypotheses. Wiley, New York, 1959(Third edition, 2005, Springer Texts in Statistics) 
2. John Pezzullo. Interactive Statistical Calculation Pages. www.Statpages.com.
3. Dean AG, Sullivan KM, Soe MM. OpenEpi: Open Source Epidemiologic Statistics for Public Health, Version. www.OpenEpi.com, updated 2013/03/21
4. Austin, Harland (Emory). Epidemiology method-II: Statistical issues for Density type follow-up studies.

Validation: WinPEPI, OpenEpi, Propective study of dietary fat and risk of prostate cancer. J national cancer institute. 1993, 85:1571-79.

Data input layout:
&O1= observed events in group-1
&PT1=person-time in group-1
&O2= observed events in group-2
&PT2= person-time in group-2

OUTPUT:

/*Original macro by Minn. The web-version has been modified by Kelly*/

%macro TWORATES(O1=,PT1=,O2=,PT2=);

*STAT FUNCTION FOR INTERVAL CALCULATION;
%macro BinP(x1=,x2=, PP=, N=);
q=&PP/(1-&PP);  k=0;  v = 1;  s=0;  tot=0;
    do while(k<=&N) ;
        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*(&N+1-k)/k;
     end;
    binP=s/tot;
%mend;

%macro BinP2(x1=,x2=, PP=, N=);
q=&PP/(1-&PP);  k=0;  v = 1;  s=0;  tot=0;
    do while(k<=&N) ;
        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*(&N+1-k)/k;
     end;
     binP2=(s/tot)*0.5;
%mend;


*input parameters;
IF &O1 ne . and &PT1 ne . and &O2 ne . and &PT2 ne . and &PT1>0 and &PT2>0 THEN DO;
vN=&O2+&O1; 
vP=&O2/vN; 


****************MID-P FOR HYPOTHESIS TESTING;
RATIO1=&O1/&PT1; RATIO2=&O2/&PT2;
O=&O1+&O2;
T=&PT1+&PT2;

*take the larger of RISKs (TESTING FOR Ha: RR>1);
IF RATIO1>=RATIO2 THEN DO;
	p1=1- probbnml(&PT1/T,O,&O1);
	p2=0.5* ( probbnml(&PT1/T,O,&O1) - probbnml(&PT1/T,O,&O1-1) );
	P3=P1+P2;
END;
ELSE DO;
	p1=1- probbnml(&PT2/T,O,&O2);
	p2=0.5* ( probbnml(&PT2/T,O,&O2) - probbnml(&PT2/T,O,&O2-1) );
	P3=P1+P2;
END;

MID_P=2*p3;
IF MID_P>1 THEN MID_P=1;
IF MID_P<0 THEN MID_P=0;


*****************RATE_RATIO AND INTERVAL ESTIMATION;
*Conditional maximum likelihood estimate of Ratio= rate2/rate1;
IF &O1 NE 0 THEN DO; *PREVENT DIVISION BY ZERO WHEN RATE RATIO1=0;

RATE_RATIO=(&O2/&PT2) / (&O1/&PT1);


if (&O2=0) then DL = 0; 
else do;  p2=vP/2; vsL=0; vsH=vP;  p=2.5/100;
        do while((vsH-vsL)>10**-5);
			%BinP(x1=&O2+1, x2=vN, PP=P2, N=VN);
			%BinP2(x1=&O2,  x2=&O2, PP=P2, N=VN);
			 if (BinP+BinP2) >p  then do;vsH=p2; p2=(vsL+p2)/2 ;end;
			 else do; 			 vsL=p2; p2=(vsH+p2)/2 ;end;
		end;
        DL = P2;
end;     

if (&O2=vN) then DU = 1; 
else do;
        p3=(1+vP)/2; vsL=vP; vsH=1; p=2.5/100;
        do while((vsH-vsL)>10**-5);
			%BinP (x1=0, x2=&O2-1, PP=P3, N=VN);
			%BinP2(x1=&O2,  x2=&O2, PP=P3, N=VN);
			if (BinP+BinP2)  <p  then do; vsH=p3; p3=(vsL+p3)/2;end;
			else do; vsL=p3; p3=(vsH+p3)/2 ;end;
		end;
        DU = P3;
end;       


ll=(dl*&PT1)/((1-dl)*&PT2);
Ul=(DU*&PT1)/((1-DU)*&PT2);

****HANDLING OF MISSING/ EXTREME VALUES;
IF &O1=. OR &PT1=. OR &PT1=0 OR &O2=. OR &PT2=. OR &PT2=0 THEN DO;
MID_P=.;RATE_RATIO=.; LL=.;UL=.;
END;


***BECAUSE RATE_RATIO=RISK2/ RISK1....
WHEN RISK GROUP2=0 then rate-ratio=0;
IF RATE_RATIO=0 THEN LL=.;

end;
END;
drop O  P3  RATIO1 RATIO2 T binP binP2 k p  p1 p2 q s tot v vN vP vsH vsL DL DU;

%mend;
quit;

********************IMPORTING AN EXTERNAL DATASET AND RUN THE MACRO*********************;

options mprint;

data CLABSIexample;/*Create a data set*/
input Observed1 PersonTime1 Observed2 PersonTime2;
cards;
4 39 5 70
;
run;

data CLABSIexample_;set CLABSIexample;/*This step calls the macro*/
%TWORATES(O1=Observed1,PT1=PersonTime1,O2=Observed2,PT2=PersonTime2);
run;

proc print;run;