Porting a WinBugs Model to Stan

Modeling
1

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
1 Like
2

what exactly don’t you understand?
what process is the BUGS program trying to model?
there are many case studies which cover different kinds of regression models.
https://mc-stan.org/users/documentation/case-studies.html

the section in the User’s Guide is pretty explicit w/r/t the differences between a BUGS program and a Stan program - in the Stan program you first declare your data inputs in the data block, your parameters in the parameters block, and you state the likelihood and priors in the model block.

the Stan wiki has a page on choice of priors:

1 Like
3

Mitzmorris,

Thanks for the pointers. I was having trouble mapping the data inputs and priors into the Stan model. I felt it was an interesting question on stackoverflow and felt it better to post here.

I.e I was trying to help the actual question owner.

My problem is that I am traveling, so had to pause on looking into the solution.

1 Like
4

you’re welcome.
it’s always good to know what different people find confusing/difficult to understand.

the priors used in a BUGS program are often extremely diffuse in order to allow the BUGS sampler to make progress. Stan’s HMC-NUTS sampler is able to avoid many of the difficulties the BUGS sampler runs into, therefore you generally don’t need or want these extremely diffuse priors (unless you really mean it).

1 Like
6

Now, following the stan documentation https://mc-stan.org/docs/2_18/stan-users-guide/stan-for-bugs-appendix.html, I am able to translate the above BUGS model to a working stan program.

Thank you in advance for any help.