[Mattice-commits] r73 - pkg/R

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

Author: andrew_hipp
Date: 2008-12-04 12:58:55 +0100 (Thu, 04 Dec 2008)
New Revision: 73

Modified:
   pkg/R/ouchSim.R
Log:
continuing to ouchtree-compatibilize ouSim

Modified: pkg/R/ouchSim.R
===================================================================
--- pkg/R/ouchSim.R	2008-12-03 21:16:14 UTC (rev 72)
+++ pkg/R/ouchSim.R	2008-12-04 11:58:55 UTC (rev 73)
@@ -1,51 +1,23 @@
-ouchSim <- function(tree, rootState = 0, shiftBranches = NULL, shiftStates = NULL, alpha = 0, variance = 1, theta = rootState, steps = 1000) {
+ouchSim.ouchtree <- function(tree, rootState = 0, shiftBranches = NULL, shiftStates = NULL, 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 tree
+##   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
-## 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
-## Andrew Hipp (ahipp at mortonarb.org), January 2008 
-## July 2008: modified to accomodate a vector of alpha and theta corresponding to branches
 
-preorderOU <- function(branchList, phy, startNode, startState, alpha, theta) {
-## Recursive function to generate the data under a Brownian motion or OU model (not needed in the Platt model)
-  startBranch = which(phy$edge[,2] == startNode)
-  if(!identical(shiftStates, NULL)) {
-    if(startBranch %in% shiftBranches) startState <- shiftStates[match(startBranch, shiftBranches)] }
-  message(paste('Working on branch',startBranch,'with starting state',startState))
-  branchList[[startBranch]][1] <- startState
-  for (i in 2:length(branchList[[startBranch]])) {
-    branchList[[startBranch]][i] <- branchList[[startBranch]][i - 1] + branchList[[startBranch]][i] + alpha[startBranch] / steps * (theta[startBranch] - branchList[[startBranch]][i - 1]) }
-  endState = branchList[[startBranch]][length(branchList[[startBranch]])]
-  daughterBranches <- phy$edge[which(phy$edge[, 1] == startNode), 2]
-  if(!identical(as.integer(daughterBranches), integer(0))) {
-    for(i in daughterBranches) branchList <- preorderOU(branchList, phy, i, endState, alpha, theta) }
-  return(branchList) }  
-
 ## 1. initialize
 
-if(length(alpha) == 1) alpha <- rep(alpha, length(phy$edge.length))
-if(length(theta) == 1) theta <- rep(theta, length(phy$edge.length))
+if(length(alpha) == 1) alpha <- rep(alpha, tree at nnodes)
+if(length(theta) == 1) theta <- rep(theta, tree at nnodes)
 ## The following creates a list of random draws from the normal distribution, with standard deviation scaled by total tree length and the number of draws for each branch equal to the number of steps in that branch. If there is a separate variance for each branch, I assume the variance is expressed in tree-length units, not branch-length units, so the scaling is the same for all branches (viz., sd = sqrt(variance / steps))
-if(model == "OU") {
-	if(length(variance) == 1) branchList <- lapply(round(phy$edge.length*steps), rnorm, sd = sqrt(variance/steps))
-	  else {
-	    branchList <- vector("list", length(phy$edge.length))
-	    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(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))
-	brLengths <- round(phy$edge.length*steps)
-	brSD <- sqrt(variance)
-	for(i in seq(length(brLengths))) branchList[[i]] <- rnorm(n = brLengths[i], mean = branchMeans[i], sd = brSD[i])
-	}
+if(length(variance) == 1) branchList <- lapply(round(phy$edge.length*steps), rnorm, sd = sqrt(variance/steps))
+else {
+  branchList <- vector("list", length(phy$edge.length))
+  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])
+  }
 timesList <- lapply(branchList, seq)
 rootNode = length(phy$tip.label) + 1
 
@@ -54,15 +26,41 @@
 for (i in 1:length(timesList)) timesList[[i]] <- timesList[[i]] + startTimes[i]
 
 ## 3. traverse
-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
+for(i in which(phy$edge[, 1] == rootNode)) {
+  branchList <- preorderOU(branchList, phy, phy$edge[i,2], rootState, alpha, theta) 
+  }
 
 value <- (list(branchList = branchList, timesList = timesList, steps = steps, parameters = list(rootState = rootState, alpha = alpha, variance = variance, theta = theta))) 
 class(value) <- "ouSim"
 return(value)}
 
+preorderOU <- function(branchList, phy, startNode, startState, alpha, theta) {
+## Recursive function to generate the data under a Brownian motion or OU model
+## modified for ouchtree (s4) Dec 08
+## branch times back from each tip are in tree at epochs, indexed by tip number
+## nodes back from each node (including tips) are in tree at lineages, indexed by node number
+## a branch length is the time of a node - the time of its ancestor
+  startBranch = which(phy$edge[,2] == startNode)
+  if(!identical(shiftStates, NULL)) {
+    if(startBranch %in% shiftBranches) startState <- shiftStates[match(startBranch, shiftBranches)] }
+  message(paste('Working on branch',startBranch,'with starting state',startState))
+  branchList[[startBranch]][1] <- startState
+  for (i in 2:length(branchList[[startBranch]])) {
+    branchList[[startBranch]][i] <- branchList[[startBranch]][i - 1] + branchList[[startBranch]][i] + alpha[startBranch] / steps * (theta[startBranch] - branchList[[startBranch]][i - 1]) }
+  endState = branchList[[startBranch]][length(branchList[[startBranch]])]
+  daughterBranches <- phy$edge[which(phy$edge[, 1] == startNode), 2]
+  if(!identical(as.integer(daughterBranches), integer(0))) {
+    for(i in daughterBranches) branchList <- preorderOU(branchList, phy, i, endState, alpha, theta) }
+  return(branchList) }  
+
+branchLength <- function(tree, endNode) {
+  endNode <- as.character(endNode)
+  ancestor <- tree at ancestors[which(tree at nodes == endNode)]
+  endTime <- ot at times[which(ot at nodes == endNode)]
+  startTime <- ot at times[which(ot at nodes == ancestor)]
+  return(endTime - startTime)
+}
+
 plot.ouSim <- function(ouSim, nodeColor = "blue", nodeDotSize = 1.4, colors = rep("black", length(ouSim$branchList)), ...) {
 ## Plot the data by calling plot(ouSimObject)
 ## To plot different clades, set the colors vector according to the branches in the original