#The data
library(arm)
library(R2WinBUGS)

 ### Prob 12-89 Montgommery & Runger - Linear regression - full Bayesian Hierarchical
 test=matrix(scan("http://civil.colorado.edu/~balajir/CVEN5454/R-sessions/sess5/prob12-93-data.txt"),ncol=8,byrow=T)
 
N=dim(test)[1]
Ystar = log(test[,2])
x3star = log(test[,5])
Xdat = cbind(x3star,test[,6],test[,7])
Xdat[,2:3]=Xdat[,2:3]/100

zz = lm(Ystar ~ Xdat)

data=list(n=N,x=Xdat,y=Ystar)
#inits = function(){list(beta0=0,beta1=1,beta2=1,taue=1,betamean0=0,betamean1=1,betamean2=1,
#tau0=1,tau1=1,tau2=1)}

## the starting values are based on the lm estimates

inits = function(){list(beta0=zz$coef[1],beta1=zz$coef[2],beta2=zz$coef[3],beta3=zz$coef[4],taue=1,
betamean0=zz$coef[1],betamean1=zz$coef[2],betamean2=zz$coef[3],betamean3=zz$coef[4],
tau0=1,tau1=1,tau2=1,tau3=1)}

parameters = c("beta0","beta1","beta2","beta3","taue","betamean0","betamean1","betamean2","betamean3","tau0","tau1","tau2","tau3")
#parameters = c("beta0","beta1","beta2","taue","betamean0","betamean1","betamean2","tau0","tau1","tau2")
 
 
 linreg.sim = bugs(data,inits=inits,parameters,"linear-reg-examp-prob12-89-full-bayes-bug.txt",n.iter=10000)
 attach.bugs(linreg.sim)

## You can see the posterior distribution of the regression paramters
## Also compute the median to obtain a point estimate
#You can look at the priors distributions as well
   
 hist(beta0)
 hist(beta1)
 hist(beta2)
 hist(beta3)
 hist(taue)