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;
... View more
