Bayesian Multinomial Model for Ordinal Data

data icecream;
   input y1-y5 brand;
   if brand = 1 then brand1 = 1;
      else brand1 = 0;
   if brand = 2 then brand2 = 1;
      else brand2 = 0;
   keep y1-y5 brand1 brand2;
   datalines;
   70 71 151 30 46 1                            
   20 36 130 74 70 2 
   50 55 140 52 50 3 
   ;

ods graphics on;
proc mcmc data=icecream nbi=10000 nmc=25000 thin=10 seed=1181
   propcov=quanew monitor=(beta1 beta2 or12 or13 or23);   
   array theta[4]; 
   array gamma[4]; 
      
   parms theta1-theta4 beta1 beta2;  
   prior beta: ~ normal(0,var=1000);
   prior theta1 ~ normal(0,var=100); 
   prior theta2 ~ normal(0,var=100,lower=theta1);
   prior theta3 ~ normal(0,var=100,lower=theta2);
   prior theta4 ~ normal(0,var=100,lower=theta3);
  
   mu = beta1*brand1 + beta2*brand2;
   do j = 1 to 4;
      gamma[j] = logistic(theta[j] + mu);
   end;
   pi1 = gamma1;
   pi2 = gamma2 - gamma1;
   pi3 = gamma3 - gamma2;
   pi4 = gamma4 - gamma3;
   pi5 = 1 - sum(of pi1-pi4);

   llike = logmpdfmultinom(of y1-y5, of pi1-pi5);
   model dgeneral(llike);

   beginnodata;;
      or12 = exp(beta1-beta2);
      or13 = exp(beta1);
      or23 = exp(beta2);
   endnodata;;         
run;
ods graphics off;