/**********************************************************************************************************************; 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: Rate Ratio= Rate in group-2 / Rate in group-1 **********************************************************************************************************************/ %macro IDRcomp(T,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; ****HANDLING OF MISSING/ EXTREME VALUES; IF (&O1=0 AND &O2=0) OR &O1=. OR &PT1=. OR &PT1=0 OR &O2=. OR &PT2=. OR &PT2=0 THEN DO; MID_P=.;RATE_RATIO=.; LL=.;UL=.; END; ELSE 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 DUn = 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)