model {
  for(i in 1:N) {
    imal[i] ~ dbern(pmal[i])
    logit(pmal[i]) <- theta[1] + theta[2]*icvs[i] + theta[3]*icdr[i] +
      theta[4]*icpr[i] + theta[5]*itri[i]
    iblk[i] ~ dbern(pblk[i])
    logit(pblk[i]) <- theta[6] + theta[7]*icvs[i] + theta[8]*ictr[i] +
      theta[9]*icpr[i] + theta[10]*irev[i] + theta[11]*itri[i]
    icjs[i] ~ dbern(pcjs[i])
    logit(pcjs[i]) <- theta[12] + theta[13]*icvm[i] +
      theta[14]*icdr[i] + theta[15]*icpr[i] + theta[16]*itri[i]
    ipri[i] ~ dbern(ppri[i])
    logit(ppri[i]) <- theta[17] + theta[18]*icvm[i] +
      theta[19]*irev[i]
    idet[i] ~ dbern(pdet[i])
    logit(pdet[i]) <- theta[20] + theta[21]*icvs[i] +
      theta[22]*icvm[i] + theta[23]*ictr[i] + theta[24]*icdr[i] +
      theta[25]*icpr[i] + theta[26]*itri[i]
    Y[i] ~ dbern(p[i])
    logit(p[i]) <- beta[county[i], 1] +
      beta[county[i], 2]*imal[i] +
      beta[county[i], 3]*iblk[i] +
      beta[county[i], 4]*icvs[i] +
      beta[county[i], 5]*icvm[i] +
      beta[county[i], 6]*ictr[i] +
      beta[county[i], 7]*icdr[i] +
      beta[county[i], 8]*icpr[i] +
      beta[county[i], 9]*icjs[i] +
      beta[county[i], 10]*ipri[i] +
      beta[county[i], 11]*idet[i] +
      beta[county[i], 12]*irev[i] +
      beta[county[i], 13]*itri[i]
      e.y[i] <- Y[i] - p[i]
  }
  for (j in 1:M) {
    for (k in 1:13) {
      beta[j, k] <- xi[k]*b[j, k]
      beta.hat[j, k] <- xi[k]*b.hat[j, k]
      e.b[j, k] <- beta[j, k] - beta.hat[j, k]
    }
    b[j, 1:13] ~ dmnorm(b.hat[j, 1:13], gamma[1:13, 1:13])
    b.hat[j, 1] <- eta.raw[1] + eta.raw[2]*cidx[j] +
      eta.raw[3]*cunr[j] + eta.raw[4]*cblp[j] + eta.raw[5]*ccon[j] +
      eta.raw[6]*csth[j] + eta.raw[7]*cgui[j]
    b.hat[j, 2] <- eta.raw[8] + eta.raw[9]*cidx[j] +
      eta.raw[10]*cunr[j] + eta.raw[11]*cblp[j] + 
      eta.raw[12]*ccon[j] + eta.raw[13]*csth[j]
    b.hat[j, 3] <- eta.raw[14] + eta.raw[15]*cidx[j] +
      eta.raw[16]*cunr[j] + eta.raw[17]*cblp[j] + 
      eta.raw[18]*ccon[j] + eta.raw[19]*csth[j]
    b.hat[j, 4] <- eta.raw[20] + eta.raw[21]*cidx[j] +
      eta.raw[22]*cunr[j] + eta.raw[23]*cblp[j] + 
      eta.raw[24]*ccon[j] + eta.raw[25]*csth[j] + eta.raw[26]*cgui[j]
    b.hat[j, 5] <- eta.raw[27] + eta.raw[28]*cidx[j] +
      eta.raw[29]*cunr[j] + eta.raw[30]*cblp[j] + 
      eta.raw[31]*ccon[j] + eta.raw[32]*csth[j] + eta.raw[33]*cgui[j]
    b.hat[j, 6] <- eta.raw[34] + eta.raw[35]*cidx[j] +
      eta.raw[36]*cunr[j] + eta.raw[37]*cblp[j] + 
      eta.raw[38]*ccon[j] + eta.raw[39]*csth[j] + eta.raw[40]*cgui[j]
    b.hat[j, 7] <- eta.raw[41] + eta.raw[42]*cidx[j] +
      eta.raw[43]*cunr[j] + eta.raw[44]*cblp[j] + 
      eta.raw[45]*ccon[j] + eta.raw[46]*csth[j] + eta.raw[47]*cgui[j]
    b.hat[j, 8] <- eta.raw[48] + eta.raw[49]*cidx[j] +
      eta.raw[50]*cunr[j] + eta.raw[51]*cblp[j] + 
      eta.raw[52]*ccon[j] + eta.raw[53]*csth[j] + eta.raw[54]*cgui[j]
    b.hat[j, 9] <- eta.raw[55] + eta.raw[56]*cidx[j] +
      eta.raw[57]*cunr[j] + eta.raw[58]*cblp[j] + 
      eta.raw[59]*ccon[j] + eta.raw[60]*csth[j]
    b.hat[j, 10] <- eta.raw[61] + eta.raw[62]*cidx[j] +
      eta.raw[63]*cunr[j] + eta.raw[64]*cblp[j] + 
      eta.raw[65]*ccon[j] + eta.raw[66]*csth[j] + eta.raw[67]*cgui[j]
    b.hat[j, 11] <- eta.raw[68] + eta.raw[69]*cidx[j] +
      eta.raw[70]*cunr[j] + eta.raw[71]*cblp[j] + 
      eta.raw[72]*ccon[j] + eta.raw[73]*csth[j]
    b.hat[j, 12] <- eta.raw[74] + eta.raw[75]*cidx[j] +
      eta.raw[76]*cunr[j] + eta.raw[77]*cblp[j] + 
      eta.raw[78]*ccon[j] + eta.raw[79]*csth[j]
    b.hat[j, 13] <- eta.raw[80] + eta.raw[81]*cidx[j] +
      eta.raw[82]*cunr[j] + eta.raw[83]*cblp[j] + 
      eta.raw[84]*ccon[j] + eta.raw[85]*csth[j] + eta.raw[86]*cgui[j]
  }
  for (k in 1:13) {
    xi[k] ~ dnorm (0, .01)
  }
# Priors for county effects:
  for (i in 1:7) {
    eta.raw[i] ~ dnorm(0.0, 0.01)
    eta[i] <- xi[1]*eta.raw[i]
  }
# Priors for individual effects:
  for (i in 1:2) {
    eta.raw[6*i+2] ~ dnorm(0.0, 0.01)
  }
  for (i in 3:8) {
    eta.raw[7*i-1] ~ dnorm(0.0, 0.01)
  }
  eta.raw[61] ~ dnorm(0.0, 0.01)
  for (i in 11:13) {
    eta.raw[6*i+2] ~ dnorm(0.0, 0.01)
  }
# Priors for interactions:
  for (i in 3:7) {
    for (j in 1:2) {
      eta.raw[6*j+i] ~ dnorm(0.0, 0.1)
    }
  eta.raw[53+i] ~ dnorm(0.0, 0.1)
  for (j in 11:12) {
    eta.raw[6*j+i] ~ dnorm(0.0, 0.1)
    }
  }
  for (i in 1:6) {
    for (j in 3:7) {
      eta.raw[7*j+i-1] ~ dnorm(0.0, 0.1)
    }
    eta.raw[61+i] ~ dnorm(0.0, 0.1)
    eta.raw[80+i] ~ dnorm(0.0, 0.1)
  }
  for (i in 8:13) {
    eta[i] <- xi[2]*eta.raw[i]
  }
  for (i in 14:19) {
    eta[i] <- xi[3]*eta.raw[i]
  }
  for (i in 20:26) {
    eta[i] <- xi[4]*eta.raw[i]
  }
  for (i in 27:33) {
    eta[i] <- xi[5]*eta.raw[i]
  }
  for (i in 34:40) {
    eta[i] <- xi[6]*eta.raw[i]
  }
  for (i in 41:47) {
    eta[i] <- xi[7]*eta.raw[i]
  }
  for (i in 48:54) {
    eta[i] <- xi[8]*eta.raw[i]
  }
  for (i in 55:60) {
    eta[i] <- xi[9]*eta.raw[i]
  }
  for (i in 61:67) {
    eta[i] <- xi[10]*eta.raw[i]
  }
  for (i in 68:73) {
    eta[i] <- xi[11]*eta.raw[i]
  }
  for (i in 74:79) {
    eta[i] <- xi[12]*eta.raw[i]
  }
  for (i in 80:86) {
    eta[i] <- xi[13]*eta.raw[i]
  }
# Hyper-priors:
  gamma[1:13, 1:13] ~ dwish(R[1:13, 1:13 ], 13)
# Priors for missing data distributions
  for(i in 1:26) {
    theta[i] ~ dnorm(0.0, 0.1)
  }
}