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

                                 Bayesvc
                                      

  DISCLAIMER:
    THIS MATERIAL IS PROVIDED BY SAS "AS IS". SAS has not necessarily
    tested this material, and use of the material is at Recipient's
    own risk. SAS is under no obligation to maintain or support this
    material. Recipient acknowledges and agrees that SAS has no
    liability for any damages whatsoever arising out of Recipient's
    use of or inability to use this material, even if SAS has been
    advised of the possibility of such damages. 

  PURPOSE:
    Examples from:
    Wolfinger, R.D. and Kass, R.E. (2000), Non-Conjugate Bayesian
      analysis of variance component models, Biometrics, 56, 768-774.
                                                                      
  REQUIRES:
    Version 7 or later of Base SAS, SAS/STAT, SAS/INSIGHT, and
    SAS/GRAPH.

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

/*-----------------------------*/
/*---Randomized Block Design---*/
/*-----------------------------*/

/*---generated data from Box and Tiao, 1973, Bayesian Inference, 
     p. 247---*/
data gen;
   input batch y;
   datalines;
1  7.298
1  3.846
1  2.434
1  9.566
1  7.990
2  5.220
2  6.556
2  0.608
2 11.788
2 -0.892
3  0.110
3 10.386
3 13.434
3  5.510
3  8.166
4  2.212
4  4.852
4  7.092
4  9.288
4  4.980
5  0.282
5  9.014
5  4.458
5  9.446
5  7.198
6  1.722
6  4.782
6  8.106
6  0.758
6  3.758
run;

/*---jeffreys prior---*/
proc mixed data=gen;
   class batch;
   model y = / s;
   random batch / s;
   prior / nsample=1e4 lognote=1e3 seed=34875770 ptrans psearch 
      out=sample;
run;

data bayes;
   set sample;
   sigma2b = covp1;
   sigma2e = covp2;
   ratio = sigma2b/sigma2e;
   icc = sigma2b/(sigma2e + sigma2b);
   beta0 = beta1;
   array uu{6} gam1-gam6;
   array gg{6} gamma1-gamma6;
   do i = 1 to 6;
      gg{i} = uu{i};
   end;
   tau1 = 5*covp1 + covp2;
   keep sigma2b sigma2e ratio icc beta0 gamma1-gamma6 tau1;
run;

/*---SAS/INSIGHT and Proc Univariate are useful for analyzing the 
     sample---*/
proc insight data=bayes;
run;

proc univariate data=bayes;
run;

/*---compute informative prior from quantiles---*/
data ig;
   input tol ptarg @@;
   array mmode{2} mode1-mode2;
   array qq{2} q1-q2;
   array cc{2,2} c1-c4;
   type = 'ig';
   do tau = 1 to 2;
      input mmode{tau} qq{tau} @@;
      do i = 1 to 2;
         input cc{tau,i} @@;
      end;
   end;
   do tau = 1 to 2;
      mode = 0;
      q = 0;
      do i = 1 to 2;
         mode = mode + cc{tau,i}*mmode{i};
         q = q + cc{tau,i}*qq{i};
      end;
      al = 1;
      au = 100;
      /*---bisection search---*/
      do while ((au - al) > tol);
         a = (au + al)/2;
         b = mode*(a+1);
         x = b/q;
         p = 1 - probgam(x,a);
         if (p < ptarg) then al = a;
         else au = a;
      end;
      parm1 = a;
      parm2 = b;
      output;
   end;
   keep type parm1 parm2;
   datalines;
1e-4 0.95  
 5  10 5 1
10  20 0 1
proc print;
run;

/*---generate sample using informative prior---*/
proc mixed data=gen;
   class batch;
   model y = / s;
   random batch / s;
   prior data=ig / nsample=1e4 lognote=1e3 seed=34875770 ptrans psearch 
      out=samplep;
run;

data bayesp;
   set samplep;
   sigma2b = covp1;
   sigma2e = covp2;
   ratio = sigma2b/sigma2e;
   icc = sigma2b/(sigma2e + sigma2b);
   beta0 = beta1;
   array uu{6} gam1-gam6;
   array gg{6} gamma1-gamma6;
   do i = 1 to 6;
      gg{i} = uu{i};
   end;
   keep sigma2b sigma2e ratio icc beta0 gamma1-gamma6;
run;



/*---kernel density estimate, just change the var statement
     to analyze other variables---*/
proc kde data=bayes out=den ng=100 method=srot bwm=1.5;
   var sigma2b;
   where sigma2b < 10;
symbol i=join;
proc gplot data=den;
   plot density*sigma2b; 
run;

/*---macro to compute kernel density estimates for 
     the random effects---*/
%macro kde(data=);
%do i = 1 %to 6;
   proc kde data=bayes&data out=g&data&i ng=100 method=srot bwm=2;
      var gamma&i;
      where gamma&i > -5 and gamma&i < 5;
   run;
   data g&data&i;
      set g&data&i;
      length parm$ 10;
      gamma = gamma&i;
      drop gamma&i;
      i = &i;
      parm = "gamma&data&i";
      drop i count;
   run;
%end;
%mend;

ods exclude all;
run;
%kde(data=)
%kde(data=p)
ods select all;
run;

data all;
   set g3 g5 g6 gp3 gp5 gp6;
   if (density >= 0);
run;

proc gplot data=all;
   plot density*gamma=parm;
run;



/*-----------------------*/
/*---Split-Plot Design---*/
/*-----------------------*/

/*---balanced data---*/
data sp;
   input block a b y;
   cards;
   1 1 1  56
   1 1 2  41
   1 2 1  50
   1 2 2  36
   1 3 1  39
   1 3 2  35
   2 1 1  30
   2 1 2  25
   2 2 1  36
   2 2 2  28
   2 3 1  33
   2 3 2  30
   3 1 1  32
   3 1 2  24
   3 2 1  31
   3 2 2  27
   3 3 1  15
   3 3 2  19
   4 1 1  30
   4 1 2  25
   4 2 1  35
   4 2 2  30
   4 3 1  17
   4 3 2  18
run;

/*---basic REML analysis---*/
proc mixed data=sp;
   class a b block;
   model y = a|b / s;
   random block a*block / s;
   lsmeans a b / diff;
run;

/*---Bayesian calculations---*/
proc mixed data=sp;
   class a b block;
   model y = a|b;
   random block a*block;
   prior / nsample=1e4 lognote=1e3 seed=34875770 psearch ptrans
      out=sample; 
run;

data bayes;
   set sample;
   sigma2b = covp1;
   sigma2w = covp2;
   sigma2e = covp3;
   sqrt2rho = sqrt(2*t3f3);
   drop covp1-covp3;
run;

proc insight data=bayes;
run;

/*---kernel density estimate---*/
%let var = sigma2e;

proc kde data=bayes out=denb gu=100 ng=1000 method=srot;
   var &var;
run;

symbol i=join;
proc gplot data=denb;
   plot density*&var;
run;


/*---unbalanced data---*/
data spu;
   input block a b y;
   cards;
   1 2 1  50
   1 2 2  36
   1 3 1  39
   1 3 2  35
   2 1 1  30
   2 1 2  25
   2 2 2  28
   2 3 1  33
   2 3 2  30
   3 1 2  24
   3 2 1  31
   3 2 2  27
   3 3 2  19
   4 1 1  30
   4 1 2  25
   4 3 2  18
run;

/*---basic REML analysis---*/
proc mixed data=spu;
   class a b block;
   model y = a|b / s;
   random block a*block / s;
   lsmeans a b / diff;
run;

/*---Bayesian calculations---*/
proc mixed data=spu;
   class a b block;
   model y = a|b / s;
   random block a*block;
   prior / nsample=1e4 lognote=1e3 seed=34875770 psearch ptrans 
      out=sample;
run;

data bayesu;
   set sample;
   sigma2b = covp1;
   sigma2w = covp2;
   sigma2e = covp3;
   sqrt2rho = sqrt(2*t3f3);
   drop covp1-covp3;
run;

proc insight data=bayesu;
run;

/*---kernel density estimate---*/
%let var = sigma2b;

proc kde data=bayesu out=denbu gu=100 ng=1000 method=srot;
   var &var;
run;

symbol i=join;
proc gplot data=denbu;
   plot density*&var;
run;



/*---------------------*/
/*---Quadratic Trend---*/
/*---------------------*/

data rb;
   input base y1-y4;
   rabbit = _n_;
   array yy{4} y1-y4;
   if (rabbit <= 11) then trt = 'A';
   else trt = 'B';
   do i = 1 to 4;
      y = yy{i};
      hour = i*1.5;
      time = hour - 3;
      output;
   end;
   drop i y1-y4;
   datalines;
 77 52 35 56 64 
 77 56 43 56 64
 77 43  . 39 64
 81 56 22 35 39
 90 73 52 60 64
103 47 22 60 99
 99  .  . 47 81
 90 35 22 26  .
 90 64 47  .  .
 90 35 12 26 68
 85 68 56 68 73
103 26 22 39 68
 85 68 56 64 81
 85 35  . 64  .
 85 35 47 77  .
 85 35  .  . 77
103 56 47 64  .
 90 52  . 35 39
 81 22  .  . 77
 85 30 22 52 94
 85 47 33 77  .
 94 39 35 52 85
run;

/*---plot profiles---*/
proc gplot data=rb;
   by trt;
   plot y*hour=rabbit / hminor=1 
      vaxis=0 to 120 by 24 vminor=1;
run;

/*---random intercept model using REML---*/
proc mixed data=rb;
   class trt rabbit;
   model y = base trt|time|time / s;
   random int / sub=rabbit;
run;

/*---random intercept and slope model using REML---*/
proc mixed data=rb;
   class trt rabbit;
   model y = base trt|time|time / s;
   random int time / type=un sub=rabbit;
run;

/*---random intercept with Bayesian calculations---*/
proc mixed data=rb;
   class trt rabbit;
   model y = base trt|hour|hour / s;
   random int / sub=rabbit;
   prior / nsample=1e4 lognote=1e3 seed=23487013 ptrans psearch
      out=sample;
run;

data bayes;
   set sample;
   s2g = covp1;
   s2e = covp2;
   base = beta2;
   aA = beta8 + beta9;
   aB = beta8;
   bA = beta5 + beta6;
   bB = beta5;
   cA = beta1 + 85*beta2 + beta3;
   cB = beta1 + 85*beta2;
   bottomA = - bA/2/aA;
   solidusA = max((bA*bA - 4*aA*(cA-85)),0);
   t85A = (-bA + sqrt(solidusA))/2/aA;
   bottomB = - bB/2/aB;
   solidusB = max((bB*bB - 4*aB*(cB-85)),0);
   t85B = (-bB + sqrt(solidusB))/2/aB;
   drop covp1 covp2 aA aB bA bB cA cB beta1-beta10 t3f1-t3f6;
run;

proc insight data=bayes;
run;

proc kde data=bayes;
   var base;
run;


/*---bivariate density estimate of variance components---*/
proc kde data=bayes out=den ng=40,40 gl=0,75 gu=400,275 bwm=2,2;
   var s2g s2e;
run;

proc g3d data=den;
   plot s2e*s2g=density; 
run;

/*---plot of predicted trends---*/
data sasuser.rbsamp;
   set sample;
run;

%macro kde(npoint,trt);

   ods exclude all;
   
   options nonotes;
   
   %do i = 1 %to &npoint;
   
      %put Computing density &i out of &npoint;
      data b;
         set sasuser.rbsamp; 
         hour = 6*(&i - 1)/(&npoint - 1);
         %if (&trt = A) %then %do;
            sugar = beta1 + 85*beta2 + beta3 + (beta5+beta6)*hour + 
                (beta8+beta9)*hour*hour; 
         %end;
         %else %do;
            sugar = beta1 + 85*beta2 + beta5*hour + beta8*hour*hour; 
         %end;
         keep sugar;
      run;
   
      proc kde data=b out=bo bwm=2 ng=100 gl=0 gu=100;
         var sugar;
         ods output percentiles=p;
      run;
   
      data bo;
         set bo;
         hour = 6*(&i - 1)/(&npoint - 1);
         keep hour sugar density;
      run;
   
      data p;
         set p;
         hour = 6*(&i - 1)/(&npoint - 1);
         keep hour sugar percent;
      run;
   
      %if (&i = 1) %then %do;
         data surf;
            set bo;
         run;
         data pct;
            set p;
         run;
      %end;
      %else %do;
         proc append base=surf data=bo;
         run;
         proc append base=pct data=p;
         run;
      %end;
   
   %end;
   
   ods select all;
   options notes;

%mend;


/*---fit densities for treatment A---*/
%kde(100,A)
data surfa;
   set surf;
   if density < 0 then density = 0;
data pcta;
   set pct;
proc sort data=pcta;
   by percent hour;
run;


/*---plot percentiles for treatment A---*/
quit;
goptions reset=all device=xcolor cback=white
   hsize   = 6 in
   vsize   = 6 in
   htext   = 3.0 pct
   htitle  = 3.5 pct
   vorigin = 0 in
   horigin = 0 in
   ftext=swissl lfactor=1;
symbol1 i=j w=5 l=2 c=black;
symbol2 i=j w=5 l=3 c=black;
symbol3 i=j w=5 l=4 c=black;
symbol4 i=j w=5 l=1 c=black;
symbol5 i=j w=5 l=4 c=black;
symbol6 i=j w=5 l=3 c=black;
symbol7 i=j w=5 l=2 c=black;
footnote1 font=italic h=1.5 j=c "(a)";
axis1 order=(0 to 6 by 2)  value=(font=swissl h=1) 
   label=(font=swissl h=1) major=(h=0.5) minor=(n=1 h=0.5);
axis2 order=(0 to 120 by 30) value=(font=swissl height=1) 
   label=(font=swissl h=1) major=(h=0.5) minor=(n=1 h=0.5);
proc gplot data=pcta;
   plot sugar*hour=percent / haxis=axis1 vaxis=axis2 noframe;
   where percent ne 1.0 and percent ne 10.0 and
         percent ne 25.0 and percent ne 75.0 and
         percent ne 90.0 and percent ne 99.0;
run;

/*
filename gsasfile "sugara.ps";
goptions device=pslepsfc gaccess=sasgaedt gsfname=gsasfile
   gsfmode=replace;
run;
*/

/*---fit densities for treatment B---*/
%kde(100,B)
data surfb;
   set surf;
   if density < 0 then density = 0;
data pctb;
   set pct;
proc sort data=pctb;
   by percent hour;
run;

/*---plot percentiles for treatment B---*/
quit;
goptions reset=all device=xcolor
   hsize   = 6 in
   vsize   = 6 in
   htext   = 3.0 pct
   htitle  = 3.5 pct
   vorigin = 0 in
   horigin = 0 in
   ftext=swissl lfactor=1;
symbol1 i=j w=5 l=2 c=black;
symbol2 i=j w=5 l=3 c=black;
symbol3 i=j w=5 l=4 c=black;
symbol4 i=j w=5 l=1 c=black;
symbol5 i=j w=5 l=4 c=black;
symbol6 i=j w=5 l=3 c=black;
symbol7 i=j w=5 l=2 c=black;
footnote1 font=italic h=1.5 j=c "(b)";
axis1 order=(0 to 6 by 2)  value=(font=swissl h=1) 
   label=(font=swissl h=1) major=(h=0.5) minor=(n=1 h=0.5);
axis2 order=(0 to 120 by 30) value=(font=swissl height=1) 
   label=(font=swissl h=1) major=(h=0.5) minor=(n=1 h=0.5);
proc gplot data=pctb;
   plot sugar*hour=percent / haxis=axis1 vaxis=axis2 noframe;
   where percent ne 1.0 and percent ne 10.0 and
         percent ne 25.0 and percent ne 75.0 and
         percent ne 90.0 and percent ne 99.0;
run;



/*---overlay plots of bottoms---*/
proc kde data=bayes out=bottomA gl=2.5 gu=4 bwm=2 method=srot ng=100;
   var bottomA;
data bottomA;
   set bottomA;
   bottom = bottomA;
   treatment = 'A';
   drop bottomA count;
proc kde data=bayes out=bottomB gl=0 gu=4 bwm=2 method=srot ng=100;
   var bottomB;
data bottomB;
   set bottomB;
   bottom = bottomB;
   treatment = 'B';
   drop bottomB count;
data bottom;
   set bottomA bottomB;
   if (density >= 0);
run;

quit;
goptions reset=all device=xcolor
   hsize   = 6 in
   vsize   = 4 in
   htext   = 3.0 pct
   htitle  = 3.5 pct
   vorigin = 0 in
   horigin = 0 in
   ftext=swiss lfactor=1;
symbol1 i=j w=5 c=black l=1;
symbol2 i=j w=5 c=black l=2;
footnote font=italic height=1.5 "(a)";
axis1 order=(0 to 4 by 1) value=(font=swissl height=1) 
   label=(font=swissl h=1) major=(h=0.5) minor=none;
axis2 order=(0 to 2.5 by 0.5) value=(font=swissl height=1) 
   label=(font=swissl h=1) major=(h=1) minor=none;
proc gplot data=bottom;
   plot density*bottom=treatment / haxis=axis1 vaxis=axis2 noframe;
run;


/*---overlay plots of t85---*/
proc kde data=bayes out=t85A bwm=3 ng=100;
   var t85A;
data t85A;
   set t85A;
   t85 = t85A;
   treatment = 'A';
   drop t85A count;
proc kde data=bayes out=t85B bwm=3 ng=100;
   var t85B;
   where t85B ne .;
data t85B;
   set t85B;
   t85 = t85B;
   treatment = 'B';
   drop t85B count;
data t85;
   set t85A t85B;
   if (density >= 0);
run;

quit;
goptions reset=all device=xcolor
   hsize   = 6 in
   vsize   = 4 in
   htext   = 3.0 pct
   htitle  = 3.5 pct
   vorigin = 0 in
   horigin = 0 in
   ftext=swiss lfactor=1;
symbol1 i=j w=5 c=black l=1;
symbol2 i=j w=5 c=black l=2;
footnote font=italic height=1.5 "(b)";
axis1 order=(5 to 8 by 1) value=(font=swissl height=1) 
   label=(font=swissl h=1) major=(h=0.5) minor=none;
axis2 order=(0 to 1.6 by 0.4) value=(font=swissl height=1) 
   label=(font=swissl h=1) major=(h=1) minor=none;
proc gplot data=t85;
   plot density*t85=treatment / haxis=axis1 vaxis=axis2 noframe;
run;