CODEWITHSUNDEEP
rbayesianjagswinbugs
Valerie Esposito
Valerie Esposito

Reputation: 1

How to convert WINBUGS code into JAGS/R code

I am trying to rerun the WinBugs code used in the paper 'A Bayesian Hierarchical Model for Risk Assessment of Methymercury'. I would like to be able to implement the same code but into JAGS (through the R Studio packages "rjags" and "coda". The WinBugs code is as follows

#Model formulation: Inverse Gaussian
model
[
for (i in 1:row)
[
y[i] <- b[i]
p.y[i] <- .1/pow(((bu[i]-b[i])/1.64),2)
p2.y[i] ~dgamma(p.y[i],.1)
y[i] ~ dnorm(mu[i], p2.y[i])
mu[i] ~ dnorm(gamma[study[i]], p.tau)
bmdso[i] <- 1/mu[i]
]
for (i in 1:nstudy) [
gamma[i] ~ dnorm(m.gam,p.gamma)
bmds[i] <- 1 / gamma[i]
]
m.gam ~ dnorm(0, .0001)
bmd <- 1 / m.gam
p.tau ~ dgamma(.001,.001)
p.gamma ~ dgamma(.001,.001)
sigo <- 1 / sqrt(p.tau)
sigs <- 1 / sqrt(p.gamma)
]
#Model formulation: Lognormal
model
[
for (i in 1:row)
[
y[i] <- log(1/b[i])
dc[i] <- .1/pow((log(y[i])-log(1/bu[i]))/1.64,2)
d[i] ~dgamma(dc[i],.1)
y[i] ~ dnorm(mu[i], d[i])
mu[i] ~ dnorm(gamma[study[i]], p.tau)
bmdso[i] <- exp(mu[i])
for (i in 1:nstudy)
[
gamma[i] ~ dnorm(m.gam,p.gamma)
bmds[i] <- exp(gamma[i])
]
m.gam ~ dnorm(3, 0.0001)
bmd <- exp( m.gam )
p.tau ~ dgamma(.001,.001)
p.gamma ~ dgamma(.001,.001)
sigo <- 1 / sqrt(p.tau)
sigs <- 1 / sqrt(p.gamma)
]
#Data input
list(b=c(
0.0081 , 0.0083 , 0.0122 , 0.0483 , 0.0091 , 0.01 ,
0.0829 , 0.0814 , 0.0786 , 0.1215 , 0.0777 ,
0.0502 , 0.0566 , 0.0353 , 0.0657 , 0.0368 ),
bu=c(
0.0402 , 0.0434 , 0.0447 , 0.0585 , 0.0426 , 0.0447 ,
0.1764 , 0.1596 , 0.165 , 0.2275 , 0.1653 ,
0.0829 , 0.0976 , 0.0677 , 0.0998 , 0.0697 ),
study=c(1,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3),
row=16, nstudy=3)

The current JAGS/R code that I have is this

# Install packages if not already installed
install.packages("rjags")
install.packages("coda")

# Load required packages
library(rjags)
library(coda)

# Data
data <- list(
  b = c(0.0081, 0.0083, 0.0122, 0.0483, 0.0091, 0.01, 0.0829, 0.0814, 0.0786, 0.1215, 0.0777, 0.0502, 0.0566, 0.0353, 0.0657, 0.0368),
  bu = c(0.0402, 0.0434, 0.0447, 0.0585, 0.0426, 0.0447, 0.1764, 0.1596, 0.165, 0.2275, 0.1653, 0.0829, 0.0976, 0.0677, 0.0998, 0.0697),
  study = c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3),
  row = 16,
  nstudy = 3
)

# JAGS model
model_string <- "
model {
  for (i in 1:row) {
    y[i] <- log(1 / b[i])
    dc[i] <- 0.1 / pow((log(y[i]) - log(1 / bu[i])) / 1.64, 2)
    d[i] ~ dgamma(dc[i], 0.1)
    y[i] ~ dnorm(mu[i], d[i])
    mu[i] ~ dnorm(gamma[study[i]], p.tau)
    bmdso[i] <- exp(mu[i])
  }

  for (i in 1:nstudy) {
    gamma[i] ~ dnorm(m.gam, p.gamma)
    bmds[i] <- exp(gamma[i])
    # No redundant definition of y[i] here
  }

  m.gam ~ dnorm(3, 0.0001)
  bmd <- exp(m.gam)
  p.tau ~ dgamma(0.001, 0.001)
  p.gamma ~ dgamma(0.001, 0.001)

  sigo <- 1 / sqrt(p.tau)
  sigs <- 1 / sqrt(p.gamma)
}
"

# Initial values
inits <- list(
  m.gam = 3,
  gamma = rnorm(data$nstudy, 3, 0.0001),
  p.tau = 0.001,
  p.gamma = 0.001
)

# Parameters to monitor
parameters <- c("bmds", "bmd", "sigo", "sigs", "mu", "gamma", "m.gam")

# Run JAGS
jags <- jags.model(file = model_string, data = data, inits = inits, n.chains = 3, n.adapt = 5000, quiet = FALSE)

# Check convergence and summarize results
summary(jags)

But it is giving me the following error:

**Error in jags.model(file = textConnection(model_string), data = data,  : 
  RUNTIME ERROR:
Compilation error on line 7.
Attempt to redefine node y[1]**

I don't know how to proceed from here is there anyone that could help me out please?

Upvotes: 0

Views: 109

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226577

The main issue, from section A4 of the JAGS manual, is that:

JAGS allows data transformations, but the syntax is different from BUGS. BUGS allows you to put a stochastic node twice on the left hand side of a relation, ...
You must put data transformations in a separate block of relations preceded by the keyword data ...

In the example below I've moved the definitions of both y[i] and dc[i] into the data block (I'm not sure this is right: I had to move the dc block to avoid an error about a cyclic path ... you might have to think about the original model more deeply to decide if this is the right way to "fix" the problem. All I can say is that I succeeded in running the model this way.)

# JAGS model
model_string <- "

data {
  for (i in 1:row) {
     y[i] <- log(1 / b[i])
     dc[i] <- 0.1 / pow((log(y[i]) - log(1 / bu[i])) / 1.64, 2)
  }
}

model {

  for (i in 1:row) {
    d[i] ~ dgamma(dc[i], 0.1)
    y[i] ~ dnorm(mu[i], d[i])
    mu[i] ~ dnorm(gamma[study[i]], p.tau)
    bmdso[i] <- exp(mu[i])
  }

  for (i in 1:nstudy) {
    gamma[i] ~ dnorm(m.gam, p.gamma)
    bmds[i] <- exp(gamma[i])
    # No redundant definition of y[i] here
  }

  m.gam ~ dnorm(3, 0.0001)
  bmd <- exp(m.gam)
  p.tau ~ dgamma(0.001, 0.001)
  p.gamma ~ dgamma(0.001, 0.001)

  sigo <- 1 / sqrt(p.tau)
  sigs <- 1 / sqrt(p.gamma)
}
}

# Initial values
inits <- list(
  m.gam = 3,
  gamma = rnorm(data$nstudy, 3, 0.0001),
  p.tau = 0.001,
  p.gamma = 0.001
)
"

# Parameters to monitor
parameters <- c("bmds", "bmd", "sigo", "sigs", "mu", "gamma", "m.gam")

I personally find the R2jags interface a little bit easier to use, so the code below is slightly modified.

fn <- tempfile()
writeLines(model_string, fn)
jags_fit <- R2jags::jags(model.file = fn, data = data,
                         inits = replicate(3, inits, simplify = FALSE),
                         n.chains = 3,
                         jags.seed = 101,
                         n.iter = 20000,
                         n.burnin = 10000, quiet = FALSE,
                         parameters.to.save = parameters)

# Check convergence and summarize results
gelman.diag(as.mcmc(jags_fit))
broom.mixed::tidy(jags_fit, conf.int = TRUE)

lattice::xyplot(as.mcmc(jags_fit), layout = c(5,6))

There are definitely still some problems with the chains, but hopefully this gets you a little farther.

Upvotes: 2

Related Questions