Skip to content

MCMCglmm priors and random effects for phylogenetic mixed model

2 messages · Alberto Gallano, Jarrod Hadfield

Original
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
4 days later
Jarrod Hadfield
# 
Hi Alberto,

1/ You can fit it at both levels if you like:

us(1+I2.M1.within.species):phylo + us(1+I2.M1.within.species):species

where the first term models between-species phylogentically correlated  
variation in intercepts and slopes  and the second term models  
between-species variation in intercepts and slopes that is not  
phylogenetically correlated. However, you have to be quite careful  
with the van der Pol & Wright method because measurement error in the  
species means (which will be high when n=3) can appear as random  
variation in slopes. The section "Within-Population Slope  
Heterogeneity" in this paper:

http://www.pnas.org/content/107/18/8292.full

discusses the problem, but unfortunately without a good resolution.

2/ I generally use parameter expanded priors of the form:

G2 = list(V = diag(2), nu = 2, alpha.mu = rep(0, 2), alpha.V= diag(500, 2, 2))

note V is an identity matrix rather than I*0.5. However, you should  
check to make sure you don't get big changes when you use other types  
of prior.

Cheers,

Jarrod







Quoting Alberto Gallano <alberto.gc8 at gmail.com> on Sun, 21 Dec 2014  
17:59:09 -0500: