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?
First let's load our packages:
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)
title(main="A) speciation rate",adj=0)
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)
title(main="B) extinction rate",adj=0)
lines(rep(mean(BD[,"d"]),2),par()$usr[3:4],lty="dashed",lwd=2,
col="red")
hist(BD[,"b"]-BD[,"d"],xlab="diversification rate",col="grey",main=NULL,
breaks=10)
title(main="C) diversification rate",adj=0)
lines(rep(mean(BD[,"b"]-BD[,"d"]),2),par()$usr[3:4],lty="dashed",lwd=2,
col="red")
Written by Liam J. Revell. Last updated 1 Jul. 2016.