neon_dynamic_disagg_multisp_sbc.txt

model {

  # Priors
  for (i in 1:4) {
    mu[i] ~ dnorm(0, 1)
    ranks_mu[i] = mu[i] > mu_true[i]
    zeros[i] = 0
  }
  
  Tau_spec[1:4, 1:4] ~ dwish(R[1:4, 1:4], 10)
  for ( i in 1:K_exp) {
    eps_spec[i, 1:4] ~ dmnorm(zeros, Tau_spec)
  }
  
  Tau_site[1:4, 1:4] ~ dwish(R[1:4, 1:4], 10)
  for (i in 1:nsite) {
    eps_site[i, 1:4] ~ dmnorm(zeros, Tau_site)
  }
  
  for (i in 1:4) {
    for (j in 1:4) {
      ranks_Tau_spec[i, j] = Tau_spec[i, j] > Tau_spec_true[i, j]
      ranks_Tau_site[i, j] = Tau_site[i, j] > Tau_site_true[i, j]
    }
    for (k in 1:K_exp) {
      ranks_eps_spec[k, i] = eps_spec[k, i] > eps_spec_true[k, i]
    }
    for (ii in 1:nsite) {
      ranks_eps_site[ii, i] = eps_site[ii, i] > eps_site_true[ii, i]
    }
  }

  for (i in 1:nsite){
    for (k in 1:K_exp) {
   
      logit_psi1[i, k]  = mu[1] + eps_site[i, 1] + eps_spec[k, 1]
      psi1[i, k]  = ilogit(logit_psi1[i, k])
      
      logit_phi[i, k]   = mu[2] + eps_site[i, 2] + eps_spec[k, 2]
      phi[i, k]  = ilogit(logit_phi[i, k])
      
      logit_gamma[i, k] = mu[3] + eps_site[i, 3] + eps_spec[k, 3]
      gamma[i, k]  = ilogit(logit_gamma[i, k])
      
      log_lambda[i, k] = mu[4] + eps_site[i, 4] + eps_spec[k, 4]
      lambda[i, k] = exp(log_lambda[i, k])
    }
  }
  
  for (k in 1:K_exp) {
    Theta[k,1:K_para] ~ ddirch(a[k,1:K_para])
    for (kprime in 1:K_para) {
      ranks_Theta[k, kprime] = Theta[k, kprime] > Theta_true[k, kprime]
    }
  }

  for (i in 1:nsite) {
    for (k in 1:K_exp) {
      z[i, k, 1] ~ dbern(psi1[i, k])
      for (t in 2:nyear) {
        z[i, k, t]  ~ dbern(z[i,k,t-1]*phi[i, k] + 
                        (1 - z[i,k,t-1])*gamma[i, k]) 
      }
      for (t in 1:nyear) {
        zlam[i, k, t] = z[i, k, t] * lambda[i, k]
      }
    }
  }
  
  for (i in 1:nsite) {
    for (j in 1:nsurv) {
      for (t in 1:nyear) {
          L[i, j, t] ~ dpois(sum(zlam[i, 1:K_exp, t]))
      }
    }
  }

  for (l in 1:Ltot) {
    pi[l, 1:K_exp] = zlam[site[l], 1:K_exp, year[l]] / sum(zlam[site[l], 1:K_exp, year[l]])
    k[l] ~ dcat(pi[l, 1:K_exp])
    y[l] ~ dcat(Theta[k[l], 1:K_para])
  }
  
  # Derived parameters
  for (i in 1:nsite) {
    for (k in 1:K_exp) {
      psi[i, k, 1] <- psi1[i, k]
      for (t in 2:nyear) {
        psi[i, k,t]        = psi[i, k, t-1]*phi[i, k] + 
                          (1 - psi[i, k, t-1])*gamma[i, k] 
        log_growth[i, k, t]     = log(psi[i, k, t]) - log(psi[i, k, t-1]) 
        turnover[i, k, t-1] = (1 - psi[i, k, t-1])*gamma[i, k]/psi[i, k, t] 
      }
    }
  }
}