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>
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
www.christophermull.weebly.com
www.earth2ocean.org
[[alternative HTML version deleted]]
_______________________________________________
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/
