[Mattice-commits] r26 - pkg/R

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

Author: andrew_hipp
Date: 2008-11-14 09:00:03 +0100 (Fri, 14 Nov 2008)
New Revision: 26

Modified:
   pkg/R/revisedBatchHansen.R
   pkg/R/revisedTreeTraversal.R
Log:
runBatchHansen now works correctly with ouch v2

Modified: pkg/R/revisedBatchHansen.R
===================================================================
--- pkg/R/revisedBatchHansen.R	2008-11-14 02:19:23 UTC (rev 25)
+++ pkg/R/revisedBatchHansen.R	2008-11-14 08:00:03 UTC (rev 26)
@@ -30,15 +30,16 @@
 # 11 nov 08: renamed to runBatchHansen
 # Runs batchHansenFit and brown.fit over a list of ouchTrees
 # Arguments:
-#  "ouchTrees" = list of OUCH-style trees; if a single tree, send as 'list(TREE)'
+#  "ouchTrees" = list of OUCH-style trees
 #  "characterStates" = vector of character states, either extracted from an ouch-style tree data.frame or a named vector
 #  REMOVED: "SEM"= standard error of the mean, vector extracted from an ouch-style tree data.frame
 #  "rescale" = factor to multiply against (times / max(times)) -- choose based on trial analyses; set at <= 0 if you don't want to rescale trees
 #  "cladeMembersList" = list of vectors containing names of the members of each clade (except for the root of the tree)
 #  "brown" = whether to analyse the data under a Brownian motion model
 #  "..." = additional arguments to pass along to hansen
-function(ouchTrees, characterStates, cladeMembersList, maxNodes = NULL, regimeTitles = NULL, logData = FALSE, brown = F, rescale = 1, ...) {
+function(ouchTrees, characterStates, cladeMembersList, maxNodes = NULL, regimeTitles = NULL, brown = F, rescale = 1, ...) {
   ## do all the objects in ouchTrees inherit ouchtree?
+  if(is(ouchTrees,'ouchtree')) ouchTrees <- list(ouchTrees)
   treeCheck <- unlist(lapply(ouchTrees, function(x) is(x,'ouchtree')))
   if(F %in% treeCheck) 
         stop(paste('This function has been rewritten to use the new S4 ', sQuote('ouchtree'), ' class.',
@@ -46,6 +47,7 @@
   
   ## Check character states to make sure that they are either named and match names in the trees, or are the same length as the tips
   for (i in 1:length(ouchTrees)) {
+    dataFlag <- NULL
     stopFlag <- F
     tree <- ouchTrees[[i]]
     terminals <- tree at nodelabels[(tree at nnodes - tree at nterm + 1):tree at nnodes]
@@ -55,27 +57,26 @@
         message("Data assumed to be in the same order as terminals")
         dataFlag <- 'sameOrderTerminals' 
         }
-      else if (length(characterStates) == tree at nnodes) {
+      if(length(characterStates) == tree at nnodes) {
         message("Data assumed to be in the same order as nodes;\nany data not associated with a terminal branch will be ignored")
         dataFlag <- 'sameOrderNodes'
         }
-      else stopFlag <- T}
+      if(identical(dataFlag, NULL)) stopFlag <- T
       message("-------------------\n")
       }
     else dataFlag <- 'named'
-  if(stopFlag) stop("Correct discrepancies between trees and data and try again!")
-  stopFlag <- F  
+    if(stopFlag) stop("Correct discrepancies between trees and data and try again!")
+    }
 
-  if(logData) characterStates <- log(characterStates)
-  if(identical(regimeTitles, NULL)) {
-    regimeTitles <- as.character(1:length(regimesList))
-    if(brown) regimeTitles <- c(regimeTitles, 'brown')
-  }
   hansenBatch = vector("list", length(ouchTrees))
   treeCounter = 0
   for (i in 1:length(ouchTrees)) {
     tree <- ouchTrees[[i]]
     regimesList = regimeVectors(tree, cladeMembersList, maxNodes)
+    if(identical(regimeTitles, NULL)) {
+      regimeTitles <- as.character(1:length(regimesList))
+      if(brown) regimeTitles <- c(regimeTitles, 'brown')
+      }
     
     ## rescale tree if requested
     if(rescale>0) tree at times <- rescale * tree at times / max(tree at times) 
@@ -85,10 +86,11 @@
 
     ## make sure data fits the tree
     dataIn <- NULL
-    if(dataFlag) == 'sameOrderTerminals' dataIn <- c(rep(NA, tree at nnodes - tree at nterm), characterStates)
-    if(dataFlag) == 'sameOrderNodes' dataIn <- characterStates
-    if(dataFlag) == 'named' dataIn <- characterStates[match(tree at nodelabels), characterStates]
+    if(dataFlag == 'sameOrderTerminals') dataIn <- c(rep(NA, tree at nnodes - tree at nterm), characterStates)
+    if(dataFlag == 'sameOrderNodes') dataIn <- characterStates
+    if(dataFlag == 'named') dataIn <- characterStates[match(tree at nodelabels), characterStates]
     if(identical(dataIn, NULL)) stop(paste("There is a problem with your data that I failed to catch at the outset of", sQuote('runBatchHansen()')))
+    else names(dataIn) <- tree at nodes
     
     ## send it off to batchHansen and just stick the results in hansenBatch... this won't work as the number of regimes gets large, 
     ##   so there should be some option here to just hang onto the coefficients for each run (i.e., hang onto 'coef(hansen(...))' rather than 'hansen(...)')
@@ -115,8 +117,10 @@
   brVars <- c("loglik", "dof", "sigma.squared", "theta")
   haVars <- c("loglik", "dof", "sigma.squared", "alpha")
   treeData <- matrix(data = NA, nrow = length(regimeTitles), ncol = length(variables), dimnames = list(regimeTitles,variables))
-  if(brown) br <- brown(data, tree)
-  treeData["brown", ] <- unlist(summary(br)[brVars])
+  if(brown) {
+    br <- brown(data, tree)
+    treeData["brown", ] <- unlist(summary(br)[brVars])
+    }
   for (i in seq(regimesList)) {
     message(paste("Running regime",i))
     ## at this point, the user has to give an initial alpha and sigma for hansen to search on... this should be relaxed

Modified: pkg/R/revisedTreeTraversal.R
===================================================================
--- pkg/R/revisedTreeTraversal.R	2008-11-14 02:19:23 UTC (rev 25)
+++ pkg/R/revisedTreeTraversal.R	2008-11-14 08:00:03 UTC (rev 26)
@@ -225,6 +225,7 @@
   regimePaintings = vector("list", length(allRegimes))
   for (i in 1:length(allRegimes)) {
     allRegimes[[i]] = c("1", allRegimes[[i]])
-    regimePaintings[[i]] = paintBranches(tree, allRegimes[[i]], as.character(allRegimes[[i]])) 
+    regimePaintings[[i]] = as.factor(paintBranches(tree, allRegimes[[i]], as.character(allRegimes[[i]])))
+    names(regimePaintings[[i]]) <- tree at nodes
     message(paste('Created regime',i))}
   return(regimePaintings) }