Challenge problem #7 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<-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.26208 11.846045 468.2897
## STATE_20000 16.59229 12.928688 477.2350
## STATE_30000 15.63759 11.302991 483.9629
## STATE_40000 18.30399 13.361793 500.4655
## STATE_50000 16.38294 12.704985 476.5338
## STATE_60000 15.84859 11.782168 480.6809
## STATE_70000 17.83694 13.863585 486.1560
## STATE_80000 17.91509 13.914363 486.8039
## STATE_90000 16.77560 13.648136 469.0235
## STATE_100000 18.85021 15.303352 483.1629
## STATE_110000 15.30197 11.364878 476.7475
## STATE_120000 13.70867 10.110107 465.1442
## STATE_130000 14.77914 11.019141 472.0973
## STATE_140000 13.75242 10.180274 464.9340
## STATE_150000 13.18526 9.241172 468.8839
## STATE_160000 14.94311 10.398498 484.8120
## STATE_170000 16.21972 11.925166 485.3762
## STATE_180000 14.24722 10.153803 475.3377
## STATE_190000 12.50726 8.417140 468.6523
## STATE_200000 13.42662 9.364049 471.9442
## STATE_210000 12.91010 8.966294 467.8987
## STATE_220000 14.06725 10.156366 471.8695
## STATE_230000 12.53388 8.363200 470.2703
## STATE_240000 14.16122 10.220275 473.0166
## STATE_250000 12.00683 7.828631 468.4835
## STATE_260000 11.47813 7.307636 465.9759
## STATE_270000 13.18344 9.124487 471.2601
## STATE_280000 14.98993 10.769537 480.4454
## STATE_290000 16.53548 11.921134 491.4277
## STATE_300000 14.08069 10.135266 472.9519
## STATE_310000 14.63411 10.691904 474.8802
## STATE_320000 14.49795 11.299339 461.5248
## STATE_330000 13.94249 11.016970 453.9769
## STATE_340000 15.30601 12.325079 460.6445
## STATE_350000 14.33163 10.297619 475.1240
## STATE_360000 14.62654 10.617516 475.8084
## STATE_370000 12.50793 9.010492 458.9456
## STATE_380000 13.36019 9.403825 470.2291
## STATE_390000 11.36661 7.445291 461.5717
## STATE_400000 11.55610 8.114760 453.5788
## STATE_410000 12.22804 8.896138 454.6970
## STATE_420000 12.34878 9.469735 446.0983
## STATE_430000 14.29176 11.375366 455.3056
## STATE_440000 15.16254 12.200130 459.6889
## STATE_450000 14.37755 11.049378 463.4753
## STATE_460000 14.23203 9.997642 478.0791
## STATE_470000 13.59571 10.471221 456.3234
## STATE_480000 12.11582 8.923821 451.3408
## STATE_490000 12.34175 9.325035 448.8539
## STATE_500000 13.47884 10.486068 453.3860
## STATE_510000 14.56118 11.129661 465.9158
## STATE_520000 13.89952 10.586776 461.3056
## STATE_530000 14.09937 10.677933 464.1186
## STATE_540000 12.69719 9.246903 458.9998
## STATE_550000 14.15961 10.749382 464.2162
## STATE_560000 13.55557 10.723818 450.4554
## STATE_570000 14.29079 11.345039 455.8633
## STATE_580000 14.67763 11.413446 463.4733
## STATE_590000 13.47627 10.062053 461.4113
## STATE_600000 13.80321 10.788418 455.2011
## STATE_610000 14.20487 11.760779 445.1899
## STATE_620000 13.15866 10.636588 442.0407
## STATE_630000 12.57560 9.831947 444.2493
## STATE_640000 12.66460 9.851352 446.1486
## STATE_650000 13.91863 10.659102 460.3444
## STATE_660000 13.86928 10.591219 460.3816
## STATE_670000 14.74383 11.698133 459.6296
## STATE_680000 14.17129 10.966764 460.3342
## STATE_690000 13.55962 10.565430 453.7866
## STATE_700000 12.82117 10.110015 444.6842
## STATE_710000 12.65642 10.112233 440.2013
## STATE_720000 14.29985 11.533687 452.3067
## STATE_730000 14.96170 12.231947 454.3350
## STATE_740000 15.15209 12.283737 457.8658
## STATE_750000 15.62232 12.572060 463.1775
## STATE_760000 14.80812 11.727531 460.5572
## STATE_770000 15.11878 12.316355 456.3709
## STATE_780000 14.60939 11.598851 458.4287
## STATE_790000 15.71920 12.447775 467.5975
## STATE_800000 15.77009 12.996607 458.3247
## STATE_810000 15.46451 12.934435 452.1122
## STATE_820000 16.70843 13.964077 461.3605
## STATE_830000 15.24398 12.663908 452.3818
## STATE_840000 14.45117 11.937799 447.5313
## STATE_850000 14.72639 12.057705 452.0391
## STATE_860000 14.42352 11.901805 447.6274
## STATE_870000 14.16073 11.419407 451.1845
## STATE_880000 13.95564 10.911933 456.3341
## STATE_890000 13.04915 10.108642 450.3984
## STATE_900000 14.58656 11.823306 453.4224
## STATE_910000 14.68679 12.038482 451.5182
## STATE_920000 15.47897 12.481414 461.6354
## STATE_930000 14.22939 11.281997 455.6389
## STATE_940000 14.01776 10.529692 464.7263
## STATE_950000 14.24416 10.991244 461.5134
## STATE_960000 15.10472 12.198426 458.4120
## STATE_970000 15.17757 12.405624 456.0521
## STATE_980000 15.92276 13.572907 450.5137
## STATE_990000 14.18896 11.309715 454.1140
## STATE_1000000 15.89839 12.284520 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.37071 11.02073 462.92919
## variance
apply(BD,2,var)
## b d logLik
## 2.051071 2.227292 140.929098
## SD
apply(BD,2,sd)
## b d logLik
## 1.432156 1.492412 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 3 August 2017.