[Mattice-commits] r45 - pkg/R

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

Author: andrew_hipp
Date: 2008-11-20 13:40:08 +0100 (Thu, 20 Nov 2008)
New Revision: 45

Modified:
   pkg/R/batchHansen.R
   pkg/R/informationCriterion.R
   pkg/R/regimes.R
Log:
fixing up multitree analysis summary

Modified: pkg/R/batchHansen.R
===================================================================
--- pkg/R/batchHansen.R	2008-11-20 06:27:51 UTC (rev 44)
+++ pkg/R/batchHansen.R	2008-11-20 12:40:08 UTC (rev 45)
@@ -75,14 +75,13 @@
     message(paste("Tree",i,"of",length(ouchTrees),"complete", "\n-----------------------------"))
   }
     
-    ## right now no summary is returned; one is needed, summarizing over trees what is summarized for each tree in batchHansen
-    ## the below returns presuppose a single tree
+    ## right now no summary is returned; perhaps one is needed, summarizing over trees what is summarized for each tree in batchHansen
   outdata <- list(hansens = hansenBatch, regList = ar$regList, regMatrix = ar$regMatrix, nodeMatrix = ar$nodeMatrix, brown = brown, N = ouchTrees[[i]]@nterm, analysisDate = date(), call = match.call())
   class(outdata) <- 'hansenBatch'
   return(outdata)}
 
 batchHansen <-
-# Runs hansen.fit and brown.fit on a tree over a batch of selective regimes
+# Runs hansen and brown on a tree over a batch of selective regimes
 # Arguments:
 #  "node" "ancestor" "times" "data" = the standard tree specification vectors of the OUCH-style tree
 #  "regimesList" = list of regime-paintings as output from regimeVectors
@@ -101,10 +100,17 @@
     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
-    ha = hansen(data, tree, regimesList[[i]], ...)
-    treeData[i, ] <- unlist(summary(ha)[haVars])
-    hansenOptima[[i]] <- summary(ha)$optima[[1]]
+    if(any(is.na(regimesList[[i]]))) {
+      message(paste("skipping regime", i))
+      treeData[i, ] <- rep(NA, dim(treeData)[[2]])
+      hansenOptima[[i]] <- NA
+      }
+    else {
+      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
+      ha = hansen(data, tree, regimesList[[i]], ...)
+      treeData[i, ] <- unlist(summary(ha)[haVars])
+      hansenOptima[[i]] <- summary(ha)$optima[[1]]
+      }
   }
   return(treeData) }
\ No newline at end of file

Modified: pkg/R/informationCriterion.R
===================================================================
--- pkg/R/informationCriterion.R	2008-11-20 06:27:51 UTC (rev 44)
+++ pkg/R/informationCriterion.R	2008-11-20 12:40:08 UTC (rev 45)
@@ -18,7 +18,7 @@
 
 informationCriterion.hansenBatch <- function(hansenBatch) {
 ## call informationCriterion for a 'hansen.batch' object
-## for right now, just returns AIC, AICc, and BIC weights for the trees analyzed in a hansenBatch object
+## Just returns AIC, AICc, and BIC weights for each of the trees analyzed in a hansenBatch object
   outdata <- list(length(hansenBatch$hansens))
   N = hansenBatch$N
   for(i in 1:length(outdata)) {

Modified: pkg/R/regimes.R
===================================================================
--- pkg/R/regimes.R	2008-11-20 06:27:51 UTC (rev 44)
+++ pkg/R/regimes.R	2008-11-20 12:40:08 UTC (rev 45)
@@ -116,19 +116,18 @@
   for(i in seq(numTrees)) {
     tree <- ouchTrees[[i]]
     regMatrices[[i]] <- regMatrix * as.numeric(matrix(nodeMatrix[i, ], dim(regMatrix)[1], dim(regMatrix)[2], byrow = T)) # multiplies regMatrix by nodes present
-      if(class(regMatrices[[i]]) != "matrix") regMatrices[[i]] <- matrix(regMatrices[[i]], nrow = 1)
-    regMatrices[[i]] <- regMatrices[[i]][which(apply(regMatrices[[i]], 1, sum) > 0), ] # subset for regimes that still have nodes
-      if(class(regMatrices[[i]]) != "matrix") regMatrices[[i]] <- matrix(regMatrices[[i]], nrow = 1)
-    regMatrices[[i]] <- regMatrices[[i]][!duplicated(apply(regMatrices[[i]], 1, as.decimal)), ] ## gets rid of non-unique regimes
-    if(class(regMatrices[[i]]) == "matrix") dimnames(regMatrices[[i]]) <- list(seq(dim(regMatrices[[i]])[1]), dimnames(regMatrices[[i]])[[2]])
-      else regMatrices[[i]] <- matrix(regMatrices[[i]], nrow = 1)
-    regMatrices[[i]] <- rbind(regMatrices[[i]], regMatrices[[i]][1,] * 0) # add one all-zero row for the OU1 model
+    regMatrices[[i]][1:(dim(regMatrices[[i]])[1] - 1), ][which(apply(regMatrices[[i]][1:(dim(regMatrices[[i]])[1] - 1), ], 1, sum) == 0), ] <- rep(NA, numNodes) # set to NA regimes that have no nodes, except for OU1 model
+    regMatrices[[i]][duplicated(apply(regMatrices[[i]], 1, as.decimal)), ] <- rep(NA, numNodes) ## set to NA non-unique regimes
+    dimnames(regMatrices[[i]]) <- list(seq(dim(regMatrices[[i]])[1]), dimnames(regMatrices[[i]])[[2]])
     numTreeRegs <- dim(regMatrices[[i]])[1]
     treeRegs <- list(numTreeRegs) # this will be assigned to regList[[i]]
     nodesVector <- unlist(lapply(nodeMembers, mrcaOUCH, tree = ouchTrees[[i]])) # as written, gets the MRCA for even invalid nodes just so indexing stays right
     for(j in seq(numTreeRegs)) {
-      treeRegs[[j]] <- as.factor(paintBranches(c("1", nodesVector[as.logical(regMatrices[[i]][j, ])]), tree))
-      names(treeRegs[[j]]) <- tree at nodes
+      if(any(is.na(regMatrices[[i]][j, ]))) treeRegs[[j]] <- NA
+      else {
+        treeRegs[[j]] <- as.factor(paintBranches(c("1", nodesVector[as.logical(regMatrices[[i]][j, ])]), tree))
+        names(treeRegs[[j]]) <- tree at nodes
+      }
     }
     regList[[i]] <- treeRegs
   }