0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ It could be the use of exp () that is causing problems. For example you use "fnb" but then take log of it - this is not necessary, just use your "lfnb" variable instead. $\endgroup$ – probabilityislogic Aug 2 '14 at 7:04
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.