Hi Chris,

OK - stick with the RAM model, the h2 is so high you will run into numerical issues otherwise. In the two-trait model you might want to add in us(at.level(trait,1)):units into the random effects (make sure it is not the last term in the random formula) in case log.dep has a h2 substantially less than 1. Having a multi-level response will help with power so I would go for it. threshBayes does handle ordinal responses but you would probably have to run it for a VERY long time to sample the posterior adequately.


Cheers,

Jarrod

On 16/12/2016 07:11, Chris Mull wrote:
Hi Jarrod,
I hadn't appreciated that the clustering of reproductive modes on the tree might limit out ability to detect some of these relationships. This is in fact a step in testing reproduction as an ordinal variable (egg-laying, lecithotrophic live-bearing, and matrotrophic live-bearing) which follows gradients of depth and latitude, and threshBayes cannot handle ordinal variables. Perhaps treating the data this way will help given more transitions. I have run the model in MCMCglmm and have appended the summary and attached the histogram of the liabilities. Thanks so much for your help with this...

summary(dep2)
Iterations = 3001:12991
 Thinning interval  = 10
 Sample size  = 1000

 DIC: 31.2585

 G-structure:  ~animal

       post.mean l-95% CI u-95% CI eff.samp
animal     82.18    35.88    140.1    6.266

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: parity ~ log.med.depth

              post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept)      0.4250 -13.5697  13.7913    28.54 0.946
log.med.depth   -0.3601  -4.4399   3.8022    16.48 0.862

On Thu, Dec 15, 2016 at 11:10 PM, Jarrod Hadfield <j.hadfi...@ed.ac.uk <mailto:j.hadfi...@ed.ac.uk>> wrote:

    Hi Chris,

    I think MCMCglmm is probably giving you the right answer. There
    are huge chunks of the phylogeny that are either egg-laying and
    live-bearing. The non-phylogenetic model shows a strong
    relationship between reproductive mode and depth, and that might
    be causal or it might just be because certain taxa tend to live at
    greater depths and *happen to have* the same reproductive mode.
    There's not much information in the phylogenetic spread of
    reproductive modes to distinguish between these hypotheses, hence
    the wide confidence intervals.   Just to be sure can you

    a) just perform independent contrasts (not really suitable for
    binary data, but probably won't give you an answer far off the
    truth and is a nice simple sanity check).

    b) using MCMCglmm (not MCMCglmmRAM) fit

    prior.dep2<-=list(R=list(V=diag(1), fix=1),
    G=list(G1=list(V=diag(1), nu=0.002)))

    dep2<-MCMCglmm(parity~log.med.depth,
                           random=~animal,
                           rcov=~units,
                           pedigree=shark.tree,
                           data=traits,
                           prior=prior.dep2,
                           pr=TRUE,
                           pl=TRUE,
                           family="threshold")

    an send me the summary and hist(dep2$Liab)

    Cheers,

    Jarrod



    On 16/12/2016 07:02, Jarrod Hadfield wrote:

    Hi Chris,

    I think ngen in threshbayes is not the number of full iterations
    (i.e. a full update of all parameters), but the number of full
    iterations multiplied by the number of nodes (2n-1). With n=600
    species this means threshbayes has only really done about 8,000
    iterations (i.e. about 1000X less than MCMCglmm). Simulations
    suggest threshbayes is about half as efficient per full iteration
    as MCMCglmm which means that it may have only collected
    0.5*1132/1200 = 0.47 effective samples from the posterior. The
    very extreme value and the surprisingly tight credible intervals
    (+/-0.007) also suggest a problem.

    **However**, the low effective sample size for the covariance in
    the MCMglmm model is also worrying given the length of chain, and
    may point to potential problems. Are egg-laying/live-bearing
    scattered throughout the tree, or do they tend to aggregate a
    lot? Can you send me the output from:

    prior.dep<-=list(R=list(V=diag(1)*1e-15, fix=1),
    G=list(G1=list(V=diag(1), fix=1)))

    dep<-MCMCglmm(parity~log.med.depth,
                           random=~animal,
                           rcov=~units,
                           pedigree=shark.tree,
                           reduced=TRUE,
                           data=traits,
                           prior=prior2,
                           pr=TRUE,
                           pl=TRUE,
                           family="threshold")

    summary(dep)

    summary(glm(parity~log.med.depth, data=traits,
    family=binomial(link=probit)))

    Cheers,

    Jarrod



    On 15/12/2016 20:59, Chris Mull wrote:
    Hi All,
    I am trying to look at the correlated evolution of traits using the
    threshold model as implemented in phytools::threshBayes (Revell 2014) and
    MCMCglmmRAM (Hadfield 2015). My understanding from Hadfield 2015 is that
    the reduced animal models should yeild equivalent results, yet having run
    both I am having trouble reconciling the results. I have a tree covering
    ~600 species of sharks. skates, and rays and am interested in testing for
    the correlated evolution between reproductive mode (egg-laying and
    live-bearing) with depth. When I test for this correlation using
    phytools:threshBayes there is a clear result with a high correlation
    coefficient (-0.986) as I would predict. When I run the same analysis using
    MCMCglmmRAM I get a very different result, with a low degree of covariation
    and CI crossing zero (-0.374; -3.638 - 3.163 95%CI). Both models are run
    for 10,000,000 generations with the same bunr-in and sampling period. I
    have looked at the trace plots for each model using coda and parameter
    estimates and likelihodd . I'm not sure how to reconcile the differences in
    these results and any advice would be greatly appreciated. I have appended
    the code and outputs below.


    #######################
    #phytools::threshBayes#
    #######################
    X<-shark.data[c("parity","log.med.depth")]
    X<-as.matrix(X)

    #mcmc paramters (also see control options)
    ngen<-10000000
    burnin<-0.2*ngen
    sample<-500

    thresh<-threshBayes(shark.tree,
                        X,
                        types=c("discrete","continuous"),
                        ngen=ngen,
                        control = list(sample=sample))
    The return correlation is -0.986 (-0.99 - -0.976  95%HPD)


    #############################
    #MCMCglmm bivariate-gaussian#
    #############################

    prior2=list(R=list(V=diag(2)*1e-15, fix=1), G=list(G1=list(V=diag(2),
    nu=0.002, fix=2)))
    ellb.log.dep<-MCMCglmm(cbind(log.med.depth,parity)~trait-1,
                           random=~us(trait):animal,
                           rcov=~us(trait):units,
                           pedigree=shark.tree,
                           reduced=TRUE,
                           data=traits,
                           prior=prior2,
                           pr=TRUE,
                           pl=TRUE,
                           family=c("gaussian", "threshold"),
                           thin=500,
                           burnin = 1000000,
                           nitt = 10000000)

    summary(ellb.log.dep)
    Iterations = 1000001:9999501
    Thinning interval  = 500
    Sample size  = 18000
    DIC: 2930.751
    G-structure:  ~us(trait):animal
                                             post.mean l-95% CI u-95% CI 
eff.samp
    traitscale.depth:traitscale.depth.animal    16.965   15.092   18.864
      18000
    traitparity:traitscale.depth.animal         -0.374   -3.638    3.163
    1132
    traitscale.depth:traitparity.animal         -0.374   -3.638    3.163
    1132
    traitparity:traitparity.animal               1.000    1.000    1.000
      0
    R-structure:  ~us(trait):units
                                            post.mean l-95% CI u-95% CI eff.samp
    traitscale.depth:traitscale.depth.units    16.965   15.092   18.864    18000
    traitparity:traitscale.depth.units         -0.374   -3.638    3.163     1132
    traitscale.depth:traitparity.units         -0.374   -3.638    3.163     1132
    traitparity:traitparity.units               1.000    1.000    1.000        0
    Location effects: cbind(scale.depth, parity) ~ trait - 1
                     post.mean l-95% CI u-95% CI eff.samp pMCMC
    traitscale.depth   0.12297 -3.63655  4.02005    18000 0.949
    traitparity       -0.02212 -1.00727  0.93387    17058 0.971

    Thanks for any and all input.

    Cheers,
    Chris

    The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
--
Christopher Mull PhD Candidate, Shark Biology and Evolutionary Neuroecology Dulvy Lab Simon Fraser University Burnaby,B.C. V5A 1S6 Canada (778) 782-3989
twitter: @mrsharkbrain e-mail:cm...@sfu.ca <mailto:e-mail%3acm...@sfu.ca>
www.christophermull.weebly.com <http://www.christophermull.weebly.com>
www.earth2ocean.org <http://www.earth2ocean.org>
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Reply via email to