[Mattice-commits] r78 - pkg/R

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

Author: andrew_hipp
Date: 2008-12-11 16:04:21 +0100 (Thu, 11 Dec 2008)
New Revision: 78

Modified:
   pkg/R/ouSim.hansenBatch.R
   pkg/R/ouSim.ouchtree.R
   pkg/R/ouSim.phylo.R
   pkg/R/ouSimHead.R
Log:
ouSim working correctly for all classes

Modified: pkg/R/ouSim.hansenBatch.R
===================================================================
--- pkg/R/ouSim.hansenBatch.R	2008-12-05 14:25:54 UTC (rev 77)
+++ pkg/R/ouSim.hansenBatch.R	2008-12-11 15:04:21 UTC (rev 78)
@@ -2,7 +2,27 @@
 ## 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, ], ...)
+  if(identical(rootState, NULL)) rootState <- analysis$thetaMatrix[treeNum, ][tree at root] # rootstate taken to be the optimum at the root
+  outdata <- ouSim(tree, rootState, alpha = analysis$modelAvgAlpha, variance = analysis$modelAvgSigmaSq, theta = analysis$thetaMatrix[treeNum, ], ...)
   return(outdata)
+}
+
+ouSim.brownHansen <- function(analysis, ...) {
+  su <- summary(analysis)
+  if(length(analysis at regimes) > 1) warning("Theta is based on analysis at regimes[[1]]")
+  if(dim(su$alpha)[1] != 1) stop("This is a one-character simulation; analysis appears to be based on > 1 character")
+  if(class(analysis) == "browntree") {
+    alpha <- 0
+    theta <- 0
+    rootState <- su$theta[[1]]
+    }
+  if(class(analysis) == "hansentree") {
+    alpha <- as.vector(su$alpha)
+    theta <- su$optima[[1]][analysis at regimes[[1]]]
+    rootState <- theta[analysis at root] # rootstate taken to be the optimum at the root
+    }
+  variance <- as.vector(su$sigma.squared)
+  tree <- ouchtree(analysis at nodes, analysis at ancestors, analysis at times) 
+  outdata <- ouSim.ouchtree(tree, rootState, alpha, variance, theta, ...)
+  return(outdata)
 }
\ No newline at end of file

Modified: pkg/R/ouSim.ouchtree.R
===================================================================
--- pkg/R/ouSim.ouchtree.R	2008-12-05 14:25:54 UTC (rev 77)
+++ pkg/R/ouSim.ouchtree.R	2008-12-11 15:04:21 UTC (rev 78)
@@ -3,7 +3,7 @@
 ## Arguments:
 ##   tree is an ouch-style (S4) tree
 ##   alpha and theta are either single values or vectors of length (length(branchList))
-
+message(paste("running sim with root =", rootState, ", alpha =", mean(alpha), ", var =", variance, "theta =", mean(theta)))
 ##embedded function---------------------
 ##can be released to the wild, but more arguments will need to be passed around
 preorderOU <- function(branchList, tree, startNode, startState, alpha, theta) {
@@ -14,7 +14,7 @@
 ## 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
-  message(paste('Working on branch',startBranch,'with starting state',startState))
+  # message(paste('Working on branch',startBranch,'with starting state',startState))
   workingBranch <- branchList[[startBranch]]
   workingBranch[1] <- startState
   for (brStep in 2:length(workingBranch)) {

Modified: pkg/R/ouSim.phylo.R
===================================================================
--- pkg/R/ouSim.phylo.R	2008-12-05 14:25:54 UTC (rev 77)
+++ pkg/R/ouSim.phylo.R	2008-12-11 15:04:21 UTC (rev 78)
@@ -7,7 +7,7 @@
 ##   shiftStates is a vector of length = length(shiftBranches) indicaing the ancestral states for the determined break points
 ## Models:
 ##  "OU" is a brownian motion or OU model 
-##  "platt" is the Platt-Molitor model, in which the only phylogenetic effect is the mean and variance for a given branch
+##  "meanVar" is a model in which the only phylogenetic effect is the mean and variance for a given branch
 ## Andrew Hipp (ahipp at mortonarb.org), January 2008 
 ## July 2008: modified to accomodate a vector of alpha and theta corresponding to branches
 ## Dec 2008: This function I'm leaving as is for the time being and just letting the phylo method be as raw as always.
@@ -40,7 +40,7 @@
 	    brLengths <- round(phy$edge.length*steps)
 	    brSD <- sqrt(variance/steps)
 	    for(i in seq(length(brLengths))) branchList[[i]] <- rnorm(n = brLengths[i], mean = 0, sd = brSD[i]) }}
-if(model == "platt") {
+if(model == "meanVar") {
 	if(length(variance) == 1) variance <- rep(variance, length(phy$edge.length))
 	branchList <- vector("list", length(phy$edge.length))
 	if(length(branchMeans) == 1) branchMeans <- rep(branchMeans, length(phy$edge))
@@ -59,7 +59,7 @@
 if(model == "OU") {
 	for(i in which(phy$edge[, 1] == rootNode)) {
 	  branchList <- preorderOU(branchList, phy, phy$edge[i,2], rootState, alpha, theta) }}
-if(model == "platt") branchList <- branchList
+if(model == "meanVar") branchList <- branchList
 
 value <- (list(branchList = branchList, timesList = timesList, steps = steps, parameters = list(rootState = rootState, alpha = alpha, variance = variance, theta = theta))) 
 class(value) <- "ouSimPhylo"

Modified: pkg/R/ouSimHead.R
===================================================================
--- pkg/R/ouSimHead.R	2008-12-05 14:25:54 UTC (rev 77)
+++ pkg/R/ouSimHead.R	2008-12-11 15:04:21 UTC (rev 78)
@@ -3,8 +3,8 @@
   switch(class(tree), 
          phylo = ouSim.phylo(tree, ...), # original function
          ouchtree = ouSim.ouchtree(tree, ...), # completed
-         brown = ouSim.brownHansen(tree, ...),
-         hansen = ouSim.brownHansen(tree, ...),
+         browntree = ouSim.brownHansen(tree, ...), # completed
+         hansentree = ouSim.brownHansen(tree, ...),
          hansenBatch = ouSim.hansenBatch(tree, ...), # completed
          hansenSummary = ouSim.hansenBatch(tree, ...), # completed
          stop("Unrecognized tree class")