Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Phylogenetic Hill Numbers
############Here's a function that calculates Hill number diversity using a phylogeny from the Leinster and Cobbold (2012) paper
#Requires a phylogeny with tip labels and a community matrix (rows = samples, columns = species)
Phylo.Z <- function(phy,samp){
#packages required for this to work:
require(reshape)
require(picante)
require(ape)
require(gtools)
phy <- prune.sample(samp,phy)
samp <- t(samp)
samp <- samp[mixedsort(rownames(samp)),]
sums <- subset(rowSums(samp), rowSums(samp) > 0)
samp <- samp[names(sums),]
samples = length(samp[1,])
taxa = length(samp[,1])
phy$tip.label <- rep(1:length(phy$tip.label))
edge <- phy$edge[,2]
Bnode = taxa + 1
Li = t(t(dist.nodes(phy)[1:taxa, Bnode]))
Lb <- cbind(phy$edge[,2], phy$edge.length)
Lb <- Lb[order(Lb[,1]),]
p=matrix(0,taxa,samples)
for (k in 1:samples){
p[,k]<-samp[,k]/sum(samp[,k])
}
That=matrix(0,1,samples)
for (k in 1:samples){
That[,k] <-sum(Li * p[,k])
}
########## This function returns the descendents of any internal node #################################
internal2tips.self = function (phy, int.node){
#require(picante); require(ape)
Ntaxa = length(phy$tip.label)
Nnode = phy$Nnode
if ((Ntaxa + Nnode - 1) != nrow(phy$edge)) {
print("tree structure error")
break
}
if (mode(int.node) == "character"){
nodes = which(phy$node.label == int.node) + Ntaxa
}else nodes = int.node
tips = c()
repeat {
nodes = phy$edge[which(phy$edge[, 1] %in% nodes), 2]
if (length(nodes) == 0)
break
tips = c(tips, nodes)
}
tips = tips[tips <= Ntaxa]
if( int.node <= Ntaxa & length(tips) == 0 ){
tips = int.node
}
tips = phy$tip.label[tips]
return(tips)
}
############################################################################################################
#####This builds the Z matrix from your phylogeny ######################################################
########################################################
S <- rep(1:length(phy$tip.label))
sp <- matrix(NA, nrow = length(S) + 1, ncol = length(edge) + 1)
sp[1,-1] <- sort(edge)
sp[-1,1] <- S
edges <- sp[1,]
for (i in (2:ncol(sp))) {
branch <- internal2tips.self(phy,edges[i])
for (j in (2:nrow(sp))){
sp[j,i] <- ifelse(sp[j,1] %in% branch, 1, 0)
}
}
hs <- sp[-1,-1]
colnames(hs) <-sort(edge)
rownames(hs) <- S
hs <- melt(hs, varnames = c("species", "branch"))
for (k in 1:samples){
hs[,k + 3] <- rep(p[,k], length(edge))
}
hs$Li <- rep(Li, length(edge))
hs <- hs[order(hs$species), ]
hs$Lb <- rep(Lb[,2], length(phy$tip.label))
hs <- (subset(hs, hs$value != 0))
hs <- hs[order(hs$species, hs$branch), ]
Z <- matrix(NA, nrow(hs), nrow(hs))
colnames(Z) <- hs$species
rownames(Z) <- hs$branch
for(i in 1:length(hs$branch)){
for(j in 1:length(hs$species)){
branch <- internal2tips.self(phy,hs$branch[i])
Z[i,j] <- ifelse(hs$species[j] %in% branch, That/Li[hs$species[j]], 0)
}
}
pi <- matrix(0, ncol(Z), samples)
for (k in 1:samples){
pi[,k] <- (hs$Lb/That[,k])*hs[,k+3]
}
lenq <- 50
qq <- seq(length = lenq, from = 0, by = .11)
# Initialise the Zp matrix to zero
Zp=matrix(0,ncol(Z), samples)
# Compute Zp
for (k in 1:samples) {
for (i in 1:ncol(Z)){
for (j in 1:ncol(Z)){
Zp[i,k]<-Zp[i,k]+Z[i,j]*pi[j,k]
}
}}
# Initialise the Diversity matrix to zero
Dqz = matrix(0, lenq ,samples)
# Loop to calculate the diversity Dqz for each value of q (iq) and each sample (k)
for (k in 1:samples) {
for (iq in 1:lenq) {
q<-qq[iq];
for (zpi in 1:length(Zp[,k])){
if (Zp[zpi,k]>0)(
Dqz[iq,k]<-Dqz[iq,k]+ pi[zpi,k]*(Zp[zpi,k])^(q-1))
}
Dqz[iq,k] <- Dqz[iq,k]^(1/(1-q));
}
}
return(Dqz)
}
####################################################### DONE!!!! ##############################################
#Here's an example:
data(phylocom)
samp = phylocom$sample
phy = prune.sample(samp, phylocom$phylo)
the.answer <- Phylo.Z(phy = phy, samp = samp)
# Plot the diversity profiles for all the samples (e.g. understorey and canopy) on the same graph.
matplot(qq,the.answer, ylim = c(min(the.answer),max(the.answer)), xlim = c(0,5),
xlab="Sensitivity parameter q",
ylab="Diversity DqZ(p)",
main="Diversity profile using a similarity matrix"
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.