Reproducible analysis code for our forthcoming Anolis roosevelti (aka. locate.yeti) Evolution paper

The following is the reproducible analysis code used in our 'in press' Evolution article in which we present a method for placing recently extinct taxa on a phylogeny using continuous character data and apply it to the case of Anolis roosevelti, a giant anole from the Puerto Rican bank Virgin Islands that is most likely extinct. The paper should be out in the first half of the year, but this will be posted to Dryad soon with our data.

------------------------------

Analyis & code for “Placing cryptic, recently extinct, or hypothesized taxa into an ultrametric phylogeny using continuous character data: A case study with the lizard Anolis roosevelti”

First load packages. Note that use of Rphylip also requires that the PHYLIP package be installed. This is used to compute the 'branch-score' distances of Kuhner & Felsenstein (1994).

set.seed(1)
library(phytools)
library(Rphylip)
library(phangorn)
library(clusterGeneration)

Next, here is our wrapper function for simulation & analysis:

foo<-function(N,m,nsim){
    cat("simulation conditions:\n")
    cat(paste("\tN =",N,"\n\tm =",m,"\n\tnsim=",nsim,"\n"))
    ## vectors for results
    rf<-bs<-vector()
    for(i in 1:nsim){
        ## simulate tree
        tt<-tree<-pbtree(n=N+1,tip.label=sample(c(paste("t",1:N,sep=""),"Yeti")),scale=1)
        if(m>1){
            ## generate a covariance matrix for simulation
            V<-genPositiveDefMat(m,covMethod="unifcorrmat")$Sigma
            X<-sim.corrs(tree,vcv=V)
        } else X<-as.matrix(fastBM(tree))
        ## prune unknown taxon "Yeti" from tree
        tree<-drop.tip(tree,"Yeti")
        ## run locate.yeti
        mltree<-locate.yeti(tree,X,plot=FALSE,search="exhaustive",quiet=TRUE)
        rf[i]<-RF.dist(tt,mltree)
        bs[i]<-Rtreedist(tt,trees2=mltree,quiet=TRUE)
        cat(".")
    }
    cat("\nDone.\n")
    return(list(rf=rf,bs=bs))
}

Now let's run our simulation showing the result for various numbers of taxa given a fixed number of traits:

m<-10 ## number of traits (for fixed number of traits)
nsim<-100 ## number of simulation
N<-seq(20,100,by=10)
objN<-lapply(N,foo,m=m,nsim=nsim)

## simulation conditions:
##  N = 20 
##  m = 10 
##  nsim= 100 
## ..............................................................

## Warning: running command 'rm outfile' had status 1

## ......................................
## Done.
## simulation conditions:
##  N = 30 
##  m = 10 
##  nsim= 100 
## ....................................................................................................
## Done.
## simulation conditions:
##  N = 40 
##  m = 10 
##  nsim= 100 
## .................................................................................

## Warning: running command 'rm outfile' had status 1

## ...................
## Done.
## simulation conditions:
##  N = 50 
##  m = 10 
##  nsim= 100 
## .............................................

## Warning: running command 'rm outfile' had status 1

## .......................................................
## Done.
## simulation conditions:
##  N = 60 
##  m = 10 
##  nsim= 100 
## ..........................................................................................

## Warning: running command 'rm intree' had status 1

## ..........
## Done.
## simulation conditions:
##  N = 70 
##  m = 10 
##  nsim= 100 
## .....................

## Warning: running command 'rm intree' had status 1

## ...............................................................................
## Done.
## simulation conditions:
##  N = 80 
##  m = 10 
##  nsim= 100 
## ....................................................................................................
## Done.
## simulation conditions:
##  N = 90 
##  m = 10 
##  nsim= 100 
## ............................

## Warning: running command 'rm intree2' had status 1

## ........................................................

## Warning: running command 'rm outfile' had status 1

## ................
## Done.
## simulation conditions:
##  N = 100 
##  m = 10 
##  nsim= 100 
## ..................................................................

## Warning: running command 'rm outfile' had status 1

## ..................................
## Done.

Next, we have a function to generate a null distribution of tree distances as if the missing taxon were attached randomly in our tree. Here is what that code looks like:

null.treedist<-function(N,nsim){
    tt<-pbtree(n=(N+1),nsim=nsim,tip.label=c(paste("t",1:N,sep=""),"Yeti"),scale=1)
    trees<-lapply(tt,drop.tip,tip="Yeti")
    trees<-lapply(trees,add.random,tips="Yeti")
    class(trees)<-"multiPhylo"
    rf<-mapply(RF.dist,tt,trees)
    bs<-Rtreedist(tt,trees2=trees,distances="corresponding",quiet=TRUE)
    list(rf=rf,bs=bs)
}

We can use it to recreate panels A & B of Figure 3, showing a comparison between locate.yeti and random placement of the missing tip:

## generate null trees
nullN<-lapply(N,null.treedist,nsim=nsim)
layout(matrix(c(1,2),1,2,byrow=TRUE))
par(mar=c(4.1,4.1,3.1,1.1))
## plot RF distances
nullRF<-sapply(nullN,function(x) x$rf)
colnames(nullRF)<-N
par(fg=grey(0.2,0.5))
boxplot(nullRF,col=grey(0.5,0.2),ylab="R-F distance",xlab="number of taxa (N-1)",ylim=c(0,40))
par(fg="black")
RF<-sapply(objN,function(x) x$rf)
colnames(RF)<-N
boxplot(RF,add=TRUE)
title("A                                                             ")
## plot branch-score distances
nullBS<-sapply(nullN,function(x) x$bs)
colnames(nullBS)<-N
par(fg=grey(0.2,0.5))
boxplot(nullBS,col=grey(0.5,0.2),ylab="B-S distance",xlab="number of taxa (N-1)",ylim=c(0,2.5))
par(fg="black")
BS<-sapply(objN,function(x) x$bs)
colnames(BS)<-N
boxplot(BS,add=TRUE)
title("B                                                             ")

Now, let's do the same analysis, but this time vary the number of traits:

set.seed(1)
N<-50 ## number of taxa (for fixed number of taxa)
m<-c(1,2,5,10,20) ## number of traits
objM<-lapply(m,foo,N=N,nsim=nsim)

## simulation conditions:
##  N = 50 
##  m = 1 
##  nsim= 100 
## ....................................................................................................
## Done.
## simulation conditions:
##  N = 50 
##  m = 2 
##  nsim= 100 
## ........................................................................

## Warning: running command 'rm intree2' had status 1

## ............................
## Done.
## simulation conditions:
##  N = 50 
##  m = 5 
##  nsim= 100 
## ....................................................................

## Warning: running command 'rm outfile' had status 1

## ................................
## Done.
## simulation conditions:
##  N = 50 
##  m = 10 
##  nsim= 100 
## ......................................................

## Warning: running command 'rm outfile' had status 1

## ..............................................
## Done.
## simulation conditions:
##  N = 50 
##  m = 20 
##  nsim= 100 
## ..............................................

## Warning: running command 'rm outfile' had status 1

## .................................

## Warning: running command 'rm outfile' had status 1

## .....................
## Done.

Now we can recreate panels C & D of Figure 3, showing a comparison between locate.yeti and random placement of the missing tip, for various numbers of traits:

layout(matrix(c(1,2),1,2,byrow=TRUE))
par(mar=c(4.1,4.1,3.1,1.1))
nullM<-nullN[[which(N==50)]]
nullRF<-matrix(rep(nullM$rf,length(m)),100,length(m),byrow=FALSE)
colnames(nullRF)<-m
par(fg=grey(0.2,0.5))
boxplot(nullRF,col=grey(0.5,0.2),ylab="R-F distance",xlab="number of characters (m)",ylim=c(0,40))
par(fg="black")
RF<-sapply(objM,function(x) x$rf)
colnames(RF)<-m
boxplot(RF,add=TRUE)
title("C                                                             ")
nullBS<-matrix(rep(nullM$bs,length(m)),100,length(m),byrow=FALSE)
colnames(nullBS)<-m
par(fg=grey(0.2,0.5))
boxplot(nullBS,col=grey(0.5,0.2),ylab="B-S distance",xlab="number of characters (m)",ylim=c(0,2.5))
par(fg="black")
BS<-sapply(objM,function(x) x$bs)
colnames(BS)<-m
boxplot(BS,add=TRUE)
title("D                                                             ")

tree<-read.tree("Revell-etal.tree.tre")
X<-read.csv("Revell-etal.data.csv",header=T,row.names=1)
X<-as.matrix(X)
X<-X[,1:20]
tip<-setdiff(rownames(X),tree$tip.label)
tip

## [1] "roosevelti"

mltree<-locate.yeti(tree,X,search="exhaustive",plot=TRUE)

## Optimizing the phylogenetic position of roosevelti using ML. Please wait....

## Done.

mltree$logL

## [1] 3392

Now we can create the contMap style plot of Figure 4.

mltree<-reorder(reorder(mltree,"pruningwise")) ## reorder hack
pca<-phyl.pca(mltree,X)
## flip PC 1 so that it loads positively (not negatively) on size
obj<-contMap(mltree,-pca$S[,1],plot=FALSE)
class(obj)<-"densityMap"
## fix tip label
obj$tree$tip.label[which(obj$tree$tip.label=="pumilis")]<-"pumilus"
## plot
plot(obj,type="fan",lwd=4,legend=0.7,
    leg.txt=c(round(obj$lims[1],3),"PC 1",round(obj$lims[2],3)),
    outline=TRUE)
## add arrows
add.arrow(tree=obj$tree,tip="roosevelti",col="red",lwd=3,hedl=0.06,angle=50)
add.arrow(tree=obj$tree,tip="cuvieri",col="blue",lwd=3,hedl=0.06,angle=50)

Finally, we can do the analyses to test the (null) hypotheses that A. roosevelti is siste to A. cuvieri, or sister-to or nested-within the majority of the other Puerto Rican anoles:

## find ML cuvieri-constraint tree
cuvieri<-which(tree$tip.label=="cuvieri")
cuvieri.tree<-locate.yeti(tree,X,search="exhaustive",constraint=cuvieri)

## Optimizing the phylogenetic position of roosevelti using ML. Please wait....
## Done.

## compute a LR to test the alternative hypothesis that the A. roosevelti
## is *not* sister
LR.cuvieri<-2*(mltree$logL-cuvieri.tree$logL)
## perform simulation for the null distribution of the LR
set.seed(1)
LR.null<-vector()
nsim<-100
obj<-phyl.vcv(X[cuvieri.tree$tip.label,],vcv(cuvieri.tree),lambda=1)
for(i in 1:nsim){
    Xsim<-sim.corrs(tree=cuvieri.tree,vcv=obj$R,anc=obj$alpha[,1])
    mltree.null<-locate.yeti(tree,Xsim,search="exhaustive",plot=FALSE,
        quiet=TRUE)
    cuvieri.tree.null<-locate.yeti(tree,Xsim,search="exhaustive",
        constraint=which(tree$tip.label=="cuvieri"),plot=FALSE,
        quiet=TRUE)
    LR.null[i]<-2*(mltree.null$logL-cuvieri.tree.null$logL)
}
P.cuvieri<-1-mean(LR.cuvieri>=LR.null)
## P-value of the LR test
P.cuvieri

## [1] 0.26

This suggests that we cannot reject the hull hypothesis that A. roosevelti is sister to A. cuvieri.

Now, let's do the same thing, but set our constraint tree to be one in which Anolis roosevelti must be found nested within or sister to the clade with most of Puerto Rican's remaining anole species:

sp<-c("cristatellus","cooki","poncensis","gundlachi","pulchellus","krugi",
    "stratulus","evermanni")
pr<-getDescendants(tree,findMRCA(tree,sp))
pr.tree<-locate.yeti(tree,X,search="exhaustive",constraint=pr)

## Optimizing the phylogenetic position of roosevelti using ML. Please wait....
## Done.

LR.pr<-2*(mltree$logL-pr.tree$logL)
LR.null<-vector()
nsim<-100
obj<-phyl.vcv(X[pr.tree$tip.label,],vcv(pr.tree),lambda=1)
for(i in 1:nsim){
    Xsim<-sim.corrs(tree=pr.tree,vcv=obj$R,anc=obj$alpha[,1])
    mltree.null<-locate.yeti(tree,Xsim,search="exhaustive",plot=FALSE,
        quiet=TRUE)
    pr.tree.null<-locate.yeti(tree,Xsim,search="exhaustive",
        constraint=getDescendants(tree,findMRCA(tree,sp)),
        plot=FALSE,quiet=TRUE)
    LR.null[i]<-2*(mltree.null$logL-pr.tree.null$logL)
}
P.pr<-1-mean(LR.pr>=LR.null)
P.pr

## [1] 0.02

Here, by contrast, we can confidently reject the hypothesis that A. roosevelti is part of the main clade of Puerto Rican anoles.

That's pretty much it.

Last updated Jan. 13, 2015.

New phytools version (phytools 0.4-45) submitted to CRAN

This afternoon I submitted a new version of phytools to CRAN. Assuming it is accepted, the new version (phytools 0.4-45), which is already available from the phytools page, should percolate through all of the CRAN mirror repositories over the next few days.

This version does not feature any significant updates over the most recent non-CRAN phytools release (although I did do some updating of the manual pages this morning before submitting); however the last CRAN update was August 26th of last year (phytools 0.4-31) so there are many changes and updates to phytools since that release.

Some updates and additions include the following (and I can't guarantee that this list is comprehensive):

(1) A faster version of the robust Newick tree reader read.newick (1, 2).

(2) A new function (ladderize.simmap) to ladderize a phylogeny with a mapped discrete character.

(3) An update to the discrete character history simulating function, sim.history, to handle errors more sensibly.

(4) A small update to the fancyTree method "phenogram95" to permit user control of shading.

(5) Replacement of the redundant phytools function repPhylo with a valid alias.

(6) A fix for a very mysterious bug in plot.contMap and plot.densityMap (not by coincidence - these two functions use the same internals).

(7) A wrapper function implementating Pagel's (1994) method to test for correlated evolution of two binary traits. (Not in phytools, but I posted code for simulating binary character correlated evolution here.)

(8) An update to locate.yeti implementing exact ML and REML estimation methods.

(9) User control of the vertical range of labels in phenogram. (Turns out this isn't particularly useful.)

(10) A fossil version of locate.yeti, called, creatively, locate.fossil.

(11) A new version of the phytools helper function setMap that can automatically invert the color map of an object of class "contMap" or "densityMap".

(12) A new function add.arrow that will add an arbitarily colored arrow to a radial or horizontal phylogram.

(13) A new option in plotSimmap (and thus plotTree) that permits arbitrary vertical spacing of tips in a right or left-facing square phylogram.

And, finally, (14) a bug fix for pbtree with user-supplied tip labels (& non-zero extinction).

As noted earlier, phytools 0.4-45 is already available from the phytools webpage, but, if all goes smoothly, should also be available from CRAN in the not-too-distant future.

Thanks for paying attention!

Bug fix for pbtree with user-supplied tip labels and extinction

I just posted a small bug fix in the phytools function pbtree.

pbtree (originally short for Pure-Birth tree) does birth-death phylogeny simulations for a variety of conditions (time-stop, taxon-stop, both together via rejection, and discrete & continuous time simulation). If the taxon-stop criterion is being used (that is, if the function is being told to return a tree with a specific number of extant terminal taxa), then the user is also allowed to supply his or her own tip labels which will be used to populated tree$tip.label in the returned tree. Since a stochastic birth-death tree with N extant tips can have an arbitrary and unknowable number of extinct tips, these are labeled X1, X2, etc. For instance:

library(phytools)
tree<-pbtree(n=26,b=1,d=0.3,tip.label=letters)

## Warning: only using labels in tip.label for extant tips.
##          extinct tips will be labeled X1, X2, etc.

plotTree(tree,ftype="i",fsize=0.9)

Unfortunately, the way this works is by first simulating the tree with arbitrary internal labels, and then assigning all N extant tips the user supplied labels after running the phytools function getExtant internally. Due to numerical precision issues, though, if the total tree length is very long (for instance, because b and d are small for a given taxon-stop, n), then (unless its tolerance argument is adjusted) getExtant will malfunction and pbtree will behave as if the entire tree consists of extinct taxa! So, for example:

set.seed(10)
tree<-pbtree(n=26,b=0.001,d=0.0003,tip.label=letters)

## Warning: only using labels in tip.label for extant tips.
##          extinct tips will be labeled X1, X2, etc.

plotTree(tree,ftype="i",fsize=0.9)

Clearly, 26 tips are extant in this tree; however the function is behaving as if they are all extinct!

The 'fix' that I've applied is to make the tolerance value in getExtant a function of the total tree length, rather the default value of getExtant.

detach("package:phytools",unload=TRUE)
install.packages("phytools_0.4-44.tar.gz",type="source")

## Installing package into 'C:/Users/Liam/Documents/R/win-library/3.1'
## (as 'lib' is unspecified)
## inferring 'repos = NULL' from 'pkgs'

library(phytools)
packageVersion("phytools")

## [1] '0.4.44'

set.seed(10)
tree<-pbtree(n=26,b=0.001,d=0.0003,tip.label=letters)

## Warning: only using labels in tip.label for extant tips.
##          extinct tips will be labeled X1, X2, etc.

plotTree(tree,ftype="i",fsize=0.9)

The updated version of phytools is here. I'm also planning to submit a new version to CRAN in the not-too-distant future. (Fingers crossed!)

That's it!

Arbitrary vertical spacing of the tips on a plotted tree

Today a R-sig-phylo reader asked the following:

I am looking for a way to manually specify varying vertical distance between the tips in a phylogeny (because i want to add several lines of text to some of the tips).

I.e. for the tree library(ape) TREE=read.tree(text=“((Tip_1:1,Tip_2:1):1,Tip_3:2);)”)

Is there a way for instance through plot phylo to print place the tips in the heights of c(1,2,6) as opposed to c(1,2,3) as they would be normally in plot.phylo.

This is not presently possible, but I realized immediately that this would be pretty straightforward to add to phytools plotSimmap because the way that functions works is by first assigning 1:N to the N tips of the tree, and then works backwards to assign the vertical position of all internal nodes via a post-order traversal.

I responded to the query, and now I have also posted code online that permits exactly this type of manipulation. The new phytools build with this update can be obtained here. Here's a demo using the code of the original query:

library(ape)
TREE=read.tree(text="((Tip_1:1,Tip_2:1):1,Tip_3:2);")
library(phytools)
packageVersion("phytools")

## [1] '0.4.43'

## normal vertical spacing
plotTree(TREE)

## modified spacing
tips<-setNames(c(1,2,6),TREE$tip.label)
tips

## Tip_1 Tip_2 Tip_3 
##     1     2     6

plotTree(TREE,tips=tips)

Note that there is nothing here to prevent us from plotting trees with line crossing! For instance:

tree<-pbtree(n=26,tip.label=LETTERS)
tips<-setNames(1:26,sample(LETTERS))
plotTree(tree,tips=tips)

Yikes!

We can also do other weird stuff, like this:

plotTree(tree,tips=setNames(log(1:26),tree$tip.label),ftype="off")

plotTree(tree,tips=setNames((1:26)^2,tree$tip.label),ftype="off")

Well, you get the idea.

That's it.

Wrapper function to optimize the λ tree transformation for a discrete trait

Today I responded to an R-sig-phylo query by posting some code to optimize Pagel's λ tree transformation for a discrete character evolving by a continuous-time Markov chain. I did this by writing a very simple wrapper around ape's ace function for ancestral character estimation, which also fits this model:

## here's the wrapper function
fitLambda<-function(tree,x,model="ER"){
    lik<-function(lambda,tree,x,model)
    logLik(ace(x,rescale(tree,model="lambda",lambda),
        type="discrete",model=model))
    obj<-optimize(lik,c(0,1),tree=tree,x=x,model=model,maximum=TRUE)
    fit<-ace(x,rescale(tree,model="lambda",lambda=obj$maximum),
        type="discrete",model=model)
    I<-fit$index.matrix
    fitted.Q=matrix(fit$rates[I],dim(I)[1],dim(I)[2],
        dimnames=list(dimnames(fit$lik.anc)[[2]],
        dimnames(fit$lik.anc)[[2]]))
    diag(fitted.Q)<--rowSums(fitted.Q,na.rm=TRUE)
    list(Q=fitted.Q,lambda=obj$maximum,logLik=logLik(fit))
}
library(geiger)
library(phytools)
## simulate some data to test it
tree<-pbtree(n=200,scale=1)
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
x<-sim.history(rescale(tree,model="lambda",lambda=0.7),Q)$states

## Done simulation(s).

fitLambda(tree,x)

## $Q
##         a       b
## a -0.7995  0.7995
## b  0.7995 -0.7995
## 
## $lambda
## [1] 0.5154
## 
## $logLik
## [1] -128.6

Note that the same model can also be fit using geiger's fitDiscrete function, but the poster reported some issues with convergence which made me want to post the second way. Here is fitDiscrete

fitDiscrete(tree,x,transform="lambda")

## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 0, R2 = 0
##  
## DINTDY-  T (=R1) illegal
## In above message, R1 = 2.8036e-222
##  
##       T not in interval TCUR - HU (= R1) to TCUR (=R2)
## In above message, R1 = 0, R2 = 0
##  
## DINTDY-  T (=R1) illegal
## In above message, R1 = 8.66551e-221
##  
##       T not in interval TCUR - HU (= R1) to TCUR (=R2)
## In above message, R1 = 0, R2 = 0
##  
## DLSODA-  Trouble in DINTDY.  ITASK = I1, TOUT = R1
## In above message, I1 = 1
##  
## In above message, R1 = 8.66551e-221
##

Here I've excluding a whole bunch of error messages....

##  
## In above message, R1 = 8.66551e-221
## 

## GEIGER-fitted comparative model of discrete data
##  fitted Q matrix:
##             a       b
##     a -0.7928  0.7928
##     b  0.7928 -0.7928
## 
##  fitted 'lambda' model parameter:
##  lambda = 0.508644
## 
##  model summary:
##  log-likelihood = -128.569446
##  AIC = 261.138891
##  AICc = 261.199805
##  free parameters = 2
## 
## Convergence diagnostics:
##  optimization iterations = 100
##  failed iterations = 56
##  frequency of best fit = NA
## 
##  object summary:
##  'lik' -- likelihood function
##  'bnd' -- bounds for likelihood search
##  'res' -- optimization iteration summary
##  'opt' -- maximum likelihood parameter estimates

Also, fitting the λ model to discrete trait data has always struck me as a somewhat peculiar enterprise, so this is posted without any prejudice towards whether this is a good idea, a bad idea, or an average idea in the first place!