How to implement Bayesian Estimation Supersedes the t Test (BEST) in SAS?

[supp]

Has anyone implemented the BEST approach to compare two group means and their differences?

I am interested in a Bayesian approach to a two sample t test. I found this paper which describes an approach referred to as BEST.

Kruschke2013JEPG.pdf (iu.edu) 

From the paper here is a visual of the model used:

Here is my attempt to apply BEST approach described in the paper to the "Behrens-Fisher Problem" described in the SAS MCMC procedure.

The data:

data behrens;
   input y ind @@;
   datalines;
121 1  94 1 119 1 122 1 142 1 168 1 116 1
172 1 155 1 107 1 180 1 119 1 157 1 101 1
145 1 148 1 120 1 147 1 125 1 126 2 125 2
130 2 130 2 122 2 118 2 118 2 111 2 123 2
126 2 127 2 111 2 112 2 121 2
;

My attempt to apply the BEST method. It seems to me the main difference is the using the t distribution as the likelihood function (as opposed to a normal distribution used in the SAS documentation)

proc sql;
	select mean(y) into :mean_y from behrens;
quit;


/**
Get pooled data
**/
proc glm data= behrens;
	class ind;
	model y = ind;
run;

* Root MSE = 19.32394 ;

%let low_pooled_std = 19.32394 / 1000;
%put &=low_pooled_std. ;


%let high_pooled_std = 19.32394 * 1000;
%put &=high_pooled_std. ;


proc mcmc data=behrens outpost=postout2 seed=123
          nmc=40000 monitor=(_parms_ mudif)
          statistics(alpha=0.01);
   ods select PostSumInt;

   parms mu1 0 mu2 0;
   parms sig21 1;
   parms sig22 1;
   parms nu 1;

   prior mu: ~ N(&mean_y., sd= &high_pooled_std.);  				* prior assumes pooled mean and normal distribution ;
   prior sig2: ~ uniform(&low_pooled_std., &high_pooled_std.);
   prior nu: ~ expon(scale= 29);									* From Kruschke paper, exponential distribution spreads prior credibility fairly evenly over nearly normal and heavy tailed data ;
      
   mudif = mu1 - mu2;

   if ind = 1 then do;
      mu = mu1;
      s2 = sig21;
   end;
   else do;
      mu = mu2;
      s2 = sig22;
   end;
   model y ~ t(mu, var=s2, nu);

/*   model y ~ n(mu, var=s2);		Use this if a normal distribution is desired. The t distribution should handle outlier better*/
run;

Here are the estimates. These match pretty close the to the example in SAS documentation:

proc sql;
select 'Probability difference of means if greater than 0', sum(mudif > 0) / count(*) as probability from postout2;
run;

The priors are subjective, but I just used the approach described in the paper. These should be modified to fit the analysis being conducted.

If anyone has attempted to implement the BEST approach I would appreciate any feedback on my approach. 

[sbxkoenk]

Hello,

I will have an in-depth look to your post tomorrow (it's 22h00 over here).

Just this quick message :

Bayesian t-test in proc MCMC
https://communities.sas.com/t5/Statistical-Procedures/Bayesian-t-test-in-proc-MCMC/m-p/709482#M34358

( This is about Paired Bayesian t test , but it might give you some more inspiration. )

Koen

[Ksharp]

https://support.sas.com/kb/23/407.html

Maybe @Rick_SAS could wrote a IML code for you .

Also could check PROC GENMOD.

[supp]

I should mention I don't have access to IML