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

[Alberto Gallano]

Hi Jarrod,

this response was incredibly helpful. I decided ultimately that I didn't
need the random slopes in the model, which simplified things.

I have a follow-up question. I read the Phillimore et al. PNAS paper you
linked to with great interest. It seems in that paper you fitted a
multivariate (i.e., multiple response) mixed model as an alternative way
(from the van der Pol & Wright method) of decomposing between- and
within-species slopes (your Equation 4). The advantage, if I understand
correctly, is that the species-mean effect and individual-level deviations
from it are estimated, rather than calculated from the sample data. I'm
having a hard time figuring out how to alter the model I specified to fit
this framework (or if this is even possible).

My question is, how would I reconfigure the formula in this model to fit a
multivariate model that provides an estimate of the between species effect
(that's all i'm interested in)?

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

Alberto


On Fri, Dec 26, 2014 at 1:50 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk>
wrote:

> 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:
>
>  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]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>>
>
> --
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
>
>
>

	[[alternative HTML version deleted]]