[R-sig-ME] Fwd: Re: [R-sig-phylo] How to incorporate intraspecific variation in MCMCglmm

[Jarrod Hadfield]

cc-ed back to the list....

Some observations seem to have no sampling variation and so have a 
residual variance of zero - which is not allowed, and doesn't make much 
sense.

Cheers,

Jarrod



-------- Forwarded Message --------
Subject: 	Re: [R-sig-phylo] How to incorporate intraspecific variation 
in MCMCglmm
Date: 	Sat, 15 Jul 2017 23:16:12 +0000
From: 	Diogo B. Provete <dbprovete at gmail.com>
To: 	Jarrod Hadfield <j.hadfield at ed.ac.uk>



Hi Jarrod,
I have tried all the suggestions for specifying the random effects. 
However, your last suggestion of making the residual covariance seems 
more appropriated in my case, since I have the standard error calculated 
already (in terms of jumped distance) and it's interpreted as 
intraspecific variance due to sampling error. But when I make my model as:

model1<-MCMCglmm(mean_distance~type_arena*type_of_stimulus,
*random=~Species, rcov=~idh(units):units*, data=df_spe, 
family="gaussian",  ginverse = list(Especie=treeAinv), nodes="ALL", 
prior=ve_priors, nitt=300000, burnin=25000, thin = 100, verbose=FALSE)

with priors:

ve_priors = list(R = list(V = diag(df_spe$se_distance^2), fix=1),
                  G=list(G1=list(V=1, nu=0.02)))

I get the following error message:

Error in priorformat(if (NOpriorG) { :
   V is not positive definite for some prior$G/prior$R elements .

I tried to google it and looked at your course notes ch. 3 and 4, but 
couldn't find anything helpful to understand it.

I'm attaching the first few lines of the prior V matrix.

Thanks once again,
Diogo

Em sex, 14 de jul de 2017 às 18:47, Diogo B. Provete 
<dbprovete at gmail.com <mailto:dbprovete at gmail.com>> escreveu:

    Dear Jarrod,
    Thanks very much for the quick reply.

    I'll try to implement the changes in the model.

    Have a nice weekend,
    Diogo


    Em Sex, 14 de jul de 2017 17:48, Jarrod Hadfield
    <j.hadfield at ed.ac.uk <mailto:j.hadfield at ed.ac.uk>> escreveu:

        Hi Diogo,

        First, your model1 is unlikely to be valid unless the residual
        variance
        happens to be 1. You should not fix it at one, and use a prior like:

        prior = list(R = list(V = 1, nu = 0.02), G=list(G1=list(V=1,
        nu=0.02)))

        Note that the residual variance (Ve) is the intra-specific variance,
        assumed constant over species, as in Lynch's h2/Pagels Lambda.
        If you
        have the standard error (se) of each observation you can also allow
        heterogeneity between observations as in random-effect
        meta-analysis:

        random=~Specie+idh(se)

        If you fix the variance associated with se at 1 in the prior

        prior = list(R = list(V = 1, nu = 0.02), G=list(G1=list(V=1,
        nu=0.02), G2=list(V=1, fix=1)))

        Then the intraspecific variance for species i is se^2_i+Ve

        If you do not fix the variance associated with se at 1, for
        example with

        prior = list(R = list(V = 1, nu = 0.02), G=list(G1=list(V=1,
        nu=0.02), G2=list(V=1, nu=0.02)))

        and estimate the associated variance (Vse). Then the
        intraspecific variance
           for species i is Vse*se^2_i+Ve. This is useful if you only
        know the standard error up to proportionality.

        Alternatively you can fit fixed-effect meta-analysis where all
        the intraspecific variance is presumed to be due to sampling error:

        random=~Specie, rcov=~idh(units):units

        with prior:

        prior = list(R = list(V = diag(se^2), fix=1),
        G=list(G1=list(V=1, nu=0.02)))

        The intraspecific variance for species i is se^2_i. This is
        quite an inefficient way of doing fixed-effect meta-analysis,
        and one day I will make such analyses easier to fit, and quicker
        to fit......

        Also, don't use nodes="TIPS" in the call to inverseA. I only
        allowed this option because its the parameterisation used by
        most other software, but its a really bad way of doing it. Stick
        with the default, node="ALL".

        Cheers,

        Jarrod









        On 14/07/2017 14:26, Diogo B. Provete wrote:
         > treeAinv<-inverseA(phylo,nodes="TIPS",scale=TRUE)$Ainv
         >
         > I included the following priors:
         >
         > prior = list(R = list(V = 1, fix = 1), G=list(G1=list(V=1,
        nu=0.02)))


        --
        The University of Edinburgh is a charitable body, registered in
        Scotland, with registration number SC005336.

    -- 
    *Diogo B. Provete, PhD.*
    FAPESP Post-doctoral fellow
    Department of Environmental Sciences
    Centre for Sciences and Technologies for Sustainability
    Federal University of São Carlos
    Sorocaba

    Gothenburg Global Biodiversity Centre
    Box 462   SE-405 30
    Göteborg, Sverige

    diogoprovete.weebly.com <http://diogoprovete.weebly.com>

    Cell phone: +5515981022137 <tel:%2815%29%2098102-2137>
    Skype: diogoprovete

    Editor: Amphibia-Reptilia | Biodiversity Data Journal

-- 
*Diogo B. Provete, PhD.*
FAPESP Post-doctoral fellow
Department of Environmental Sciences
Centre for Sciences and Technologies for Sustainability
Federal University of São Carlos
Sorocaba

Gothenburg Global Biodiversity Centre
Box 462   SE-405 30
Göteborg, Sverige

diogoprovete.weebly.com <http://diogoprovete.weebly.com>

Cell phone: +5515981022137 <tel:%2815%29%2098102-2137>
Skype: diogoprovete

Editor: Amphibia-Reptilia | Biodiversity Data Journal
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20170716/475ee5f3/attachment.ksh>