Undefined real result in a zero-inflated negative binomial

Question

I have run a zero-inflated Poisson model in WinBUGS without problems, and now I am trying to run its equivalent negative binomial. However, I get an "undefined real result" trap message over and over. Here is the data set http://www.megafileupload.com/en/file/474809/data1.html

#Data handling
library(R2WinBUGS)
setwd('') # working directory
dget('data1')
data1 <- data1[order(data1$C), ]
    n = nrow(data1)
    n0 = length(which(data1$C == 0))

sink("ZINB.txt")
cat("
model {

# Pr(Y=0|x,z) # probability of zeros
K <- 10000
for(i in 1:n0){
zeros0[i] <- 0
zeros0[i] ~ dpois(mu0[i])
mu0[i] <- -log(p00[i] + (1-p00[i])*pow(overd[i], r)) + K

overd[i] <- r/(r + lambda0[i])
p00[i] <- min( 0.999999, max(.000001,p0[i]) ) # avoid stackoverflow
logit(p0[i]) <- a0 + b0.A*A[i] + b0.M*M[i] + b0.cat*cat[i] + b0.d*d[i]
log(lambda0[i])<- a + logh[i] + b.A*A[i] + b.M*M[i] + b.cebo*cebo[i] + b.descart*descart[i] + b.d*d[i]
}

# Pr(Y>0|x,z) # probability of positive counts
for(j in (n0+1):n){
zeros[j] <- 0
zeros[j] ~ dpois(mu[j])
mu[j] <- -log((1-p1[j])*fnb[j]) + K

lfnb[j] <-  r*log(overd[j]) + C[j]*log(1 - overd[j]) + loggam(C[j] + r) - loggam(r) - loggam(C[j] + 1) # - log(1 - pow(overd[j], r))
fnb[j] <- exp( lfnb[j] )
overd[j] <- r/(r + exp(lambda[j]))
p1[j] <- min( 0.999999, max(.000001,p[j]) )
logit(p[j]) <- a0 + b0.A*A[j] + b0.M*M[j] + b0.cat*cat[j] + b0.d*d[j] + b0.descart*descart[j]
log(lambda[j])<- a + logh[j] + b.A*A[j] + b.M*M[j] + b.cebo*cebo[j] + b.descart*descart[j] + b.d*d[j] + b.dcos*dcos[j]
pred[j] <- (1 - p1[j])*lambda[j] # predicted values (see Zuur et al 2009)
var[j] <- (1 - p1[j])*(lambda[j] + lambda[j]*lambda[j]/r) + lambda[j]*lambda[j]*(p1[j]*p1[j] + p1[j]) # variance (see Zuur et al 2009)
res[j] <- (C[j] - lambda[j]*(1 - p[j]))/sqrt( var[j] ) # Pearson residuals
disp[j] <- pred[j]/var[j] # dispersion
}

# Vague priors for model coefficients
r ~ dunif(0,10)
a0 ~ dnorm(0,0.001)
a ~ dnorm(0,0.001)
b.A ~ dnorm(0,0.001)
b0.A ~ dnorm(0,0.001)
b.M ~ dnorm(0,0.001)
b0.M ~ dnorm(0,0.001)
b.d ~ dnorm(0,0.001)
b0.d ~ dnorm(0,0.001)
b.descart ~ dnorm(.0,0.001)
b0.descart ~ dnorm(0,0.001)
b0.cat ~ dnorm(0,0.001)
b.dcos ~ dnorm(0,0.001)   
b.cebo ~ dnorm(0,0.001)
}
",fill=TRUE)
sink()

# Data passed to WinBUGS
win.data <- list( C = data1$C, A = as.numeric(data1$A), M = as.numeric(data1$M), d = as.numeric(scale(data1$dcol)), dcos = as.numeric(scale(data1$dcos)), cebo = as.numeric(data1$cebo), descart = as.numeric(data1$descart), cat = as.numeric(data1$cat), logh = log(data1$h), n = n, n0 = n0 )

# Initial values
inits <- function(){ list( r = runif(1, 1, 5), a0 = rnorm(1), a = rnorm(1), b.A = rnorm(1), b0.A = rnorm(1), b.M = rnorm(1), b0.M = rnorm(1), b.d = rnorm(1), b0.d = rnorm(1), b.descart = rnorm(1), b0.descart = rnorm(1), b0.cat = rnorm(1), b.dcos = rnorm(1), b.cebo = rnorm(1) ) }

# Parameters
params <- c( "r", "mu", "lambda", "a", "a0", "b.A", "b0.A", "b.M", "b0.M", "b.d", "b0.d", "b.descart", "b0.descart", "b0.cat", "b.dcos", "b.cebo", "res", "pred", "var" )

# MCMC settings
nc <- 3
ni <- 30000
nb <- 3000
nt <- 3

# WinBUGS execution
out <- bugs( data = win.data, inits = inits, parameters.to.save = params, model.file = "ZINB.txt", n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni, debug = TRUE )

print(out, dig = 2)

The trap arises from the Pr(Y > 0) part of the model since the Pr(Y = 0) runs without problems. I've tried different inits and priors to no avail.

Answer 1

Ok, I found a missing zero in one of the priors:

b.descart ~ dnorm(.0,0.001)  

Apparently:

1. BUGS doesn't take .0 as a 0.0, and
2. WinBUGS doesn't take .0 as a syntax error

Another point: the model neither runs with more than three covariates nor more than two categorical covariates. Any clue?

As an aside, the value I gave to the constant K makes the deviance and DIC to get wildly high values.