[R-sig-ME] prior specification for multinomial model MCMCglmm

[Jarrod Hadfield]

Hi Vicenzo,

This is a difficult problem.

First, you should drop the us(trait):Family.Genus.Species term: each 
species is only represented once and so this is confounded with your 
rcov specification.

My feeling is that the priors for the two random effect covariance 
matrices should be of the form V = (1/k)*(I + J) too as this would 
represent equal between-group variance in the three outcomes and also 
independence among the three outcomes (after accounting for the outcomes 
being negatively correlated because if you are a cavity nester you can't 
be an open cup nester). As far as I am aware there is still no treatment 
of parameter expanded priors for covariance matrices, but my guess is 
that they don't work well in this instance because they would shrink the 
covariances to zero rather than 2/k.

However, even with reasonable priors MCMCglmm fails on this data set 
because of numerical under/overflow. The between-Genus and 
between-Family variances are so large you end up with extreme category 
problems that are not well controlled by the variance components. You 
can see this by having pl=TRUE in the call to MCMCglmm and then 
obtaining the probability of the two non-base-line categories for each 
observation. If these probabilities are less than 1e-11 or greater than 
1-1e-11 (this is probably not stringent enough) then numerical problems 
are likely. With nu=1.002, probabilities this extreme occur often:

prob<-exp(mod$Liab[,1:(ncol(mod$Liab)/2)])/(exp(mod$Liab[,1:(ncol(mod$Liab)/2)])+exp(mod$Liab[,ncol(mod$Liab)/2+1:(ncol(mod$Liab)/2)]))

prop.table(table(prob<1e-11 | prob>(1-1e-11)))

In summary, I don't think MCMCglmm is going to give robust answers on 
this data-set. Sorry!

Cheers,

Jarrod








On 06/04/2017 20:44, Vincenzo Ellis wrote:
> library(gsheet)
> dat <- gsheet2tbl('
> https://docs.google.com/spreadsheets/d/1qVABAjNSBIlsoG3FYjREejBKyGkJvjaPiIJ_T_bzO_Y/edit?usp=sharing
> ')
>
> ## group taxonomic categories into unique categories for MCMCglmm for use
> as nested random effects
> dat$Family.Genus <- paste0(dat$Family, dat$Genus)
> dat$Family.Genus.Species <- paste0(dat$Family, dat$Genus, dat$Species)
>
> ## make variables factors...they come out as characters in the download
> dat <- as.data.frame(unclass(dat))
>
> ## set up model
> library(MCMCglmm)
>
> ## k = number of categories in response variable
> k <- length(levels(dat$Nest.Type))
>
> ## I and J are matrices that will set up constraints on the residuals of
> the model
> I <- diag(k-1)
> J <- matrix(1, k-1, k-1)
>
> ## set up prior; how to set up the G structure?
> prior1 <- list(R = list(V = (1/k)*(I + J), fix = 1))
>
> ## model
> mod <- MCMCglmm(Nest.Type ~ trait - 1, rcov = ~us(trait):units,
>                  random = ~us(trait):Family + us(trait):Family.Genus +
> us(trait):Family.Genus.Species,
>                  prior = prior1, data = dat, family = "categorical")


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.