V. Ancestral Reconstruction with Maximum Likelihood on Discrete Traits in R

Primary Contact(s)	Created	Required Software
Rich Glor	9 March 2009	R
Example Datafile	Prerequisites
See Introduction	R for Phylogenetics parts I, II, III, and IV

Introduction|I. Getting Started|II. Tree Basics|III. Loading Character Data|IV. Testing Phylogenetic Signal|V. Ancestral Reconstruction|VI. Testing Patterns

One of the simplest things one can do with a phylogeny and character data is to reconstruct the history of character evolution by inferring ancestral character states. We can use the ace function of ape to reconstruct ancestral character states using maximum likelihood:

ace(micro, anolisComparativeTree, type="discrete") -> anolisMicro

When we try this, however, we get an error message:

Error in ace(micro, anolisComparative, type = "discrete") : 
  "phy" is not rooted AND fully dichotomous.

This error message is fairly explicit. Many models of character evolution do not permit polytomies (i.e., trees that are not comprised exclusively of dichotomous branching events). To confirm that our tree is not fully dichotomous, we can use the following command:

is.binary.tree(anolisComparativeTree)

The response - FALSE - confirms that are tree is not fully dichotomous (a fact that should have been obvious earlier given the fact that the number of internal nodes in our tree is less than taxon number - 1). To render our tree entirely dichotomous, we can use another R function from the ape package:

multi2di(anolisComparativeTree) -> anolisComparativeTree_resolved

Now we can retry our effort to calculate ancestral character states:

ace(micro, anolisComparativeTree, type="discrete") -> anolisMicro

If this function successfully runs to completion, we should be able to investigate the stored output by typing anolisMicro (see below). This object contains several important elements, including the overall likelihood score ($loglik), the inferred rate of evolution ($rates), the standard error ($se), a matrix indicating the probabilities of change among all possible states ($index.matrix), and large matrix comprised of the marginal likelihood of each character state at each node ($lik.anc).

$loglik 
[1] -78.05017

$rates
[1] 0.4646661

$se
[1] 0.0840395

$index.matrix
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   NA    1    1    1    1    1
[2,]    1   NA    1    1    1    1
[3,]    1    1   NA    1    1    1
[4,]    1    1    1   NA    1    1
[5,]    1    1    1    1   NA    1
[6,]    1    1    1    1    1   NA

$lik.anc
                 0            1
 [1,] 2.059730e-01 1.303059e-01
 [2,] 2.014837e-01 1.769795e-01
 [3,] 1.257452e-01 2.538911e-01
 [4,] 1.193377e-02 9.344931e-01
 [5,] 1.491693e-01 2.683564e-01

We can view the final element of this object graphically by plotting the marginal likelihoods of each state at each node using pie diagrams:

plot(anolisComparativeTree, show.node.label=FALSE, cex=0.4) 
nodelabels(pie=anolisMicro$lik.anc, cex=0.5)

NOTE: An important warning is necessary here. Although it has improved dramatically since first introduced, the ace function sometimes fails to properly reconstruct ancestral character states. This is certainly true when you observe negative marginal likelihood values in the lik.anc portion of your ace output. Some other error messages may be ignored, but negative marginal likelihoods should not be ignored.