Analyzing Trait Evolution with Approximate Bayesian Computation in TreEvo

For this vignette, we will use the tree of ammonite genera, and the continuous character data for those ammonites, taken from Raia et al. (2015, American Naturalist). This data is taken from their supplemental appendix, available at the American Naturalist website, and is available as an example dataset in R package paleotree.

library(paleotree)

## Loading required package: ape

data(RaiaCopesRule)

ls()

##  [1] "ammoniteTraitsRaia"   "ammoniteTreeRaia"     "ceratopsianTreeRaia" 
##  [4] "cervidTreeRaia"       "charData"             "coreLimit"           
##  [7] "generation.time"      "multicore"            "nInitialSimsPerParam"
## [10] "nRuns"                "nStepsPRC"            "numParticles"        
## [13] "resultsBM"            "resultsBound"         "shellSize"           
## [16] "sutureComplexity"     "tree"

The trees for the three groups examined in this paper all appear to be trees dated to the last appearance times (as opposed to the first appearance time) and specifically the end-boundary of the interval containing the last appearance.

First, let’s plot the tree:

plot(ladderize(ammoniteTreeRaia), show.tip.label=FALSE)
axisPhylo()

This tree has polytomies - we will have to deal with that. Let’s randomly resolve tree using multi2di from ape, and then use function addTermBranchLength from package paleotree

tree<-multi2di(ammoniteTreeRaia)

# let's apply ATBL
tree<-addTermBranchLength(tree,0.01)

# let's try plotting it again
plot(ladderize(tree), show.tip.label=FALSE)
axisPhylo()

Because we aren’t adjust the internal edges at all, multi2di has just inserted a number of zero-length internal edges between internodes within the tree. This means that alternative resolutions of the tree (as multi2di is stochastic) won’t create meaningfully different variance-covariance matrices, so we likely don’t need to worry about resolving multiple trees in this case (because we’re taking the trees used by Raia et al. at total face-value, which perhaps is the truly unsupported assumption).

Let’s plot a traitgram of the data.

maybe use phytools rather than mine?

## Plotting ancestral reconstructions for nodes.
## Plotting ancestral reconstructions for nodes.

how small is the smallest branch length, even after extending terminal branches?

brlen<-tree$edge.length
min(brlen)

## [1] 0

# zero-length branches!

# what about not zero-length branches?
min(brlen[brlen!=0])

## [1] 0.199982

# in years
min(brlen[brlen!=0])*10^6

## [1] 199982

# let's look at the distribution of these edge lengths, ignoring zero-length
hist(brlen[brlen!=0], main="",xlab="edge length")

# total number of zero-length branches
sum(brlen==0)

## [1] 113

# proportion of edges that are zero-length
sum(brlen==0)/length(brlen)

## [1] 0.2354167

Yeesh - nearly a quarter of branches on this tree are zero-length branches. Not looking so good.

# number of zero-length internal branches
sum(brlen==0 & tree$edge[,2]>Ntip(tree))

## [1] 113

# number of zero-length terminal branches 
  # (shouldn't be any because of our application of ATBL)
sum(brlen==0 & tree$edge[,2]<=Ntip(tree))

## [1] 0

Let’s use TreEvo!

library(TreEvo)

## Warning: replacing previous import 'phytools::scores' by 'pls::scores' when
## loading 'TreEvo'

packageVersion("TreEvo")

## [1] '0.22.1'

# character data for doRun must be in matrix form
  # with rows labeled with taxon names

charData<-matrix(sutureComplexity,ncol=1)
rownames(charData)<-names(sutureComplexity)

Let’s analyze it with ordinary BM

resultsBM<-doRun_prc(
  phy = tree,
  traits = charData,
  intrinsicFn=brownianIntrinsic,
  extrinsicFn=nullExtrinsic,
  startingPriorsFns="normal",
  startingPriorsValues=matrix(c(mean(charData[,1]), sd(charData[,1]))),
  intrinsicPriorsFns=c("exponential"),
  intrinsicPriorsValues=matrix(c(10, 10), nrow=2, byrow=FALSE),
  extrinsicPriorsFns=c("fixed"),
  extrinsicPriorsValues=matrix(c(0, 0), nrow=2, byrow=FALSE),
  generation.time=generation.time,
  standardDevFactor=0.2,
  plot=FALSE,
  StartSims=10,
  epsilonProportion=0.7,
  epsilonMultiplier=0.7,
  nStepsPRC=3,
  numParticles=20,
  jobName="typicalBMrun",
  stopRule=FALSE,
  multicore=multicore,
  coreLimit=coreLimit,
  verboseParticles=TRUE
  )

Let’s analyze it with ordinary BM with a bound

resultsBound<-doRun_prc(
  phy = tree,
  traits = charData,
  intrinsicFn=boundaryMinIntrinsic,
  extrinsicFn=nullExtrinsic,
  startingPriorsFns="normal",
  startingPriorsValues=matrix(c(mean(charData[,1]), sd(charData[,1]))),
  intrinsicPriorsFns=c("exponential","normal"),
  intrinsicPriorsValues=matrix(c(10, 10, -10, 1), nrow=2, byrow=FALSE),
  extrinsicPriorsFns=c("fixed"),
  extrinsicPriorsValues=matrix(c(0, 0), nrow=2, byrow=FALSE),
  generation.time=generation.time,
  standardDevFactor=0.2,
  plot=FALSE,
  StartSims=10,
  epsilonProportion=0.7,
  epsilonMultiplier=0.7,
  nStepsPRC=3,
  numParticles=20,
  jobName="BMwithBoundRun",
  stopRule=FALSE,
  multicore=multicore,
  coreLimit=coreLimit,
  verboseParticles=FALSE
  )

The function summarizePosterior provides us with the mean, standard deviation, and highest posterior density (HPD) for each free parameter in the posterior sample, for each individual run. We can specify the probability density we would like to receive the bounds for using the alpha argument.

For example, for each Brownian Motion run:

summarizePosterior(resultsBM[[1]]$particleDataFrame, alpha = 0.8)
summarizePosterior(resultsBM[[2]]$particleDataFrame, alpha = 0.8)
summarizePosterior(resultsBM[[3]]$particleDataFrame, alpha = 0.8)

Or just to look at the first run of the Bounded BM run:

summarizePosterior(resultsBound[[1]]$particleDataFrame, alpha = 0.8)

# This function calculates Effective Sample Size (ESS) on results.  
# Performs the best when results are from multiple runs.

# ESS for single runs
pairwiseESS(resultsBM)

pairwiseESS(resultsBound)

# ESS for parameters shared between BM and bound
#pairwiseESS(list(resultsBM$particleDataFrame,resultsBound$particleDataFrame))

# plotPosteriors
    # for each free parameter in the posterior, a plot is made of the distribution of values estimate in the last generation

# collect the particleDataFrames into lists

resultsBMpart<-lapply(resultsBM,function(x) x$particleDataFrame)

resultsBoundpart<-lapply(resultsBound,function(x) x$particleDataFrame)

plotPosteriors(particleDataFrame=resultsBMpart,
   priorsMat=resultsBM[[1]]$PriorMatrix)

plotPosteriors(particleDataFrame=resultsBoundpart,
   priorsMat=resultsBound[[1]]$PriorMatrix)