Bym and besag models in brms

Interfaces brms
1

Hi, I was wondering if there’s currently a way to easily fit besag and bym CAR models in brms.

Thanks,
Paul

2

To be honest, I am lost in the various names of CAR models and always forget which of those is named how. The CAR models supported in brms are explained in two case studies, which are linked to under ?cor_car. It may very well be that one of those (perhaps the one by @mitzimorris?) also goes under one of the names you mentioned.

3

the functions for the CAR models are those supplied in the Max Joseph case study:

https://mc-stan.org/users/documentation/case-studies/mbjoseph-CARStan.html

the Stan code itself is here:

I wrote a case study on the ICAR models: https://mc-stan.org/users/documentation/case-studies/icar_stan.html - the critical implementation detail there is the use of a soft sum-to-zero constraint for the ICAR component which really speeds up computation.

I thought that RStanArm implemented the R package car-bayes - maybe @bgoodri knows the status of this?

1 Like
4

Imad had a branch but we haven’t merged it. I’m not sure if it is ready but we really can’t add any more models to rstanarm until the compilation time goes way down.

5

Thanks, I’ve gone through your case study as well as a recent publication that used a similar approach. It is hard for me to tell if the ICAR specification in brms is equivalent to a bym model, though.

6

the devil’s in the definition of ‘BYM model’. I think the answer is no - the BYM model is a Poisson GLM with both a multi-variate normal random effects component as well as an ICAR component and you can’t just do that. there are lots of flavors of BYM models.

I would love to see the BYM2 model go into BRMS - @paul.buerkner is this feasible? how can I help?

1 Like
7

Thanks @mitzimorris for offering to help out. I am not an expert for these kinds of models so first I need to see the mathematical formulation of the BYM2 model, before I can judge if is feasible in brms. Is such a formulation available somewhere?

8

mathematical formulation (thanks @anon75146577!) in this paper: https://arxiv.org/pdf/1403.4630.pdf, equation 6.1

and more discussion of BYM models in Riebler et al: https://arxiv.org/pdf/1601.01180.pdf

and it would be super-cool to add in stuff I worked out for maps which have disconnected sub-regions - see this paper - https://www.sciencedirect.com/science/article/abs/pii/S1877584517301600?via%3Dihub

2 Likes
9

Thanks! This is very helpful. It is only now that I realize you are discussing the BYM2 model in your case study already. Please excuse my question for a reference earlier. I should have looked at your case study beforehand.

Just for clarity, what you discuss in your case study as ICAR model (and what I have implemented in brms via cor_icar), is some flavor of the classical BYM model, right? More precisely, when me model

brm(y ~ ... + (1 | obs), family = poisson(), autocor = cor_icar(...))

where obs is an indicator of the rows (i.e. unique per observation), we end up with the BYM model you discuss in Section “Example: disease mapping using the Besag York Mollié model”. Is my understanding correct?

Provided that the above is roughly correct, The BYM2 model (where we have a joint prior over the spatial and non spatial residuals) should be relatively straightforward to implemented. From the user perspective it would simply look like

brm(y ~ ..., family = poisson(), autocor = cor_icar_het(...))

(or a better name than cor_icar_het for the icar structure with depending non-spatial residuals) as the non-spatial residuals are now part of the autocorrelation structure so (1|obs) is no longer needed.

I will take a look at the disconneted sub-regions paper later and then respond separately.

10

not really - it’s an ICAR model as introduced in Besag 1974. the innovation of the BYM model was to add a 2nd component for extra-spatial variation. Max Joseph implemented both the CAR and ICAR models in Stan. my understanding is that the epidemiologists use some version of the BYM model, and that they generally use INLA to fit it.

11

That makes sense. Thanks for taking the time to help me understand this topic. I have opened an issue for brms (https://github.com/paul-buerkner/brms/issues/710) so that I don’t forget to work on this topic.

3 Likes
12

Thanks to both of you! Yes it’s something that I have to use INLA for at the moment, and while INLA has its difficulties to learn it is straightforward to bring in a shapefile from GIS, build an adjacency matrix, and plug that into a besag or bym model in INLA.

13

what format, exactly? the Stan model is written in terms of the set of <i,j> edges in the graph - we can handle the conversion - it’s just a matter of deciding on the set of accepted input formats and documenting them properly.

1 Like
14

I believe I’ve been using the nb2adj function in spdep, starting with a shapefile.

15

Hi, i’m transitioning from INLA into Stan, and using these models, since cor_car() takes a matrix argument, shouldn’t nb2mat() be what’s used?

1 Like
16

Perhaps - though since it seems besag/bym models aren’t implemented in brms I haven’t yet tried to run a spatial cor_car() model.

I’m finishing up a paper now where instead of using a bym model i’ve run GAMs using the polygon centroid coordinate pairs. It would be interesting to see how that differs from a CAR type model, as it’s certainly very easy to implement GAMs.

17

I’m going back to things i’ve done in winbugs and INLA and using Stan, so i’m trying to recover all the estimates from before, most of mine were CAR spatial models. Once i’m happy with it i’ll publish to rpubs or something

18

Bym/bym2 is supported as of the latest brms version via a specific type in the cor_car Function.

1 Like
19

(Just in case) you need to throw massive warnings if anyone tries to use bayes factors on a model with a bym/ICAR because the normalizing constant is almost definitely wrong!

20

I think I have added such a warning at some point but I need to check if the check is still valid.

1 Like