**************************************************************************
**************************************************************************
THE FOLLOWING MACRO CAN BE USED TO COMPARE AN SIR TO 1 AND IS CONSISTENT
WITH THE METHODS USED BY NHSN FOR SUCH COMPARISONS.
(JUNE 2011)
**************************************************************************
**************************************************************************;

%macro SIRcomp(nobs,nexp);
    format sir_l sir_u 8.3 SIR_pval 6.4;
	SIR = round(&nobs/&nexp,.001) ;
	if 0 <= &nobs < 15 then do; 
		alpha = .95 ;                                                                                                                
		a = .025;
		z =    1.96;                                                                                                                              
		if &nobs > 0 then SIR_L = round(.5*cinv(  a,(2*&nobs+0))  /&nexp,.001);   /* a(alpha) = .025 here  */                                                                       
		SIR_U = round(.5*cinv(1-a,(2*&nobs+2))/&nexp,.001);   /* and .975 here  */
		if &nobs > &nexp then SIR_pval = round(cdf('chisq',2*&nexp,2*&nobs),.0001);
			else if &nobs < &nexp then SIR_pval = round(sdf('chisq',2*&nexp,2*&nobs+2),.0001);
		if &nobs = 0 then sir_l = .;
	end;
	else if &nobs >=15 then do;                                                                                                             
		alpha=.95;                                                                                                                            
		z=1.96;                                                                                                                               
		SIR_L = round( &nobs    * (1 - ((9* &nobs+0)   **-1) - (z*((9* &nobs)   **(-1/2))))**3 / &nexp,.001) ; 
		SIR_U = round((&nobs+1) * (1 - ((9*(&nobs+1))**-1) + (z*((9*(&nobs+1))**(-1/2))))**3 / &nexp,.001) ;
		if &nobs > &nexp then SIR_pval = round(1 - cdf('normal',(3*(&nobs**(1/2))) * (1 - (9*&nobs)**-1 - (&nobs/&nexp)**-(1/3)),0,1),.01);
			else if &nobs < &nexp then SIR_pval = round(1 - cdf('normal', (3*((&nobs+1)**(1/2))) * ((((&nobs+1)/&nexp)**(-1/3)) + ((9*(&nobs+1))**-1) - 1),0,1),.01);
	end;
	if &nexp < 1 then sir = .;
    drop a alpha z ;
%mend SIRcomp;