/**************************************************************************************************
Bayesian Ideas and Data Analysis
Chapter 11: Count Data
SAS code for the FMD example in the ZIP subsection
*****************************************************************************************************/

proc import out=FMDdata1
  datafile='G:\MyDocuments\BAYES\PoissonRegression\FMDdata\FMDData.xls'
  dbms=excel replace;
  getnames=yes;
  sheet="Sheet1";
run;
proc print data=FMDdata1;
run;

proc means data=FMDdata1;   
 var Cattle;
run;


data FMDdata;
  set FMDdata1;
  stdcattle=(Cattle-162667.88)/ 112105.31;
run;
proc print data=FMDdata;
run;



/* 1. MLEs for ordinary Poisson regression of the FMD 1998 data.  
      Response: yi = FMD count in province i
      Covariates:  EasternTurkey (1=east; 0=west);  standardized cattle population 
      AIC=231.2    BIC=237.8
*/
proc genmod data=FMDdata;
  model FMD1998=stdcattle EasternTurkey / dist=Poisson;
run;




/* 2. Bayesian ordinary Poisson regression of the FMD 1998 data.  
      Response: yi = FMD count in province i
      Covariates:  EasternTurkey (1=east; 0=west);  standardized cattle population  
      DIC=231.2.  Estimates similar to MLEs
*/
proc mcmc data=FMDdata nbi=15000 nmc=1000000 dic propcov=quanew;
parms beta1 0 beta2 0 beta3 0;
prior beta: ~ normal(0,var=1000);
mu = exp(beta1 + beta2*stdcattle + beta3*EasternTurkey);
model FMD1998 ~ poisson(mu);
run;





/* 3. Bayesian ZIP regression of the FMD 1998 data.  
      Response: yi = FMD count in province i
      Covariates:  EasternTurkey (1=east; 0=west);  standardized cattle population 
      DIC=172.3
*/
proc mcmc data=FMDdata nbi=15000 nmc=1000000 dic propcov=quanew;
parms beta1 0 beta2 0 beta3 0 pi 0.5;
prior beta: ~ normal(0,var=1000);
mu = exp(beta1 + beta2*stdcattle + beta3*EasternTurkey);
prior pi ~ uniform(0,1);
llike=log(pi*(FMD1998 eq 0) + (1-pi)*pdf("poisson",FMD1998,mu));
model general(llike);
run;