Dear list,

I am using MCMCglmm to fit a binomial model to metadata on learning ability
in animals, where the main question is "what affects learning?", and I am
seeking reassurance that I have been sensible about model fitting before I
go on to assess and interpret the results. The response is the proportion
of animals in a group that solve a task, the random effects are the species
(provided as a phylogenetic tree of 9 primates), and the study (25), and I
have 7 fixed effect covariates: 5 continuous, 1 ordinal with four levels
(task_difficulty), and one categorical with three levels
(as.factor(task_opaque). The dataset is quite small, at 103 observations.

I tried to use start=list(QUASI=FALSE) to get overdispersed starting values
but the model would not converge. I have used three different priors for
the random effects to look at the effect of the prior on the on the
posterior and it appears, based on looking at plots and HPDinterval() for
the model$VCV and model$Sol, that all three priors give similar posteriors.
I have not been using MCMCglmm very long and would like some advice on
whether a) the prior specification is correct, b) the priors are varied
enough to check for sensitivity.

## prior specifications

# prior 1: low degree of belief


b <- 1e+08

B1_V <- diag(8)*b; B1_mu <- rep(0,8); B1_V; B1_mu # an 8x8 matrix for the
fixed effects prior because I have 6 continuous variables and one factor
with 3 levels: 6+(3-1)=8


prior1 <- list(R=list(V=diag(1), n=1,
nu=0.002), G=list(G1=list(V=diag(1), n=Gn1,
nu=0.002), G2=list(V=diag(1), n=Gn2, nu=0.002)), B=list(V=B1_V, mu=B1_mu))


# prior2: moderate degree of belief


prior2 <- list(R=list(V=diag(1), n=1, nu=0.01), G=list(G1=list(V=diag(1),
n=Gn1, nu=0.01), G2=list(V=diag(1), n=Gn2, nu=0.01)), B=list(V=B1_V,
mu=B1_mu))


# prior3: high degree of belief


prior3 <- list(R=list(V=diag(1), n=1, nu=0.01), G=list(G1=list(V=diag(1), n=
Gn1, nu=0.01), G2=list(V=diag(1), n=Gn2, nu=0.01)), B=list(V=B1_V, mu=B1_mu
))

## model specification

model <- MCMCglmm(cbind(n_solve, group_size-n_solve) ~ n_demo + trials +
task_difficulty + as.factor(task_opaque) + baseline + prop_male + prop_lifespan
- 1,

family="multinomial2", random= ~animal + study,

start=list(QUASI=TRUE),

                         prior=prior, data=subdata1,
pedigree=myreducedtree_1,

  thin=500, nitt=1000000, burnin=500000, verbose=F, pr=T, pl=T,
nodes="TIPS")


All three priors give very similar results. The pMCMC values change, as
expected, but the same fixed effects appear to be important. For the random
effects, the way the variance is divided up into the three components
(species-effect, study-effect, residual-variance-effect) also changes very
little under the three different priors. Is this something I should be
worried about, does it mean the priors I tried were not varied enough? Or
can I trust it as a genuine signal in the data? I have explored the shape
of these priors using dinvgamma() with the appropriate parameters but I am
not very familiar with the inverse Wishart and I want to make sure I am
providing proper priors so I have not experimented much with changing the
value of V to be anything but 1.


I was not sure if it is appropriate in this case to use this prior for the
fixed effects:

B = list(mu = c(0, 0), V = diag(1) * (1+pi^2/3))


I have not standardised the variance components. Do I need to add pi^2/3 to
the sum of the posterior variance estimates in order to work out the
contribution of each component to the total variance? When is this
recommended? I could only find a fixed effect model example in the MCMCglmm
course notes.

species.effect <- model$VCV[, "animal"]/(rowSums(model$VCV) + pi^2/3)

So far as I have been working it out as follows:

species.effect <- model$VCV[, "animal"]/(rowSums(model$VCV))

Can I interpret the proportion directly as the proportional contribution of
each component to the total variance?

While on the topic of the random effects, I have been wondering if I can do
better than modelling the effect of species and study separately, which is
what, to my understanding, is done when the random formula is ~ animal +
study. Because a single species will feature in more than one study, I
wondered if there was a nested structure that might be better at
representing the structure of the data. I have read up about the idh() and
us() structures, but I am still not sure whether either is appropriate in
my case.


While I am aware of the fact that the predict function is incomplete and
needs further testing, I would still like to make predictions for each
species on the scale of the response.

The following code does not seem to work for me:

predict(model, marginal= ~study, type = "response", interval="confidence")

I get the error message:

Error in M[, which(rm.v), ] <- 0 : incorrect number of subscripts

Is this because it is a binomial model or have I got something coded
wrongly?


Any help and/or advice will be very much appreciated. I hope I have
provided all the necessary information.

Thank you in advance for your time.


Kind regards,

Theoni



-- 

Theoni Photopoulou, Research Fellow
 Centre for Social Learning and Cognitive Evolution
Bute Medical Building, University of St Andrews, Scotland KY16 9TS, UK
tel. +44 1334 467230 / skype. theoni_p / email. tp14@st-andrews.ac.uk

The University of St Andrews is a charity registered in Scotland : No
SC013532
" Be silly. Be honest. Be kind " Ralph Waldo Emerson

	[[alternative HTML version deleted]]

