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

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

Author: regetz
Date: 2009-11-09 19:10:33 +0100 (Mon, 09 Nov 2009)
New Revision: 61

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


Modified: branches/single-tree/R/indices.R
===================================================================
--- branches/single-tree/R/indices.R	2009-11-09 18:09:14 UTC (rev 60)
+++ branches/single-tree/R/indices.R	2009-11-09 18:10:33 UTC (rev 61)
@@ -1,19 +1,3 @@
-PAE <- function(tree) {
-
-  # Calculate PD
-  PD <- pd(tree)
-
-  # Extract tip lengths and abundances of those taxa
-  TL <- tipLength(tree)
-  N <- abundance(tree)
-
-  # Calculate PAE
-  numer <- PD + sum(TL * (N-1))
-  denom <- PD + (mean(N)-1) * sum(TL)
-  return(numer/denom)
-
-}
-
 IAC <- function(tree) {
 
   # Count number of lineages originating at each internal node (i.e.

Added: branches/single-tree/R/pae.R
===================================================================
--- branches/single-tree/R/pae.R	                        (rev 0)
+++ branches/single-tree/R/pae.R	2009-11-09 18:10:33 UTC (rev 61)
@@ -0,0 +1,41 @@
+##
+## PAE methods
+##
+
+setGeneric("pae", function(x, na.rm=TRUE) {
+    standardGeneric("pae")
+})
+
+setMethod("pae", signature(x="phylo4d"), function(x, na.rm=TRUE) {
+    phyc <- phylo4com(x)
+    pae(phyc, na.rm=na.rm)
+})
+
+setMethod("pae", signature(x="phylo4com"), function(x, na.rm=TRUE) {
+
+    N <- abundance(x)
+    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))]
+
+    # now for each subtree...
+    #PD <- sapply(subtrees, pd)[as.character(comms)]
+    PD <- pd(x)
+    tmp <- setNames(rep(0, nTips(x)), tipLabels(x))
+    TL <- lapply(subtrees, function(tree) {
+        res <- tipLength(tree)
+        tmp[match(names(res), names(tmp))] <- res
+        tmp
+    })
+    TL <- do.call("cbind", TL[as.character(comms)])
+
+    numer <- PD + colSums(TL * (N - 1))
+    denom <- PD + (colSums(N, na.rm = na.rm) / colSums(presence(x),
+        na.rm=na.rm) - 1) * colSums(TL)
+    res <- numer/denom
+    names(res) <- names(comms)
+    return(res)
+
+})


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