In PROC GENMOD (v9.4 TS 1M3, x64 for Windows), why is there a discrepancy between the exponentiated estimates as output using the EXP option in the LSMEANS statement and the calculated values (i.e. doing the exponentiation calculation in a data step) when the BAYES statement is used, but not when a frequentist analysis is done? The values differ very slightly in the Bayesian analysis. Am I missing something?
Here is some code, run on an example in the book by Stokes, Davis and Koch (p. 354).
data respire;
do region=1,2;
do agegp=2,3,4,5,6,1;
input cases total;
ltotal=log(total);
output;
end;
end;
datalines;
75 220407
68 198119
63 134084
45 70708
27 34233
64 1074246
76 564535
98 592983
104 450740
63 270908
80 161850
61 2880262
run;
proc format;
value agfmt 1='<35'
2='35-44'
3='45-54'
4='55-64'
5='65-74'
6='>75'
;
value regfmt 1='South'
2='North'
;
run;
proc print data=respire;
format agegp agfmt. region regfmt.;
run;
/* Frequentist*/
ods output diffs=frediffs_all ;
proc genmod data=respire;
class agegp region;
model cases=agegp region/dist=poisson link=log offset=ltotal;
lsmeans agegp /diff=control("1") exp;
lsmeans region /diff=control("1") exp;
title 'Koch p 354 Frequentist Poisson analysis';
run;
ods output close;
/* Bayesian*/
ods output diffs=baydiffs_all ;
proc genmod data=respire;
class agegp region;
model cases=agegp region/dist=poisson link=log offset=ltotal;
bayes seed=555 STATS(alpha=0.05 percent=2.5 25 50 75 97.5 )=all nbi=2000 nmc=300000 diagnostics=(autocorr ess gelman(nchain=3)) ;
lsmeans agegp /diff=control("1") exp;
lsmeans region /diff=control("1") exp;
title 'Koch p 354 Bayesian Poisson analysis';
run;
quit;
ods output close;
proc sort data=frediffs_all;
by agegp region;
run;
proc sort data=baydiffs_all;
by agegp region;
run;
/* Compare exponentiated (as calculated) vs. printed*/
data allestdiff;
set frediffs_all(in=f)
baydiffs_all(in=b);
by agegp region;
if f then antype='Frequentist';
else if b then antype='Bayesian';
expest=exp(estimate);
label expest='Expomnentiated estimate (calculated)';
run;
options ls=132;
proc compare data=allestdiff ;
id antype agegp region;
var expestimate ;
with expest;
title 'Comparison of exponentiated estimates as output, and calculated';
run;
proc print data=allestdiff label;
format region regfmt. agegp agfmt. estimate expest expestimate 8.6;
id antype;
var agegp region estimate expest expestimate;
title 'Estimates from GENMOD using frequentist and Bayesian analyses : exponentiated as output and as calculated';
run;
