CODEWITHSUNDEEP
rglmbayesianmixed-modelsphylogeny
Mila
Mila

Reputation: 51

Set G in prior using MCMCglmm, with categorical response and phylogeny

I am new to the MCMCglmm package in R, and rather new to glm models in general. I have a dataset of species traits and whether or not they have been introduced outside of their native range.

I would like to test whether being introduced (as a binary 0/1 response variable) can be explained by any of the species traits. I would also like to correct for phylogeny between species.

I was told that for a binary response I could use family =“threshold” and I should fix the residual variance at 1. But I am having some trouble with the other parameters needed for the prior.

I've specified the R value for the random effects, but if I specify R I must also specify G and it is not clear to me how to decide the values for this parameter. I've tried putting default values but I get error messages:

Error in MCMCglmm(fixed, random = ~species, data = data2, family = "threshold",  : 
prior$G has the wrong number of structures

I have read the help vignettes and course but have not found an example with a binary response, and it is not clear to me how to decide the values for the priors. This is what I have so far:

fixed=Intro_binary ~ Trait1+ Trait2 + Trait3 
Ainv=inverseA(redTree1)$Ainv

binary_model = MCMCglmm(fixed, random=~species, data = data, family = "threshold", ginverse=list(species=Ainv),
 prior = list( 
    G = list(),    #not sure about the parameters for random effects.
    R = list(V = 1, fix = 1)),  #to fix the residual variance at one
  nitt = 60000, burnin = 10000)

Any help or feedback would be greatly appreciated!

Upvotes: 5

Views: 2640

Answers (1)

Thomas Guillerme
Thomas Guillerme

Reputation: 1857

This one is a bit tricky with the information you provide. I'd say you can define G as a "weak" prior using:

priors <- list(R = list(V = 1, nu = 0.002),
               G = list(V = 1, fix = 1)))

binary_model <- MCMCglmm(fixed, random = ~species, data = data,
                         family = "threshold",
                         ginverse = list(species = Ainv),
                         prior = priors,
                         nitt = 60000, burnin = 10000)

However, without more information on your analysis, I strongly suggest you plot your posteriors to have a look at the results and see if anything looks wrong. Have a look at the MCMCglmm package Course Notes for more info on how to set these priors (especially on what not to do in section 1.5 - you can also find more specific info on how to tune it to your model if it fits in the categories of the tutorial).

Upvotes: 0

Related Questions