[Mattice-commits] r74 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 4 19:41:30 CET 2008
Author: andrew_hipp
Date: 2008-12-04 19:41:30 +0100 (Thu, 04 Dec 2008)
New Revision: 74
Added:
pkg/R/ouSim.ouchtree.R
pkg/R/ouSim.phylo.R
pkg/R/plot.ouSim.R
Removed:
pkg/R/ouSim.R
pkg/R/ouchSim.R
Log:
ouSim modifications for ouchtree + cleaning up files
Deleted: pkg/R/ouSim.R
===================================================================
--- pkg/R/ouSim.R 2008-12-04 11:58:55 UTC (rev 73)
+++ pkg/R/ouSim.R 2008-12-04 18:41:30 UTC (rev 74)
@@ -1,75 +0,0 @@
-ouSim.phylo <- function(phy, rootState = 0, shiftBranches = NULL, shiftStates = NULL, alpha = 0, variance = 1, theta = rootState, model = "OU", branchMeans = NULL, steps = 1000) {
-## function to plot a simulated dataset under brownian motion or Ornstein-Uhlenbeck (OU) model
-## Arguments:
-## phy is an ape-style 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
-## Dec 2008: This function I'm leaving as is for the time being and just letting the phylo method be as raw as always.
-## New developments will be in the ouchtree, brown, hansen, and hansenBatch methods
-
-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))
-## 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])
- }
-timesList <- lapply(branchList, seq)
-rootNode = length(phy$tip.label) + 1
-
-startTimes <- max(branching.times(phy)) - branching.times(phy) ## WORKS, BUT ASSUMES ULTRAMETRICITY
-startTimes <- round(startTimes[match(phy$edge[,1], names(startTimes))] * steps)
-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
-
-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) }
Copied: pkg/R/ouSim.ouchtree.R (from rev 73, pkg/R/ouchSim.R)
===================================================================
--- pkg/R/ouSim.ouchtree.R (rev 0)
+++ pkg/R/ouSim.ouchtree.R 2008-12-04 18:41:30 UTC (rev 74)
@@ -0,0 +1,75 @@
+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 (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
+
+ ## 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
+
+ ## 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.
+ ## If there is a separate variance for each branch, I assume the variance is expressed in terms of total tree-length,
+ ## so the scaling is the same for all branches
+ sdStep <- sqrt(variance / steps)
+ if(length(variance) == 1) branchList <- lapply(round(brLengths*steps), rnorm, sd = sdStep)
+ else {
+ branchList <- vector("list", length(tree at nnodes))
+ for(i in seq(tree at nnodes)) branchList[[i]] <- rnorm(n = brLengths[i], mean = 0, sd = sdStep[i])
+ }
+ names(branchList) <- tree at nodes # branches are indexed by their end nodes
+ timesList <- lapply(branchList, seq) # creates a list of sequential numbers corresponding to the rnorm vectors
+ startTimes <- round(unlist(sapply(tree at nodes, ancestorTime, tree = tree)) * steps) ## ASSUMES ULTRAMETRICITY
+ for (i in 1:length(timesList)) timesList[[i]] <- timesList[[i]] + startTimes[i]
+
+ ## 3. traverse
+ for(i in which(tree at ancestors == tree at root)) {
+ 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
+## 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)])
+}
+
+branchLength <- function(endNode, tree) {
+ return(nodeTime(endNode, tree) - ancestorTime(endNode, tree))
+}
+
+ancestorTime <- function(endNode, tree) {
+ endNode <- as.character(endNode)
+ ancestor <- tree at ancestors[which(tree at nodes == endNode)]
+ if(is.na(ancestor)) startTime <- NA
+ else startTime <- tree at times[which(tree at nodes == ancestor)]
+ return(startTime)
+}
Copied: pkg/R/ouSim.phylo.R (from rev 72, pkg/R/ouSim.R)
===================================================================
--- pkg/R/ouSim.phylo.R (rev 0)
+++ pkg/R/ouSim.phylo.R 2008-12-04 18:41:30 UTC (rev 74)
@@ -0,0 +1,75 @@
+ouSim.phylo <- function(phy, rootState = 0, shiftBranches = NULL, shiftStates = NULL, alpha = 0, variance = 1, theta = rootState, model = "OU", branchMeans = NULL, steps = 1000) {
+## function to plot a simulated dataset under brownian motion or Ornstein-Uhlenbeck (OU) model
+## Arguments:
+## phy is an ape-style 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
+## Dec 2008: This function I'm leaving as is for the time being and just letting the phylo method be as raw as always.
+## New developments will be in the ouchtree, brown, hansen, and hansenBatch methods
+
+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))
+## 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])
+ }
+timesList <- lapply(branchList, seq)
+rootNode = length(phy$tip.label) + 1
+
+startTimes <- max(branching.times(phy)) - branching.times(phy) ## WORKS, BUT ASSUMES ULTRAMETRICITY
+startTimes <- round(startTimes[match(phy$edge[,1], names(startTimes))] * steps)
+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
+
+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) }
Deleted: pkg/R/ouchSim.R
===================================================================
--- pkg/R/ouchSim.R 2008-12-04 11:58:55 UTC (rev 73)
+++ pkg/R/ouchSim.R 2008-12-04 18:41:30 UTC (rev 74)
@@ -1,71 +0,0 @@
-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 (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
-
-## 1. initialize
-
-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(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
-
-startTimes <- max(branching.times(phy)) - branching.times(phy) ## WORKS, BUT ASSUMES ULTRAMETRICITY
-startTimes <- round(startTimes[match(phy$edge[,1], names(startTimes))] * steps)
-for (i in 1:length(timesList)) timesList[[i]] <- timesList[[i]] + startTimes[i]
-
-## 3. traverse
-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
-## 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) }
Added: pkg/R/plot.ouSim.R
===================================================================
--- pkg/R/plot.ouSim.R (rev 0)
+++ pkg/R/plot.ouSim.R 2008-12-04 18:41:30 UTC (rev 74)
@@ -0,0 +1,8 @@
+plot.ouSim <- 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)
+ 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)
+}
\ No newline at end of file
More information about the Mattice-commits
mailing list