/* ********************************************************************************************************************** Used in output for calculating rate CIs (called by StandardRateTable.sas). rate_l and rate_u are the output parameters. Mid-P exact CI for Incidence Density Rate Developed by Minn M. Soe, MPH, MBBS, 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. Rothman KJ, Greenland S. Modern Epidemiology, 2nd Edition. Lippincott-Raven Publishers, Philadelphia, 1998. **********************************************************************************************************************; */ *-------------------------------------mid-P exact CI for Incidence Density Rate---------------------------------------; %macro fish(z=,x1=,x2=); q=1; tot=0; s=0; k=0; do while(k<&z or q>(tot*10**-10)); tot=tot+q; if(k>=&x1 and k<=&x2) then s=s+q ; if(tot>10**30)then do; s=s/10**30; tot=tot/10**30; q=q/10**30; end; k=k+1; q=q*&Z/k; end; poisP=s/tot; %mend; %macro fish2(z=,x1=,x2=); q=1; tot=0; s=0; k=0; do while(k<&z or q>(tot*10**-10)); tot=tot+q; if(k>=&x1 and k<=&x2) then s=s+q ; if(tot>10**30)then do; s=s/10**30; tot=tot/10**30; q=q/10**30; end; k=k+1; q=q*&Z/k; end; poisP2=0.5*(s/tot); %mend; %macro rateCIComp(numer,denom); format rate_l rate_u 8.3 ; rate=(&numer/&denom)*1000; *Lower tail; if &numer=0 then rate_l=.; else do; v=0.5; dv=0.5; p=2.5/100; do while(dv>10**-5); dv=dv/2; %fish(z=(1+&numer)*v/(1-v), x1=(&numer+1), x2=10**10); %fish2(z=(1+&numer)*v/(1-v), x1=&numer, x2=&numer); if( (PoisP+PoisP2)>p) then v=v-dv; else v=v+dv; end; if &denom > 0 then rate_l=round((((1+&numer)*v/(1-v))/&denom)*1000,.001); end; *Upper tail; v=0.5; dv=0.5; p=2.5/100; do while(dv>10**-5); dv=dv/2; %fish(z=(1+&numer)*v/(1-v), x1=0, x2=(&numer-1) ); %fish2(z=(1+&numer)*v/(1-v), x1=&numer, x2=&numer); if( (PoisP+PoisP2) 0 then rate_u=round((((1+&numer)*v/(1-v))/&denom)*1000,.001); drop v dv p q tot k s poisP poisP2; %mend; ****************IMPORT AN EXTERNAL DATASET AND RUN THE MACRO*********************; options mprint; data ClabsiExample; input Numerator Denominator; cards; 12 35 ; run; data ClabsiExample_;set ClabsiExample; %rateCIcomp(numer=Numerator, denom=Denominator); run; proc print;run;