model {            
            for (i in 1:N2) {
                    y[i] ~ dnorm(mean.y[i],tau.y)
                    mean.y[i]<-b0+b1*x1.mar[i]+b2*x2.mar[i]+b3*x3.mar[i]
                    logit(mean.x1[i])<-mu0.x1+mu1.x1*x2.mar[i]
                    x1.mar[i]~dbern(mean.x1[i])
                    mean.x2[i]<-mu0.x2
                    mean.x3[i]<-mu0.x3+mu1.x3*x1.mar[i]
                    x2.mar[i]~dnorm(mean.x2[i],tau.x2)
                    x3.mar[i]~dnorm(mean.x3[i],tau.x3)                            
            } 
            tau.y<-pow(sigma.y,-2)
            sigma.y~dunif(.00001,100)
            tau.x2<-pow(sigma.x2,-2)
            sigma.x2~dunif(.00001,100)
            tau.x3<-pow(sigma.x3,-2)
            sigma.x3~dunif(.00001,100)
            b0 ~ dnorm(0,.01)
            b1 ~ dnorm(0,.01)
            b2 ~ dnorm(0,.01)
            b3 ~ dnorm(0,.01)
            mu0.x1~dnorm(0,.01)
            mu1.x1~dnorm(0,.01)
            mu0.x2~dnorm(0,.01)
            mu0.x3~dnorm(0,.01)
            mu1.x3~dnorm(0,.01)
        }