[Ecopd-commits] r85 - in branches/single-tree: . R

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

Author: regetz
Date: 2009-11-18 20:14:46 +0100 (Wed, 18 Nov 2009)
New Revision: 85

Modified:
   branches/single-tree/NAMESPACE
   branches/single-tree/R/pd.R
Log:
implemented abundance-weighted pd calculation for phylo4com (but only
using minimum tip lengths automatically inferred from Supertree, hence
only valid for species represented therein)


Modified: branches/single-tree/NAMESPACE
===================================================================
--- branches/single-tree/NAMESPACE	2009-11-18 18:33:57 UTC (rev 84)
+++ branches/single-tree/NAMESPACE	2009-11-18 19:14:46 UTC (rev 85)
@@ -34,8 +34,6 @@
 # pd-related functions
 export(minTL)
 export(`minTL<-`)
-export(getMinTL)
-export(weightByAbund)
 
 # miscellaneous
 export(caterpillar)
@@ -46,3 +44,8 @@
 export(Supertree)
 
 #----------------------------------------------------------------------
+
+# unexported functions
+#export(getSubtrees)
+#export(getMinTL)
+#export(weightByAbund)

Modified: branches/single-tree/R/pd.R
===================================================================
--- branches/single-tree/R/pd.R	2009-11-18 18:33:57 UTC (rev 84)
+++ branches/single-tree/R/pd.R	2009-11-18 19:14:46 UTC (rev 85)
@@ -30,33 +30,22 @@
 
     method <- match.arg(method)
 
-    ##TODO fix abundance-weighted PD!
     if (method %in% c("polytomy", "yule")) {
-        stop("Non-trad PD not yet implemented for phylo4com")
-    } 
-
-    .pd <- function(phy, method) {
-        ## if using non-traditional pd, adjust tip lengths based on abundance
-        if (method %in% c("polytomy", "yule")) {
-            phy <- weightByAbund(phy, method)
-        } 
-        # exclude root edge from calculation (if it exists)
-        if (isRooted(phy)) {
-            nonroot.nodes <- setdiff(nodeId(phy), rootNode(phy))
-            tot.length <- sum(edgeLength(phy, nonroot.nodes))
-        } else {
-            tot.length <- sum(edgeLength(phy))
+        message("estimating minTL values from Supertree")
+        res <- sapply(communities(x), function(community) {
+            phyd <- phylo4d(x, community)
+            n <- tipData(phyd)[[community]]
+            pd(extractTree(weightByAbund(phyd, n, method)))
+        })
+    } else if (method == "traditional") {
+        comms <- x at metadata$comms
+        if (is.null(comms)) {
+            return(.pd(extractTree(x), method))
         }
-        return(tot.length)
+        subtrees <- x at subtrees[as.character(comms)]
+        res <- sapply(subtrees, pd)[as.character(comms)]
+        names(res) <- names(comms)
     }
-
-    comms <- x at metadata$comms
-    if (is.null(comms)) {
-        return(.pd(extractTree(x), method))
-    }
-    subtrees <- x at subtrees[as.character(comms)]
-    res <- sapply(subtrees, .pd, method)[as.character(comms)]
-    names(res) <- names(comms)
     return(res)
 
 })
@@ -120,26 +109,17 @@
 }
 
 # function to weight tip length based on abundance
-weightByAbund <- function(tree, method=c("polytomy", "yule")) {
+weightByAbund <- function(tree, n, method=c("polytomy", "yule")) {
 
   method <- match.arg(method)
 
-  if (is.null(abundance(tree))) {
-    stop("tree contains no abundance information")
-  }
-
   if (is.null(minTL(tree))) {
-    stop("tree contains no minTL information")
+    minLength <- getMinTL(tree, genera(tree))
+  } else {
+    minLength <- minTL(tree)
   }
 
-  n <- abundance(tree)
-  minLength <- minTL(tree)
   tipLen <- tipLength(tree)
- 
-  # Test statement:
-  if (!identical(names(n), names(minLength))) {
-    stop("mismatch between abundance and minTL vectors")
-  }
 
   if (method=="polytomy") {
     newLength <- tipLen + (n-1) * minLength