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
--
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:[email protected]
www.christophermull.weebly.com
www.earth2ocean.org
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list - [email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/[email protected]/
