Sampling from a posterior distribution in SAS

Question

Lets say I have a dataset where I want to estimate the relative risk of outcome X based on a binary treatment level Y, using PROC GENMOD to fit a logistic regression model. I can use the BAYES statement to have SAS 9.2 produce nice, MCMC-based prior distributions for my data.

My question is, what if I then want to sample from that posterior distribution? For example, if I estimate a OR and its 95% credible interval for the effect of Y, and subsequently want to use that estimate in another analysis. If I wanted to perform a sensitivity analysis based on the possibility of error in my estimate of the effect of Y, it makes sense to just run a Monte Carlo simulation based on the distribution of Y. But in this case, rather than being a really straightforward random number generation, I need to pull from that prior.

Is the best way to do this using the OUTPOST option to get the data set of generated posterior samples and then just randomly sample from that new data set? Or is there some other, more clever mechanism I'm not seeing?

Answer 1

Here is my proposed answer to my own question, in SAS, using the method I alluded to previously. I can't necessarily think of another way to go about doing this, and I had a good long boring drive to ponder it. Again, if someone has a better way, I'd be happy to hear it.

Generally, what this code does is generate a simple cohort-type data set with a fixed relative risk, then run a GLM on said cohort using the inbuilt MCMC routines to get a posterior distribution. Because this is a proof of concept, I'm using uninformative priors. This distribution is then sampled without replacement into a new data set.

data work.tinkering;
do id=1 to 20000;
        disroll = ranuni(68273);
yran = ranuni(102842);
if yran <= 0.333 then y = 1; else if yran > 0.333 then y = 0;
pdis = exp(-2.30258+(1.0986*y));
if disroll<= pdis then case=1; else if disroll>pdis then case=0;
output;
end;

ods graphics on;
proc genmod data=work.tinkering descending;
    model case = y / dist=binomial link=log;
    bayes nbi=2000 nmc=10000 plots=all outpost=work.out;       
ods graphics off;

proc sgplot data=work.out;
    density y;
    density y / type=kernel;

proc surveyselect data=work.out
        method=urs outhits n=1000 out=work.SampledPrior;
run;

Graphical comparison of the two: