Dear all,

I would need to run Bayesian regression models with a negative binomial distribution on zero-inflated data.

I am fitting a brm zero_inflated_negbinomial model with no specific prior distribution.

My data are several matrices, each one is a rectangular matrix with biological species in columns and locations in rows.

The entry of the matrix is count data, with many zeros.

You can find an example of such matrix in the link below

Now, my problem is that for each matrix I want to fit several models, as many as the number of biological species (i.e., columns).

Each model consists in a single species j as response and all other species as predictors, such as

mod <- glm.nb(mat1[,j] ~mat1[,-j], data = data.frame(mat1))

This analysis consists of several thousand models to fit.

The problem I encounter using the brm function is that it requires to specify the name of the response as is the dataframe as well as the entire formula in full with all the names of the predictors.

Indeed, it does not work if I run the model as

brm(mat1[,j] ~mat1[,-j], family = zero_inflated_negbinomial(), data = data.frame(mat1), control = list(max_treedepth = 12), cores = numCores, seed=2504)

but I have to specify each species name like

mod<-brm(Androsace.chamaejasme ~ Arenaria.ciliata + Botrychium.lunaria, family = zero_inflated_negbinomial(), data = mat1, control = list(max_treedepth = 12), cores = numCores, seed=2504)

However, this makes the analyses quite complicated, as I cannot materially type-in all the names of the species – the formula – each time.

Do you might have an idea how to solve that?

Thank you very much,

All the best,

Gianalberto