4
$\begingroup$

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?

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ 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. $\endgroup$
    – Björn
    Jan 10, 2016 at 7:06
  • $\begingroup$ 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. $\endgroup$
    – user14241
    Jan 10, 2016 at 16:44
  • $\begingroup$ 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 $\endgroup$
    – Melissa
    Nov 11, 2020 at 12:52

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.