Skip to content

Instantly share code, notes, and snippets.

@dwinter
Created July 24, 2012 23:36
Show Gist options
  • Save dwinter/3173395 to your computer and use it in GitHub Desktop.
Save dwinter/3173395 to your computer and use it in GitHub Desktop.
calculate gsi in 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)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment