**********************************************************************************************************************; ****************************compare 2 SIRs by EXACT BINOMIAL TEST and MID-P 95%CI********************************************** February 1, 2014 Used in statistics calculator to compare two SIRs (called by SIRComp.sas) Output parameters: ll and ul Function: compare 2 SIRs (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. WHO. Statistical Methods in Cancer Research - Volume II - The Design and Analysis of Cohort Studies: Comparison of standardized mortality ratios. http://www.iarc.fr/en/publications/pdfs-online/stat/sp82/SP82_vol2-3.pdf 2. Erich L. Lehmann, Joseph P. Romano. Testing Statistical Hypotheses. Wiley, New York, 1959 (Third edition, 2005, Springer Texts in Statistics) 3. John Pezzullo. Interactive Statistical Calculation Pages. www.Statpages.com. 4. Dean AG, Sullivan KM, Soe MM. OpenEpi: Open Source Epidemiologic Statistics for Public Health, Version. www.OpenEpi.com, updated 2013/03/21 Validation: with WHO. Statistical Methods in Cancer Research Data input layout: &O1= observed events in current year &E1= expected events in current year &O2= observed events in past year &E2= expected events in past year Output: TWOTAILP=MID-P. RATIO= SIR in the past year/ SIR in the current year = SIR2/SIR1 = (&O2/&E2)/(&O1/&E1) assumption: reduction of HAI incidence over the years LL & UL =LOWER AND UPPER LIMITS OF 95%CI. **********************************************************************************************************************; %macro binom(O1,E1,O2,E2); %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=0 AND &O2=0) OR &O1=. OR &E1=. OR &O2=. OR &E2=. OR &E1<1 OR &E2<1 THEN DO; midP=.; RATIO=.; LL=.; UL=.; END; ELSE DO; vN=&O2+&O1; vP=&O2/vN; *******************MID-P VALUE FOR HYPOTHESIS TESTING; RATIO1=&O1/&E1; RATIO2=&O2/&E2; O=&O1+&O2; T=&E1+&E2; *take the larger of SIRs (TESTING FOR Ha: RR>1); IF RATIO1>=RATIO2 THEN DO; p1=1- probbnml(&E1/T,O,&O1); p2=0.5* ( probbnml(&E1/T,O,&O1) - probbnml(&E1/T,O,&O1-1) ); P3=P1+P2; END; ELSE DO; p1=1- probbnml(&E2/T,O,&O2); p2=0.5* ( probbnml(&E2/T,O,&O2) - probbnml(&E2/T,O,&O2-1) ); P3=P1+P2; END; midP=2*p3; IF midP>1 THEN midP=1; IF midP<0 THEN midP=0; *******************RATIO AND INTERVAL ESTIMATION (SIR2/SIR1, apriori assumption: SIR2>SIR1); *Conditional maximum likelihood estimate of Ratio; if &O1 ne 0 then do; *<-----------to prevent division by zero when ratio1=zero ; RATIO=(&O2/&E2) / (&O1/&E1); 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)