calculate gsi in R

gsi.R

#calculate Cummings et al (2008) _gsi_ - a meaure of the exlusivity of a 
#predefiend group of leaves in a phylogenetic tree


#example
#
# tr <- rtree(10)
# grp <- paste("t", 1:5, sep="")
# gsi(tr, grp)

library(ape)
library(geiger) #for node.leaves()


gsi <- function(tr, grp){
    n <- length(grp) - 1
    #only consider internal nodes (tips get index 1:Ntip(tr))
    internal.nodes <- seq(Ntip(tr)+1, Ntip(tr) + Nnode(tr))
    # For each internal node, what are the tips
    descendants <- lapply(internal.nodes, function(n) node.leaves(tr, n))
    # Which nodes have descendants in the group being considered?
    required <- sapply(descendants, function(x) any(grp %in% x) )
    # How many connections to those nodes have? (tree is not necessarily 
    # dichotomous)
    degree <- table(tr$edge)[ internal.nodes[required] ]
    #Ape takes one connection off the root node, so if d=2 then treat it as d=3
    # (not other nodes can have d=2)
    obs.gs <- n / (sum( degree - 2 ) + sum(degree==2))
    #minGS (basically same procedure but for  whole tree)
    degree.total <- table(tr$edge)[seq(Ntip(tr)+1, Ntip(tr) + Nnode(tr))]
    min.gs <- n / (sum( degree.total - 2 ) + sum(degree.total==2))
    gsi <- (obs.gs - min.gs) / (1 - min.gs)
    return(gsi)
}