我有下面的 WinBUGS 模型，我需要在 R 中使用它们。但是，由于我不熟悉 WinBUGS 并且为了与分析的其他部分保持一致，我将不胜感激将它们转换为 rstan的任何帮助。
以下是模型和数据（完整脚本）：
错误模型
# Check your OS Version
if ( Sys.info()['sysname'] == "Darwin" ) { WINDOWS <- FALSE } else { WINDOWS <- TRUE }
if (!WINDOWS) {
# You can run WINBugs from MacOSX using wine
# https://www.r-bloggers.com/running-r2winbugs-on-a-mac-running-osx/
#
setwd(path.expand("~/dev/stackoverflow/53809468"))
bugs.directory = path.expand("~/.wine/drive_c/Program Files/WinBUGS14")
} else {
# You have to run RStudio in administrator mode to
# correctly launch and execute WINBugs
setwd(path.expand("~/dev/stackoverflow/53809468"))
bugs.directory <- 'C:/Program Files/WinBUGS14'
}
library(R2WinBUGS)
dat <-
data.frame(
A = c(1, 1, 0, 0),
B = c(1, 0, 1, 0),
Pass = c(278,
100, 153, 79),
Fail = c(743, 581, 1232, 1731))
N <- length(dat$Pass)
case <- dat$Pass
nn <- dat$Fail + dat$Pass
A <- dat$A
B <- dat$B
data <- list(N = N, case = case, A = A, B = B, nn = nn)
regi1.bug <- "model {
for (i in 1:N){
pp[i] <- exp(a0) * (1 + a1*A[i] + a2*B[i] + a3*A[i]*B[i])
case[i] ~ dbin(pp[i], nn[i])}
a0 ~ dnorm(0, 5)
a1 ~ dnorm(0, 5)
a2 ~ dnorm(0, 5)
a3 ~ dnorm(0, 5)
ones <- 1
ones ~ dbern(C1)
C1 <- step(0 - a0)*step(1 + a1)*step(exp(- a0) - 1 - a1)*step(1 - a2)*
step(exp(-a0) -1 - a2)*step(1 + a1 + a2 + a3)*step(exp(- a0) - 1 - a1 - a2 - a3)}"
writeLines(regi1.bug, "regi1.bug")
inits <- function() {
list(a0 = -4, a1 = 0, a2 = 0, a3 = 0)
}
sim1 <- bugs(
data,
inits = inits,
model.file = "regi1.bug",
parameters.to.save = c("a0", "a1", "a2", "a3"),
n.chains = 4,
n.burnin = 5e+05,
n.iter = 1e+06,
bugs.directory = bugs.directory,
debug = FALSE)
print(sim1)
# mean sd 2.5% 25% 50% 75% 97.5% Rhat
# a0 -2.8 0.1 -2.9 -2.8 -2.8 -2.7 -2.6 1
# a1 1.3 0.2 0.9 1.2 1.3 1.5 1.8 1
# a2 0.8 0.1 0.5 0.7 0.8 0.9 1.0 1
# a3 1.0 0.3 0.5 0.9 1.0 1.2 1.5 1
regi2.bug <- "model {
for (i in 1:N){
odds[i] <- exp(b0) * (1 + b1*A[i] + b2*B[i] + b3*A[i]*B[i])
pp[i] <- odds[i]/(1 + odds[i])
case[i] ~ dbin(pp[i], nn[i])
}
b0 ~ dnorm(0, 5)
b1 ~ dnorm(0, 5)
b2 ~ dnorm(0, 5)
b3 ~ dnorm(0, 5)
ones <- 1
ones ~ dbern(C1)
C1 <- step(1 + b1)*step(1 + b2)*step(1 + b1 + b2 + b3)
}"
writeLines(regi2.bug, "regi2.bug")
inits <- function() {
list(b0 = 0, b1 = 1, b2 = 1, b3 = 1)
}
sim2 <-
bugs(
data,
inits = inits,
model.file = "regi2.bug",
parameters.to.save = c("b0", "b1", "b2", "b3"),
n.chains = 4,
n.burnin = 5e+05,
n.iter = 1e+06,
bugs.directory = bugs.directory,
debug = FALSE
)
print(sim2)
# mean sd 2.5% 25% 50% 75% 97.5% Rhat
# b0 -2.7 0.1 -2.8 -2.7 -2.7 -2.6 -2.5 1
# b1 1.4 0.2 1.0 1.3 1.4 1.6 2.0 1
# b2 0.9 0.2 0.5 0.7 0.9 1.0 1.3 1
# b3 1.3 0.3 0.7 1.1 1.3 1.5 1.9 1
控制台输出：
> source('.../dev/stackoverflow/53809468/53809468.R')
Loading required package: coda
Loading required package: boot
Inference for Bugs model at "regi1.bug", fit using WinBUGS,
4 chains, each with 1e+06 iterations (first 5e+05 discarded), n.thin = 2000
n.sims = 1000 iterations saved
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
a0 -2.8 0.1 -2.9 -2.8 -2.8 -2.7 -2.6 1 560
a1 1.3 0.2 0.9 1.2 1.3 1.5 1.8 1 400
a2 0.8 0.1 0.5 0.7 0.8 0.9 1.0 1 520
a3 1.0 0.3 0.5 0.9 1.0 1.2 1.5 1 1000
deviance 42.6 5.6 33.7 38.5 42.0 46.0 55.1 1 290
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = Dbar-Dhat)
pD = 2.9 and DIC = 45.5
DIC is an estimate of expected predictive error (lower deviance is better).
Inference for Bugs model at "regi2.bug", fit using WinBUGS,
4 chains, each with 1e+06 iterations (first 5e+05 discarded), n.thin = 2000
n.sims = 1000 iterations saved
mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
b0 -2.7 0.1 -2.8 -2.7 -2.7 -2.6 -2.5 1 1000
b1 1.4 0.2 1.0 1.3 1.4 1.6 2.0 1 1000
b2 0.9 0.2 0.5 0.7 0.9 1.0 1.3 1 420
b3 1.3 0.3 0.7 1.1 1.3 1.5 1.9 1 1000
deviance 49.3 6.8 38.4 44.4 48.6 53.7 64.8 1 1000
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = Dbar-Dhat)
pD = 3.0 and DIC = 52.3
DIC is an estimate of expected predictive error (lower deviance is better).
>
预先感谢您的任何帮助。