Translating WinBUGS models to rstan
.everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty,.everyoneloves__bot-mid-leaderboard:empty{ height:90px;width:728px;box-sizing:border-box;
}
4
I have WinBUGS models below which I need to use in R. However because I am not familiar with WinBUGS and for consistency with other parts of the analysis, I would appreciate any help to convert them to rstan.
Here are the models and data (full script):
BUGS Model
# 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
Console Output:
> 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).
>
Thank you in advance for any help.
r winbugs stan rstan r2winbugs
asked Dec 17 '18 at 5:35
KrantzKrantz
3881416
|
show 4 more comments
4
I have WinBUGS models below which I need to use in R. However because I am not familiar with WinBUGS and for consistency with other parts of the analysis, I would appreciate any help to convert them to rstan.
Here are the models and data (full script):
BUGS Model
# 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
Console Output:
> 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).
>
Thank you in advance for any help.
r winbugs stan rstan r2winbugs
asked Dec 17 '18 at 5:35
KrantzKrantz
3881416
3
There's a chapter at the end of the Stan User's Guide specifically aimed at converting BUGS models to Stan.
– Bob Carpenter
Dec 17 '18 at 20:59
4
There's also an R2WinBUGS interface on CRAN and also wrappers for OpenBUGS, which is just a port of WinBUGS and should run the same programs.
– Bob Carpenter
Dec 19 '18 at 16:40
1
@Krantz Do you have the full WinBUGS model? The data you provide is not in the right format for WinBUGS. It should be a list (i.e.list ( A = c( 1, 1, 0, 0) , B = c( 1, 0, 1, 0) , case = c( 278, 100, 153, 79) , negative = c( 743, 581, 1232, 1731) )having the full working WinBUGS model would make it easier to convert. Are you using R2WinBUGS - do you have the complete model script?
– Technophobe01
Dec 25 '18 at 23:57
2
Here is a link to the porting documentation for Stan mc-stan.org/docs/2_18/stan-users-guide/…
– Technophobe01
Dec 26 '18 at 0:03
1
@Krantz I updated the code example to be self-contained so people can copy and run if they have WinBUGs installed.
– Technophobe01
Dec 26 '18 at 15:50
|
show 4 more comments
4
4
4
1
I have WinBUGS models below which I need to use in R. However because I am not familiar with WinBUGS and for consistency with other parts of the analysis, I would appreciate any help to convert them to rstan.
Here are the models and data (full script):
BUGS Model
# 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
Console Output:
> 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).
>
Thank you in advance for any help.
r winbugs stan rstan r2winbugs
asked Dec 17 '18 at 5:35
KrantzKrantz
3881416
I have WinBUGS models below which I need to use in R. However because I am not familiar with WinBUGS and for consistency with other parts of the analysis, I would appreciate any help to convert them to rstan.
Here are the models and data (full script):
BUGS Model
# 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
Console Output:
> 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).
>
Thank you in advance for any help.
r winbugs stan rstan r2winbugs
r winbugs stan rstan r2winbugs
asked Dec 17 '18 at 5:35
KrantzKrantz
3881416
asked Dec 17 '18 at 5:35
KrantzKrantz
3881416
edited Jan 11 at 13:25
Krantz
asked Dec 17 '18 at 5:35
KrantzKrantz
3881416
asked Dec 17 '18 at 5:35
KrantzKrantz
3881416
asked Dec 17 '18 at 5:35
KrantzKrantz
3881416
3881416
3
There's a chapter at the end of the Stan User's Guide specifically aimed at converting BUGS models to Stan.
– Bob Carpenter
Dec 17 '18 at 20:59
4
There's also an R2WinBUGS interface on CRAN and also wrappers for OpenBUGS, which is just a port of WinBUGS and should run the same programs.
– Bob Carpenter
Dec 19 '18 at 16:40
1
@Krantz Do you have the full WinBUGS model? The data you provide is not in the right format for WinBUGS. It should be a list (i.e.list ( A = c( 1, 1, 0, 0) , B = c( 1, 0, 1, 0) , case = c( 278, 100, 153, 79) , negative = c( 743, 581, 1232, 1731) )having the full working WinBUGS model would make it easier to convert. Are you using R2WinBUGS - do you have the complete model script?
– Technophobe01
Dec 25 '18 at 23:57
2
Here is a link to the porting documentation for Stan mc-stan.org/docs/2_18/stan-users-guide/…
– Technophobe01
Dec 26 '18 at 0:03
1
@Krantz I updated the code example to be self-contained so people can copy and run if they have WinBUGs installed.
– Technophobe01
Dec 26 '18 at 15:50
|
show 4 more comments
3
There's a chapter at the end of the Stan User's Guide specifically aimed at converting BUGS models to Stan.
– Bob Carpenter
Dec 17 '18 at 20:59
4
There's also an R2WinBUGS interface on CRAN and also wrappers for OpenBUGS, which is just a port of WinBUGS and should run the same programs.
– Bob Carpenter
Dec 19 '18 at 16:40
1
@Krantz Do you have the full WinBUGS model? The data you provide is not in the right format for WinBUGS. It should be a list (i.e.list ( A = c( 1, 1, 0, 0) , B = c( 1, 0, 1, 0) , case = c( 278, 100, 153, 79) , negative = c( 743, 581, 1232, 1731) )having the full working WinBUGS model would make it easier to convert. Are you using R2WinBUGS - do you have the complete model script?
– Technophobe01
Dec 25 '18 at 23:57
2
Here is a link to the porting documentation for Stan mc-stan.org/docs/2_18/stan-users-guide/…
– Technophobe01
Dec 26 '18 at 0:03
1
@Krantz I updated the code example to be self-contained so people can copy and run if they have WinBUGs installed.
– Technophobe01
Dec 26 '18 at 15:50
3
3
There's a chapter at the end of the Stan User's Guide specifically aimed at converting BUGS models to Stan.
– Bob Carpenter
Dec 17 '18 at 20:59
There's a chapter at the end of the Stan User's Guide specifically aimed at converting BUGS models to Stan.
– Bob Carpenter
Dec 17 '18 at 20:59
4
4
There's also an R2WinBUGS interface on CRAN and also wrappers for OpenBUGS, which is just a port of WinBUGS and should run the same programs.
– Bob Carpenter
Dec 19 '18 at 16:40
There's also an R2WinBUGS interface on CRAN and also wrappers for OpenBUGS, which is just a port of WinBUGS and should run the same programs.
– Bob Carpenter
Dec 19 '18 at 16:40
1
1
@Krantz Do you have the full WinBUGS model? The data you provide is not in the right format for WinBUGS. It should be a list (i.e.
list ( A = c( 1, 1, 0, 0) , B = c( 1, 0, 1, 0) , case = c( 278, 100, 153, 79) , negative = c( 743, 581, 1232, 1731) ) having the full working WinBUGS model would make it easier to convert. Are you using R2WinBUGS - do you have the complete model script?– Technophobe01
Dec 25 '18 at 23:57
@Krantz Do you have the full WinBUGS model? The data you provide is not in the right format for WinBUGS. It should be a list (i.e.
list ( A = c( 1, 1, 0, 0) , B = c( 1, 0, 1, 0) , case = c( 278, 100, 153, 79) , negative = c( 743, 581, 1232, 1731) ) having the full working WinBUGS model would make it easier to convert. Are you using R2WinBUGS - do you have the complete model script?– Technophobe01
Dec 25 '18 at 23:57
2
2
Here is a link to the porting documentation for Stan mc-stan.org/docs/2_18/stan-users-guide/…
– Technophobe01
Dec 26 '18 at 0:03
Here is a link to the porting documentation for Stan mc-stan.org/docs/2_18/stan-users-guide/…
– Technophobe01
Dec 26 '18 at 0:03
1
1
@Krantz I updated the code example to be self-contained so people can copy and run if they have WinBUGs installed.
– Technophobe01
Dec 26 '18 at 15:50
@Krantz I updated the code example to be self-contained so people can copy and run if they have WinBUGs installed.
– Technophobe01
Dec 26 '18 at 15:50
|
show 4 more comments
0
active
oldest
votes
Your Answer
StackExchange.ifUsing("editor", function () {
StackExchange.using("externalEditor", function () {
StackExchange.using("snippets", function () {
StackExchange.snippets.init();
});
});
}, "code-snippets");
StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "1"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});
function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});
}
});
draft saved
draft discarded
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53809468%2ftranslating-winbugs-models-to-rstan%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
0
active
oldest
votes
0
active
oldest
votes
active
oldest
votes
active
oldest
votes
draft saved
draft discarded
Thanks for contributing an answer to Stack Overflow!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
To learn more, see our tips on writing great answers.
draft saved
draft discarded
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53809468%2ftranslating-winbugs-models-to-rstan%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
3
There's a chapter at the end of the Stan User's Guide specifically aimed at converting BUGS models to Stan.
– Bob Carpenter
Dec 17 '18 at 20:59
4
There's also an R2WinBUGS interface on CRAN and also wrappers for OpenBUGS, which is just a port of WinBUGS and should run the same programs.
– Bob Carpenter
Dec 19 '18 at 16:40
1
@Krantz Do you have the full WinBUGS model? The data you provide is not in the right format for WinBUGS. It should be a list (i.e.
list ( A = c( 1, 1, 0, 0) , B = c( 1, 0, 1, 0) , case = c( 278, 100, 153, 79) , negative = c( 743, 581, 1232, 1731) )having the full working WinBUGS model would make it easier to convert. Are you using R2WinBUGS - do you have the complete model script?– Technophobe01
Dec 25 '18 at 23:57
2
Here is a link to the porting documentation for Stan mc-stan.org/docs/2_18/stan-users-guide/…
– Technophobe01
Dec 26 '18 at 0:03
1
@Krantz I updated the code example to be self-contained so people can copy and run if they have WinBUGs installed.
– Technophobe01
Dec 26 '18 at 15:50