/**********************************************************************************************************************

							Two-sided P-value for SIR compared to 1 or other nominal value

***********************************************************************************************************************

2 SIDED P-VALUE 

Exact p-value with mid-p method based on the Poisson distribution will be computed only when observed numbers of events are <=100 (point estimate of mu <=100).
Otherwise, Byar approximation method (based on the large sample method) is used for computational efficiency, 
assuming similar results with mid-P method. The one-tailed P-values thus obtained are doubled to provide two-tailed P-values 
in the final results.


REFERENCE: 
1. JH ABRAMSON. WINPEPI PROGRAMS. DESCRIBE MANUAL (VERSION 2.42), PAGE-52. Available at 'http://www.brixtonhealth.com/pepi4windows.html'
2. Dean AG, Sullivan KM, Soe MM. OpenEpi: Open Source Epidemiologic Statistics for Public Health, Version. www.OpenEpi.com.
3. 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.
4. Rothman KJ, Greenland S.  Modern Epidemiology, 2nd Edition.  Lippincott-Raven Publishers, Philadelphia, 1998.
5. GEOFFREY RC, SHU-YING Y. MID-P CONFIDENCE INTERVALS FOR THE POISSON EXPECTATION. STATISTICS IN MEDICINE, 13,2189-2203 (1994)

How to use the macro?
A facility’s SIR can be compared to a single nominal value, such as a “goal” or “target” SIR, using this macro. 
To do this, you will need to re-calibrate the number of predicted events for each hospital based on the target SIR. 

For example, facility A has '2' observed events and '2.04' predicted events, resulting in the SIR of 0.98 for a given time period,
and we would like to compare this SIR value (0.98) to the target SIR of 0.80. 
* Fist, we must multiply the number of predicted events (2.04) by the target SIR (0.80), and obtain the newly calculated
number of predicted events (2.04 x 0.80 = 1.63). 
* In the macro parameter values, the number of observed events is '2' for the macro variable 'OBS' and number of predicted events
is '1.63' in 'EXP'. 
* Run the macro to obtain the p-value. If the two-sided p-value is <= 0.05, we can conclude that facility A’s SIR of 0.98 
is significantly different from the target SIR of 0.80. 

**********************************************************************************************************************/;


options mprint mlogic symbolgen;

%macro SIR(OBS,EXP);

e=2.718281828459045235; 

**********missing data handling and when &EXP<1;
IF  &OBS = . or &EXP = . or &EXP<1 THEN DO; 
midp=.;
end;
else do;

****when observed numbers of events 'mu' are <=100;
if (&OBS<=100)  THEN DO;

if (&OBS<&EXP) then do;
        k=&OBS-1; 
        do while (k>=0); 
    	numerator=(e**-&EXP)* (&EXP**k);
	
	    kk=k;
	    denom=1;
	    do while (kk>0);
			denom=denom*kk; kk=kk-1; *calculate the value of a factorial;
	    end;
	    subtotal=numerator/denom;	   
	    total=sum(total,subtotal);
    	k=k-1;
		end;

*first part of equation;
	num=0;  deno=0;  i=0;  aa=0;
	num=(e**-&EXP)* (&EXP**&OBS);
	deno=1;
	i=&OBS;
	do while (i>0); deno=deno*i; i=i-1;end; 
	
	aa=(num/deno)*0.5;
	
	*combine both parts of equation;
	MidP=2*(sum(aa,total));

end;


if (&OBS>=&EXP) then do;
        k=&OBS-1; 
        do while (k>=0); 
    	numerator=(e**-&EXP)* (&EXP**k);
	
	    kk=k;
	    denom=1;
	    do while (kk>0);
			denom=denom*kk; kk=kk-1; *calculate the value of a factorial;
	    end;
	    subtotal=numerator/denom;	   
	    total=sum(total,subtotal);
    	k=k-1;
		end;

*first part of equation;
	num=0;  deno=0;  i=0;  aa=0;
	num=(e**-&EXP)* (&EXP**&OBS);
	deno=1;
	i=&OBS;
	do while (i>0); deno=deno*i; i=i-1;end; 
	
	aa=(num/deno)*0.5;
	
	*combine both parts of equation;
	MidP=2*(1-(sum(aa+total)));

end;
if midp>1 then MidP=1;
if midp<0 then MidP=0;
END;


****when observed numbers of events are >100, use Byar poisson approximation (Rothman KJ, Boice JD);
ELSE DO;

Byar_p=0; z=0; OBS_=0;
OBS_=&OBS;
if (&OBS<=&EXP) then OBS_=&OBS+1;
z=  sqrt(9*OBS_) * ( 1 - 1/(9*OBS_) - (&EXP/OBS_)**(1/3) );

if &OBS<=&EXP then Byar_p=2*probnorm(z);
else Byar_p=2*(1-probnorm(z));


if Byar_p>1 then Byar_p=1;
if Byar_p<0 then Byar_p=0;
midp=Byar_p;

END;
END;

drop Byar_p aa num deno denom  e  i  k   kk   numerator   subtotal   total z OBS_;          

%mend;



********************IMPORTING AN EXTERNAL DATASET AND RUN THE MACRO*********************
Assume a nominal value of 0.80;

data ClabsiExample;
input OBS EXP;
cards;
101 112
10 2
13 11
1 9
;
run;

data ClabsiExample_;set ClabsiExample;
EXP=EXP*0.8;*<----Assuming a nominal value of 0.80;
%sir(OBS,EXP); 
RUN;
title;

proc print;run;