我正在尝试重新运行论文“甲基汞风险评估的贝叶斯分层模型”中使用的 WinBugs 代码。我希望能够实现相同的代码，但进入 JAGS（通过 R Studio 包“rjags”和“coda”）。WinBugs 代码如下
#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)
我当前的 JAGS/R 代码是这个
# 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)
但它给了我以下错误：
**Error in jags.model(file = textConnection(model_string), data = data, :
RUNTIME ERROR:
Compilation error on line 7.
Attempt to redefine node y[1]**
我不知道如何从这里开始，有人可以帮助我吗？
JAGS 手册 A4 部分的主要问题是：
JAGS 允许数据转换，但语法与 BUGS 不同。 BUGS 允许您在关系的左侧放置两次随机节点，...在下面的示例中，我已将您必须将数据转换放在单独的关系块中，前面带有关键字
data...
y[i]
和
dc[i]的定义移至
data块中（我不确定这是否正确：我必须移动
dc块以避免出现以下错误）循环路径...)
# 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")
我个人觉得R2jags
界面比较好用一点，所以下面的代码稍作修改。
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))
链条肯定还存在一些问题，但希望这能让您走得更远。