-
-
Save pitakakariki/3173954 to your computer and use it in GitHub Desktop.
calculate gsi in R
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #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) | |
| } | |
| gsi2 <- 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)) | |
| # Which nodes have descendants in the group being considered? | |
| grp.index <- match(grp, tr$tip.label) | |
| required <- ancestors(tr, grp.index) | |
| # 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) | |
| } | |
| # Returns a logical vector over the internal nodes. | |
| # | |
| # rval[i] == TRUE iff internal node i has a descendant in grp | |
| ancestors <- function(tr, grp.index) { | |
| rval <- logical(Nnode(tr)) # initially all FALSE | |
| parents <- tr$edge[,1][match(grp.index, tr$edge[,2])] | |
| rval[parents-Ntip(tr)] <- TRUE | |
| if(all(is.na(parents))) rval else rval | ancestors(tr, parents) | |
| } |
Author
pitakakariki
commented
Jul 25, 2012
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment