[Mattice-commits] r76 - pkg/R

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

Author: andrew_hipp
Date: 2008-12-04 23:56:32 +0100 (Thu, 04 Dec 2008)
New Revision: 76

Added:
   pkg/R/ouSim.hansenBatch.R
Modified:
   pkg/R/ouSim.ouchtree.R
   pkg/R/ouSimHead.R
Log:
ouSim.hansenBatch now takes a hansenBatch object or hansenSummary object + tree and simulates on the model-averaged parameters

Added: pkg/R/ouSim.hansenBatch.R
===================================================================
--- pkg/R/ouSim.hansenBatch.R	                        (rev 0)
+++ pkg/R/ouSim.hansenBatch.R	2008-12-04 22:56:32 UTC (rev 76)
@@ -0,0 +1,8 @@
+ouSim.hansenBatch <- function(analysis, tree, treeNum = 1, rootState = NULL, ...) {
+## runs ouSim.ouchtree for a hansenBatch or hansenSummary object, using the model-averaged alpha, sigma.squared, and theta vector from one tree
+## tree, rootState = 0, alpha = 0, variance = 1, theta = rootState, steps = 1000
+  if(class(analysis) == "hansenBatch") analysis <- summary(analysis)
+  if(identical(rootState, NULL)) rootState <- mean(analysis$thetaMatrix[treeNum, ])
+  outdata <- ouSim(tree, rootState, alpha = analysis$modelAvgAlpha, theta = analysis$thetaMatrix[treeNum, ], ...)
+  return(outdata)
+}
\ No newline at end of file

Modified: pkg/R/ouSim.ouchtree.R
===================================================================
--- pkg/R/ouSim.ouchtree.R	2008-12-04 22:24:53 UTC (rev 75)
+++ pkg/R/ouSim.ouchtree.R	2008-12-04 22:56:32 UTC (rev 76)
@@ -1,10 +1,8 @@
-ouSim.ouchtree <- function(tree, rootState = 0, shiftBranches = NULL, shiftStates = NULL, alpha = 0, variance = 1, theta = rootState, steps = 1000) {
+ouSim.ouchtree <- function(tree, rootState = 0, alpha = 0, variance = 1, theta = rootState, steps = 1000) {
 ## function to plot a simulated dataset under brownian motion or Ornstein-Uhlenbeck (OU) model
 ## Arguments:
 ##   tree is an ouch-style (S4) tree
 ##   alpha and theta are either single values or vectors of length (length(branchList))
-##   shiftBranches is a vector indicating any branches at which an OU or brownian motion model has a determined shift in ancestral state
-##   shiftStates is a vector of length = length(shiftBranches) indicaing the ancestral states for the determined break points
 
 ##embedded function---------------------
 ##can be released to the wild, but more arguments will need to be passed around
@@ -16,15 +14,12 @@
 ## a branch length is the time of a node - the time of its ancestor
 ## not fixed yet for ouchtree
   startBranch <- startNode ## startNode really means start branch... it's the end node of hte branch starting this process
-  if(!identical(shiftStates, NULL)) {
-    if(startBranch %in% shiftBranches) startState <- shiftStates[match(startBranch, shiftBranches)] 
-    }
   message(paste('Working on branch',startBranch,'with starting state',startState))
   workingBranch <- branchList[[startBranch]]
   workingBranch[1] <- startState
   for (brStep in 2:length(workingBranch)) {
     workingBranch[brStep] <- 
-      workingBranch[brStep - 1] + workingBranch[brStep] + alpha[startBranch] / steps * (theta[startBranch] - workingBranch[brStep - 1]) 
+      workingBranch[brStep - 1] + workingBranch[brStep] + alpha[startBranch] / steps * (theta[startBranch] - workingBranch[brStep - 1]) # denom was mult'd by steps... should be? 
     }
   branchList[[startBranch]] <- workingBranch
   endState <- branchList[[startBranch]][length(branchList[[startBranch]])]

Modified: pkg/R/ouSimHead.R
===================================================================
--- pkg/R/ouSimHead.R	2008-12-04 22:24:53 UTC (rev 75)
+++ pkg/R/ouSimHead.R	2008-12-04 22:56:32 UTC (rev 76)
@@ -1,11 +1,12 @@
 ouSim <- function(tree, ...) {
 # right now this is just a switcher, but eventually these should be turned into proper methods of a generic ouSim
   switch(class(tree), 
-         phylo = ouSim.phylo(tree, ...),
-         ouchtree = ouSim.ouchtree(tree, ...),
+         phylo = ouSim.phylo(tree, ...), # original function
+         ouchtree = ouSim.ouchtree(tree, ...), # completed
          brown = ouSim.brownHansen(tree, ...),
          hansen = ouSim.brownHansen(tree, ...),
-         batchHansen = ouSim.batchHansen(tree, ...),
-         "Unrecognized tree class"
+         hansenBatch = ouSim.hansenBatch(tree, ...),
+         hansenSummary = ouSim.hansenBatch(tree, ...),
+         stop("Unrecognized tree class")
          )
 }
\ No newline at end of file