I contacted the author of the mentioned paper and he was kind enough to help me with providing a code for the exact problem described above. I will copy the code here so others with the same question may also use it.
I agree with Doc@Duke's comments -- there's some sense in Bayesian statistics that with moderate to large sample sizes, the sample should overwhelm the prior when you're making inferences, as long as you've chosen reasonable, relatively non-informative priors. A Bayesian analysis still has "merit" -- you can interpret the posterior estimates as true probabilities, and you have the ability to incorporate your prior beliefs into your estimates. That said, this model will probably take a long time to run, and will probably take a large number of iterations to reach convergence, so you'll have to weight whether the properties of a Bayesian analysis are worth the headache.
Good luck!
Jake
data niam; input id choice $ age female ses1 - ses4; /* Recode choice into dummy variables */ choice_y = 0; choice_e = 0; choice_n = 0; if choice eq "Y" then choice_y = 1; else if choice eq "E" then choice_e = 1; else if choice eq "N" then choice_n = 1; cards; 1 Y 29 1 1 0 0 0 2 E 39 0 0 1 0 0 3 Y 18 1 0 0 0 1 4 N 30 1 1 0 0 0 ; run; proc mcmc data = niam nmc = 1000 nbi = 100 outpost = posterior; /* Data */ array Y[3] choice_y choice_e choice_n; /* Parameters */ array p[3]; array beta0[2]; array beta1[2]; array beta2[2]; array beta3[2]; array beta4[2]; array beta5[2]; parms beta:; /* You may want to divide these up into separate "parms" statements, because I think if you use a single statement, all the parameters are proposed together, which might lead to slower convergence */ /* Priors */ prior beta: ~ norm(0, var = 100); /* Independent, normal priors for the beta coefficients in the model. You should change the priors to represent your prior beliefs */ /* linear predictor */ array eta[2]; do i = 1 to 2; eta = beta0 + beta1 * age + beta2 * female + beta3 * ses2 + beta4 * ses3 + beta5 * ses4; end; denominator = 1 + exp(eta[1]) + exp(eta[2]); /* Inverse logit transformation */ do i = 1 to 2; p = exp(eta) / denominator; end; p[3] = 1/denominator; /* Sampling model */ model y ~ multinom(p); run;
... View more