I have a dataset consisting of two groups of participants (25 participants in one group and 26 in the other), with each participating in two conditions, and resulting a single binary outcome for each condition. I'm trying to analyze this model using JAGS, by combining some of Kruschke's suggestions in his book but I'm getting strange results. This is the model (B represents between, and W within. There are two levels for each. BxW represents the interaction) The problem is that I get posterior predictives distributions that are more extreme than the actual data (proportion of 1's in each group/condition; See plots below). This is really my first time building a Bayesian model for data analysis so any help would be appreciated. Isaac.

Code: 
model{
for (i in 1:nTotal){
clicks[i]~dbern(theta[i])
theta[i]<-ilogit(a0+aB[grp[i]]+aW[cond[i]]+aS[subj[i]]+aBxW[grp[i],cond[i]])
}

a0 ~ dnorm(0, 1/2^2) 
aB[1]~dnorm( 0 , 1/sigmaB^2 )
aB[2]~dnorm( 0 , 1/sigmaB^2 )
sigmaB ~ dgamma(1.64,0.32) ###all sigma's of effects are set to this prior following Kruschke's suggesting of using a gamma distribution with a mode of 2 and sd of 4, thus allowing for the maximum possible deflections  (logit(0.999)-log(0.001)
aW[1]~dnorm( 0 , 1/sigmaW^2 )
aW[2]~dnorm( 0 , 1/sigmaW^2 )
sigmaW ~ dgamma(1.64,0.32)
for (s in 1:nSubj) {aS[s] ~ dnorm( 0 , 1/sigmaS^2 ) }
sigmaS ~ dgamma(1.64,0.32)
for ( i in 1:2) { 
  for ( j in 1:2) {
    aBxW[i,j] ~ dnorm( 0 , 1/sigmaBxW^2 )
  } 
}
sigmaBxW ~ dgamma(1.64,0.32)

####Sum to zero

for (s in 1:nSubj) {
  for (j in 1:2) {
   mSxW[s,j]<-(a0+aB[grpPerSubj[s]]+aW[j]+aBxW[grpPerSubj[s],j]+aS[s]) 

   ####grpPerSubj is a vector of thegroup of each subject. It's length is   #51(25+26). It is needed because the length of the dataset if 51*2 (two within-#subject conditions)

    for (i in 1:2) {
      mSxBxW[s,i,j] <- ifelse(grpPerSubj[s]==i,
                            (a0 + aB[grpPerSubj[s]] + aW[j] +   aBxW[grpPerSubj[s],j] + aS[s]),0)
  } } }

for ( i in 1:2) {
  for ( j in 1:2) { 
    mBxW[i,j] <- sum(mSxBxW[1:nSubj,i,j])/nSubjGrpCond[i,j]
    bBxW[i,j] <- mBxW[i,j] - mB[i] - mW[j] + m
} }

for (s in 1:nSubj) {
  mS[s] <- mean(mSxW[s,1:2] )
  bS[s] <- mS[s] - mB[grpPerSubj[s]]
}
for (i in 1:2){
  mB[i] <- mean(mBxW[i,1:2])
  bB[i] <- mB[i] - m
}
for (j in 1:2){
  mW[j] <- mean(mBxW[1:2,j])
  bW[j] <- mW[j] - m
}
m <- mean(mBxW[1:2,1:2] )
b0 <- m

}