Dear Liam,
thank you very much for the response! This is really helpful for me, especially the part related to the appropriateness of using this approach. But it introduces more questions than solutions :) Regardless of that, I very appreciate your efforts to prepare and distribute posts on the phytools blog! Thank you for that.

You raised caveats regarding the approach I asked, but now I am puzzled what to do with the uncertainty of tips in the phylogeny...

My main aim is to compare the evolution of a continuous trait in regard to two regimes (aquatic and terrestrial). However, defining these two states is not a straightforward task since some species tend to be aquatic and terrestrial at the same time with tendencies towards one or the other state. All of them can be considered as "semi-aquatic" (but "semi-terrestrial" is ok too I'd say...). Likewise, few species lack information about habitat preferences. So, I wanted to include such uncertainty into the analysis. I am not able to make e.g. 100 maps and use each map as an input in model fitting (we did that way in the course in Barcelona in 2016 where you was the instructor) because I am doing an exhausting multivariate model fitting and computation is super long (+ failed approaches I performed as I am beginner)... Nevertheless, I came to the idea to include uncertainty of tip states into stochastic mapping, obtain a single tree and use it for further analysis. Though there is no explanation in the phytools package, I suspected that if one does following will get PP of tips being in particular state:
mtrees <- [set of stochastic maps]
XX <- describe.simmap (mtrees)
XX$tips ## this will give PP of each tip being in a particular state

If this is correct, then the stochastic mapping can be used not just to estimate ancestral states, but also PP of tips given a model and data. This involves some circularity, but otherwise, I do not know what are the values of XX$tips and how we can use/interpret them...Please, can you provide some comment for this?

Note that I set priors for 3 states, but obj <- densityMap (mtrees) returns a tree where tips are assigned either as "aquatic" or "terrestrial". So obj tree has two states what I actually needed, though I am not sure if this is the correct way of thinking... The distribution of two states is extremely similar to that if one would assign "semiaquatic" as "aquatic". However, there are few species estimated as "terrestrial" even I set the greatest prior to "semi-aquatic" state. I believe that this approach (if the rate of character evolution is relatively low as you pointed out) can be used as justification why particular tip are assigned to certain states. Imagine that one aims to test a bunch of models of continuous trait evolution allowing change of mode and tempo across two discrete ecological states, but ecological states for some tips are not certain or are unknown. Should we exclude species with uncertain states? Or, should we still guess a value of their states?

Please, can you provide some reproducible code that compares if marginal posterior probabilities and the joint likelihood of states are more-or-less similar (case when the rate of character evolution is relatively low as you pointed out) and if they are, how to obtain binary-colored tree after densityMap function? I tried later but I failed to accomplish it...

I am sorry for the long email.

On 2017-02-25 01:03, Liam J. Revell wrote:
Hi Marko.

This is possible to do, of course. We could just traverse the tree,
find the state with the highest posterior probability at each node,
and paint the edge by that state using phytools::paintBranches. We
could even (with a bit more difficult) traverse our "densityMap"
object and find the state with the highest PP at the midpoint of each

However, I might counsel against it for the following reasons.

Firstly, doing so (that is, identifying a set of values at nodes with
highest posterior probability at each node) ignores the substantial
uncertainty that can be associated with ancestral character
estimation. For instance, if we imagine that a set of internal nodes
all have posterior probability of state "A" more or less equal to but
slightly greater than 0.5 (for a binary character), we would
'reconstruct' all of these nodes to be in state "A" - but this would
actually be virtually meaningless with respect to our confidence that
these nodes are indeed in state "A".

Secondly, to the extent that we *do* want one set of states for all
internal nodes in the tree, a reasonable choice for these states would
be those that maximized the probability of the data (that is, the
Maximum Likelihood states). Unfortunately, these are not guaranteed to
be the states that have the highest marginal posterior probability at
all the nodes of the tree. (This is explained in various places,
including in the book 'Computational Molecular Evolution' by Yang.)

On the other hand, if the rate of character evolution for our
character is relatively low then our marginal posterior probabilities
and joint likelihood states will likely more-or-less coincide, and
these are also probably quite similar to our MP states. In that case
it would perhaps be reasonable to summarize our reconstructed
character histories in the way you've proposed.

I hope these comments are helpful - and thank you for using phytools!

All the best, Liam

Liam J. Revell, Associate Professor of Biology
University of Massachusetts Boston

On 2/24/2017 6:38 PM, wrote:
Dear r-sig-phylo participants,

I have data set where tips are assigned to 3 discrete states (aquatic,
semi-aquatic and terrestrial). Because the definition of "semi-aquatic" is quite arbitrary and some species in the tree lack field observations
of their ecology, I decided to use a matrix of state priors for
stochastic mapping in the phytools package. That approach will allow me
to account for uncertainties/lack of information for some tips in the
phylogeny. I fitted 3 models (ER, SYM, and ARD) where Q was empirically
estimated and nsim was set to 1000. According to AIC value, the SYM
model was the best-fitted one. Describe.simmap showed that mean total
time spent in the state "semi-aquatic" was 0. Thus, all mapped trees
were actually binary and I was able to employ the densityMap function to obtain an object, let's say "obj", which contains a single tree with the
posterior density for the "aquatic" and "terrestrial" states from 1000
stochastic maps.

Here is the question:
Is it straightforward idea to paint back the obj$tree with just two
colors where colors are determined by a threshold value of posterior
probability (indicated as a legend bar in the bottom left part of the
For instance, is it appropriate to paint a tree edges with a color A if
PP is lower or equal than 0.5 and color B if PP is greater than 0.5?

Some R packages for model fitting allow simmap tree as input, but if my question makes sense, it would be better to provide consensus tree from
n stochastic maps instead to use one stochastic map as input.

Thank you for your time.
Kind regards,

R-sig-phylo mailing list -
Searchable archive at

R-sig-phylo mailing list -
Searchable archive at

Reply via email to