model{
  
  # prior distributions
  alpha ~ dnorm(0.37, 1536.292)
  for(k in 1:K){
    b[k] ~ dlogis(0, 1)
  }
  
  lambda[1] ~ dnorm(1.01, 1536.292)
  lambda[2] ~ dnorm(1.01, 1536.292)
  
  ### CALCULATING STABLE STAGE DISTRIBUTION
  for(i in 1:G){  # looping over all groups
    
    # survival probability
    logit(p[i, 1]) <- b[1] + b[2] * year[i]
    g[i, 1] <- p[i, 1]  # survive to next stage class
    z[i, 1] <- 0  # remain in current stage class
    w[i, 1] <- 1  # proportion to stable stage
    
    for(s in 2:(S - 1)){  # loop through all stage classes
      logit(p[i, s]) <- b[1] + b[2] * year[i]
      g[i, s] <- p[i, s]
      z[i, s] <- 0
      w[i, s] <- g[i, s - 1] / (lambda[year[i] + 1] - z[i, s]) * w[i, s - 1]
    }
    
    logit(p[i, S]) <- b[1] + b[2] * year[i]
    g[i, S] <- 0
    z[i, S] <- p[i, S]
    w[i, S] <- g[i, S - 1] / (lambda[year[i] + 1] - z[i, S]) * w[i, S - 1]
    
    # stable stage distribution
    C[i, 1:S] <- w[i, 1:S] / sum(w[i, 1:S])
    
    ### EXPECTED COUNT IN EACH STAGE CLASS
    ey[i, 2] <- p[i, 1] * C[i, 1]
    ey[i, 3] <- p[i, 2] * C[i, 2]
    ey[i, 4] <- p[i, 3] * C[i, 3]
    ey[i, 5] <- p[i, 4] * C[i, 4]
    ey[i, 6] <- p[i, 5] * C[i, 5] + p[i, 6] * C[i, 6]
    ey[i, 1] <- alpha * (lambda[year[i] + 1] - sum(ey[i, 2:6]))
    
    ### LIKELIHOOD
    y[i, 1:S] ~ dmulti(ey[i, 1:S], N[i])
    
    ### MODEL CHECKING
    y_new[i, 1:S] ~ dmulti(ey[i, 1:S], N[i])
    
    for(n in 1:S){
    
      g_obs[i, n] <- y[i, n] *
        log(y[i, n] / (N[i] * ey[i, n] / sum(ey[i, 1:S])))
      
      g_new[i, n] <- y_new[i, n] *
        log(y_new[i, n] / (N[i] * ey[i, n] / sum(ey[i, 1:S])))
      
    }
    
  }
  
  cs_obs <- 2 * sum(g_obs)
  cs_new <- 2 * sum(g_new)
  bp <- step(cs_new - cs_obs)
  
}