model{
   #Likelihood:
   for(i in 1:N){
     Y[i]~dnorm(mu[i],tau)
     mu[i]<-int+inprod(x.c[i,],beta[])
   }
   tau~dgamma(0.1,0.1)
   int~dnorm(0,0.001)

   #SSVS prior:
   for(j in 1:p){
      delta[j]~dbern(0.05)
      alpha[j]~dnorm(0,0.1)
      beta[j]<-delta[j]*alpha[j]
   }

   #keep track of individual model probabilities, e.g. 
   #model[2,1,2] indicates the model is x1+x3:
   for(j1 in 1:2){for(j2 in 1:2){for(j3 in 1:2){
     model[j1,j2,j3]<-equals(delta[1],j1-1)*equals(delta[2],j2-1)*equals(delta[3],j3-1)
   }}}

   #center the data:
   for(i in 1:N){for(j in 1:3){
      x.c[i,j]<- (x[i,j]-mean(x[,j]))/sd(x[,j])
   }}
}

#Initial values:
list(delta=c(1,1,1),alpha=c(0,0,0),tau=1,int=15)

#Stacks data:
list(p = 3,   N = 21,  
	Y = c(42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15),
	x = structure(.Data = c(80, 27, 89,
				80, 27, 88,
				75, 25, 90,
				62, 24, 87,
				62, 22, 87,
				62, 23, 87,
				62, 24, 93,
				62, 24, 93,
				58, 23, 87,
				58, 18, 80,
				58, 18, 89,
				58, 17, 88,
				58, 18, 82,
				58, 19, 93,
				50, 18, 89,
				50, 18, 86,
				50, 19, 72,
				50, 19, 79,
				50, 20, 80,
				56, 20, 82,
				70, 20, 91),
 .Dim = c(21, 3))