Working similarity matrix code sensu Leinster & Cobbold 2012

WorkingZCode

require(reshape)
require(picante)
require(gtools)

#Function that lists all terminal tip labels descended from ancestral 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)
}

#Simulated ultrametric phylogeny with 3 species
phy <- read.tree(text="(3:3,(1:2,2:2):1);")
phy$tip.label <- as.numeric(phy$tip.label)
plot(phy)
nodelabels()
edge <- phy$edge[,2]

taxa <- length(phy$tip.label)
Bnode = taxa + 1
Lb = t(t(dist.nodes(phy)[1:taxa, Bnode]))
That <- sum(Lb * p)
p = rep((1/taxa), taxa)

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)
  }
}

#Number of historic species
Nhist <- sum(sp[-1,-1])
histsp <- sp[-1,-1]
colnames(histsp) <-sort(edge)
rownames(histsp) <- S
hs <- melt(histsp, varnames = c("species", "branch"))
hs$abund <- rep(p, length(edge))
hs$Lb <- rep(Lb, length(edge))
hs <- (subset(hs, hs$value != 0))
hs <- hs[order(hs$species, -hs$branch), ]

#Initialize Z similarity matrix
Z <- matrix(NA, Nhist, Nhist)

for(i in 1:Nhist){
  
  branch <- internal2tips.self(phy,hs$branch[i])
  
  for(j in 1:Nhist){
    
    Z[i,j] <- ifelse(hs$species[j] %in% branch, That/Lb[hs$species[j]], 0)
  }
}

#Print Similarity matrix
Z

pi = matrix((hs$Lb / That)*hs$abund)

Z %*% pi