model {
  ## priors
  for(k in 2:ncat) {
    ## priors for model coefficients
    for(j in 1:2) {
      beta[j,k] ~ dnorm(0, 0.04)
    }
    ## hyperpriors for random effects
    sigmaind[k] ~ dt(0,1,2)I(0,)
    sigmafam[k] ~ dt(0,1,2)I(0,)
  }
  sigmaind[1] <- 0
  sigmafam[1] <- 0

  ## random effects 
  for(i in 1:nind) {rind[i, 1] <- 0}         # zero for baseline category
  for(f in 1:nfam) {rfam[f, 1] <- 0}         # zero for baseline category
  for(k in 2:ncat) {
    for(i in 1:nind) {
      rind[i,k]~dnorm(0, 1)                  # random individual effect
    }
    for(f in 1:nfam) {
      rfam[f,k]~dnorm(0, 1)                  # random family effect
    }
  }

  ## likelihood
  for(i in 1:n) {
    roost[i]~dcat(p[i,1:ncat])
    for(k in 1:ncat) {
      ## linear predictors including the random factors
      z[i,k]<-beta[1,k] + beta[2,k]*temp[i] + sigmafam[k]*rfam[fam[i],k] + sigmaind[k]*rind[ind[i],k] 
      expz[i,k] <- exp(z[i,k])
      p[i,k] <- expz[i,k] / sum(expz[i,1:ncat])# logit link
    }
  }

  ## constrain coefficients of the baseline category to zero
  for(j in 1:2) {
    beta[j,1] <- 0
  }

  ## predict site selection probabilities for different temperatures 
  for(k in 1:ncat) {
    for(i in 1:nnew) {
      pnew[i,k] <- expznew[i,k] / sum(expznew[i,1:ncat])
      znew[i,k] <- beta[1,k] + beta[2,k]*newtemp[i] 
      expznew[i,k] <- exp(znew[i,k])
    }
  }
}