model {

  # reconstructing the data from what is passed from R...
  # nn <- c(n) # lays n out as a sequence of columns
  # dd <- c(d) # lays d out at a sequence of columns

  for (i in 1:N) {
    for (k in 1:K) {
#     n[i,k] <- nn[i + (k-1)*N]
      n[k,i] <- nn[i + (k-1)*N]
    }
    for (j in 1:J) {
      d[i,j] <- dd[i + (j-1)*N]
    }
  }
  
  for (k in 1:K) {
#   ntot[k] <- sum(n[,k])
    ntot[k] <- sum(n[k,])
  }
  
  # Level 1
  for (k in 1:K) {
    for (i in 1:N) {
      for (j in 1:J) {
        kern[i,j,k] <- p[i,j,k]*d[i,j] + (1-p[i,j,k])*(1-d[i,j])
        logit(p[i,j,k]) <- alpha[k] + beta[j]
      }
#     pp[i,k] <- prod(kern[i,,k])
      pp[k,i] <- prod(kern[i,,k])
    }
#   n[1:N,k] ~ dmulti(pp[1:N,k],ntot[k])
    n[k,1:N] ~ dmulti(pp[k,1:N],ntot[k])
  }

  # Level 2
  for (k in 1:K) {
    alpha[k] <- 0     #    alpha[k] ~ dnorm(mu.a,tau.a)
  }
  for (j in 1:J) {    #  beta[1] <- 0
  beta[j] ~ dnorm(mu.b,tau.b)
  }

  # Level 3

  mu.b ~ dnorm(0,1.0E-5)
  tau.b ~ dgamma(0.001,0.001) # or dgamma(1,0.001)

  sig.b <- pow(tau.b,-0.5)

}