Account for relationships between species in a multivariate brms model

Modeling
, , ,
1

Hi all. I am currently working with a dataset containing 30 years of fish abundance for 12 fish taxa across 5 river systems, and I would like to construct joint species distribution models to determine where artificial barriers (e.g. dams, weirs, culverts) are spatially structuring fish communities. I have devoted a lot of time to this in HMSC (a terrific package for the record), but unfortunately the long run-times and lack of support for offset terms made it unsuibtable for my study with the resources available to me. I am now giving the same analysis a go in brms, but am struggling to implement some of the features I was using in HMSC.

This is how my model is currently specified. Distance, DateNum and Season are used to capture “natural” spatiotemporal patterns in fish abundance, while up_bar is a factor which is there to capture discrete “jumps” in community structure and fish abundance across barriers. SampleID is included as a random effect because there are multiple observations nested in samples.

nam_bars_bform_ts <- bf(mvbind(Bony_herring, Common_carp, Golden_perch, Murray_Darling_rainbowfish, Carp_gudgeon_species_complex, Spangled_perch, Murray_cod, Eastern_Gambusia, Australian_smelt, Goldfish, Unspecked_hardyhead, Freshwater_catfish) 
                       ~ s(Distance, bs = "ts") + s(DateNum, bs = "ts") + Season + 
                          up_bar + 
                          (1|d|SampleID) + 
                          offset(log(ElectrofishingDuration)))

nam_bars_fit <- brm(nam_bars_bform_ts, data = nam_wide,
                    family = zero_inflated_negbinomial(),
                    chains = 4, cores = 4,
                    iter = 2000, thin = 10,
                    backend = "cmdstanr", threads = threading(4))

A key feature of JSDMs is accounting for relationships and covariance between species, and even allowing some variance to be shared among species for more robust estimates for rarer species. It seems if I had a Gaussian model I could do a version of this with set_rescor = TRUE, but since my data are zero-inflated overdispersed counts, this won’t work for me. I’ve done a bit of reading around on this forum and in the GitHub issues, and it looks like what I am after is currently in the works for brms 3.0, but isn’t ready yet. Since my honours thesis is due in November, I’ll definitely need to find another solution.

Previously in my project I worked with glmmTMB (but had to go Bayesian due to singular fits), which uses a latent variables approach to account for species intercorrelation, via a reduced rank covariance structure: Covariance structures with glmmTMB. Is this something I could implement in a multivariate brms model?

Another option would be to put my data in long-format, with “Count” as my response variable and “Species” as a factor, and add a random slope by species for all of my fixed terms. Then I could use fcor to add a covariance matrix which accounts for intercorrelation between species, perhaps based on bray-curtis similarity or something similar.

nam_bars_fit <- brm(Count ~ s(Distance, bs = "ts") + s(DateNum, bs = "ts") + Season + up_bar + 
                    (s(Distance, bs = "ts") + s(DateNum, bs = "ts") + Season + up_bar | Species) +
                    fcor(fish_covariance_matrix) +
                    data = nam_long,
                    family = zero_inflated_negbinomial())

The two problems here are that 1. I don’t know if this would actually do what I want it to achieve, and 2. this creates a much more complex model than the multivariate fit, and puts my effects of interest as random effects, making then harder to interpret. I would therefore much prefer to incorporate this in a multivariate model.

Any help here would be greatly appreciated. My undergraduate thesis has taken me on a tour of just about every statistical package R has to offer, but I’m incredibly keen to finally get some results!

2

You can do this with Stan using the mvgam 📦. There is now a dedicated wrapper called jsdgam() for adding latent factors to approximate unmodelled correlations. Similarly to gllvm and glmmtmb, you can include a wide range of effects in the factor formulae (i.e. to capture spatial, temporal or other sources of autocorrelation).

1 Like