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.
_______________________________________________
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