Challenge problem #4 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)
## Loading required package: ape
library(phytools) ## just in case
## Loading required package: maps
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<-matrix(NA,length(eel.trees),3)
rownames(BD)<-paste("tree",1:length(eel.trees))
colnames(BD)<-c("b","d","logLik")
mlfit<-list()
for(i in 1:length(eel.trees)){
obj<-make.bd(eel.trees[[i]],
sampling.f=0.2)
x.init<-runif(n=2,min=0,max=2)*
rep(log(Ntip(eel.trees[[i]])/0.2)-
log(2),2)/max(nodeHeights(eel.trees[[i]]))
mlfit[[i]]<-find.mle(obj,x.init,method="optim",
lower=0)
BD[i,]<-c(mlfit[[i]]$par,mlfit[[i]]$lnLik)
}
BD
## b d logLik
## tree 1 15.26208 11.846037 468.2897
## tree 2 16.59229 12.928689 477.2350
## tree 3 15.63759 11.302984 483.9629
## tree 4 18.30388 13.361542 500.4655
## tree 5 16.38296 12.705509 476.5338
## tree 6 15.84857 11.782158 480.6809
## tree 7 17.83674 13.863329 486.1560
## tree 8 17.91509 13.914364 486.8039
## tree 9 16.77584 13.648413 469.0235
## tree 10 18.85012 15.303184 483.1629
## tree 11 15.30195 11.365484 476.7475
## tree 12 13.70867 10.110108 465.1442
## tree 13 14.77915 11.019149 472.0973
## tree 14 13.75241 10.180266 464.9340
## tree 15 13.18529 9.241219 468.8839
## tree 16 14.94225 10.397712 484.8120
## tree 17 16.22020 11.925648 485.3762
## tree 18 14.24723 10.153835 475.3377
## tree 19 12.50734 8.417234 468.6523
## tree 20 13.42651 9.363886 471.9442
## tree 21 12.90997 8.966075 467.8987
## tree 22 14.06721 10.156288 471.8695
## tree 23 12.53383 8.363140 470.2703
## tree 24 14.16125 10.220319 473.0166
## tree 25 12.00683 7.828630 468.4835
## tree 26 11.47813 7.307637 465.9759
## tree 27 13.18341 9.124497 471.2601
## tree 28 14.98991 10.769536 480.4454
## tree 29 16.53549 11.921158 491.4277
## tree 30 14.08013 10.134994 472.9519
## tree 31 14.63412 10.691930 474.8802
## tree 32 14.49794 11.299332 461.5248
## tree 33 13.94248 11.016959 453.9769
## tree 34 15.30591 12.324879 460.6445
## tree 35 14.33187 10.297897 475.1240
## tree 36 14.62745 10.618576 475.8084
## tree 37 12.50766 9.010174 458.9456
## tree 38 13.36128 9.405607 470.2291
## tree 39 11.36668 7.445377 461.5717
## tree 40 11.55606 8.114705 453.5788
## tree 41 12.22832 8.896533 454.6970
## tree 42 12.34867 9.469616 446.0983
## tree 43 14.29178 11.375385 455.3056
## tree 44 15.16254 12.200130 459.6889
## tree 45 14.37755 11.049379 463.4753
## tree 46 14.23207 9.997679 478.0791
## tree 47 13.59571 10.471218 456.3234
## tree 48 12.11582 8.923806 451.3408
## tree 49 12.34174 9.325029 448.8539
## tree 50 13.47884 10.486086 453.3860
## tree 51 14.56118 11.129663 465.9158
## tree 52 13.89964 10.586902 461.3056
## tree 53 14.09962 10.678220 464.1186
## tree 54 12.69701 9.246740 458.9998
## tree 55 14.15944 10.749201 464.2162
## tree 56 13.55554 10.723819 450.4554
## tree 57 14.29091 11.345113 455.8633
## tree 58 14.67784 11.413730 463.4733
## tree 59 13.47627 10.062051 461.4113
## tree 60 13.80328 10.788509 455.2011
## tree 61 14.20488 11.760790 445.1899
## tree 62 13.15860 10.636550 442.0407
## tree 63 12.57560 9.831948 444.2493
## tree 64 12.66460 9.851353 446.1486
## tree 65 13.91865 10.659143 460.3444
## tree 66 13.86934 10.591263 460.3816
## tree 67 14.74381 11.698129 459.6296
## tree 68 14.17129 10.966765 460.3342
## tree 69 13.55680 10.562385 453.7866
## tree 70 12.82115 10.109992 444.6842
## tree 71 12.65641 10.112225 440.2013
## tree 72 14.29950 11.533309 452.3067
## tree 73 14.96171 12.231974 454.3350
## tree 74 15.15209 12.283737 457.8658
## tree 75 15.62231 12.572040 463.1775
## tree 76 14.80812 11.727534 460.5572
## tree 77 15.11848 12.315858 456.3709
## tree 78 14.60939 11.598846 458.4287
## tree 79 15.71846 12.446941 467.5975
## tree 80 15.77012 12.996635 458.3247
## tree 81 15.46454 12.934472 452.1122
## tree 82 16.70808 13.963679 461.3605
## tree 83 15.24396 12.663894 452.3818
## tree 84 14.45112 11.937659 447.5313
## tree 85 14.72667 12.057961 452.0391
## tree 86 14.42357 11.901850 447.6274
## tree 87 14.16073 11.419407 451.1845
## tree 88 13.95595 10.912310 456.3341
## tree 89 13.04915 10.108650 450.3984
## tree 90 14.58657 11.823311 453.4224
## tree 91 14.68675 12.038414 451.5182
## tree 92 15.47899 12.481450 461.6354
## tree 93 14.22939 11.282001 455.6389
## tree 94 14.01845 10.530017 464.7263
## tree 95 14.24414 10.991226 461.5134
## tree 96 15.10475 12.198475 458.4120
## tree 97 15.17756 12.405614 456.0521
## tree 98 15.92272 13.572831 450.5137
## tree 99 14.18897 11.309723 454.1140
## tree 100 15.89839 12.284521 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.37069 11.02072 462.92919
## variance
apply(BD,2,var)
## b d logLik
## 2.051057 2.227211 140.929098
## SD
apply(BD,2,sd)
## b d logLik
## 1.432151 1.492384 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 6 December 2016.