The AnthroTree Website & "The Comparative Approach in Evolutionary Anthropology and Biology" by C. Nunn
HistoryEdited Feb 21, 2012 2:44 PM by admin…

4.1 Maximum Likelihood Reconstruction for Continuous Characters in R

Table of Contents - Chapter Overview

Previous Section - Next Section


ML Reconstruction for Continuous Characters in R


In this section, we will use the 'ape' package in R to reconstruct ancestral states using the example of the intermembral index in apes from the book. The procedure is similar to that in Section 3.4. The data can be found in these files:

Download file "anthrotree41_tree.txt"

Download file "anthrotree41_imi.csv"


If you haven't done so already, download R. Open R and navigate to the folder that contains the datafile from above (e.g., on a Mac, go to the Misc menu and change the working directory, or use the setwd function). Load the 'ape' library (see Section 1.1.2 for installation instructions):

library("ape")


Now, let's load the data and phylogeny:

imi_data = read.csv("anthrotree41_imi.csv")
imi_tree = read.nexus("anthrotree41_tree.txt")


Check out the data by typing "imi_data". You should see the following output on the command line:

                 Species   IMI
1        Gorilla_gorilla 115.8
2           Homo_sapiens  72.0
3       Hylobates_agilis 121.0
4     Hylobates_concolor 140.0
5      Hylobates_hoolock 129.0
6      Hylobates_klossii 126.0
7          Hylobates_lar 129.7
8   Hylobates_leucogenys 140.0
9       Hylobates_moloch 127.0
10    Hylobates_muelleri 129.0
11 Hylobates_syndactylus 147.0
12          Pan_paniscus 102.7
13       Pan_troglodytes 104.5
14        Pongo_pygmaeus 139.0

Similarly, take a look at the tree by typing:

plot(imi_tree)

Notice that some of the branches represent polytomies (set to be zero branch length in this file), with the phylogenetic information taken from Bininda-Emonds et al. (2007, 2008). Let's also be sure that the species in the tree and data files are identical by typing "imi_tree$tip.label". You should see something like this (depending on the width of your R console screen):

 [1] "Gorilla_gorilla"       "Homo_sapiens"         
 [3] "Hylobates_agilis"      "Hylobates_concolor"   
 [5] "Hylobates_hoolock"     "Hylobates_klossii"    
 [7] "Hylobates_lar"         "Hylobates_leucogenys" 
 [9] "Hylobates_moloch"      "Hylobates_muelleri"   
[11] "Hylobates_syndactylus" "Pan_paniscus"         
[13] "Pan_troglodytes"       "Pongo_pygmaeus"


Reading these in order (left to right, starting with Gorilla), you should see that they are in the same order as the list of taxa in 'imi_data.csv'. So, we are good to go with the analysis.


Now, we can use the ace function to obtain estimates of ancestral states, using maximum likelihood, and plot those on the tree:

ace_call = ace(imi_data$IMI, imi_tree, type="continuous")
nodelabels(round(ace_call$ace, digits=1))

You should see the following:


Let's compare that to Figure 4.1 in the book, which was generated using ANCML (Schluter et al. 1997)...


It appears that the root node reconstruction (and others) differs. ANCML gave a root node value of 122.4, while ace is giving 124.6. Not a huge difference, but disturbing nonetheless.


Let's check out another package to see what results it gives. Specifically, the package 'geiger' (Harmon et al. 2008) has a function called 'getAncStates', which can also produce ancestral state estimates. The help function says the following about the function:

This function uses an algorithm to calculate Maximum-likelihood estimates of ancestral character states for continuous characters.

The ape function ace carries out a similar function but finds estimates by actually maximizing the likelihood across all ancestors simultaneously.

I have implemented an algorithm that is faster and sometimes more reliable especially for large trees.


Let's give it a try: 
library("geiger")
gas_call = getAncStates(imi_data$IMI, imi_tree)
plot(imi_tree)
nodelabels(round(gas_call, digits=1))


You should see the following in your graphics window:

These results are identical to what is found with ANCML and shown in the book in Figure 4.1, as reproduced above. (If you want to try ANCML, go to this page).

Clearly, different programs are producing slightly different answers, and it is worth carefully checking and cross-checking the results.


References


Bininda-Emonds, O. R. P., M. Cardillo, K. E. Jones, R. D. E. MacPhee, R. M. D. Beck, R. Grenyer, S. A. Price, R. A. Vos, J. L. Gittleman, and A. Purvis. 2007. The delayed rise of present-day mammals. Nature 446:507-512.

Bininda-Emonds, O. R. P., M. Cardillo, K. E. Jones, R. D. E. MacPhee, R. M. D. Beck, R. Grenyer, S. A. Price, R. A. Vos, J. L. Gittleman, and A. Purvis. . 2008. Corrigendum: The delayed rise of present-day mammals. Nature 456:274.

Harmon, L., J. Weir, C. Brock, R. Glor, W. Challenger, and G. Hunt. 2008. geiger: analysis of evolutionary diversification. R package version:1.2-13.

Schluter, D., T. Price, A. O. Mooers, and D. Ludwig. 1997. Likelihood of ancestor states in adaptive radiation. Evolution 51:1699-1711.


Contributed by Charles Nunn



Previous Section - Next Section

Table of Contents - Chapter Overview



Comments

/groups/pica/search/index.rss?tag=hotlist/groups/pica/search/?tag=hotWhat’s HotHotListHot!?tag=hot0/groups/pica/sidebar/HotListNo items tagged with hot.hot/groups/pica/search/index.rss?sort=modifiedDate&kind=all&sortDirection=reverse&excludePages=wiki/welcomelist/groups/pica/search/?sort=modifiedDate&kind=all&sortDirection=reverse&excludePages=wiki/welcomeRecent ChangesRecentChangesListUpdates?sort=modifiedDate&kind=all&sortDirection=reverse&excludePages=wiki/welcome0/groups/pica/sidebar/RecentChangesListmodifiedDateallRecent ChangesRecentChangesListUpdateswiki/welcomeNo recent changes.reverse5search