[Mattice-commits] r160 - in pkg: inst/doc misc

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

Author: andrew_hipp
Date: 2009-01-15 19:42:09 +0100 (Thu, 15 Jan 2009)
New Revision: 160

Added:
   pkg/misc/multiModels.R
Modified:
   pkg/inst/doc/maticce.Rnw
Log:
updated vignette, adding multiModels.R

Modified: pkg/inst/doc/maticce.Rnw
===================================================================
--- pkg/inst/doc/maticce.Rnw	2009-01-14 23:26:14 UTC (rev 159)
+++ pkg/inst/doc/maticce.Rnw	2009-01-15 18:42:09 UTC (rev 160)
@@ -115,13 +115,17 @@
 
 <<runBatch, fig=FALSE, echo=TRUE>>=
 # First, analyze with maxNodes set to 2
-ha.4.2 <- runBatchHansen(ovales.tree, ovales.data, ovales.nodes[1:4], maxNodes = 2, brown = T)
+ha.4.2 <- runBatchHansen(ovales.tree, ovales.data, 
+          ovales.nodes[1:4], maxNodes = 2, brown = T)
 print(summary(ha.4.2))
 # Then, analyze with maxNodes set to 4
-ha.4.4 <- runBatchHansen(ovales.tree, ovales.data, ovales.nodes[1:4], maxNodes = 4, brown = T)
+ha.4.4 <- runBatchHansen(ovales.tree, ovales.data, 
+          ovales.nodes[1:4], maxNodes = 4, brown = T)
 print(summary(ha.4.4))
 @
 
+In the examples above, support for node two is relatively little affected by the value of \code{maxNodes}.
+
 \begin{figure}[h]
 \centering
 <<ouSim, fig=TRUE, width=30, height =15, echo=TRUE>>=

Added: pkg/misc/multiModels.R
===================================================================
--- pkg/misc/multiModels.R	                        (rev 0)
+++ pkg/misc/multiModels.R	2009-01-15 18:42:09 UTC (rev 160)
@@ -0,0 +1,59 @@
+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')
+  pSum <- c('loglik', 'dof') # parameters to sum for evaluating partitioned trees
+  aboveNodeTaxa <- phy$tip.label[-which(phy$tip.label %in% node)]
+  wholeTree <- ape2ouch(phy)
+  upTree <- ape2ouch(drop.tip(phy, node))
+  downTree <- ape2ouch(drop.tip(phy, aboveNodeTaxa))
+  partialTree <- list(upTree = upTree, downTree = downTree)
+  outMatrix <- matrix(NA, nrow = length(models), ncol = length(parameters), dimnames = list(models, parameters))
+  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
+   outMatrix[model,  ] <- ifelse(i[1] == 'whole', 
+                                 wholeModel(wholeTree, dat, i[2], node, parameters)$params,
+                                 partialModel(partialTree, dat, i[2], c('uptree', 'downtree'), parameters)$params
+                                 )
+   }
+  MAKE A WEIGHTS COLUMN using informationCriterion
+  return(outMatrix)
+  }
+
+wholeModel <- function(phy, dat, model, nodes, 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)
+  params <- unlist(summary(analysis)[parameterVector])[parameterVector]
+  out <- list(analysis = analysis, params = params)
+  return(out)
+  }
+
+partialModel <- function(phyList, dat, model, treeNames, parameterVector = NULL, pSum = NULL) {
+# sums a subset of parameters (indicated by pSum) and leaves the others separate
+  allParams <- pSum
+  for(i in parameterVector) {
+   if(i %in% pSum) next
+   else allParams <- c(allParams, paste(treeNames, i, sep = ""))
+   }
+  analysis <- ifelse(model == "brown",
+                     lapply(phyList, brown, data = dat),
+                     lapply(phyList, hansen, data = dat, 
+                            regimes = structure(rep(1, phyList[[1]]@nnodes), 
+                                                names = phyList[[1]]@nodes, 
+                                                levels = 1, class = 'factor'),       ## won't work... multiple subtrees, different regimes!
+                            sigma = 1, alpha = 1)
+                     )
+  params <- lapply(analysis, function(x) {unlist(summary(x)[pV])[pV]}, pV = allParams)
+  rawMat <- matrix(unlist(params), nrow = length(params), ncol = length(allParams), dimnames = list(treeNames, allParams))
+  params <- colSums(rawMat, na.rm = T)
+  out <- list(analysis = analysis, rawMat = rawMat, params = params)
+  return(out)
+  }
+  
\ No newline at end of file