# Solved – WinBUGS error with zero values in binomial distribution: value of order of binomial must be greater than zero

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:

Contents

## 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"), 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"), 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"), 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"), 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) ``

@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.

Rate this post