diversity dependent diversification tutorial

November 8, 2015

``library("TreeSim")``
``````## Loading required package: ape
``library("TreePar")``
``````## Loading required package: Matrix
##
## 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()[]``````

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)[][]``

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``````

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?