Failed JAGS: Node inconsistant with parents

Question

Hi I am new to jags and I encountered the dreaded 'node inconsistent with parents' problem. I was trying to run my data with a Bayesian integrated SDM following the instructions on Fidino's blog post in R, with an adaptation of changing the occupancy model to an abundance model.

I got the following error:

All chains have finished Note: the model did not require adaptation Error in runjags.readin(directory = startinfo$directory, silent.jags = silent.jags, : Simulation numbers 2, 3, 4 appear to have crashed; check the output in failed.jags() for clues Note: Either one or more simulation(s) failed, or there was an error in processing the results. You may be able to retrieve any successful simulations using: results.jags("C:/Users/jc817859/AppData/Local/Temp/RtmpInTLt3/runjagsfiles4ea43b04d2b", recover.chains=TRUE) See the help file for that function for possible options. To remove failed simulation folders use cleanup.jags() - this will be run automatically when the runjags package is unloaded

When look into the failed.jags():

chain 2 . Initializing model Error in node z[274969] Node inconsistent with parents

chain 3 . Initializing model Error in node ones[131] Invalid parent values Deleting model

---

The model (which is 95% the same as Fidino's):

model{
  # Bayesian version of the Koshkina (2017) model.
  #
  # The latent-state model
  for(pixel in 1:npixel){
    # latent state linear predictor
    #
    # x_s  = covariates for latent state
    # beta = latent state model regression coefficients
    # cell_area = log area of grid cell 
    #
    log(lambda[pixel]) <-inprod(x_s[pixel,], beta) + cell_area
    # Species presence in a gridcell as a Bernoulli trial
    n[pixel] ~ dpois(lambda[pixel])
    # presence only thinning prob linear predictor
    #
    # h_s = covariates for thinning probability
    # cc  = presence-only data model regression coefficients
    #
    logit(b[pixel]) <-  inprod(h_s[pixel,] , cc)
  }
  # The presence only data model.
  #
  # This part of the model just uses the
  #  what we have calculated above (lambda
  #  and b). The denominator of this likelihood
  #  is actually a scalar so we can calculate it
  #  outside of a for loop. Let's do that first.
  #
  # The presence_only data model denominator, which
  #  is the thinned poisson process across the
  #  whole region (divided by the total number of 
  #  data points because it has to be 
  #  evaluated for each data point).
  po_denominator <- inprod(lambda[1:npixel], b[1:npixel] ) / m
  #
  # Loop through each presence-only datapoint
  #  using Bernoulli one's trick. The numerator
  #  is just the thinned poisson process for
  #  the ith data point.
  for(po in 1:m){
    ones[po] ~ dbern(
      exp(
        log(lambda[po_pixel[po]]*b[po_pixel[po]]) -
          log(po_denominator)
      ) / CONSTANT
    )
  } 
#
# Detection / non-detection data model
for(site in 1:nsite){
  # detection/non-detection data model linear predictor
  #
  #  v = detection covariates for the entire region B
  #  a = det/non-det data model regression coefficients
  #
  logit(rho[site]) <-inprod(v[pa_pixel[site], ],a)
  # The number of detections for site is a binomial
  #  process with Pr(rho[site]) | z = 1 with
  #  W sampling occasions per site.
  y[site] ~ dbin(
    z[pa_pixel[site]] * rho[site],
    W
  )
}
#
# Priors for latent state model
for(latent in 1:nlatent){
  beta[latent] ~ dnorm(0, 0.01)
}
# Priors for presence-only data model
for(po in 1:npar_po){
  cc[po] ~ dlogis(0, 1)
}
# Priors for det/non-det data model
for(pa in 1:npar_pa){
  a[pa] ~ dlogis(0, 1)
}
# Derived parameter, the number of cells occupied
nsum <- sum(n)
}

m1 <- run.jags(
  model = "./JAG models/integrated_model.R", 
  data = input, 
  n.chains = 10, 
  inits = inits, 
  monitor = c("beta", "cc", "a", "zsum"), 
  adapt = 1000, 
  burnin = 2000, 
  sample = 30000,
  thin = 2,
  modules = "glm",
  method = 'parallel'
)

---

The input data: https://shorturl.at/pvGL3 (It's a large RDS file)

I looked into the old posts, and it seems that the priors, inits, and out of distribution boundary can all be the causes of the error.

For the inits, I have manually set all the inits to 1, the same as the original model. I also tried to increase the CONSTANT, but it didn't work.

I also tried the functions and distribution boundary outside of the model, but couldn't find anything weird. Maybe it's because I don't have the mathmatical brain?

Any hints and modifications are highly appreciated!