## Solution to Challenge Problem 8: Estimating diversification rates

Challenge problem #8 involves using the file anguila.nex which contains a sample of 100 trees from the Bayesian posterior probability distribution of the phylogeny of elopomorph eels. Only about 20% of exant eels are represented in this tree!

The challenge problem was as follows:

Use diversitree to estimate the speciation and extinction rate for each tree. What is the average speciation and extinction rate across all trees? In what terms are these rates measured? What is the variance and speciation & extinction rates among trees?

``````library(diversitree)
library(phytools) ## just in case
``````

Now, let's try & read our trees from file using the ape function `read.nexus`:

``````eel.trees<-read.nexus("anguila.nex")
eel.trees
``````
``````## 100 phylogenetic trees
``````

Next, let's estimate speciation and extinction rates across all trees. Here I will use an `sapply` call, but I could equally well use a `for` loop and can demonstrate that in class:

``````bd<-function(tree,sampling.f){
obj<-make.bd(tree,sampling.f=sampling.f)
x.init<-runif(n=2,min=0,max=2)*rep(log(Ntip(tree)/sampling.f)-
log(2),2)/max(nodeHeights(tree))
fit<-find.mle(obj,x.init,method="optim",lower=0)
setNames(c(fit\$par,fit\$lnLik),c("b","d","logLik"))
}
BD<-t(sapply(eel.trees,bd,sampling.f=0.2))
BD
``````
``````##                      b         d   logLik
## STATE_10000   15.26269 11.846735 468.2897
## STATE_20000   16.59230 12.928702 477.2350
## STATE_30000   15.63759 11.302990 483.9629
## STATE_40000   18.30401 13.361823 500.4655
## STATE_50000   16.38296 12.705502 476.5338
## STATE_60000   15.84848 11.781954 480.6809
## STATE_70000   17.83664 13.863213 486.1560
## STATE_80000   17.91509 13.914375 486.8039
## STATE_90000   16.77584 13.648413 469.0235
## STATE_100000  18.84874 15.301599 483.1629
## STATE_110000  15.30194 11.365480 476.7475
## STATE_120000  13.70869 10.110131 465.1442
## STATE_130000  14.77919 11.018581 472.0973
## STATE_140000  13.75249 10.180369 464.9340
## STATE_150000  13.18522  9.241129 468.8839
## STATE_160000  14.94317 10.398624 484.8120
## STATE_170000  16.21990 11.925326 485.3762
## STATE_180000  14.24723 10.153815 475.3377
## STATE_190000  12.50754  8.417461 468.6523
## STATE_200000  13.42662  9.364053 471.9442
## STATE_210000  12.91002  8.966185 467.8987
## STATE_220000  14.06749 10.156582 471.8695
## STATE_230000  12.53384  8.363129 470.2703
## STATE_240000  14.16118 10.220225 473.0166
## STATE_250000  12.00673  7.828437 468.4835
## STATE_260000  11.47813  7.307635 465.9759
## STATE_270000  13.18337  9.124375 471.2601
## STATE_280000  14.98991 10.769544 480.4454
## STATE_290000  16.53561 11.921299 491.4277
## STATE_300000  14.07994 10.134729 472.9519
## STATE_310000  14.63413 10.691931 474.8802
## STATE_320000  14.49795 11.299337 461.5248
## STATE_330000  13.94248 11.016959 453.9769
## STATE_340000  15.30612 12.325263 460.6445
## STATE_350000  14.33185 10.297860 475.1240
## STATE_360000  14.62654 10.617513 475.8084
## STATE_370000  12.50773  9.010250 458.9456
## STATE_380000  13.36017  9.403798 470.2291
## STATE_390000  11.36658  7.445226 461.5717
## STATE_400000  11.55751  8.116437 453.5788
## STATE_410000  12.22822  8.896327 454.6970
## STATE_420000  12.34867  9.469615 446.0983
## STATE_430000  14.29181 11.375461 455.3056
## STATE_440000  15.16599 12.204050 459.6889
## STATE_450000  14.37759 11.049424 463.4753
## STATE_460000  14.23208  9.997690 478.0791
## STATE_470000  13.59571 10.471220 456.3234
## STATE_480000  12.11584  8.923852 451.3408
## STATE_490000  12.34174  9.325029 448.8539
## STATE_500000  13.47861 10.485827 453.3860
## STATE_510000  14.56117 11.129660 465.9158
## STATE_520000  13.89972 10.586977 461.3056
## STATE_530000  14.09962 10.678221 464.1186
## STATE_540000  12.69716  9.246909 458.9998
## STATE_550000  14.15961 10.749382 464.2162
## STATE_560000  13.55540 10.723697 450.4554
## STATE_570000  14.29099 11.345199 455.8633
## STATE_580000  14.67784 11.413737 463.4733
## STATE_590000  13.47627 10.062054 461.4113
## STATE_600000  13.80331 10.788566 455.2011
## STATE_610000  14.20488 11.760788 445.1899
## STATE_620000  13.15866 10.636592 442.0407
## STATE_630000  12.57567  9.832053 444.2493
## STATE_640000  12.66460  9.851358 446.1486
## STATE_650000  13.91861 10.659093 460.3444
## STATE_660000  13.86928 10.591220 460.3816
## STATE_670000  14.74382 11.698139 459.6296
## STATE_680000  14.17131 10.966787 460.3342
## STATE_690000  13.55960 10.565392 453.7866
## STATE_700000  12.82117 10.110016 444.6842
## STATE_710000  12.65640 10.112216 440.2013
## STATE_720000  14.29985 11.533689 452.3067
## STATE_730000  14.96167 12.231898 454.3350
## STATE_740000  15.15209 12.283737 457.8658
## STATE_750000  15.62229 12.572005 463.1775
## STATE_760000  14.80796 11.727377 460.5572
## STATE_770000  15.11849 12.315864 456.3709
## STATE_780000  14.60941 11.598906 458.4287
## STATE_790000  15.71849 12.446987 467.5975
## STATE_800000  15.77012 12.996633 458.3247
## STATE_810000  15.46454 12.934472 452.1122
## STATE_820000  16.70841 13.964050 461.3605
## STATE_830000  15.24398 12.663908 452.3818
## STATE_840000  14.45108 11.937610 447.5313
## STATE_850000  14.72667 12.057961 452.0391
## STATE_860000  14.42356 11.901848 447.6274
## STATE_870000  14.16073 11.419408 451.1845
## STATE_880000  13.95597 10.912319 456.3341
## STATE_890000  13.04915 10.108642 450.3984
## STATE_900000  14.58657 11.823312 453.4224
## STATE_910000  14.68675 12.038409 451.5182
## STATE_920000  15.47906 12.481548 461.6354
## STATE_930000  14.22939 11.282001 455.6389
## STATE_940000  14.01706 10.528837 464.7263
## STATE_950000  14.24414 10.991224 461.5134
## STATE_960000  15.10475 12.198475 458.4120
## STATE_970000  15.17756 12.405610 456.0521
## STATE_980000  15.92276 13.572912 450.5137
## STATE_990000  14.18897 11.309724 454.1140
## STATE_1000000 15.89838 12.284519 474.1133
``````

How are the rates measured? Rates are in terms of the units of branch length in the tree.

What is the mean & variance of estimated speciation & extinction rates?

``````## mean
apply(BD,2,mean)
``````
``````##         b         d    logLik
##  14.37075  11.02077 462.92919
``````
``````## variance
apply(BD,2,var)
``````
``````##          b          d     logLik
##   2.050921   2.227166 140.929098
``````
``````## SD
apply(BD,2,sd)
``````
``````##         b         d    logLik
##  1.432104  1.492369 11.871356
``````
``````## visualize
par(mfrow=c(3,1))
hist(BD[,"b"],xlab="speciation rate",col="grey",main=NULL,breaks=10)
lines(rep(mean(BD[,"b"]),2),par()\$usr[3:4],lty="dashed",lwd=2,
col="red")
hist(BD[,"d"],xlab="extinction rate",col="grey",main=NULL,breaks=10)