0

I made a model using R2jags. I like the jags syntax but I find the output produced by R2jags not easy to use. I recently read about the rstanarm package. It has many useful functions and is well supported by the tidybayes and bayesplot packages for easy model diagnostics and visualisation. However, I'm not a fan of the syntax used to write a model in rstanarm. Ideally, I would like to get the best of the two worlds, that is writing the model in R2jags and convert the output into a Stanreg object to use rstanarm functions.

Is that possible? If so, how?

Junitar
  • 905
  • 6
  • 13

1 Answers1

1

I think then question isn't necessarily whether or not it's possible - I suspect it probably is. The question really is how much time you're prepared to spend doing it. All you'd have to do is try to replicate in structure the object that gets created by rstanarm, to the extent that it's possible with the R2jags output. That would make it so that some post-processing tasks would probably work.

If I might be so bold, I suspect a better use of your time would be to turn the R2jags object into something that could be used with the post-processing functions you want to use. For example, it only takes a small modification to the JAGS output to make all of the mcmc_*() plotting functions from bayesplot work. Here's an example. Below is the example model from the jags() function help.

# An example model file is given in:
model.file <- system.file(package="R2jags", "model", "schools.txt")

# data
J <- 8.0
y <- c(28.4,7.9,-2.8,6.8,-0.6,0.6,18.0,12.2)
sd <- c(14.9,10.2,16.3,11.0,9.4,11.4,10.4,17.6)


jags.data <- list("y","sd","J")
jags.params <- c("mu","sigma","theta")
jags.inits <- function(){
  list("mu"=rnorm(1),"sigma"=runif(1),"theta"=rnorm(J))
}

jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
                n.iter=5000, model.file=model.file, n.chains = 2)

Now, what the mcmc_*() plotting functions from bayesplot expect is a list of matrices of MCMC draws where the column names give the name of the parameter. By default, jags() puts all of them into a single matrix. In the above case, there are 5000 iterations in total, with 2500 as burnin (leaving 2500 sampled) and the n.thin is set to 2 in this case (jags() has an algorithm for identifying the thinning parameter), but in any case, the jagsfit$BUGSoutput$n.keep element identifies how many iterations are kept. In this case, it's 1250. So you could use that to make a list of two matrices from the output.

jflist <- list(jagsfit$BUGSoutput$sims.matrix[1:jagsfit$BUGSoutput$n.keep, ], 
               jagsfit$BUGSoutput$sims.matrix[(jagsfit$BUGSoutput$n.keep+1):(2*jagsfit$BUGSoutput$n.keep), ])

Now, you'd just have to call some of the plotting functions:

mcmc_trace(jflist, regex_pars="theta")

enter image description here

or

mcmc_areas(jflist, regex_pars="theta")

enter image description here

So, instead of trying to replicate all of the output that rstanarm produces, it might be a better use of your time to try to bend the jags output into a format that would be amenable to the post-processing functions you want to use.

EDIT - added possibility for pp_check() from bayesplot.

The posterior draws of y in this case are in the theta parameters. So, we make an object that has elements y and yrep and make it of class foo

x <- list(y = y, yrep = jagsfit$BUGSoutput$sims.list$theta)
class(x) <- "foo"

We can then write a pp_check method for objects of class foo. This come straight out of the help file for bayesplot::pp_check().

pp_check.foo <- function(object, ..., type = c("multiple", "overlaid")) {
  y <- object[["y"]]
  yrep <- object[["yrep"]]
  switch(match.arg(type),
         multiple = ppc_hist(y, yrep[1:min(8, nrow(yrep)),, drop = FALSE]),
         overlaid = ppc_dens_overlay(y, yrep[1:min(8, nrow(yrep)),, drop = FALSE]))
}

Then, just call the function:

pp_check(x, type="overlaid")

enter image description here

DaveArmstrong
  • 18,377
  • 2
  • 13
  • 25
  • Thanks for your input. "[…] try to bend the `jags` output into a format that would be amenable to the post-processing functions you want to use": my question was probably not clear enough but that's what I'm trying to do. I can actually already use the `mcmc_trace()` and `mcmc_areas()` functions after converting the `jags` output using `as.mcmc()`. There're however other functions that I would like to use that require a `rstanarm` object. E.g. the `pp_check()` function from the `rstanarm` package or the `add_fitted_draws()` from the `tidybayes` package. – Junitar Sep 09 '20 at 14:25
  • For `add_fitted_draws()` I think it would be pretty difficult. It relies on all of the underlying infrastructure of an `rstanarm` model. It seems like it would take a lot of time to try to replicate that. You might be better off trying to consistently use the same coding conventions in your own code (i.e., `mu[i] <- inprod(X[i, ], b[])` such that it would be easy to generate predictions for new data, which is what `add_fitted_draws()` does. – DaveArmstrong Sep 09 '20 at 15:06
  • As for `pp_check()`, the documentation in `bayesplot` suggests that it's supposed to be generic - allowing other package authors to generate methods for `pp_check` for their own purposes. Answer above edited to show one possible way. – DaveArmstrong Sep 09 '20 at 15:09
  • Thank you for the additional info and for the edit. That's useful even if it only solves half of the problem. – Junitar Sep 10 '20 at 07:17