I am trying to port a model from WinBugs to rstan, the model is self-contained yet I appear to be missing something in the porting process.
I have ported my own models but this is not my own. I am working through the port documentation but all pointers or recommendations will be gratefully received. As well as attribution combined with good karma.
WinBug 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)
# data
(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
pp <- 0
A <- dat$A
B <- dat$B
data <- list(N = N, case = case, A = A, B = B, nn = nn, pp = pp)
# BUGS model regi1.bug
# It is easier to contain the model in the WIN2Bugs Script and write to a file
# this way you can edit your model from a single R file.
#
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)
# ensuring appropriate constraints for a0, a1, a2 and a3
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
)
# regi1.bug based estimates
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
# BUGS model regi2.bug
# BUGS model regi2.bug
# It is easier to contain the model in the WIN2Bugs Script and write to a file
# this way you can edit your model from a single R file.
#
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)
# ensuring appropriate constraints for b1, b2 and b3
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
)
# regi2.bug based estimates
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