**********************************************************************************************
This code was created 090914 by Anna Törner, anna.torner@smi.se, Swedish Institute for 
Infectious Disease Control.

The purpose of the code is to illustrate how calculations for a cumulative SIR graph may be 
set up.

In addition to the code below, a macro for time splitting  is needed. The macro 
used below is the lexis macro written by Bendix Carstensen and this macro may be downloaded from:
(http://staff.pubhealth.ku.dk/~bxc/Lexis/). The code given below includes the use of this lexis macro, 
but the actual macro itself must be downloaded form the referenced page.

The macro takes one parameter, i, size of the time window. This value is currently set to 10, for 
trial runs, but should be adjusted to an appropriate value, e.g. 730 (2 years).
Two data sets, COHORT and INC, are needed as input data sets for the code and macros below:

COHORT is the data set holding information on the cohort. The dataset COHORT should contain the following variables:

SEX			Man/Woman (coded as "1" and "2")
BIRTHDATE 	Date of birth (date format)
ENTRY 		Date of entry into cohort, i.e. when person is identified as member of cohort (date format)
EXIT		Date of studied outcome, death, emigration or end of follow-up, whichever occurs first (date format)
EXITSTATUS	Is set to "1" if the individual experiences the outcome of interest, otherwise "0".
DIADATE		The date of diagnosis of the outcome studied. If patient does not experience this outcome 
			this is set to "." (missing) (date format).


INC is the data set holding incidence rate of the studied outcome in the reference population.
The INC data set should contain the following variables:

SEX			Man/Woman (coded as "1" and "2")
YEAR		Calender year for which the incidence rate is valid
AGE			Age group for which the incidence rate is given
INC			Incidence rate per 100 000 individuals

 
*********************************************************************************************;
*STEP 1 - SETTING UP PROCEDURE FOR CALCULATING SIR
*********************************************************************************************;
*The macro SIR calculates the observed and expected number of events from day i until day 730;

%Macro SIR(i);

Data COHORT_MAC ;
	set COHORT;
	Exit=min((entry+730),exit);
	entry=entry+&i;
	if diadate ne . then do;
	if diadate=>entry and diadate=<exit then fail=exitstatus;
	else fail=0;
	end;
	else fail=0;
run; 

/* split by attained age *//*The use of the lexis macro is described at http://staff.pubhealth.ku.dk/~bxc/Lexis/   */
%Lexis ( origin = birthdate,
         entry = entry,
		 exit = exit,
          scale = 365.25,
		  data = COHORT_MAC,
		  out = c,
           breaks = %str(0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,120 ), /*split by same age groups as in INC data*/
           left = age
        ) ;

/* split by calendar year */
/* In Sweden the first year of the cancer register was 1958. The date given in macro below should be adjusted to give orgin as first year of cancer register*/
%Lexis ( origin = %str('01jan58'd),
         entry = entry,
		 exit = exit,
          scale = 365.25,
		  data = c,
		  out = d,
           breaks = %str(0 to 51 by 1),
           left = year
        ) ;


*Summarizing risk time divided by calender year, sex and age group;
Proc sort data=d;
		by  year sex age;
run;
Proc summary data=d;
		var risk fail;
		by year sex age;
		output out=RISKTIME sum=TIME EVENT;
run;

*Transforming the timescale into calender time;
Data RISKTIME2 (drop= _TYPE_ _FREQ_ );
		set RISKTIME;
		year=year+1958;
run;
********************************************************************************************************;
*Preparing to merge data set with risk time and the data set INC holding the incidence rates in the 
reference population.
*******************************************************************************************************;

Proc sort DATA=INC;
	by age year sex;
run;
Proc sort data=RISKTIME2;
	by age year sex;
run;
Data EXPECTED;
	merge INC(in=a) RISKTIME2(in=b);
	by age year sex;
	if b;
	exp=inc*time/100000;
run;

proc summary data=EXPECTED;
	var exp event;
	output out=OBSEXP (drop= _FREQ_) sum=EXP OBS;
run;


*Adding the observed and expected number of cases for follow-up time from day 730 and onwards;
**********************************************************************************************;
DATA OBSEXP2;
	set OBSEXP OBSEXP730;
run;
Proc summary data=OBSEXP2;
	var OBS EXP;
	output out= OBSEXP4 (drop= _FREQ_) sum=OBS EXP;
run;

**********************************************************************************************;
*Calculating SIR and adding confidence interval to the calculated SIR;
*Confidence intervals calculated using the Max method ref PK Kulkami et al, AJE 1998
**********************************************************************************************;
Data sir;
	set OBSEXP4;
	SIR=round(OBS/EXP,0.01);
	if OBS>0 then do;
		cl=0.5*cinv(0.025,2*obs)/exp;
		cu=0.5*cinv(0.975,2*(obs+1))/exp;
	end;
	else do;
		cl=0;
		cu=0.5*cinv(0.975,2*(obs+1))/exp;
	end;
EXP=round(exp,0.01);
CL=round(cl,0.01);
CU=round(cu,0.01);
Time=&i;
run;
*Adding calculated SIR for window i to dataset ALLSIR;
Data ALLSIR; 
set SIR ALLSIR;
run;
%MEND;
***************************************************************************************************
STEP 2 - REPEATED CALCULATIONS FOR DIFFERENT WINDOW SIZES
First part of code calculates the expected and observed number of events from day 731 and onwards. 
This "bulk" calculation is then added to the expected and observed number of events calculated for 
the first 730 days where the window size is increased one day at the time.
This calculation is performed once and then added to the calculations were one day is subtracted at 
the time. 

***************************************************************************************************;
%Macro doit;
***************************************************************************************************;

Data FROM2YEARS;
	set COHORT;
	entry=entry+730;
	if diadate ne . then do;
	if diadate=>entry and diadate=<exit then fail=exitstatus;
	else fail=0;
	end;
	else fail=0;
run; 

/* split by attained age */
%Lexis ( origin = birthdate,
         entry = entry,
		 exit = exit,
          scale = 365.25,
		  data = FROM2YEARS,
		  out = c,
           breaks = %str(0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,120 ), 
           left = age
        ) ;

/* split by calendar year */

%Lexis ( origin = %str('01jan58'd),
         entry = entry,
		 exit = exit,
          scale = 365.25,
		  data = c,
		  out = d,
           breaks = %str(0 to 51 by 1),
           left = year
        ) ;


*Summarizing risk time divided by follow-up after entry into the cohort;
Proc sort data=d;
		by  year sex age;
run;
Proc summary data=d;
		var risk fail;
		by year sex age;
		output out=RISKTIME (drop= _TYPE_ _FREQ_ ) sum=TIME EVENT;
run;

*Transforming the timescale into normal calender time (cancer register started in 1958 in Sweden);
Data RISKTIME2 ;
		set RISKTIME;
		year=year+1958;
run;
************************************************************************************;
*Preparing to merge with data file containing incidence rates in general population
************************************************************************************;
Proc sort DATA=INC;
	by age year sex;
run;
Proc sort data=RISKTIME2;
	by age year sex;
run;

Data EXPECTED;
	merge INC(in=a) RISKTIME2(in=b);
	by age year sex;
	if b;
	exp=inc*time/100000;
run;


proc summary data=EXPECTED;
	var exp event;
	output out=OBSEXP730 (drop= _FREQ_) sum=EXP OBS;
run;

	%do i=0 %to 10;
		%SIR(&i);
	%end;
%Mend;

Data ALLSIR;*Creating an empty data set to which calculated SIRs are added;
run;

%Doit; *Running all code


**************************************************************;
*STEP 3 - PLOTTING OF CUMULATIVE SIR 
**************************************************************;

title 'Graph of cumulative SIR';

symbol1 i=Join c=black line=2 repeat=2;
symbol2 i=Join c=black line=1;


proc gplot data=ALLSIR;
	plot cu*time cl*time SIR*time/overlay;
run;