diversity dependent diversification tutorial
November 8, 2015
## Loading required package: ape ## Loading required package: geiger ## Loading required package: laser
## 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
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
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
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
TreePar there are faster options in
expoTree (not very intuitive).
What are the input choices for
bd.densdep.optim? What would you specify in order to find MLEs?
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
#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
Now, fit the density dependent model and the shift dependent model, which one should fit better?