## Comparing the Statistical Fit of Different Evolutionary Models in R

library("geiger")

tree_primates = read.nexus("57.Tree.txt")

data_primates = read.table("57.MassData.txt",header=T,sep="\t")

We are now ready to estimate the likelihood of the data on this tree given a variety of evolutionary models. For example, we can test for a directional trend in body size by comparing the likelihood of the "trend" model to the likelihood of the "BM" (Brownian motion) model. Brownian motion is the simplest (fewest parameter) model of trait evolution in which variance in a continuous trait accumulates at a constant rate in random directions (addition to or subtraction from the trait mean). Because of this constant variance assumption, it is common to logarithmically transform body size data so that proportional rather than absolute changes in size are modeled on the tree. You can very easily transform your data in R as follows (note that lines starting with # are ignored by R).

#convert kilograms to gramsdata_primates$Mass_gm = data_primates$Mass_kg*1000#take the natural logarithm of gramsdata_primates$Log_Mass_gm = log(data_primates$Mass_gm)

Take a look at the two new columns of data you added to your dataframe:

data_primateshist(data_primates$Log_Mass_gm)

You can see that the structure of the data table is preserved but the weights are all now the natural logarithm of grams. Now calculate the likelihood of the data under Brownian motion by typing:

brownian = fitContinuous(tree_primates, data_primates$Log_Mass_gm,model="BM")

Then call "brownian". R will print the following

$Trait1$Trait1$lnl[1] -30.98123$Trait1$beta[1] 0.02371125$Trait1$aic[1] 65.96246$Trait1$aicc[1] 66.27825$Trait1$k[1] 2

R is first saying these are the values from the table for trait 1 (there is only one trait in this case). After the second $ sign is says which parameter or result is printed on the next line. So, the "$Trait1$lnl" means the next line shows the log likelihood of the data under Brownian motion, which is -30.98. We can calculate the comparable likelihood under model of Brownian motion with a directional trend by typing:

trend = fitContinuous(tree_primates, data_primates$Log_Mass_gm,model="trend")

If you type "trend," R will return to the command line:

$Trait1 $Trait1$lnl [1] -30.86872 $Trait1$mean [1] 9.867361 $Trait1$beta [1] 0.02358446 $Trait1$mu [1] -0.05444671 $Trait1$aic [1] 67.73744 $Trait1$aicc [1] 68.38609 $Trait1$k [1] 3

pchisq(0.22,1,lower.tail=F)

[1] 0.6390399

lambda = fitContinuous(tree_primates, data_primates$Log_Mass_gm,model="lambda")

$Trait1$lambda [1] 0.811984

OU = fitContinuous(tree_primates, data_primates$Log_Mass_gm,model="OU")

## References

Harmon, L., J. Weir, C. Brock, R. Glor, W. Challenger, and G. Hunt (2009). geiger: Analysis of evolutionary diversification. R package version 1.3-1. http://CRAN.R-project.org/package=geiger

R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for

Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.