This function will calculate diversity indices for groups of species through time, given a phylogenetic tree and groupings for the species. The supporting data are: function code, phylogeny and species groups.
# Required packages
library(geiger)
library(phytools)
library(vegan)
source("phylo.diversity.R")
phy <- read.nexus("example.phy.70.spp.nex")
data <- read.csv("example.group.data.csv")
output <- phylo.diversity(phy=phy,groups=data)
# Look at the produced data file
tail(output)## Age shannon simpson Family_A Family_B Family_C Family_D
## 96 -0.22626883 1.625287 0.7880726 0 41 24 4
## 97 -0.18101506 1.623028 0.7876276 0 42 25 4
## 98 -0.13576130 1.618528 0.7867381 0 44 27 4
## 99 -0.09050753 1.621603 0.7870270 0 46 30 5
## 100 -0.04525377 1.620854 0.7869741 0 48 31 5
## 101 0.00000000 1.623278 0.7876538 0 50 31 5
## Family_E Family_F Family_G Family_H
## 96 39 0 15 41
## 97 40 0 15 42
## 98 42 0 15 44
## 99 44 0 14 46
## 100 46 0 15 48
## 101 48 0 17 50To plot an example, we first need to generate overall lineage through time data for the phylogeny
attach(output)
ltt <- as.data.frame(ltt.plot.coords(phy))
root.time <- ltt[1,1]We can now plot this using the results generated above
par(mar=c(5,12,4,4)+0.1)
plot(ltt$time, ltt$N, axes=F, ylim=c(0,max(ltt$N)), xlab="", ylab="",type="l",col="black", main="",xlim=c(root.time,0))
# points(ltt$time, ltt$N,pch=1,col="black")
axis(2, ylim=c(0,max(ltt$N)),col="black",lwd=2)
mtext(2, text="Number of Species", line=2)
par(new=T)
col <- grep("simpson", colnames(output))
r <- range(output[,col])
r[1] <- r[1]-0.1
r[2] <- r[2]+0.1
plot(Age,output[,col], axes=F,ylim=r,
xlab="", ylab="", type="l",lty=2, main="",xlim=c(root.time,0),lwd=2)
axis(2,ylim=c(0,max(output[,col])),lwd=2,line=3.5)
# points(output[,"Age"], output[,col],pch=2)
mtext(2,text="Simpson's Diversity Index",line=5.5)
# add in the x axis
axis(1, pretty(range(ltt$time),10))
mtext("Age (Mya)",side=1,col="black",line=2)
legend("top",legend=c("Number of Species","Simpson Diversity Index"),lty=c(1,2))