[R-sig-ME] prior specification for multinomial model MCMCglmm
Jarrod Hadfield
j.hadfield at ed.ac.uk
Tue Apr 11 05:30:49 CEST 2017
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.
More information about the R-sig-mixed-models
mailing list