[Ecopd-commits] r63 - branches/single-tree/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 9 19:15:32 CET 2009
Author: regetz
Date: 2009-11-09 19:15:32 +0100 (Mon, 09 Nov 2009)
New Revision: 63
Added:
branches/single-tree/R/ed.R
Modified:
branches/single-tree/R/indices.R
Log:
rewrote original set of ED (and related) functions as set of 'ed' (and
related) S4 methods
Added: branches/single-tree/R/ed.R
===================================================================
--- branches/single-tree/R/ed.R (rev 0)
+++ branches/single-tree/R/ed.R 2009-11-09 18:15:32 UTC (rev 63)
@@ -0,0 +1,221 @@
+##
+## ED and related methods
+##
+
+setGeneric("ed", function(x, na.rm=TRUE) {
+ standardGeneric("ed")
+})
+
+setMethod("ed", signature(x="phylo4d"), function(x, na.rm=TRUE) {
+ phyc <- phylo4com(x)
+ ed(phyc, na.rm=na.rm)
+})
+
+# TODO: This function includes its own code for not counting root edge
+# length. Maybe this should maybe be done at a higher level?
+setMethod("ed", 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))]
+ .edi <- function(tree) {
+ # set length of root edge to zero
+ edgeLength(tree)[edgeId(tree, "root")] <- 0
+
+ all.nodes <- nodeId(tree, type = "all")
+ des <- descendants(tree, all.nodes, type="tips")
+ nv <- edgeLength(tree, all.nodes) / sapply(des, length)
+ names(nv) <- all.nodes
+
+ tip.nodes <- nodeId(tree, "tip")
+ anc <- ancestors(tree, tip.nodes, "ALL")
+
+ res <- sapply(anc, function(n) sum(nv[as.character(n)], na.rm=TRUE))
+ names(res) <- tipLabels(tree)
+ res
+ }
+ res <- lapply(subtrees, .edi)[as.character(comms)]
+ names(res) <- names(comms)
+ return(res)
+
+})
+
+setGeneric("hed", function(x, na.rm=TRUE) {
+ standardGeneric("hed")
+})
+
+setMethod("hed", signature(x="phylo4d"), function(x, na.rm=TRUE) {
+ phyc <- phylo4com(x)
+ hed(phyc, na.rm=na.rm)
+})
+
+setMethod("hed", signature(x="phylo4com"), function(x, na.rm=TRUE) {
+ ED <- ed(x)
+ PD <- pd(x)
+ res <- sapply(seq_len(length(PD)), function(i) {
+ scaledED <- ED[[i]] / PD[[i]]
+ -sum(scaledED * log(scaledED))
+ })
+ names(res) <- names(PD)
+ return(res)
+})
+
+setGeneric("eed", function(x, na.rm=TRUE) {
+ standardGeneric("eed")
+})
+
+setMethod("eed", signature(x="phylo4d"), function(x, na.rm=TRUE) {
+ phyc <- phylo4com(x)
+ eed(phyc, na.rm=na.rm)
+})
+
+setMethod("eed", 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))]
+ hed(x) / sapply(subtrees, nTips)[as.character(comms)]
+})
+
+setGeneric("aed", function(x, na.rm=TRUE) {
+ standardGeneric("aed")
+})
+
+setMethod("aed", signature(x="phylo4d"), function(x, na.rm=TRUE) {
+ phyc <- phylo4com(x)
+ aed(phyc, na.rm=na.rm)
+})
+
+# TODO: This function includes its own code for not counting root edge
+# length. Maybe this should maybe be done at a higher level?
+setMethod("aed", 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))]
+ .isD <- function(tree) {
+ # get all node IDs, but excluding root node
+ nonroot.nodes <- setdiff(nodeId(tree), rootNode(tree))
+ # Create logical matrix indicating which tips (in columns) are
+ # descendants of each node (in rows), self-inclusive
+ t(sapply(ancestors(tree, tipLabels(tree), "ALL"),
+ function(n) nonroot.nodes %in% n))
+ }
+
+ .elen <- function(tree) {
+ nonroot.nodes <- setdiff(nodeId(tree), rootNode(tree))
+ edgeLength(tree, nonroot.nodes)
+ }
+ isDescendant <- lapply(subtrees, .isD)[as.character(comms)]
+ edge.length <- lapply(subtrees, .elen)[as.character(comms)]
+ N <- abundance(x)
+
+ res <- lapply(seq_along(N), function(i) {
+ spp <- row.names(isDescendant[[i]])
+ dAbund <- N[spp, i] * isDescendant[[i]]
+
+ # Calculate individual-based AED of each species
+ AED <- colSums(edge.length[[i]] * t(prop.table(dAbund, margin=2)))
+ AED/N[spp, i]
+ })
+ names(res) <- names(comms)
+ return(res)
+
+})
+
+# .isD <- function(tree, template) {
+# # get all node IDs, but excluding root node
+# nonroot.nodes <- setdiff(nodeId(tree), rootNode(tree))
+# # Create logical matrix indicating which tips (in columns) are
+# # descendants of each node (in rows), self-inclusive
+# res <- t(sapply(ancestors(tree, tipLabels(tree), "ALL"),
+# function(n) nonroot.nodes %in% n))
+# template[match(rownames(res), rownames(template)), ] <- res
+# template
+# }
+#
+# .el <- function(tree) {
+# nonroot.nodes <- setdiff(nodeId(tree), rootNode(tree))
+# edgeLength(tree, nonroot.nodes)
+# }
+# tmp <- setNames(rep(NA, nTips(x)), tipLabels(x))
+# tmp <- matrix(NA, nrow=nTips(x), ncol=nNodes(x)+nTips(x)-1)
+# rownames(tmp) <- tipLabels(x)
+# isDescendant <- lapply(subtrees, .isD, tmp)
+# isDescendant <- array(unlist(isDescendant), dim=c(nTips(x),
+# nTips(x)+nNodes(x)-1, length(subtrees)))
+#
+# # Create vector of ancestral edge lengths
+# edge.length <- sapply(subtrees, .elen)
+#
+# # Create matrix containing number of individuals of each species
+# # descending from each interior node
+# N <- as.matrix(abundance(x))
+# dAbund <- sweep(isDescendant, c(1,3), N, "*")
+#
+# # Calculate individual-based AED of each species
+# pt <- prop.table(dAbund, margin=c(2,3))
+# AED <- sweep(pt, c(2,3), edge.length, "*")
+# AED <- apply(AED, 3, rowSums) / N
+#
+# return(AED)
+
+setGeneric("haed", function(x, na.rm=TRUE) {
+ standardGeneric("haed")
+})
+
+setMethod("haed", signature(x="phylo4d"), function(x, na.rm=TRUE) {
+ phyc <- phylo4com(x)
+ haed(phyc, na.rm=na.rm)
+})
+
+setMethod("haed", signature(x="phylo4com"), function(x, na.rm=TRUE) {
+ # Recast AED in terms of individuals
+ AED <- aed(x)
+ PD <- pd(x)
+ N <- abundance(x)
+ scaled.AED <- lapply(seq_along(N), function(i) {
+ spp <- names(AED[[i]])
+ rep(unname(AED[[i]]), N[spp, i]) / PD[[i]]
+ })
+ res <- sapply(scaled.AED, function(x) -sum(x * log(x)))
+ names(res) <- names(AED)
+ return(res)
+})
+
+setGeneric("eaed", function(x, na.rm=TRUE) {
+ standardGeneric("eaed")
+})
+
+setMethod("eaed", signature(x="phylo4d"), function(x, na.rm=TRUE) {
+ phyc <- phylo4com(x)
+ eaed(phyc, na.rm=na.rm)
+})
+
+setMethod("eaed", signature(x="phylo4com"), function(x, na.rm=TRUE) {
+ haed(x) / log(colSums(abundance(x)))
+})
+
+setGeneric("value",
+ function(x, na.rm=TRUE) {
+ standardGeneric("value")
+})
+
+setMethod("value", signature(x="phylo4d"),
+ function(x, na.rm=TRUE) {
+ phyc <- phylo4com(x)
+ value(phyc, na.rm=na.rm)
+})
+
+setMethod("value", signature(x="phylo4com"),
+ function(x, na.rm=TRUE) {
+ AED <- aed(x)
+ lapply(AED, function(x) x/sum(x))
+})
+
Property changes on: branches/single-tree/R/ed.R
___________________________________________________________________
Name: svn:eol-style
+ native
Modified: branches/single-tree/R/indices.R
===================================================================
--- branches/single-tree/R/indices.R 2009-11-09 18:11:58 UTC (rev 62)
+++ branches/single-tree/R/indices.R 2009-11-09 18:15:32 UTC (rev 63)
@@ -1,72 +0,0 @@
-# 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) {
-
- # set length of root edge to zero
- edgeLength(tree)[edgeId(tree, "root")] <- 0
-
- all.nodes <- nodeId(tree, type = "all")
- des <- descendants(tree, all.nodes, type="tips")
- nv <- edgeLength(tree, all.nodes) / sapply(des, length)
- names(nv) <- all.nodes
-
- tip.nodes <- nodeId(tree, "tip")
- anc <- ancestors(tree, tip.nodes, "ALL")
- EDI <- sapply(anc, function(n) sum(nv[as.character(n)], na.rm=TRUE))
- names(EDI) <- tipLabels(tree)
-
- return(EDI)
-
-}
-
-HED <- function(tree) {
- scaledED <- ED(tree) / pd(tree)
- return(-sum(scaledED * log(scaledED)))
-}
-
-EED <- function(tree) {
- HED(tree) / log(nTips(tree))
-}
-
-AED <- function(tree) {
-
- # get all node IDs, but excluding root node
- nonroot.nodes <- setdiff(nodeId(tree), rootNode(tree))
- # Create logical matrix indicating which tips (in columns) are
- # descendants of each node (in rows), self-inclusive
- isDescendant <- sapply(ancestors(tree, tipLabels(tree), "ALL"),
- function(n) nonroot.nodes %in% n)
-
- # Create vector of ancestral edge lengths
- edge.length <- edgeLength(tree, nonroot.nodes)
-
- # Create matrix containing number of individuals of each species
- # descending from each interior node
- abundance <- abundance(tree)
- dAbund <- abundance * t(isDescendant)
-
- # Calculate individual-based AED of each species
- AED <- colSums(edge.length * t(prop.table(dAbund, margin=2)))
- AED <- AED/abundance
-
- return(AED)
-
-}
-
-HAED <- function(tree) {
- # Recast AED in terms of individuals
- AED <- AED(tree)
- abundance <- abundance(tree)
- scaledIndivAED <- rep(AED, abundance) / pd(tree)
- return(-sum(scaledIndivAED * log(scaledIndivAED)))
-}
-
-EAED <- function(tree) {
- HAED(tree) / log(sum(abundance(tree)))
-}
-
-value <- function(tree) {
- aed <- AED(tree)
- aed/sum(aed)
-}
-
More information about the Ecopd-commits
mailing list