Phylogenetic Hill Numbers

Phylo.Z.code

############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"
          )