Computing 2 independent models at once leads to wrong results in WinBUGS

Question

I wanted to compare how two different (independent) models perform under winbugs, So I created code like:

# likelihood
for (i in 1:sites) {
    M[i] ~ dpois(lambda) # model 1
    M2[i] ~ dpois(lambda2) # model 2

    for (j in 1:sample) {
        obs[i, j] <- sum(y[i, j,])

        # model 1
        y[i, j, 1] ~ dbin(p, M[i])
        obs[i, j] ~ dbin(p_komb, M[i])

        # model 2
        obs[i, j] ~ dbin(p2, M2[i])
    }
} 

But it leads to wrong results in model 2! But if I separated the model 2:

# likelihood
for (i in 1:sites) {
    M2[i] ~ dpois(lambda2)
    for (j in 1:sample) {
        obs[i, j] <- sum(y[i, j,])
        obs[i, j] ~ dbin(p2, M2[i])
    }
}

then it performs well and returns good result! I can't really find out why the first example should return different result. The two models combined in the example seems to be completely independent. WinBUGS bug? Any ideas?

(It doesn't matter if I separate the models in two for loops - result is the same.)

I have WinBUGS 1.4.3 (August 2007) with immortality patch installed. Here is complete reproducible code for R and package R2WinBUGS (with data generation):

################## 
# data generation
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)

############################################
# two models together - simple model with complex model together 


sink("tmp_bugs/model.txt")
cat("

model {

# likelihood
for (i in 1:sites) {
    M[i] ~ dpois(lambda) # model 1
    M2[i] ~ dpois(lambda2) # model 2

    for (j in 1:sample) {
        obs[i, j] <- sum(y[i, j,])

        # model 1
        y[i, j, 1] ~ dbin(p, M[i])
        obs[i, j] ~ dbin(p_komb, M[i])

        # model 2 (simple)
        obs[i, j] ~ dbin(p2, M2[i])
    }
}

# derived parameters
Mtot <- sum(M[])
M2tot <- sum(M2[])

# priors

tau <- 1/(4 * 4)

p <- 1/(1+exp(-logit_p))
logit_p ~ dnorm(0, tau)

p2 <- 1/(1+exp(-logit_p2))
logit_p2 ~ dnorm(0, tau)

p_komb <- p + (1 - p) * p


lambda ~ dunif(0, 100)
lambda2 ~ dunif(0, 100)

}


")
sink()


win.data = list(y = y, sample = sample, sites = sites)

inits = function () { list(
    M = apply(y_sample_total, 1, max), 
    M2 = apply(y_sample_total, 1, max), 
    logit_p = rnorm(1, 0, 4),
    logit_p2 = rnorm(1, 0, 4)
) }

params = c("M", "M2", "p", "p2", "Mtot", "M2tot", "lambda", "lambda2")
#params = c("M2", "p2", "M2tot", "lambda2")

ni <- 2500
nt <- 16
nb <- 1000
nc <- 3

date()
out1 <- 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()


#############################
# simpler model itself

sink("tmp_bugs/model.txt")
cat("

model {

# likelihood
for (i in 1:sites) {
    M2[i] ~ dpois(lambda2)
    for (j in 1:sample) {
        obs[i, j] <- sum(y[i, j,])
        obs[i, j] ~ dbin(p2, M2[i])
    }
}

# derived parameters
M2tot <- sum(M2[])

# priors

tau <- 1/(4 * 4)

p2 <- 1/(1+exp(-logit_p2))
logit_p2 ~ dnorm(0, tau)

lambda2 ~ dunif(0, 100)

}


")
sink()


win.data = list(y = y, sample = sample, sites = sites)

inits = function () { list(
#   M = apply(y_sample_total, 1, max), 
    M2 = apply(y_sample_total, 1, max), 
#   N = y_sample_total, 
#   after_removal = y[,,2],
#   logit_p = rnorm(1, 0, 4),
#   logit_q = rnorm(1, 0, 4), 
    logit_p2 = rnorm(1, 0, 4)
) }


#params = c("M", "M2", "p", "p2", "Mtot", "M2tot", "lambda", "lambda2")
params = c("M2", "p2", "M2tot", "lambda2")

ni <- 2500
nt <- 16
nb <- 1000
nc <- 3

date()
out2 <- 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()



#############################
# summary - results of simple model (M2, M2tot) differ depending on whether the 
# model was evaluated along with the more complex model or alone!!!
# Shouldn't be!!!

#print(out, dig = 3)

par(mfrow = c(2, 2))

out <- out1

hist(out$sims.list$M2tot, breaks = 100)
abline(v = out$mean$M2tot, col = "red", lwd = 2)
abline(v = sum(M), col = "green", lwd = 2)
#lines(quantile(out$sims.list$M2tot, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2)
lines(quantile(out$sims.list$M2tot, 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$lambda2)

out <- out2

hist(out$sims.list$M2tot, breaks = 100)
abline(v = out$mean$M2tot, col = "red", lwd = 2)
abline(v = sum(M), col = "green", lwd = 2)
#lines(quantile(out$sims.list$M2tot, c(0.025, 0.975)), rep(sum(par("usr")[3:4]*c(0.9,0.1)), 2), lwd = 2)
lines(quantile(out$sims.list$M2tot, 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$lambda2)

Answer 1

Basically, you have two lines giving a distribution to obs[i, j], and WinBUGS isn't clever enough to figure out that you're trying to do two things at once. The simple solution would be to define distinct variables:

    obsa[i, j] <- sum(y[i, j,])
    obsb[i, j] <- sum(y[i, j,])

    # model 1
    obsa[i, j] ~ dbin(p_komb, M[i])

    # model 2
    obsb[i, j] ~ dbin(p2, M2[i])