Ken Kleinman

Here is SAS code to replicate example 8.18 from the blog, but collecting CI or equivalents instead of just means.  Note that different output data sets are grabbed with the ods output statements and that the renaming in the merge statement is more complicated.

data rlog;
do trial = 1 to 100;
x = 0; events = ranbin(0,100,.001); n = 100; output;
x = 1; events = ranbin(0,100,.05);    n = 100; output;
end;
run;

ods select none;
ods output clparmpl = glm;
proc logist data = rlog;
by trial;
model  events/n = x /  clparm=pl;
run;

ods output clparmpl = firth;
proc logist data = rlog;
by trial;
model  events/n = x / firth clparm=pl;
run;

ods output exactparmest = exact;
proc logist data = rlog;
by trial;
model  events/n = x;
exact x / estimate;
run;

data prior;
input _type_ $ Intercept x;
datalines;
Var 25 25
Mean 0 0
;
run;

ods output postsummaries = mcmcmean postintervals = mcmcints;
proc genmod data = rlog;
by trial;
model  events/n = x / dist=bin;
bayes nbi=1000 nmc=6000
coeffprior=normal(input=prior) diagnostics=none
statistics=(summary interval);
run;

data lregsep;
merge
glm (where = (parameter = "x") rename = (estimate = glm lowercl=gl uppercl=gu))
firth (where = (parameter = "x") rename = (estimate = firth lowercl=fl uppercl=fu))
exact (rename = (estimate = exact  lowercl=el uppercl=eu))
mcmcmean (where = (parameter = "x") rename = (mean=mcmc))
mcmcints (where = (parameter = "x"))
rlog (where = (x = 1) rename = (events = events1))
rlog (where = (x = 0) rename = (events = events0));
by trial;
events0 = events0 + (uniform(0) * 0.05);
events1 = events1 + (uniform(0) * 0.05);
run;

ods select all;
proc print data = lregsep; var events1 events0 gl glm gu fl firth fu el exact eu crediblelower mcmc credibleupper; run;