/*
  Original simulated example: Bayesian logistic model for two binomial groups.
  Run this file from the folder that contains example.csv.
*/

options nodate nonumber;

proc import datafile="example.csv"
            out=bayes_counts
            dbms=csv replace;
  guessingrows=max;
run;

ods graphics on;
proc mcmc data=bayes_counts
          outpost=posterior
          seed=20260821
          nbi=5000 nmc=40000 thin=5
          diag=(mcse ess)
          plots=trace
          monitor=(alpha beta p_control p_druga risk_diff odds_ratio);
  parms alpha 0 beta 0;
  prior alpha beta ~ normal(0, var=100);

  p=logistic(alpha + beta*arm);
  p_control=logistic(alpha);
  p_druga=logistic(alpha + beta);
  risk_diff=p_druga - p_control;
  odds_ratio=exp(beta);

  model responders ~ binomial(total, p);
run;
ods graphics off;

/* 2.5% and 97.5% limits are already printed in MCMC Posterior Intervals. */
proc means data=posterior n mean std p50 maxdec=4;
  var p_control p_druga risk_diff odds_ratio;
run;

data posterior_decision;
  set posterior;
  any_benefit=(risk_diff > 0);
  clinically_relevant=(risk_diff > 0.10);
run;

proc means data=posterior_decision mean maxdec=4;
  var any_benefit clinically_relevant;
  label any_benefit='Posterior Pr(RD > 0)'
        clinically_relevant='Posterior Pr(RD > 0.10)';
run;