Error in checkForRemoteErrors(val) : 3 nodes produced errors; first error: Error in node Y[3,5,3] Node inconsistent with parents

Question

I am working with the R-package R2jags. After running the code I attach below, R produced the error message: "Error in checkForRemoteErrors(val) : 3 nodes produced errors; first error: Error in node Y[3,5,3] Node inconsistent with parents".

I tried to solve it. However, the error message persists.

loading libraries
library(reshape)
library (vegan)
library(plotrix)
library(ggplot2)
library(tidyr)
library(tidyverse)
library(tibble)
library(snow)
library(rjags)
library(dclone)
library(abind)
library(R2jags)

# BUGS model
modelFilename = "dmsom.txt"

cat("
model {
# priors
omega ~ dunif(0,1)
psiMean ~ dunif(0,1)
for (t in 1:T) {
    pMean[t] ~ dunif(0,1)
}
for (t in 1:(T-1)) {
    phiMean[t] ~ dunif(0,1)
    gamMean[t] ~ dunif(0,1)
}
lpsiMean <- log(psiMean) - log(1-psiMean)
for (t in 1:T) {
    lpMean[t] <- log(pMean[t]) - log(1-pMean[t])
}
for (t in 1:(T-1)) {
    lphiMean[t] <- log(phiMean[t]) - log(1-phiMean[t])
    lgamMean[t] <- log(gamMean[t]) - log(1-gamMean[t])
}
lpsiSD ~ dunif(0,10)
lpsiPrec <- pow(lpsiSD,-2)
lpSD ~ dunif(0,10)
lphiSD ~ dunif(0,10)
lgamSD ~ dunif(0,10)
for (t in 1:T) {
    lpPrec[t] <- pow(lpSD,-2)
}
for (t in 1:(T-1)) {
    lphiPrec[t] <- pow(lphiSD,-2)
    lgamPrec[t] <- pow(lgamSD,-2)
}

# likelihood
for (i in 1:M) {
    # initial occupancy state at t=1
    w[i] ~ dbern(omega)
    b0[i] ~ dnorm(lpsiMean, lpsiPrec)T(-12,12)
    lp[i,1] ~ dnorm(lpMean[1], lpPrec[1])T(-12,12)
    p[i,1] <- 1/(1+exp(-lp[i,1]))
    for (j in 1:n.site) {
        lpsi[i,j,1] <- b0[i]
        psi[i,j,1] <- 1/(1 + exp(-lpsi[i,j,1]))
        mu.z[i,j,1] <- w[i] * psi[i,j,1]
        Z[i,j,1] ~ dbern(mu.z[i,j,1])
        mu.y[i,j,1] <- p[i,1]*Z[i,j,1]
        Y[i,j,1] ~ dbin(mu.y[i,j,1], K_tot[j,1])
    }
    # model of changes in occupancy state for t=2, ..., T
    for (t in 1:(T-1)) {
        lp[i,t+1] ~ dnorm(lpMean[t+1], lpPrec[t+1])T(-12,12)
        p[i,t+1] <- 1/(1+exp(-lp[i,t+1]))
        c0[i,t] ~ dnorm(lgamMean[t], lgamPrec[t])T(-12,12)
        d0[i,t] ~ dnorm(lphiMean[t], lphiPrec[t])T(-12,12)
        for (j in 1:n.site) {
            lgam[i,j,t] <- c0[i,t]
            gam[i,j,t] <- 1/(1+exp(-lgam[i,j,t]))
            lphi[i,j,t] <- d0[i,t]
            phi[i,j,t] <- 1/(1+exp(-lphi[i,j,t]))
            psi[i,j,t+1] <- phi[i,j,t]*psi[i,j,t] + gam[i,j,t]*(1-psi[i,j,t])
            mu.z[i,j,t+1] <- w[i] * (phi[i,j,t]*Z[i,j,t] + gam[i,j,t]*(1-Z[i,j,t]))
            Z[i,j,t+1] ~ dbern(mu.z[i,j,t+1])
            mu.y[i,j,t+1] <- p[i,t+1]*Z[i,j,t+1]
            Y[i,j,t+1] ~ dbin(mu.y[i,j,t+1], K_tot[j,t+1])
        }
    }
}

# Derive total species richness for the metacommunity
N_tot <- sum(w[])

# Derive yearly species richness
for (i in 1:M) {
    for (t in 1:T) {
        tmp[i,t] <- sum(Z[i,,t])
        tmp2[i,t] <- ifelse(tmp[i,t]==0,0,1)
    }
}
for (t in 1:T) {
    N[t] <- sum(tmp2[,t])
}

# Derive WPI
for (i in 1:nspeciesWPI) {
    for (t in 1:T) {
        psi_ratio[i,t] <- psi[id_speciesWPI[i],1,t]/psi[id_speciesWPI[i],1,1] 
        log_psi_ratio[i,t] <- log(psi_ratio[i,t])
    }
}
for (t in 1:T) {
    WPI[t] <-  exp((1/(nspeciesWPI))*sum(log_psi_ratio[1:nspeciesWPI,t]))
}

# rate of change in WPI
for (t in 1:(T-1)) {
    lambda_WPI[t] <- WPI[t+1]/WPI[t]    
}

}
", fill=TRUE, file=modelFilename)


# getting value of number of traps
nsites <- dim(Yaug_tot)[2]

# getting value of number of primary periods (years)
T <- dim(Yaug_tot)[3]

# latent states
w <- c(rep(1,length(allnames3)),rep(NA,nzeros))

# creating value of parameters monitored
parameters <- c("omega","psiMean","pMean","phiMean","gamMean",
                "lpsiSD","lpSD","lphiSD","lgamSD",
                "N","N_tot","WPI","lambda_WPI")

# creating data frame of each species at each site with WPI 
bugs.data <- list(M=dim(Yaug_tot)[1],n.site=nsites,K_tot=K_tot,Y=Yaug_tot,T=T,
                  nspeciesWPI=nspeciesWPI,id_speciesWPI=id_speciesWPI)

# creating function of initial values
inits <- function() { list(omega=runif(1),Z=(Yaug_tot>0)*1,w=w,
                           psiMean=runif(1),pMean=runif(T),phiMean=runif(T-1),gamMean=runif(T-1),
                           lpsiSD=runif(1,0,4),lpSD=runif(1,0,4),lphiSD=runif(1,0,4),lgamSD=runif(1,0,4))}

# # Check dimensions
# dim(Yaug_tot)
# dim(K_tot)
# Yaug_tot[3, 5, 3]
# cat(readLines("dmsom.txt"), sep = "\n")
# any(K_tot == 0)
# inits()

#markov chain settings
#creating summaries of the parameter posterior distribution calculated from three markov chains 
#initialised with random starting values run 30,000 times, after a 10,000 burn in and thinned every 10 draws
#pre burn in value
n.adapt <- 5000 

#burn in value
n.update <- 10000

#post-burnin value
n.iter <- 30000     

#thin value
thin <- 10

#markov chain value
chains<-3

# create data frame of cluster of markov chains
cl <- makeCluster(chains, type = "SOCK")

#record start time model being run
start.time = Sys.time()

out <- jags.parfit(cl, data = bugs.data, 
                   params = parameters, 
                   model = "dmsom.txt",
                   inits = inits,
                   n.adapt = n.adapt,    
                   n.update = n.update,
                   n.iter = n.iter,
                   thin = thin, n.chains = chains)

end.time = Sys.time()
elapsed.time = difftime(end.time, start.time, units='mins')
cat(paste(paste('Posterior computed in ', elapsed.time, sep=''), ' minutes\n', sep=''))
stopCluster(cl)