Bayesian Multivariate Prior for Multiple Linear Regression

/*-----------------------------------------------------------------
  Example: Fitting a Bayesian Linear Regression model with a 
           Multivariate Prior     
  Requires: SAS/STAT
  Version: 9.2                                                                       
  ------------------------------------------------------------------*/

data brainweight;
   input body gestation litter @@;
   log_brain = log(brain); 
   log_body = log(body);
   log_gestation = log(gestation); 
   datalines;
   17.5 3.5 26 1
   3.5 0.93 34 4.6
   3.15 0.15 46 3
   1.14 0.049 51 1.5
   1.37 0.064 46 1.5
   22 2.1 135 1
   12.8 1.2 90 1.2
   9.9 0.7 135 1
   54 7.7 139 1
   73 3.7 180 1
   114 9.1 140 1
   109 7.7 140 1
   7.8 0.22 145 2
   84.6 6 175 1
   107 8.7 165 1.1
   183 21 180 1
   179 32 180 1
   67 4.6 195 1
   65.5 5.8 168 1
   102 5.5 210 1
   343 37 270 1
   360 45 230 1
   406 140 265 1
   1300 65 270 1
   12 3.7 120 4
   9.6 2.2 31 5
   13.3 2.9 41 2.5
   6.23 0.33 38 3
   1.89 0.052 40 3.1
   40 20 128 2.9
   45 25 128 4
   0.68 0.027 23 3.7
   0.63 0.026 23 5
   0.52 0.017 24 5
   0.69 0.024 24 5
   0.67 0.036 21 4.6
   1.12 0.13 16 6.3
   1.04 0.065 21 4
   0.72 0.05 23 7.3
   2.38 0.34 21 8
   0.45 0.024 19 5
   1.18 0.15 27 5.6
   37 14 112 1.2
   24 6.6 113 1
   4.28 0.97 67 2.6
   76 30 123 3
   20.3 2.8 104 1.3
   9.9 0.78 98 1.2
   5.25 0.43 110 2
   23 5 132 5.5
   1600 160 360 1
   537 56 270 1
   70.2 8.5 63 4
   48 6 52 4
   37.3 3.8 63 3.7
   28.5 3.2 65 4
   400 250 219 2.3
   500 250 240 1.8
   41.6 5.3 63 3.5
   31.2 2 77 1.1
   53 6 60 2.2
   28.4 2.5 63 4
   75 12 60 2.5
   157 46 92 2.5
   260 180 108 3
   302 210 104 3
   355 250 254 1
   363 100 343 1
   442 110 240 1
   550 400 310 1
   20.5 3.8 225 2.4
   712 480 330 1
   250 230 390 1
   185 150 120 4
   180 190 115 8
   590 1400 240 1
   260 150 205 1
   225 93 330 1
   198 45 300 1.1
   124 16 183 1.1
   223 80 240 1
   219 89 218 1
   435 200 255 1
   365 120 235 1
   383 120 246 1.1
   288 110 225 1
   480 560 255 1
   334 250 255 1
   456 520 280 1
   93 13 120 1
   200 39 180 1
   210 66 158 1.2
   125 49 150 2.4
   106 30 151 2
   ;

proc template;
   define statgraph brainweight;
      begingraph;

     layout lattice / rows=2 columns=2;

    layout overlay;
         histogram brain;
      endlayout;

    layout overlay;
         histogram body;
      endlayout;

    layout overlay;
         histogram gestation;
      endlayout;

    layout overlay;
         histogram litter;
      endlayout;

     endlayout;
     endgraph;
end;
run;

ods graphics on;
proc sgrender data=brainweight template=brainweight;
run;
ods graphics off;

ods graphics on;
proc mcmc data=brainweight nbi=5000 nmc=25000 thin=5 seed=1181
   propcov=quanew monitor=(_parms_ ExpBeta2 ExpBeta3);
   
   array mu_pr[4] mu_pr0-mu_pr3;
   array beta[4] beta0-beta3;
   array data[4] 1 log_body log_gestation litter;

   begincnst;
      call zeromatrix(mu_pr);     
      rc = logmpdfset('lp',16,0,16,0,3,2.25,0,-2,0,1);
   endcnst;

   parms beta: 0;                                 
   parms sig2 1;
   
   ExpBeta2=exp(beta2);
   ExpBeta3=exp(beta3);

   lmvn = logmpdfnormal(of beta0-beta3, of mu_pr0-mu_pr3, 'lp');
   prior beta: ~ general(lmvn);
   prior sig2 ~ igamma(shape = 2.001, scale = 1.001);

   call mult(beta, data, mu);
   model log_brain ~ n(mu, var = sig2);
run;

data _null_;
   rc = logmpdffree();
   put rc=;
run;
ods graphics off;