3
$\begingroup$

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)
Improve this question
$\endgroup$
1
  • $\begingroup$ Isn't this somehow related to joint probability? But I would suppose that WinBUGS recognizes that the two models are not "jointed" in any way and will not "joint" them... :) $\endgroup$
    – Tomas
    Commented Feb 15, 2012 at 7:49

1 Answer 1

Reset to default
1
$\begingroup$

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])
Improve this answer
$\endgroup$
3
  • $\begingroup$ you are right!!! It is stupid behaviour, because obs[] is only derived from data, so it is a fixed number... but he probably takes it as a parameter which then "joins" the two models together! Thanks, your fix works, the model 2 is then computed the same way as if separate! $\endgroup$
    – Tomas
    Commented Feb 15, 2012 at 14:44
  • $\begingroup$ It also works if you just add obs2[i, j] <- obs[i, j] and use obs in model 1 and obs2 in model 2. This is interesting, that pure copy of data makes him think that the data is different. I'd expect the contrary. So suprisingly no need to run the sum twice. $\endgroup$
    – Tomas
    Commented Feb 15, 2012 at 14:48
  • $\begingroup$ I also tried to supply the obs as input data, already sum-ed in R. If I supply only one copy of this data, according to expectation, the behaviour is buggy until I do the same as above, i.e. clone the data by obs2[i, j] <- obs[i, j] and use obs in model 1 and obs2 in model 2. $\endgroup$
    – Tomas
    Commented Feb 15, 2012 at 14:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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