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