[Mattice-commits] r193 - in pkg: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 29 22:18:09 CEST 2009
Author: andrew_hipp
Date: 2009-05-29 22:18:08 +0200 (Fri, 29 May 2009)
New Revision: 193
Modified:
pkg/R/ouSim.hansenBatch.R
pkg/R/ouSim.ouchtree.R
pkg/R/plot.ouSim.R
pkg/man/carex.Rd
pkg/man/ouSim.Rd
Log:
modifications to deal with multiple trees (new data set) and very short branches, which turn out to cause difficulties in ouSim.
Modified: pkg/R/ouSim.hansenBatch.R
===================================================================
--- pkg/R/ouSim.hansenBatch.R 2009-05-29 18:50:13 UTC (rev 192)
+++ pkg/R/ouSim.hansenBatch.R 2009-05-29 20:18:08 UTC (rev 193)
@@ -5,7 +5,7 @@
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, ], ...)
class(outdata) <- "ouSim"
- return(outdata)
+ return(outdata)t
}
ouSim.hansenBatch <- function(object, ...) ouSim(summary(object), ...)
Modified: pkg/R/ouSim.ouchtree.R
===================================================================
--- pkg/R/ouSim.ouchtree.R 2009-05-29 18:50:13 UTC (rev 192)
+++ pkg/R/ouSim.ouchtree.R 2009-05-29 20:18:08 UTC (rev 193)
@@ -5,31 +5,32 @@
## alpha and theta are either single values or vectors of length (length(branchList))
tree <- object
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) {
-## 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
- # 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]) # denom was mult'd by steps... should be?
- }
- 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)
-}
-##--------------------------------------
+
+ ##embedded function---------------------
+ ##could 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 the branch starting this process
+ 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]) # denom was mult'd by steps... should be?
+ }
+ 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
@@ -45,7 +46,7 @@
## 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)
+ if(length(variance) == 1) branchList <- lapply(sapply(round(brLengths*steps), max, 1), rnorm, sd = sdStep) # makes each branch a minimum of 2 steps long
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])
@@ -60,7 +61,7 @@
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)))
+ value <- list(branchList = branchList, timesList = timesList, steps = steps, parameters = list(rootState = rootState, alpha = alpha, variance = variance, theta = theta))
class(value) <- "ouSim"
return(value)
}
Modified: pkg/R/plot.ouSim.R
===================================================================
--- pkg/R/plot.ouSim.R 2009-05-29 18:50:13 UTC (rev 192)
+++ pkg/R/plot.ouSim.R 2009-05-29 20:18:08 UTC (rev 193)
@@ -1,4 +1,4 @@
-plot.ouSim <- function(x, nodeColor = "blue", nodeDotSize = 1.4, colors = NULL, ...) {
+trplot.ouSim <- function(x, nodeColor = "blue", nodeDotSize = 1.4, colors = NULL, ...) {
## To plot different clades, set the colors vector according to the branches in the original
## only passes the ... along to lines
branches = length(x$branchList)
Modified: pkg/man/carex.Rd
===================================================================
--- pkg/man/carex.Rd 2009-05-29 18:50:13 UTC (rev 192)
+++ pkg/man/carex.Rd 2009-05-29 20:18:08 UTC (rev 193)
@@ -8,7 +8,7 @@
}
\usage{data(carex)}
\format{
- A list with three items:
+ A list with four items:
\item{\code{ovales.tree}}{
An ultrametric tree in \code{phylo} format with 53 tips
}
@@ -18,6 +18,10 @@
\item{\code{ovales.nodes}}{
A list of eight taxon vectors defining the eight nodes studied.
}
+ \item{\code{ovales.bayesTrees}}{
+ A list of 100 ultrametric trees in \code{phylo} format, subsampled from the MCMC analysis underlying
+ \code{carex$ovales.bayesTrees}.
+ }
}
\details{
Phylogeny (\code{ovales.tree}) was estimated for approximately 80 species (Hipp 2006), branch lengths optimized using
@@ -43,14 +47,14 @@
\examples{
library(maticce)
data(carex)
- ovales.tree <- ape2ouch(carex$ovales.tree) # tree comes in in phylo format, but we need an ouchtree object
- trial <- runBatchHansen(ovales.tree, carex$ovales.data, carex$ovales.nodes[1:4], maxNodes = 2) # for expedience, only tests for changes at up to 2 of the first 4 nodes
+ trees <- lapply(carex$ovales.bayesTrees[1:10], ape2ouch) # tree comes in in phylo format, but we need an ouchtree object
+ trial <- runBatchHansen(trees, carex$ovales.data, carex$ovales.nodes[1:4], maxNodes = 2) # for expedience, only tests for changes at up to 2 of the first 4 nodes
summary(trial) # summarizes results using summary.hansenBatch and displays the results using print.hansenSummary
multiModel(carex$ovales.tree, carex$ovales.data, carex$ovales.nodes[[2]]) # compares five different models of character change at node 2
- trialSim <- ouSim(trial, ovales.tree) # simulates the evolution of the chromosome number under the model-averaged values
+ trialSim <- ouSim(trial, trees[[1]], treeNum = 1) # simulates the evolution of the chromosome number under the model-averaged values, using tree # 1 as the guide tree
plot(trialSim) # plots the character simulation, with all branches black
- plot(trialSim, colors = paintBranches(carex$ovales.nodes[1:4], ovales.tree)) # plots the character simulation, with branch colors changing at all 8 nodes
- plot(trialSim, colors = paintBranches(list(carex$ovales.nodes[[2]]), ovales.tree)) # plots the character simulation, with branch colors changing only at node 2
+ plot(trialSim, colors = paintBranches(carex$ovales.nodes[1:4], trees[[1]])) # plots the character simulation, with branch colors changing at all 8 nodes
+ plot(trialSim, colors = paintBranches(list(carex$ovales.nodes[[2]]), trees[[1]])) # plots the character simulation, with branch colors changing only at node 2
}
\author{Andrew L. Hipp <ahipp at mortonarb.org>}
\keyword{datasets}
\ No newline at end of file
Modified: pkg/man/ouSim.Rd
===================================================================
--- pkg/man/ouSim.Rd 2009-05-29 18:50:13 UTC (rev 192)
+++ pkg/man/ouSim.Rd 2009-05-29 20:18:08 UTC (rev 193)
@@ -42,7 +42,7 @@
}
\item{\code{hansenBatch}}{
Model-averaged parameters from the \code{hansenBatch} object are used for analysis.
- One of the trees used for analysis must be provided, and a corrsponding tree number must be provided so that
+ One of the trees used for analysis must be provided, and a corresponding tree number must be provided so that
branches are indexed correctly.
}
@@ -51,10 +51,10 @@
}
\item{\code{ouSim.phylo}}{
- A very basic simulation engine, but also the most flexible. As written, the user has to specify the
- model using two vectors that correspond to the branches in an \code{ape}-format tree. It is important to note that
- this simulation method is really a heuristic device, not appropriate for estimating parameter distributions. For
- analysis purposes, you should utilize the \code{simulate} and \code{bootstrap} methods in \pkg{ouch}.
+ A very basic simulation engine, but also the most flexible. As written, the user has to specify the
+ model using two vectors that correspond to the branches in an \code{ape}-format tree. It is important to note that
+ this simulation method is really a heuristic device, not appropriate for estimating parameter distributions. For
+ analysis purposes, you should utilize the \code{simulate} and \code{bootstrap} methods in \pkg{ouch}.
}
}
\arguments{
More information about the Mattice-commits
mailing list