One final comment on this. Since the conditional scaled likelihoods of the subtree and the marginal ancestral state reconstruction are equivalent at the root node of the tree, a computationally simple (but slow) method for getting the marginal ancestral state reconstructions is to move the root to each node of the tree and recompute the conditional likelihoods using ace. (Described in Felsenstein 2004, p. 259.) This works only for a reversible model of evolution. I posted a function for this on my blog (although it is straightforward), here: http://blog.phytools.org/2013/03/a-little-more-on-ancestral-state.html.

There are much more efficient algorithms for this, and I believe that this is what is implemented in diversitree by Rich Fitzjohn: http://www.zoology.ubc.ca/prog/diversitree/.

And, to reiterate, you should not be using the conditional scaled likelihoods from ace as ancestral state estimates (except at the root).

All the best, Liam

Liam J. Revell, Assistant Professor of Biology
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://blog.phytools.org

On 3/20/2013 6:02 PM, Liam J. Revell wrote:
Just to follow up on ace(...,type="discrete"), which I believe is not
doing what many people think it is doing:
http://blog.phytools.org/2013/03/conditional-scaled-likelihoods-in-ace.html.


Please don't take this as a criticism of ace - but I believe that the
difference between conditional, marginal, and joint reconstructions is
not appreciated by many phylogeny method users.

One option to get what I believe will asymptotically converge to the
joint reconstructions would be to use the posterior frequencies from
stochastic mapping under the model. This can be done using make.simmap
as described in my blog post - but please make sure that you have the
latest phytools release (phytools>=0.2-27) if you want to do this, as
there was a bug in earlier versions of make.simmap.

Finally, something that I neglected to note in my post is that I believe
that Rich Fitzjohn's package diversitree can be used to obtain the joint
reconstructions under likelihood.

All the best, Liam

Liam J. Revell, Assistant Professor of Biology
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://blog.phytools.org

On 3/20/2013 4:50 PM, Liam J. Revell wrote:
Hi Andres.

To answer your question at its face value, the line of code you want is:

anc<-setNames(apply(ANC$lik.anc,1,function(x)
   names(which.max(x))),length(tree$tip.label)+1:tree$Nnode)

(The setNames ensures that your resultant vector gets the node numbers
for names.)

This message really deserves a more lengthy response for two reasons:

1) The scaled likelihoods returned by ace(...,type="discrete") are
actually conditional likelihoods of the subtree obtained via the pruning
algorithm of Felsenstein. These should generally not be used for
ancestral state estimation (we should use the marginal or joint
reconstructions instead).

2) Extracting only the most likely state as *the* ancestral state at the
node ignores the uncertainty of our estimate. Really we have a
probability distribution for the states at each node.

- Liam

Liam J. Revell, Assistant Professor of Biology
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://blog.phytools.org

On 3/20/2013 2:23 PM, Andrés Parada wrote:
Sorry if this was addressed before.

I reconstructed a binary discrete character with* ace* (coded as 0-1).
As suggested in APER 2nd edition with the* Sylvia *case I can obtain the
"lik.anc" values but I wanted to derive the ancestral states from this
matrix.

*First*

*ANC <- ace (data$character, tree, type="d", method="ML", model="ARD")*

ANC$lik.anc look like this:

                  0            1
   [1,] 9.999185e-01 8.147752e-05

*So I applied the same strategy as in the example mentioned above:*

*anc <- apply(ANC$lik.anc, 1, which.max)*


This returned a matrix of with values of 1 or 2.

Is there a way to modify the formula and sort this so it returns only
0 or
1 depending on which state has higher probability?

Thank you very much for your help.

    [[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/


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


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


_______________________________________________
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