Beta - Binomial Bayesian Model

[kc]

I am trying to fit a logistic regression model to a dataset that where the outcome is binary (0, 1), and the independent variable is treatment group (2 values, 0 and 1). I would like a beta (1,1) / Uniform(0, 1) prior on the probability of success, p.

I am using the following code to model the outcome. There are missing values in the outcome variable which are excluded from the analysis using proc genmod.

Question: Is it possible to impute missing values in proc genmod, so all of the data can be used in the analysis?

/***************************************************************************************/

proc genmod data=ds;
model s = trt /d=bin link=logit;
bayes seed=17 coeffprior=uniform nmc=20000 thin=2 outpost=outgen;
run;

/***************************************************************************************/

As an alternative, I am using proc mcmc which imputes the missing outcome values. But, I am a bit confused on how to assign a beta prior to the parameter p in the binomial distribution.

/***************************************************************************************/

proc mcmc data=ds seed=123 outpost=logit_bayes nbi=1000 nmc=10000 nthin=1 DIC statistics=all diagnostics=all;
parms (beta0 beta1 /*p*/);
prior beta0 ~ normal(0, var=10000);
prior beta1 ~ normal(0, var=10000);
/*prior p ~ beta(1,1);*/
p=logistic(beta0 + beta1*trt);
model s ~ binomial(1,p);
run;

/***************************************************************************************/

Question: Is there a way to modify/change the mcmc code above to specify a beta prior so inferences about the parameter p can be made?

[StatDave]

That all seems like overkill for what is really just a 2x2 table. If your goal is to estimate the probability of the event in each treatment and to compare them, then a model-based approach is not needed - you can do that in PROC FREQ using the RISKDIFF option. There are even options to do an exact analysis if desired. Imputation can be done using PROC MI, but with no other variables available, I can't see how any meaningful imputation can be done.

[kc]

Agreed. But, the exercise is to run a bayesian analysis. If you were to still use proc mcmc, how do you specify a beta prior on the binomial parameter p, instead of specifying priors on the regression parameters beta0 and beta1. Apologies if the question itself isn't clear or doesn't make sense.

[SteveDenham]

The documentation for PROC FMM and HPFMM have the equations for the likelihood for a beta-binomial. Those could be incorporated into MCMC.

SteveDenham

[kc]

When using the proc mcmc code below , I get the following error:

ERROR: Parameters cannot have values assigned to them in the program. Change the assignment statement for the parameter prob.

proc mcmc data=ds outpost=PosteriorSample nmc=20000

  propcov=quanew ntu=1000 nthin=2 seed=7 diag=all stats=all dic;

  parms b0 b1 prob;

  prior b0 ~ normal(0.01, var=10000);

  prior b1 ~ normal(0.01, var=10000);

  prior prob ~ beta(1,1);

  xb = b0 + b1*rand_arm;

  prob = exp(xb)/(1+exp(xb));

  ll = log((prob**m6)*((1-prob)**(1-m6)));

  model m6 ~ general(ll);

  ods select PostSummaries PostIntervals TADPanel;

run;

[SAS_Rob]

Since Prob is just the predicted value after applying a link function to the linear predictor, it should not be listed as a parameter in the model.  In other words, you should remove it from the PARMS and PRIOR statements.

[kc]

But, the purpose is to make inferences about the binomial probability 'prob' as mentioned in the original post. This code was just to show that I tried writing out the likelihood function for the model.

[SAS_Rob]

You can still get posterior intervals, plots, summaries etc. for Prob.  But it is just a function of the other parameters and thus is not truly a parameter.

proc mcmc data=mylib.ds outpost=PosteriorSample nmc=20000

propcov=quanew ntu=1000 nthin=2 seed=7 diag=all stats=all dic monitor=(prob);

parms b0 b1;

prior b0 ~ normal(0.01, var=10000);

prior b1 ~ normal(0.01, var=10000);

xb = b0 + b1*rand_arm;

prob = exp(xb)/(1+exp(xb));

ll = log((prob**m6)*((1-prob)**(1-m6)));

model m6 ~ general(ll);

ods select PostSummaries PostIntervals TADPanel;

run;