It seems that WinBUGS has problems **if it has only zero draws from one binomial distribution**:

**1. case – simple model**

`for (i in 1:sites) { N[i] ~ dpois(lambda) for (j in 1:sample) { y[i, j] ~ dbin(p, N[i]) } } `

If all values `y[i,]`

are zero, then the following error message appears:

value of order of binomial y[53,1] must be greater than zero

**EDIT:** note that `N[i]`

was not zero in generated data upon failure.

**2. case – advanced removal model**

`for (i in 1:sites) { M[i] ~ dpois(lambda) for (j in 1:sample) { y[i, j, 1] ~ dbin(p, M[i]) after_removal[i, j] <- M[i] - y[i, j, 1] y[i, j, 2] ~ dbin(p, after_removal[i, j]) } } `

Here it seems that if some `y[i, j, 2]`

is zero then the error message (same as in case 1) appears!

Note that in this case, I found a workaround: I reparametrized the model to different one which **I hope** is equivalent:

`for (i in 1:sites) { M[i] ~ dpois(lambda) for (j in 1:sample) { y[i, j, 1] ~ dbin(p, M[i]) obs[i, j] <- y[i, j, 1] + y[i, j, 2] obs[i, j] ~ dbin(p_komb, M[i]) } } p_komb <- p + (1 - p) * p `

… anyway an explanation, clean solution and solution for case 1 wanted:

## Why this happens? Is it a WinBUGS bug? How can it be overriden?

I have WinBUGS 1.4.3 (August 2007) with immortality patch installed.

Below is complete **reproducible code for R and package R2WinBUGS (with data generation)**:

**1. case**

`require(vcd) sites <- 120 # 60 mean_N <- 16 N <- rpois(sites, mean_N) p <- 0.4 sample <- 3 # 3 y = matrix(nrow = sites, ncol = sample) for (i in 1:sites) { y[i,] = rbinom(sample, N[i], p) } y[20,] = 0 sink("tmp_bugs/model.txt") cat(" model { # likelihood for (i in 1:sites) { N[i] ~ dpois(lambda) for (j in 1:sample) { y[i, j] ~ dbin(p, N[i]) } } # derived parameters Ntot <- sum(N[]) # priors p <- 1/(1+exp(-logit_p)) tau <- 1/(4 * 4) logit_p ~ dnorm(0, tau) lambda ~ dunif(0, 100) } ") sink() win.data = list(y = y, sample = sample, sites = sites) #win.data = list(y = y) #inits = function () { list(N = apply(y, 1, max), p = mean(apply(y, 1, mean)/apply(y, 1, max))) } inits = function () { list(N = apply(y, 1, max), logit_p = rnorm(1, 0, 4)) } params = c("N", "p", "Ntot", "lambda") ni <- 500 nt <- 12 nb <- 200 nc <- 3 out <- bugs(win.data, inits, params, "model.txt", nc, ni, nb, nt, bugs.directory = bugs.dir, working.directory = paste(getwd(), "/tmp_bugs/", sep = ""), debug = TRUE ) print(out, dig = 3) par(mfrow = c(2, 2)) hist(out$sims.list$p, breaks = 100) abline(v = out$mean$p, col = "red", lwd = 2) abline(v = p, col = "green", lwd = 2) #lines(quantile(out$sims.list$p, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2) lines(quantile(out$sims.list$p, c(0.025, 0.975)), rep(par("usr")[3], 2), lwd = 4) legend("topleft", c("estimated", "real", "95% cred. int."), col = c("red", "green", "black"), lty = 1, box.lty = 0, cex = 0.7) #hist(mean(y) / out$sims.list$p, breaks = 1000, xlim = c(0, 50)) hist(out$sims.list$Ntot, breaks = 100) abline(v = out$mean$Ntot, col = "red", lwd = 2) abline(v = sum(N), col = "green", lwd = 2) #lines(quantile(out$sims.list$Ntot, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2) lines(quantile(out$sims.list$Ntot, c(0.025, 0.975)), rep(par("usr")[3], 2), lwd = 4) legend("topright", c("estimated total N", "real total N", "95% credible int."), col = c("red", "green", "black"), lty = 1, box.lty = 0, cex = 0.7) `

**2. case**

(one of the `y[i, j, 2]`

will likely be zero; if not, it can be assigned directly)

`require(vcd) sites <- 120 # 60 mean_M <- 16 M <- rpois(sites, mean_M) p <- 0.4 #0.64 sample <- 2 # 3 y = rep(NA, sites * sample * 2) dim(y) = c(sites, sample, 2) for (i in 1:sites) { # obs[i,] = rbinom(sample, M[i], p) for (j in 1:sample) { y[i,j,1] = rbinom(1, M[i], p) y[i,j,2] = rbinom(1, M[i] - y[i,j,1], p) } } y_sample_total = apply(y, c(1, 2), sum) obs = y_sample_total ############################################ sink("tmp_bugs/model.txt") cat(" model { # likelihood for (i in 1:sites) { M[i] ~ dpois(lambda) for (j in 1:sample) { y[i, j, 1] ~ dbin(p, M[i]) after_removal[i, j] <- M[i] - y[i, j, 1] y[i, j, 2] ~ dbin(p, after_removal[i, j]) } } # derived parameters Mtot <- sum(M[]) # priors tau <- 1/(4 * 4) p <- 1/(1+exp(-logit_p)) logit_p ~ dnorm(0, tau) lambda ~ dunif(0, 100) } ") sink() win.data = list(y = y, #obs = y_sample_total, sample = sample, sites = sites) inits = function () { list( M = apply(y_sample_total, 1, max), logit_p = rnorm(1, 0, 4) ) } #params = c("M", "M2", "p", "p2", "Mtot", "M2tot", "lambda", "lambda2") params = c("M", "p", "Mtot", "lambda") ni <- 2500 nt <- 16 nb <- 1000 nc <- 3 date() out <- bugs(win.data, inits, params, "model.txt", nc, ni, nb, nt, bugs.directory = "C:/Program Files/WinBUGS14/", working.directory = paste(getwd(), "/tmp_bugs/", sep = ""), debug = TRUE ) date() ############################# print(out, dig = 3) par(mfrow = c(2, 2)) hist(out$sims.list$p, breaks = 100) abline(v = out$mean$p, col = "red", lwd = 2) abline(v = p, col = "green", lwd = 2) #lines(quantile(out$sims.list$p, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2) lines(quantile(out$sims.list$p, c(0.025, 0.975)), rep(par("usr")[3], 2), lwd = 4) legend("topleft", c("estimated", "real", "95% cred. int."), col = c("red", "green", "black"), lty = 1, box.lty = 0, cex = 0.7) #hist(mean(y) / out$sims.list$p, breaks = 1000, xlim = c(0, 50)) hist(out$sims.list$Mtot, breaks = 100) abline(v = out$mean$Mtot, col = "red", lwd = 2) abline(v = sum(M), col = "green", lwd = 2) #lines(quantile(out$sims.list$Mtot, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2) lines(quantile(out$sims.list$Mtot, c(0.025, 0.975)), rep(par("usr")[3], 2), lwd = 4) legend("topright", c("estimated total M", "real total M", "95% credible int."), col = c("red", "green", "black"), lty = 1, box.lty = 0, cex = 0.7) hist(out$sims.list$lambda) `

#### Best Answer

**@Stephane appears to be correct:**

For your 1st case, the problem is when N[i]=0 !

This problem is easy to reproduce and it's straightforward to test his conclusion. E.g., compile this model

`model { for (i in 1:sites) { N[i] ~ dpois(lambda) for (j in 1:sample) { y[i, j] ~ dbin(p, N[i]) } } } `

with a simple set of data such as

`list( y = structure( .Data = c(0,0,0,0,0,0,0,0), .Dim = c(4,2) ), lambda = 1, sites = 4, sample = 2, p = 0.5 ) `

(Note that all data values are zero.) When you try to generate initial values for the chains, your error is raised:

value of order of binomial y[1,1] must be greater than zero

Now change lambda from a small value to a large one (for which N[i]==0 will essentially never happen):

`... lambda = 100, ... `

The model runs fine.

You can also vary the data `c(0,0,0,0,0,0,0,0)`

to confirm that their values do not affect the error. *Ergo*, the error *must* be due to either to an invalid value of N[i] or p. Clearly p (= 0.5) is valid, QED.

### Similar Posts:

- Solved – representing quantile like quartile in form of normal distribution curve
- Solved – How to incorporate uncertainty of actual historical data into forecast prediction intervals
- Solved – What graphs / plots are best suited to visualise percentiles
- Solved – Runtime error on binomial model/node error on exponential model
- Solved – How does boxplot in R calculate quantiles