Translating WinBUGS model to JAGS (using R)

Question

I am trying to implement the following model written for WinBUGS in JAGS:

model {
  for (i in 1:N) {                         
    wtp[i] ~ dweib(r[G[i]], mu[i])I(lower[i], upper[i])
    mu[i] <- exp(beta[G[i]])
    G[i] ~ dcat(P[])       
  }                                       
  P[1] ~ dunif(0.01, 0.99)
  P[2] <- 1 - P[1]
  r[1] ~ dunif(1, 10)
  r[2] ~ dunif(0.1, 10)
  beta[1] ~ dunif(0, 1000)      
  beta[2] ~ dunif(-1000, 0)                
  weibmed[1] <- pow(log(2) * exp(-beta[1]), 1 / r[1]) 
  weibmed[2] <- pow(log(2) * exp(-beta[2]), 1 / r[2])
  weibmed[3] <- pow(log(1 / (1 - 0.5 + P[1])) * exp(-beta[2]), 1 / r[2])
  weibmean[1] <- pow(exp(-beta[1]), 1 / r[1]) * exp(loggam((1 + r[1]) / r[1]))
  weibmean[2] <- pow(exp(-beta[2]), 1 / r[2]) * exp(loggam((1 + r[2]) / r[2]))
  weibmean[3] <- P[1] * weibmean[1] + P[2] * weibmean[2]
}

I thought it would be straightforward to get it going in JAGS with:

library(rjags)

txt <- 'model {
  for (i in 1:N) {                         
    wtp[i] ~ dweib(r[G[i]], mu[i])T(lower[i], upper[i])
    mu[i] <- exp(beta[G[i]])
    G[i] ~ dcat(P[])       
  }                                       
  P[1] ~ dunif(0.01, 0.99)
  P[2] <- 1 - P[1]
  r[1] ~ dunif(1, 10)
  r[2] ~ dunif(0.1, 10)
  beta[1] ~ dunif(0, 1000)      
  beta[2] ~ dunif(-1000, 0)                
  weibmed[1] <- pow(log(2) * exp(-beta[1]), 1 / r[1]) 
  weibmed[2] <- pow(log(2) * exp(-beta[2]), 1 / r[2])
  weibmed[3] <- pow(log(1 / (1 - 0.5 + P[1])) * exp(-beta[2]), 1 / r[2])
  weibmean[1] <- pow(exp(-beta[1]), 1 / r[1]) * exp(loggam((1 + r[1]) / r[1]))
  weibmean[2] <- pow(exp(-beta[2]), 1 / r[2]) * exp(loggam((1 + r[2]) / r[2]))
  weibmean[3] <- P[1] * weibmean[1] + P[2] * weibmean[2]
}'

set.seed(3.14159)
dat <- list(N = 1000, lower = rep(0, 1000), upper = runif(1000, 5, 200000))
ini <- list(P = c(0.4, NA), r = c(8.2, 1.2), beta = c(3.8, -6.5))

mod <- jags.model(
  file = textConnection(txt), 
  data = dat,
  inits = c(ini, .RNG.name = 'base::Mersenne-Twister', .RNG.seed = 314159),
  n.chains = 1,
  n.adapt  = 100
)

sam.jags <- coda.samples(
  model = mod,
  variable.names = c('P', 'r', 'beta', 'weibmed', 'weibmean'),
  n.iter = 400,
  n.thin = 1
)

by just replacing the I() with a T(). This produces the coda.samples() error:

Error: Error in node weibmed[3]
Invalid parent values

If I ignore monitoring of weibmed and weibmean then coda.samples() works but the parameter estimates:

                Mean          SD    Naive SE Time-series SE
P[1]       0.4840704   0.2769491  0.01384746     0.01384746
P[2]       0.5159296   0.2769491  0.01384746     0.01384746
beta[1]  509.3614647 295.0860473 14.75430237    14.75430237
beta[2] -487.5362940 285.4126899 14.27063449    14.27063449
r[1]       5.2054730   2.6330434  0.13165217     0.13165217
r[2]       5.0478143   2.9480476  0.14740238     0.14740238

are not comparable to those that I get when using WinBUGS:

library(R2WinBUGS)

sam.bugs <- bugs(
  model.file = 'model.bug',
  data = dat,
  inits = list(ini),
  parameters.to.save = c('P', 'r', 'beta'), #, 'weibmed', 'weibmean'),
  n.chains = 1,
  n.burnin = 100,
  n.iter = 500,
  n.thin = 1,
  debug = F,
  DIC = F,
  bugs.seed = 314159
)

 Inference for Bugs model at "3mixout2.bug", fit using WinBUGS,
 1 chains, each with 500 iterations (first 100 discarded)
 n.sims = 400 iterations saved
        mean  sd 2.5%  25%  50%  75% 97.5%
P[1]     0.4 0.0  0.3  0.4  0.4  0.4   0.4
P[2]     0.6 0.0  0.6  0.6  0.6  0.6   0.7
r[1]     7.2 0.8  6.2  6.6  6.9  7.7   9.3
r[2]     1.5 0.0  1.4  1.4  1.5  1.5   1.6
beta[1]  5.5 0.5  4.6  5.0  5.4  5.8   6.6
beta[2] -7.2 0.2 -7.5 -7.3 -7.2 -7.1  -6.9

Any thoughts or suggestions?

cuttlefish44 · Answer 1 · 2016-11-25T07:54:56.920

1

JAGS manual says;

pow(x, z) || Power function || Real || If x < 0 then z is integer

When 0.5 < P[1], log(1 / (1 - 0.5 + P[1])) * exp(-beta[2]) < 0, so weibmed[3], pow(negative, non_integer), becomes NaN. So far as I know, coda.samples() doesn't permit monitored variables to take NaN.

If you use P[1] ~ dunif(0.01, 0.5) or except weibmed[3] from the monitered variable list by changing its name, such as weibmed3 <- pow(log(1 / (1 - 0.5 + P[1])) * exp(-beta[2]), 1 / r[2]), your code runs.

edited Nov 25 '16 at 07:54

answered Nov 25 '16 at 07:39

cuttlefish44

6,586
2
17
34

Thanks for pointing that out. Excepting both `weibmed` and `weibmean` from the monitored variable list do you have any idea why WinBUGS and JAGS would produce such different estimates for `P`, `r`, and `beta`? – johnson-shuffle Nov 25 '16 at 18:22
@feats-by-jake; Sorry, I've thrown BUGS environment with my old PC away. Could you show me the result ? (and whether it changes by excepting both `weibmed` and `weibmean` from `parameters.to.save` ?) – cuttlefish44 Nov 26 '16 at 01:41
I excepted both `weibmed` and `webmean` in both cases and added output from both `coda.samples()` and `bugs()`. The WinBUGS output does not change if `weibmed` and `webmean` are added back. – johnson-shuffle Nov 26 '16 at 01:58

score 1 · Answer 2 · answered Nov 28 '16 at 02:14

It appears that this is a case where dinterval should be used instead of T():

txt2 <- '
data {
  x <- rep(1, N)
}
model {
  for (i in 1:N) {     
    x[i] ~ dinterval(wtp[i], B[i, ])                 
    wtp[i] ~ dweib(r[G[i]], mu[i])
    mu[i] <- exp(beta[G[i]])
    G[i] ~ dcat(P[])
  }                                     
  P[1] ~ dunif(0, 1)
  P[2] <- 1 - P[1]
  r[1] ~ dunif(0, 10)
  r[2] ~ dunif(0, 10)
  beta[1] ~ dunif(1, 1000)      
  beta[2] ~ dunif(-1000, 0)                
  weibmed[1] <- pow(log(2) * exp(-beta[1]), 1 / r[1]) 
  weibmed[2] <- pow(log(2) * exp(-beta[2]), 1 / r[2])
  # weibmed[3] <- pow(log(1 / (1 - 0.5 + P[1])) * exp(-beta[2]), 1 / r[2])
  weibmean[1] <- pow(exp(-beta[1]), 1 / r[1]) * exp(loggam((1 + r[1]) / r[1]))
  weibmean[2] <- pow(exp(-beta[2]), 1 / r[2]) * exp(loggam((1 + r[2]) / r[2]))
  # weibmean[3] <- P[1] * weibmean[1] + P[2] * weibmean[2]
}'

mod <- jags.model(
  file = textConnection(txt2), 
  data = list(N = dat[[1]], B = cbind(dat$lower, dat$upper)),
  inits = c(ini, .RNG.name = 'base::Mersenne-Twister', .RNG.seed = 314159),
  n.chains = 1,
  n.adapt  = 100
)

I am sorry that my reply to you is delayed. The estimations by this code are are not comparable to those both JAGS and BUGS makes. For example, `P[1]` is about 0 and `P[2]` is about 1, are they correct ? — cuttlefish44, Nov 28 '16 at 05:21
You're right! I think this is because my reproducible example is probably not well designed. I've been using more carefully simulated data and JAGS and BUGS are producing similar results. — johnson-shuffle, Nov 29 '16 at 16:46

Translating WinBUGS model to JAGS (using R)

2 Answers2