options nofmterr;
libname in 'c:\';

data one; set in.sub;
if birth_yr='9999' then birth_yr=' ';
b_year=birth_yr+0;   ***change to numeric for birth year;
survey_sy=year(survey_start); ***only use the survey calendar year in calculations;
survey_ey=year(survey_end);   ***all survey ended in 1997!!!;
proc freq;
  tables survey_sy survey_ey;
run;

%macro py(a);
data two; set one;
if survey_sy ge 1986 and survey_sy le &a and survey_ey ge &a; ***py for the calendar year;
if dead=1 then delete; **** dead were excluded from p-y calculation;
age_py=&a-b_year;      **** age was repeatly calculated by calendar year;
if age_py ge 18 and age_py lt 25 then agecat=1;
if age_py ge 25 and age_py lt 30 then agecat=2;
if age_py ge 30 and age_py lt 35 then agecat=3;
if age_py ge 35 and age_py lt 40 then agecat=4;
if age_py ge 40 and age_py lt 45 then agecat=5;
if age_py ge 45 and age_py lt 50 then agecat=6;
if age_py ge 50 and age_py lt 55 then agecat=7;
if age_py ge 55 and age_py lt 60 then agecat=8;
if age_py ge 60 and age_py lt 65 then agecat=9;
if age_py ge 65 and age_py lt 70 then agecat=10;
if age_py ge 70 and age_py lt 75 then agecat=11;
if age_py ge 75 and age_py lt 80 then agecat=12;
if age_py ge 80 and age_py lt 85 then agecat=13;
if age_py ge 85 then agecat=14;
/*proc univariate normal plot;
  var age_py;
*/;
proc freq;
  tables agecat*sex/norow nopercent nocol noprint missing out=py&a;
run;

data py&a; set py&a;
year=&a;
rename count=p_at_risk;
drop percent;
proc sort; by agecat year sex;
run;
%mend py;

%py(1986);
%py(1987);
%py(1988);
%py(1989);
%py(1990);
%py(1991);
%py(1992);
%py(1993);
%py(1994);
%py(1995);
%py(1996);
%py(1997);

data p_at_risk; set py1986 py1987 py1988 py1989 py1990 py1991
                    py1992 py1993 py1994 py1995 py1996 py1997;
 by agecat sex;
proc sort; by sex year agecat;
proc print;
run;

data in.p_at_risk; set p_at_risk; if agecat ne .;  *** a few with age missing were excluded;
proc sort; by agecat sex;
proc means noprint;
  var p_at_risk;
  by agecat sex;
  output out=age_sex_py sum=py;   ***py was added for each age group;
run;

data age_sex_py; set age_sex_py;
keep agecat sex py;
proc sort; by agecat;
proc print;
run;

data mortality;    **** mortality rate for us population in 2000;
input agecat sex $ mrate@@;
datalines;
1 1 142.0 1 2 48.2
2 1 141.9 2 2 56.5
3 1 157.3 3 2 76.0
4 1 209.8 4 2 115.1
5 1 303.3 5 2 172.2
6 1 461.7 6 2 254.0
7 1 658.3 7 2 386.3
8 1 1007.5 8 2 611.8
9 1 1565.5 9 2 982.0
10 1 2399.3 10 2 1527.5
11 1 3705.4 11 2 2381.8
12 1 5591.2 12 2 3812.6
13 1 8956.9 13 2 6444.8
14 1 16605.4 14 2 14768.6
;
proc sort; by agecat;
run;

data smr; merge age_sex_py mortality; by agecat;
exp_death=py*mrate/100000;
proc sort; by sex;
proc means noprint;
  var exp_death;
  output out=exp sum=exp_death;
  by sex;          ***subtotals of expected deaths by sex;
run;

data exp_death; set exp;
keep sex exp_death;
proc sort; by sex;
run;

data four; set one; if dead=1;  ***is there other way to define observed death?;
proc freq;
  tables sex/out=obs;   ***observed deaths by sex;
run;

data obs_death; set obs;
rename count=obs_death;
drop percent;
proc sort; by sex;
run;

data smr1; merge exp_death obs_death; by sex;
smr=round(obs_death/exp_death,.01);
ll=obs_death*(1-1/(9*obs_death)-(1.96/3)*sqrt(1/obs_death))**(3)/exp_death;
ul=(obs_death+1)*(1-1/(9*(obs_death+1))+(1.96/3)*sqrt(1/(obs_death+1)))**(3)/exp_death;
ll=round(ll,.01);
ul=round(ul,.01);
if (ll ge 0 and ul lt 1) or ll gt 1 then sig='*';
proc print;
run;