D1.education1[i]<-equals(education[i],1)

logit(p[i,j]) <- (b1cons + b1age*x1[i] + b1education[1]*D1.education1[i] +  b1education[2]*D1.education2[i]+ b1education[3]*D1.education3[i]+b1residence*D1.residence[i] + b1working*D1.working[i] +b1breast*D1.breast[i]+ group[i] +b2cons + b2age*x1[j] + b2education[1]*D2.education1[j] +  b2education[2]*D2.education2[j]+ b2education[3]*D2.education3[j]+b2residence*D2.residence[j] + b2working*D2.working[j] +b2breast*D2.breast[j]+ group[j])

    }
    }
     for(i in 1:N) {
    for(j in 1:M){

    Y_ppd[i,j] ~ dbern(p[i,j])
    }
    }
    # M groups
    for (i in 2:K) {
    for (j in 2:K){

        group1[i]~ dnorm(0,tau1)
        group2[j]~dnorm(0,tau2)
    }
    }
    b1education[4] <- 0.0
    b2education[4] <- 0.0
    # Priors
    b1cons ~ dnorm(0.0,1.0E-6)
    b1age ~ dnorm(0.0,1.0E-6)
    b1residence ~ dnorm(0.0,1.0E-6)
    b1working  ~ dnorm(0.0,1.0E-6)
    b1breast ~ dnorm(0.0,1.0E-6)
    for(l in 1:3){
    b1education[l]~dnorm(0.0,1.0E-6)
    }

    b2cons ~ dnorm(0.0,1.0E-6)
    b2age ~ dnorm(0.0,1.0E-6)
    b2residence ~ dnorm(0.0,1.0E-6)
    b2working  ~ dnorm(0.0,1.0E-6)
    b2breast ~ dnorm(0.0,1.0E-6)
    for(l in 1:3){
    b2education[l]~dnorm(0.0,1.0E-6)

    }
    # Hyperprior
    tau1 ~ dgamma(0.001,0.001)
    tau2 ~ dgamma(0.001,0.001)