[R-sig-ME] MCMCglmm priors and random effects for phylogenetic mixed model

[Alberto Gallano]

I have a question about prior and random effects specification for a
phylogenetic mixed model. I am fitting a linear mixed model using MCMCglmm,
accounting for phylogenetic dependence in the residuals. My fixed effects
are two continuous variables (ratios of central and lateral incisor width
to first molar length) in various species of animals (n = 106). I have
measurements for multiple individuals within each species (ranging from n=3
to n=110). I am mainly interested in the among-species slope and intercept.
Therefore, I have used the methods outlined in (van der Pol & Wright, 2009)
to split the explanatory variable into two: one with species means, the
other within-species centered. This disentangles effects for within- and
between species. For random effects, I am fitting random intercepts for the
species level phylogenetic random effect (called "phylo"), random
intercepts for the species level non-phylogenetic random effect (called
"species"), and random slopes for the within-species centered variable.

I have two questions:

1) I want to a model that allows random slopes to vary for each species,
but i'm not sure if I have this specified correctly? Also, should the
random slopes be coupled with the "species" or "phylo" random effect?

2) What is a good prior for the random slope / random intercept (here G2)?
I'm not sure whether this is specified correctly (e.g., should I used
parameter expanded priors)?

Here are the priors and model:

priors1 <- list(
    B = list(mu = rep(0, 3), V = diag(9, 3)),
    G = list(G1 = list(V = 1, nu = 0.002), G2 = list(V = diag(2)/2, nu =
0.002)),
    R = list(V = 1, nu = 0.002)
    )

# parameter expanded version
priors2 <- list(
    B = list(mu = rep(0, 3), V = diag(9, 3)),
    G = list(G1 = list(V = 1, nu = 0.002),
             G2 = list(V = diag(2)/2, nu = 2, alpha.mu = rep(0, 2), alpha.V
= diag(500, 2, 2))),
    R = list(V = 1, nu = 0.002)
    )

# inverse of sigma matrix of phylogenetic correlation
inv_phylo_mat <- inverseA(tree, nodes = "TIPS", scale = TRUE)

fit <- MCMCglmm(
    fixed = I1.M1 ~ I2.M1.species.mean + I2.M1.within.species,
    rcov = ~ units,
    random = ~ phylo + idh(1+I2.M1.within.species):species,
    data = incisor.dat,
    family = "gaussian",
    ginverse = list(phylo = inv_phylo_mat$Ainv),
    prior = priors1,
    pr = TRUE,
    pl = TRUE,
    nitt = 1.1e+6, thin = 10, burnin = 1e+5,
    verbose = FALSE
    )


best,
Alberto

	[[alternative HTML version deleted]]