Created
September 2, 2012 15:35
-
-
Save sckott/3600585 to your computer and use it in GitHub Desktop.
Comparison of functions for comparative phylogenetics
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
> # Load | |
require(motmot); require(geiger); require(picante) | |
source("http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/phylosig/v0.3/phylosig.R") | |
source("http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/fastBM/v0.4/fastBM.R") | |
# Make tree | |
tree <- rcoal(10) | |
# Transform branch lengths | |
> system.time( replicate(1000, transformPhylo(tree, model = "lambda", lambda = 0.5)) ) # motmot | |
user system elapsed | |
1.757 0.004 1.762 | |
> system.time( replicate(1000, lambdaTree(tree, 0.9)) ) # geiger | |
user system elapsed | |
3.708 0.008 3.716 | |
> # motmot wins!!! | |
# Simulate trait evolution | |
system.time( replicate(1000, transformPhylo.sim(tree, model = "bm")) ) # motmot | |
user system elapsed | |
3.732 0.007 3.741 | |
> system.time( replicate(1000, rTraitCont(tree, model = "BM")) ) # ape | |
user system elapsed | |
0.312 0.009 0.321 | |
> system.time( replicate(1000, fastBM(tree)) ) # Revell | |
user system elapsed | |
1.315 0.005 1.320 | |
> # ape wins!!! | |
# Phylogenetically independent contrasts | |
trait <- rnorm(10) | |
names(trait) <- tree$tip.label | |
> system.time( replicate(10000, pic.motmot(trait, tree)$contr[,1]) ) # motmot | |
user system elapsed | |
3.062 0.007 3.070 | |
> system.time( replicate(10000, pic(trait, tree)) ) # ape | |
user system elapsed | |
2.846 0.007 2.853 | |
> # ape wins!!! | |
# Phylogenetic signal, Blomberg's K | |
> system.time( replicate(100, Kcalc(trait, tree)) ) # picante | |
user system elapsed | |
1.311 0.005 1.316 | |
> system.time( replicate(100, phylosig(tree, trait, method = "K")) ) # Revell | |
user system elapsed | |
0.201 0.000 0.202 | |
> # Liam Revell wins!!! | |
# Ancestral character state estimation | |
> system.time( replicate(100, ace(trait, tree)$ace) ) # ape | |
user system elapsed | |
4.988 0.018 5.007 | |
> system.time( replicate(100, getAncStates(trait, tree)) ) # geiger | |
user system elapsed | |
2.253 0.005 2.258 | |
> # geiger wins!!! |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment