# Specifying the priors for multivariate MCMCglmm mixed model in R (Poisson distribution)

I am trying to build a model using MCMCglmm. Ideally, I would use a negative binomial distribution for my response; however, this is not an option in MCMCglmm. I don't know of any open-source alternatives for multivariate mixed models, so I am attempting to use a Poisson distribution instead. Suppose we wanted to look at diet in sharks based on water temperature and salinity. Multiple sharks were sampled at different sites, so we want to include site as a random effect to control for spatial dependence. Prey was in 3 groups: bony fishes, crustaceans, and rays (count data). In MCMCglmm, this model would be represented as:

formula = cbind(bony_fish, crustaceans, rays) ~ temp + sal


with the random effect specified as: random = ~site. This model will run if I fit it using family = rep("gaussian", 3). However, when I run it with family = rep("poisson", 3), R crashes. I can change family = rep("zipoisson", 3) and run:

model= MCMCglmm(formula,
random=~site, # G-structure
data=sharks,
family=rep("zipoisson", 3), # for 3 prey groups
rcov = ~us(trait):units, # R-structure
#prior = prior)


the crash stops, but I get the error:

Error in MCMCglmm(formula, random = ~site, data = sharks, family = rep("zipoisson",  :
Mixed model equations singular: use a (stronger) prior


So I read through these instructions, but I still haven't been able to grasp how to specify my priors. The simple form of the prior is:

prior=list(R=R, G=G)


That is simple enough. The trouble I am having is specifying R- and G-structures correctly. According to the package author here, the G-structure is the covariance matrix of the random effects, and the R-structure is the covariance matrix of the residuals. I have tried setting G as:

G=list(G1=list(V=diag(1), nu=0.002))


and R as:

R=list(V=diag(3), nu=0.002)


but it still won't work. Does anyone know the proper priors to a multivariate poisson model?

• If you can fit a random effects Poisson model with a gamma random effect on the rate for each individual (=shark), then you have effectively a negative binomial regression. One alternative you have is using Stan via the rstan package. Jan 10, 2016 at 7:06
• Thanks for the input. I am new to Bayesian, but do have some modeling experience. STAN sounds like a very interesting alternative. I did not realize that it was capable of multivariate analysis. I see there is also a Python interface now. This is something I look forward to learning. Right now it would be a problem for me to learn such a complex interface in the time allotted for me for this analysis. Jan 10, 2016 at 16:44
• Did you try this for the random effects prior: G=list(G1=list(V=diag(3), nu=0.002)) given that you have 3 response variables Nov 11, 2020 at 12:52