Long live BUGS !

[Kery Marc]

Hi all,

looking at some post in the BUGS mailing list, I just realized that WinBUGS/OpenBUGS is not dead yet :) It seems to me that these projects have grown into another one called 'multibugs', see here: https://www.multibugs.org/

This is really great because it means that we continue to have three MCMC engines at our avail in which complex hierarchical models can be specified in a hierarchical manner, that is, without having to sum out from the likelihood any discrete latent variables such as abundance (N) or presence/absence (z):  JAGS, Win/Open/MultiBUGS and Nimble

They seem to have a paper in press: Goudie, R. J. B., Turner, R. M., De Angelis, D., Thomas, A. (in press) MultiBUGS: A parallel implementation of the BUGS modelling framework for faster Bayesian inference. Journal of Statistical Software. arXiv:1704.03216 (slides pdf)

Have to check out soon and I'd be interested if anyone has any experience with multibugs for the types of models we talk about on this list ?

Best regard  --- Marc

[Kery Marc]

Dear all,

I just got a question "And what about greta ?". greta is a new MCMC-based statistical modeling software developed by Nick Golding and I have heard (and a little also seen) many good things about it. But I don't know it well enough to say more about it.

My main point in my earlier email was this:

- in some MCMC engines you can specify the likelihoods of hierarchical models 'directly', i.e., by a series of linked distributional statements, e.g., for an occupancy model:

   z_i ~ Bernoulli(psi)

   y_ij ~ Bernoulli(z_i * p)

- in this hierarchical construction of the likelihood, you have the z's which are latent variables or discrete random effects

- in other HM software, you cannot write the likelihood in this 'hierarchical' manner, but you have to get rid of the discrete latent variables (e.g., for abundance N or presence/absence z) by summation.

- for trained statisticians, summing out the N's and the z's does not seem to be a big deal, at all, but for us ecologists it is. I am convinced that the ability to specify hierarchical models with discrete latent variables in this 'direct' way is absolutely key for a wide adoption of such a software by ecologists: in my experience with hundreds of colleagues, this is a manner in which fairly many ecologists can learn to actually fit their hierarchical models without having to pester a statistician for help every time they want to fit a new model or they want to change anything a little more substantial in an existing model.

- In JAGS, Nimble, Win/Openbugs and presumably now also in multibugs, you can specify your HMs in that way, and this is crucial selling point for these software pieces with us ecologists.

- I am sure that you can fit all the models that you can fit in BUGS et al also in greta. But I don't know about the details of the model-definition language in greta and how ecologist-friendly greta is.

Best regards  --- Marc

---

From: hmec...@googlegroups.com [hmec...@googlegroups.com] on behalf of Kery Marc [marc...@vogelwarte.ch]
Sent: 15 May 2019 11:19
To: hmec...@googlegroups.com
Subject: Long live BUGS !

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2019)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To post to this group, send email to hmec...@googlegroups.com.
Visit this group at https://groups.google.com/group/hmecology.
To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/F801018864E940468323F405DF95D2DD76005DB9%40MailNT.vogelwarte.ch.
For more options, visit https://groups.google.com/d/optout.

[Kery Marc]

Dear all,

here's a 'short' message about 'greta' from its developer Nick Golding, to whom I forwarded the couple messages in this thread. He wrote "I didn't have the time to write a short message" --- love that !

Best regards  ---- Marc

# --------------------------------  Start of Nick's message ---------------------------------------------------

Like Stan, greta can't currently handle discrete parameters since it has until recently been primarily focussed on HMC, which doesn't work with discrete parameters. That's changing though. There are other samplers in there now, and support for discrete parameters is on the horizon (it's the thing my colleagues bug me about the most). Probably won't be released until next year though.

On Marc's point about statistician vs. ecologist users, and whether greta is intended for ecologists, I think greta is quite different to the other packages. greta is designed to be extended with other R packages, written by people other than me and the core developers, that add or streamline functionality. So statistically-minded people who understand these issues can write and nicely wrap up code that others can use.

An example of this is this prototype extension package (*which is not tested and I think does not yet work correctly, so please don't use it!*) which includes the exact marginalisation thing that Marc describes, implemented as a new probability distribution. This is the source code: https://github.com/goldingn/greta.multivariate/blob/master/R/conditional_bernoulli.R but that code never needs to be seen by users. All they would need to do is load the extension package and define the following the likelihood alongside their greta model for p and psi:

  distribution(y) <- conditional_bernoulli(p, psi)

There are other examples of prototype extension packages, such as one for Gaussian process modelling (spatial and temporal random effects) and for adding GAM smoothers, via the exact same formula & smoother syntax as mgcv. The idea is that people can write handy packages for specific problems (think SCR, distance sampling, JSDMs, etc.) with components that can be used together.

The crucial point here is that getting experts in these methods to write the code once, well, and then distribute to users is a much more sustainable and bug-resistant approach than requiring users to copy-paste snippets of code from stackoverflow and try to adapt it to their problems.

One of my favourite counterexamples is this Stan code for sampling the multivariate probit distirbution in Stan: https://github.com/stan-dev/example-models/blob/master/misc/multivariate-probit/probit-multi-good.stan It's very clever and very useful, but distributing it in this way, rather than as a distribution makes it both very hard to use and prone to errors in people copy-pasting it into their model!

So unlike (Win/Open/Multi)BUGS, JAGS and Stan which have a single interface which is used by both beginner and advanced users, greta has a user API, a extension-developer API (which requires using R6 classes, but I'm working on making it easier), and a package management system (it's just R packages, so they can be distributed on CRAN or GitHub!). Hopefully this will make it easier for the statistical ecology community to develop and share new models and clever ideas than is possible with existing MCMC software.

One final point. It looks like multibugs has the nice feature of being able to parallelise computations *within* chains. Ie. in addition to running each MCMC chain on a separate CPU, each one of those can use multiple CPUs. If I understand their paper correctly (from a very cursory skim read) this still works with the BUGS ideas of nodes, but can use multiple cores for bits of the computational graph that are in parallel. I think Stan has ben working on something similar, but it's quite limited and the interface is not amenable to most users. greta automatically enables parallel computation within a single chain of MCMC, but does it in a different way. greta does it by vectorising operations (rather than having to loop through indices), and using very highly optimised parallel linear algebra routines that can run on multiple CPUs and GPUs. That's all handled by Tensorflow, so it'll run much much faster than the multiBUGS approach when the model can be expressed in a vectorised way or with linear algebra. But greta's approach cannot parallelise when the user doesn't vectorise, but writes for-loops over elements of vectors/matrices. So it takes a different approach to writing the code than BUGS, and is less suited to things like iterative processes (some time series models that cannot be expressed with vectorised operations), where I'd expect the Stan/multiBUGS approach to be faster.

# --------------------------------  End of Nick's message ---------------------------------------------------

---

From: hmec...@googlegroups.com [hmec...@googlegroups.com] on behalf of Kery Marc [marc...@vogelwarte.ch]

Sent: 15 May 2019 15:11
To: hmec...@googlegroups.com
Subject: RE: Long live BUGS !

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/F801018864E940468323F405DF95D2DD76005F7A%40MailNT.vogelwarte.ch.

[Perry de Valpine]

Hi Marc,

It looks like you opened up more than you intended with this thread :)

To add to Nick's comments on Greta, I'll share some brief comments on NIMBLE (R package nimble):

* (Nearly) compatible with the BUGS/JAGS model language (but we've made the language extensible).  There are notes on converting from BUGS/JAGS here (https://r-nimble.org/quick-guide-for-converting-from-jags-or-bugs-to-nimble).  The examples from Kéry and Royle's Applied Hierarchical Modeling Vol. I are given for NIMBLE, with some comparisons of MCMC performance (without customization of any kind) here (https://github.com/nimble-dev/AHMnimble).

* All models are written from R and automatically compiled to C++ along with algorithms such as MCMC.  Hence, the workflow is in R.

* It is easy to write new distributions and new functions to use in model code and have them automatically compiled to C++ (no need to know C++).  A new "distribution" is a somewhat general concept that can include, for example, integrating over discrete latent states.  Examples with multi-state / multi-event capture-recapture showing very large MCMC efficiency gains compared to JAGS and/or simpler approaches also in nimble are here (https://link.springer.com/article/10.1007/s10651-016-0353-z).  

* (Very technical comment: NIMBLE's generated C++ uses the Eigen package, which also underlies CPU computations in TensorFlow.)

* NIMBLE includes a variety of MCMC samplers.  One can use a default configuration of samplers or easily modify the configuration.

* It is easy to write new MCMC samplers.  Samplers that work well with particular models can give very large boosts in MCMC efficiency.

* An example of writing a new sampler is a reversible-jump sampler for variable selection that is given on our examples page here (https://r-nimble.org/examples).  This is something we're developing more, so please let me know if you are interested.

* One fairly specialized set of methods supported in nimble are Bayesian non-parametric distributions.  See here (https://r-nimble.org/bayesian-nonparametric-models-in-nimble-part-1-density-estimation).

* Like Nick's point about Greta, we intend for statistical developers to be able to do the careful work of writing such extensions, which they can then provide as their own R packages that rely on nimble.  Others should then be able to use such packages without specialized statistical knowledge.

* We are working on a "nimbleEcology" package that will provide some basic tools for capture-recapture, occupancy, and N-mixture models as a start.  Please let me know if you are interested in helping to put this together.

* In collaboration with multiple groups (this is a casual email, so I won't go into details), we have very good results boosting MCMC efficiency for spatial capture-recapture models (see tweet from ISEC 2018 on Daniel Turek's 41x speedup for a SCR model here: https://twitter.com/MethodsEcolEvol/status/1014823959148269568).

* We have HMC based on automatic differentiation (similar sampling methodology as Stan; similar AD methodology as TMB, based on CppAD) working in development, but not yet released.  It will be possible to use HMC in conjunction with other methods.

* NIMBLE is not just an MCMC engine.  It is a hierarchical model and model-generic algorithm programming system (see here: https://www.tandfonline.com/doi/abs/10.1080/10618600.2016.1172487).

* Materials from a workshop for ecologists last year are available here (https://github.com/perrydv/Vogelwarte_NIMBLE_workshop).

* For anyone attending The Wildlife Society conference in Reno in September/October, I will be giving a one-day nimble workshop there.

* We do our best to respond quickly to user questions.

* You can learn more at r-nimble.org.

Thanks for reading!

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/F801018864E940468323F405DF95D2DD76006183%40MailNT.vogelwarte.ch.

[Kery Marc]

Dear Perry,

 

not at all ! Thank you very much for your excellent description of your exciting software; no doubt this is very informative and helpful for the readers of this list. And of course thank you very much (to you and your colleagues) for creating Nimble and to work so hard to make it progress all the time.

 

Thank you also for translating all the examples in the AHM1 book. You may soon want to consider translating a whole new set of models. AHM2 will soon be done, finally (argh)....

 

As you know, two of the things I really love about Nimble are:

      (1) that you have adopted almost the original BUGS language for model definition and

      (2) that Nimble can do discrete latent variables.

In my opinion, these two things alone are absolutely crucial for making a command-line-driven modeling software really widely usable for ecologists who are interested in fitting the types of hierarchical models that we talk about on this list.

 

As opposed to statisticians, for whom these two things are not such a big deal. Or, statisticians who may not have encountered cases where in such mixture models we do have a real interest in the values of these latent variables, because they mean something eminently important in biology (e.g., abundance, presence/absence, being alive or dead etc.).

 

At the same time, I understand that one of the goals of Nimble is to make it readily extensible, e.g. by writing new distributions, for people who know how to do this. Nick's goal for greta, to have people write the own packages for his software, sounded very similar in my ears.

 

Currently, I get the impression that perhaps Nimble looks best among the Bayesian programming languages to combine the needs and interests of ecologists (and scientists and managers with similar stats/comp education) on the one hand and of research statisticians (and similarly quantitative folks) on the other. Congratulations !

 

Best regards  ---  Marc

 

 

 

---

From: Perry de Valpine [pdeva...@berkeley.edu]
Sent: 16 May 2019 17:43
To: Kery Marc
Cc: hmec...@googlegroups.com
Subject: Re: Long live BUGS !