data diabetes;

input diabetes age py accident;
py000 = py/1000;
ln_py = log(py000);  diab_py = diabetes*py000;

lines;
1 20    65.2    3
0 20 26657.9 2177
1 30    81.2    6
0 30 27145.3 1326
1 40   136.2    9
0 40 18500.9  830
1 50   306.1   14
0 50 11620.0  456
1 60   502.1   24
0 60 10515.1  336
1 70   717.7   32
0 70 10625.3  340
;
run;


options ls=75 ps=55 nodate ;
*   **** for mantel-haenszel  **** ;
proc format;
   value diabfmt 0='not'
                 1='diabetic';
title1 'tabular... ' ;
data for_mh; keep age diabetes crash fr;
format diabetes diabfmt.;
set diabetes;
crash = 1; fr = accident; output;
crash = 0; fr = py - accident; output;
run;


title1 'tabular... ' ;
title2 'crude... ' ;
proc freq data=for_mh order=formatted;
 tables diabetes*crash / riskdiff relrisk chisq cmh;
 weight fr;
run;

title1 'tabular... ' ;
title2 'mantel-haenszel... ' ;
proc freq data=for_mh order=formatted;
 tables age*diabetes*crash / riskdiff relrisk chisq cmh;
 weight fr;
run;



title1 'by regression ' ;
title2 'crude rate difference ' ;
proc genmod data=diabetes;
model accident = py000 diab_py/
      link=identity dist = poisson noint ;
run;

title2 'crude rate ratio ' ;
proc genmod data=diabetes;
model accident = diabetes/
      link=log dist = poisson offset = ln_py;
run;

title2 'adjusted rate difference ' ;
proc genmod data=diabetes;
class age;
model accident = age*py000 py000 diab_py/
      link=identity dist = poisson noint ;
run;

title2 'adjusted rate ratio ' ;
proc genmod data=diabetes;
class age;
model accident = age diabetes  /
      link=log dist = poisson offset = ln_py;
run;