[Ecopd-commits] r62 - branches/single-tree/R

[noreply at r-forge.r-project.org]

Author: regetz
Date: 2009-11-09 19:11:58 +0100 (Mon, 09 Nov 2009)
New Revision: 62

Added:
   branches/single-tree/R/iac.R
Modified:
   branches/single-tree/R/indices.R
Log:
rewrote original IAC function as 'iac' S4 methods


Added: branches/single-tree/R/iac.R
===================================================================
--- branches/single-tree/R/iac.R	                        (rev 0)
+++ branches/single-tree/R/iac.R	2009-11-09 18:11:58 UTC (rev 62)
@@ -0,0 +1,52 @@
+##
+## IAC methods
+##
+
+setGeneric("iac", function(x, na.rm=TRUE) {
+    standardGeneric("iac")
+})
+
+setMethod("iac", signature(x="phylo4d"), function(x, na.rm=TRUE) {
+    phyc <- phylo4com(x)
+    iac(phyc, na.rm=na.rm)
+})
+
+setMethod("iac", signature(x="phylo4com"), function(x, na.rm=TRUE) {
+
+    comms <- x at metadata$comms
+    if (is.null(comms)) {
+        stop("no community data specified in phylo4com object")
+    }
+    subtrees <- x at subtrees[unique(as.character(comms))]
+    .denom <-  function(tree, template) {
+        # Count number of lineages originating at each internal node
+        # (i.e. number of splits)
+        int.nodes <- nodeId(tree, "internal")
+        nSplits <- sapply(int.nodes, function(x) length(children(tree,
+            x)))
+        names(nSplits) <- int.nodes
+        # For each tip, take the product of the number of splits across
+        # all of its ancestral nodes
+        res <- sapply(ancestors(tree, tipLabels(tree)), function(x)
+            prod(nSplits[as.character(x)]))
+        template[match(names(res), names(template))] <- res
+        template
+    }
+
+    # now for each subtree...
+    tmp <- setNames(rep(NA, nTips(x)), tipLabels(x))
+    denom <- lapply(subtrees, .denom, tmp)
+    denom <- do.call("cbind", denom[as.character(comms)])
+    nnodes <- sapply(x at subtrees, nNodes)[as.character(comms)]
+
+    # Calculate expected number of individuals under null hypothesis
+    # of equal allocation to each lineage at each (node) split 
+    N <- abundance(x)
+    expected <- t(colSums(N, na.rm=na.rm) / t(denom))
+
+    # IAC: summed absolute difference between expected and observed
+    # abundances, divided by number of nodes
+    colSums(abs(expected - N), na.rm=na.rm) / nnodes
+
+})
+


Property changes on: branches/single-tree/R/iac.R
___________________________________________________________________
Name: svn:eol-style
   + native

Modified: branches/single-tree/R/indices.R
===================================================================
--- branches/single-tree/R/indices.R	2009-11-09 18:10:33 UTC (rev 61)
+++ branches/single-tree/R/indices.R	2009-11-09 18:11:58 UTC (rev 62)
@@ -1,28 +1,3 @@
-IAC <- function(tree) {
-
-  # Count number of lineages originating at each internal node (i.e.
-  # number of splits)
-  int.nodes <- nodeId(tree, "internal")
-  nSplits <- sapply(int.nodes, function(x) length(children(tree, x)))
-  names(nSplits) <- int.nodes
-
-  # For each tip, take the product of the number of splits across all of
-  # its ancestral nodes
-  denom <- sapply(ancestors(tree, nodeId(tree, "tip")), function(x)
-    prod(nSplits[as.character(x)]))
-
-  abundance <- abundance(tree)
-
-  # Calculate expected number of individuals under null hypothesis of
-  # equal allocation to each lineage at each (node) split 
-  expected <- sum(abundance) / denom
-
-  # IAC: summed absolute difference between expected and observed
-  # abundances, divide by number of nodes
-  sum(abs(expected-abundance)) / nNodes(tree)
-
-}
-
 # TODO: This function includes its own code for not counting root edge
 # length. Maybe this should maybe be done at a higher level?
 ED <- function(tree) {