Update #1 in phytools 0.7-47: Estimating speciation & extinction from trees

I've been working furiously for the past few days (and less intensely for a few weeks) on prepping a new version of phytools for submission to CRAN. Today I submitted it, and this new version is already availabe on CRAN. Note that the most recent version of phytools can always be installed directly from the phytools GitHub page using devtools. (If you don't know how to do this, just search this blog!)

There are many new features in this phytools version. Unfortunately, I haven't had time of late to keep up with blogging about these changes as they have been added to the package. Over the next few days I will try to make a handful of posts highlighting some of these phytools updates.

In the current phytools version phytools 0.7-47 I have included a number of new datasets not previously packaged with phytools. One of these is a phylogeny of Plethodon salamanders from Highton & Larson (1979). This phylogeny is kind of fun, because it was one of the trees used in the original literature by Sean Nee and colleagues (e.g., Nee et al. 1994) on the estimation of speciation & extinction rates from reconstructed phylogenies.

Let's load this phylogeny & plot it:

packageVersion("phytools")

## [1] '0.7.47'

library(phytools)
data(salamanders)
plotTree(salamanders,ftype="i")

Nee et al. (1994) showed an LTT plot for this tree. Let's duplicate their plot:

plethodon.ltt<-ltt(salamanders,plot=FALSE)
plethodon.ltt

## Object of class "ltt" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal
##     nodes.
## 
## (2) Vectors containing the number of lineages (ltt) and
##     branching times (times) on the tree.
## 
## (3) A value for Pybus & Harvey's "gamma" statistic of
##     gamma = 0.7076, p-value = 0.4792.

plot(plethodon.ltt,log.lineages=FALSE,bty="l",lwd=2,
    main=paste("lineage through time plot for",
    "\nPlethodon salamanders"),font.main=3,
    log="y")

Ok, now let's fit our birth-death model. According to Wikipedia the genus Plethodon has about 55 species. Nee et al. were not able to take incomplete taxon sampling into account; however, subsequent methodological developments now readily permit this. Let's fit a birth-death model assuming that the true species richness of Plethodon is 55.

plethodon.bd<-fit.bd(salamanders,
    rho=Ntip(salamanders)/55)
plethodon.bd

## 
## Fitted birth-death model:
## 
## ML(b/lambda) = 0.01974 
## ML(d/mu) = 0.01492 
## log(L) = -84.419 
## 
## Assumed sampling fraction (rho) = 0.4727 
## 
## R thinks it has converged.

This analysis reveals a rather high rate of extinction relative to speciation.

Finally, fit.bd in phytools now exports the likelihood function which allows is to do things like evaluate alternative values of λ or μ, or compute and plot the likelihood surface.

Let's do the latter:

ngrid<-100
b<-seq(0.01,0.03,length.out=ngrid)
d<-seq(0,0.025,length.out=ngrid)
logL<-sapply(d,function(d,b) sapply(b,function(b,d) 
    plethodon.bd$lik(c(b,d)),d=d),b=b)
contour(x=b,y=d,logL,nlevels=100,
    xlab=expression(lambda),
    ylab=expression(mu),bty="l")
title(main="Likelihood surface for Plethodon diversification",
    font.main=3)
points(plethodon.bd$b,plethodon.bd$d,cex=1.5,pch=4,
    col="blue",lwd=2)
legend("bottomright","ML solution",pch=4,col="blue",
    bg="white",pt.cex=1.5,pt.lwd=2)

Ooh. That's kind of neat.