Phyloseminar 2015: Tutorial for density dependent diversification

diversitydependent.pdf diversitydependent.Rmd

diversity dependent diversification tutorial

Rosana Zenil-Ferguson

November 8, 2015

library("TreeSim")
## Loading required package: ape
## Loading required package: geiger
## Loading required package: laser
library("TreePar")
## Loading required package: Matrix
## Loading required package: subplex
## Loading required package: deSolve
##
## Attaching package: 'deSolve'
##
## The following object is masked from 'package:graphics':
##
##     matplot

If you have not installed them you should use the function

#install.packages("nameofthepackage")

Tutorial modified from Tanja Stadler’s PhyloDynamics seminar.

Part 1: Density Dependent Diversification Estimates

Simulate a tree using TreeSim package. We are assuming that the change in shape of the tree is given by a shift in the diversification rates (like a mass extinction) and not because of diversity dependence (we will test this in the next section). Input for the simulation will be:

  • 10 extant species
  • lambdavect=c(2,0.1), this represents a shift in birthrate
  • muvect=c(0,0.05), this represents a shift in extinction rate
  • All species are sampled before and after the change in diversification rates
  • time=c(0,1) the shift in rate occure 1 unit of time ago (i.e. mya)
  • We do not see lineages extinct
set.seed(1)
#tree1<-sim.rateshift.taxa()[[1]]

Plot the tree and the lineage through time. Do you notice the “pull of the present”?
Let’s calculate speciation times using the following instruction

#x<-sort(getx(tree1),decreasing=TRUE)

So now we estimate our density dependent diversification rates. Three parameters are of importance (K, lambda_n, mu) mu is fixed but not zero, lambda_n depends on the number of lineages and K is the carrying capacity. For this part we will use function bd.densdep.optim in TreePar there are faster options in DDD and expoTree (not very intuitive).
What are the input choices for bd.densdep.optim? What would you specify in order to find MLEs?

#res.DDD<-bd.densdep.optim()

Part 2: Comparing against a model with a shift in rates

If we want to compare to the model with shifts we can optimize using bd.shifts.optim.

#res.shifts <- bd.shifts.optim(x,c(rho,1),0.1,0.1,2.1)[[2]][[2]]

How are your estimates for the parameters? What are the potential problems?
To make model comparisons we can calculate AIC values for both models

#aic.ddd <- 2*3+2*res.DDD$value
#aic.shifts <- 2*5+2*res.shifts[1]

What should be happening given what you simulated?

Part 3: Challenge

Can you simulate a tree with density dependent diversification? Check the notes of function sim.rateshift.taxa()
Now, fit the density dependent model and the shift dependent model, which one should fit better?