[Mattice-commits] r75 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 4 23:24:53 CET 2008
Author: andrew_hipp
Date: 2008-12-04 23:24:53 +0100 (Thu, 04 Dec 2008)
New Revision: 75
Modified:
pkg/R/ouSim.ouchtree.R
pkg/R/ouSim.phylo.R
pkg/R/plot.ouSim.R
Log:
ouSim.ouchtree works for an ouchtree object, but regime must be specified manually
Modified: pkg/R/ouSim.ouchtree.R
===================================================================
--- pkg/R/ouSim.ouchtree.R 2008-12-04 18:41:30 UTC (rev 74)
+++ pkg/R/ouSim.ouchtree.R 2008-12-04 22:24:53 UTC (rev 75)
@@ -1,4 +1,4 @@
-ouchSim.ouchtree <- function(tree, rootState = 0, shiftBranches = NULL, shiftStates = NULL, alpha = 0, variance = 1, theta = rootState, steps = 1000) {
+ouSim.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 (S4) tree
@@ -6,11 +6,43 @@
## 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
+preorderOU <- function(branchList, tree, 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
+## 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])
+ }
+ branchList[[startBranch]] <- workingBranch
+ endState <- branchList[[startBranch]][length(branchList[[startBranch]])]
+ daughterBranches <- tree at nodes[which(tree at ancestors == startNode)]
+ if(!identical(as.integer(daughterBranches), integer(0))) {
+ for(i in daughterBranches) branchList <- preorderOU(branchList, tree, i, endState, alpha, theta) }
+ return(branchList)
+}
+##--------------------------------------
+ #debug(preorderOU)
+
## 1. initialize
if(length(alpha) == 1) alpha <- rep(alpha, tree at nnodes)
if(length(theta) == 1) theta <- rep(theta, tree at nnodes)
brLengths <- c(0, unlist(lapply(2:tree at nnodes, branchLength, tree = tree))) # assumes first node is root; this should be relaxed
names(brLengths) <- tree at nodes # branches are indexed by end node
+ names(alpha) <- tree at nodes
+ names(theta) <- tree at nodes
## 2. 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.
@@ -29,35 +61,14 @@
## 3. traverse
for(i in which(tree at ancestors == tree at root)) {
- branchList <- preorderOU(branchList, phy, phy$edge[i,2], rootState, alpha, theta)
+ branchList <- preorderOU(branchList, tree, tree at nodes[i], 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
-## not fixed yet for ouchtree
-
- 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)
-}
-
nodeTime <- function(node, tree) {
return(tree at times[which(tree at nodes == node)])
}
Modified: pkg/R/ouSim.phylo.R
===================================================================
--- pkg/R/ouSim.phylo.R 2008-12-04 18:41:30 UTC (rev 74)
+++ pkg/R/ouSim.phylo.R 2008-12-04 22:24:53 UTC (rev 75)
@@ -62,14 +62,6 @@
if(model == "platt") branchList <- branchList
value <- (list(branchList = branchList, timesList = timesList, steps = steps, parameters = list(rootState = rootState, alpha = alpha, variance = variance, theta = theta)))
-class(value) <- "ouSim"
-return(value)}
-
-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
-## only passes the ... along to lines
- branches = length(ouSim$branchList)
- plot(1:ouSim$steps, ylim = range(unlist(ouSim$branchList)), type = "n", ylab = "Trait value", xlab = "Time")
- for(i in 1:branches) lines(ouSim$timesList[[i]], ouSim$branchList[[i]], col = colors[i], ...)
- for(i in 1:branches) points(ouSim$timesList[[i]][1], ouSim$branchList[[i]][1], pch = 19, col = nodeColor, cex = nodeDotSize) }
+class(value) <- "ouSimPhylo"
+return(value)
+}
\ No newline at end of file
Modified: pkg/R/plot.ouSim.R
===================================================================
--- pkg/R/plot.ouSim.R 2008-12-04 18:41:30 UTC (rev 74)
+++ pkg/R/plot.ouSim.R 2008-12-04 22:24:53 UTC (rev 75)
@@ -1,4 +1,4 @@
-plot.ouSim <- function(ouSim, nodeColor = "blue", nodeDotSize = 1.4, colors = rep("black", length(ouSim$branchList)), ...) {
+plot.ouSimPhylo <- function(ouSim, nodeColor = "blue", nodeDotSize = 1.4, colors = rep("black", length(ouSim$branchList)), ...) {
## To plot different clades, set the colors vector according to the branches in the original
## only passes the ... along to lines
branches = length(ouSim$branchList)
More information about the Mattice-commits
mailing list