[Mattice-commits] r165 - pkg/misc

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

Author: andrew_hipp
Date: 2009-02-09 16:46:04 +0100 (Mon, 09 Feb 2009)
New Revision: 165

Modified:
   pkg/misc/multiModels.R
Log:
fixing multiModels to format data correctly

Modified: pkg/misc/multiModels.R
===================================================================
--- pkg/misc/multiModels.R	2009-01-27 18:43:02 UTC (rev 164)
+++ pkg/misc/multiModels.R	2009-02-09 15:46:04 UTC (rev 165)
@@ -1,7 +1,7 @@
 multiModels <- function(phy, dat, node, models = c('whole.brown', 'whole.ou1', 'whole.ou2','part.brown', 'part.ou'), ic = "BIC") {
 # test the support for alternative models on simple and partitioned trees
 # currently only works on one tree; fix so it runs on a set of trees
-  parameters <- c('loglik', 'dof', 'alpha', 'sigma.squared', 'theta', 'optima.dat.1', 'optima.dat.n')
+  parameters <- c('loglik', 'dof', 'alpha', 'sigma.squared', 'theta', 'optima.nodeAncestor', 'optima.nodeDescendent')
   pSum <- c('loglik', 'dof') # parameters to sum for evaluating partitioned trees
   aboveNodeTaxa <- phy$tip.label[-which(phy$tip.label %in% node)]
   wholeTree <- ape2ouch(phy)
@@ -12,7 +12,7 @@
   modelsSplit <- strsplit(models, ".", fixed = T)
   for(i in modelsSplit) {
    model <= paste(i[1], ".", i[2], sep = "")
-   # as written, the following won't correctly chunk parameters into the matrix... mismatch between partitioned- and whole-tree tests
+   # as written, the following won't correctly chunk parameters into the matrix... mismatch between partitioned- and whole-tree tests -- NOW MAY WORK... DOUBLE CHECK RELATIVE TO NAMES SPIT OUT IN WHOLE MODEL, AND CHEKC PARTIALMODEL NAMES
    outMatrix[model,  ] <- ifelse(i[1] == 'whole', 
                                  wholeModel(wholeTree, dat, i[2], node, parameters)$params,
                                  partialModel(partialTree, dat, i[2], c('uptree', 'downtree'), parameters)$params
@@ -22,15 +22,22 @@
   return(outMatrix)
   }
 
-wholeModel <- function(phy, dat, model, nodes, parameterVector) {
+wholeModel <- function(phy, dat, model, node, parameterVector) {
   if(model == "brown") analysis <- brown(dat, phy)
   if(model == "ou1") analysis <- hansen(dat, phy, 
                                         regimes = structure(rep(1, phy at nnodes), 
                                         names = phy at nodes, levels = 1, 
                                         class = 'factor'),
                                         alpha = 1, sigma = 1)
-  if(model = "ou2") analysis <- hansen(dat, phy, regimes = paintBranches(nodes, phy), alpha = 1, sigma = 1)
+  if(model = "ou2") {
+    regime <- paintBranches(node, phy)
+    ancNum <- as.character(tree at root)
+    descNum <- as.character(unique(regime))[unique(regime) != tree at root]
+    analysis <- hansen(dat, phy, regime, alpha = 1, sigma = 1)
+    } # close if
   params <- unlist(summary(analysis)[parameterVector])[parameterVector]
+  params['optima.nodeAncestor'] <- summary(analysis)$optima[[1]][ancNum]
+  params['optima.nodeDescendent'] <- summary(analysis)$optima[[1]][descNum]
   out <- list(analysis = analysis, params = params)
   return(out)
   }