[Mattice-commits] r251 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Mar 22 17:53:41 CET 2011
Author: andrew_hipp
Date: 2011-03-22 17:53:41 +0100 (Tue, 22 Mar 2011)
New Revision: 251
Modified:
pkg/R/ouSim.hansenBatch.R
pkg/R/ouSim.ouchtree.R
pkg/R/ouSim.phylo.R
pkg/R/summarizingAnalyses.R
Log:
Fixing ouSim so that it uses alpha, not sqrt.alpha; this was causing an error in the vignette
Modified: pkg/R/ouSim.hansenBatch.R
===================================================================
--- pkg/R/ouSim.hansenBatch.R 2011-03-15 10:14:51 UTC (rev 250)
+++ pkg/R/ouSim.hansenBatch.R 2011-03-22 16:53:41 UTC (rev 251)
@@ -1,9 +1,9 @@
ouSim.hansenSummary <- function(object, tree, treeNum = 1, rootState = NULL, ...) {
-## runs ouSim.ouchtree for a hansenBatch or hansenSummary object, using the model-averaged sqrt.alpha, sigma.squared, and theta vector from one tree
+## runs ouSim.ouchtree for a hansenBatch or hansenSummary object, using the model-averaged alpha, sigma.squared, and theta vector from one tree
analysis <- object
# if(class(analysis) == "hansenBatch") analysis <- summary(analysis)
if(identical(rootState, NULL)) rootState <- analysis$thetaMatrix[treeNum, ][tree at root] # rootstate taken to be the optimum at the root
- outdata <- ouSim(tree, rootState, sqrt.alpha = analysis$modelAvgAlpha, variance = analysis$modelAvgSigmaSq, theta = analysis$thetaMatrix[treeNum, ], ...)
+ outdata <- ouSim(tree, rootState, alpha = analysis$modelAvgAlpha, variance = analysis$modelAvgSigmaSq, theta = analysis$thetaMatrix[treeNum, ], ...)
class(outdata) <- "ouSim"
return(outdata)
}
@@ -14,13 +14,13 @@
analysis <- object
su <- summary(analysis)
if(length(analysis at regimes) > 1) warning("Theta is based on analysis at regimes[[1]]")
- if(dim(su$sqrt.alpha)[1] != 1) stop("This is a one-character simulation; analysis appears to be based on > 1 character")
- sqrt.alpha <- as.vector(su$sqrt.alpha)
+ if(dim(su$alpha)[1] != 1) stop("This is a one-character simulation; analysis appears to be based on > 1 character")
+ 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, sqrt.alpha, variance, theta, ...)
+ outdata <- ouSim.ouchtree(tree, rootState, alpha, variance, theta, ...)
outdata$colors <- analysis at regimes[[1]]
class(outdata) <- "ouSim"
return(outdata)
@@ -30,13 +30,13 @@
analysis <- object
su <- summary(analysis)
if(length(analysis at regimes) > 1) warning("Theta is based on analysis at regimes[[1]]")
- if(dim(su$sqrt.alpha)[1] != 1) stop("This is a one-character simulation; analysis appears to be based on > 1 character")
- sqrt.alpha <- 0
+ if(dim(su$alpha)[1] != 1) stop("This is a one-character simulation; analysis appears to be based on > 1 character")
+ alpha <- 0
theta <- 0
rootState <- su$theta[[1]]
variance <- as.vector(su$sigma.squared)
tree <- ouchtree(analysis at nodes, analysis at ancestors, analysis at times)
- outdata <- ouSim.ouchtree(tree, rootState, sqrt.alpha, variance, theta, ...)
+ outdata <- ouSim.ouchtree(tree, rootState, alpha, variance, theta, ...)
outdata$colors <- analysis at regimes[[1]]
class(outdata) <- "ouSim"
return(outdata)
Modified: pkg/R/ouSim.ouchtree.R
===================================================================
--- pkg/R/ouSim.ouchtree.R 2011-03-15 10:14:51 UTC (rev 250)
+++ pkg/R/ouSim.ouchtree.R 2011-03-22 16:53:41 UTC (rev 251)
@@ -1,14 +1,14 @@
-ouSim.ouchtree <- function(object, rootState = 0, sqrt.alpha = 0, variance = 1, theta = rootState, steps = 1000, ...) {
+ouSim.ouchtree <- function(object, rootState = 0, alpha = 0, variance = 1, theta = rootState, steps = 1000, ...) {
## function to plot a simulated dataset under brownian motion or Ornstein-Uhlenbeck (OU) model
## Arguments:
## object is an ouch-style (S4) tree
-## sqrt.alpha and theta are either single values or vectors of length (length(branchList))
+## alpha and theta are either single values or vectors of length (length(branchList))
tree <- object
-message(paste("running sim with root =", rootState, ", sqrt.alpha =", mean(sqrt.alpha), ", var =", variance, "theta =", mean(theta)))
+message(paste("running sim with root =", rootState, ", alpha =", mean(alpha), ", var =", variance, "theta =", mean(theta)))
## embedded function---------------------
## could be released to the wild, but more arguments would need to be passed around
- preorderOU <- function(branchList, tree, startNode, startState, sqrt.alpha, theta) {
+ 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
@@ -21,24 +21,24 @@
else {
for (brStep in 1:length(workingBranch)) {
workingBranch[brStep] <-
- startState + workingBranch[brStep] + sqrt.alpha[startBranch] / steps * (theta[startBranch] - startState) # denom was mult'd by steps... should be?
+ startState + workingBranch[brStep] + alpha[startBranch] / steps * (theta[startBranch] - startState) # denom was mult'd by steps... should be?
startState <- workingBranch[brStep]
}
branchList[[startBranch]] <- workingBranch
endState <- branchList[[startBranch]][length(branchList[[startBranch]])]
}
if(!identical(as.integer(daughterBranches), integer(0))) {
- for(i in daughterBranches) branchList <- preorderOU(branchList, tree, i, endState, sqrt.alpha, theta) }
+ for(i in daughterBranches) branchList <- preorderOU(branchList, tree, i, endState, alpha, theta) }
return(branchList)
}
## --------------------------------------
## 1. initialize
- if(length(sqrt.alpha) == 1) sqrt.alpha <- rep(sqrt.alpha, tree at nnodes)
+ 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(sqrt.alpha) <- tree at nodes
+ 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
@@ -60,10 +60,10 @@
## 3. traverse
for(i in which(tree at ancestors == tree at root)) { ## calls preorderOU for each descendent from the root.
- branchList <- preorderOU(branchList, tree, tree at nodes[i], rootState, sqrt.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, sqrt.alpha = sqrt.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/ouSim.phylo.R
===================================================================
--- pkg/R/ouSim.phylo.R 2011-03-15 10:14:51 UTC (rev 250)
+++ pkg/R/ouSim.phylo.R 2011-03-22 16:53:41 UTC (rev 251)
@@ -1,19 +1,19 @@
-ouSim.phylo <- function(object, rootState = 0, shiftBranches = NULL, shiftStates = NULL, sqrt.alpha = 0, variance = 1, theta = rootState, model = "OU", branchMeans = NULL, steps = 1000, ...) {
+ouSim.phylo <- function(object, 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
-## sqrt.alpha and theta are either single values or vectors of length (length(branchList))
+## 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
## "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 sqrt.alpha and theta corresponding to branches
+## 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
phy <- object
-preorderOU <- function(branchList, phy, startNode, startState, sqrt.alpha, theta) {
+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)) {
@@ -21,16 +21,16 @@
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] + sqrt.alpha[startBranch] / steps * (theta[startBranch] - branchList[[startBranch]][i - 1]) }
+ 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, sqrt.alpha, theta) }
+ for(i in daughterBranches) branchList <- preorderOU(branchList, phy, i, endState, alpha, theta) }
return(branchList) }
## 1. initialize
-if(length(sqrt.alpha) == 1) sqrt.alpha <- rep(sqrt.alpha, length(phy$edge.length))
+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") {
@@ -58,10 +58,10 @@
## 3. traverse
if(model == "OU") {
for(i in which(phy$edge[, 1] == rootNode)) {
- branchList <- preorderOU(branchList, phy, phy$edge[i,2], rootState, sqrt.alpha, theta) }}
+ branchList <- preorderOU(branchList, phy, phy$edge[i,2], rootState, alpha, theta) }}
if(model == "meanVar") branchList <- branchList
-value <- (list(branchList = branchList, timesList = timesList, steps = steps, parameters = list(rootState = rootState, sqrt.alpha = sqrt.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)
}
\ No newline at end of file
Modified: pkg/R/summarizingAnalyses.R
===================================================================
--- pkg/R/summarizingAnalyses.R 2011-03-15 10:14:51 UTC (rev 250)
+++ pkg/R/summarizingAnalyses.R 2011-03-22 16:53:41 UTC (rev 251)
@@ -17,7 +17,7 @@
nnodes <- length(nodeSums) # number of nodes being studied
nodes <- dimnames(hansenBatch$regMatrix$overall)[[2]] # grab the overall regMatrix, which includes all possible nodes
sigmaSqVector <- numeric(ntrees) # vector to capture model-averaged sigma^2 for each tree
- alphaVector <- numeric(ntrees) # vector to capture model-averaged sqrt.alpha for each tree # s/b sqrt.alpha
+ alphaVector <- numeric(ntrees) # vector to capture model-averaged alpha for each tree
modelsMatrix <- vector('list', ntrees) # list of matrices, indexed by tree, holding the weight for each model
matrixRows <- c('AICwi', 'AICcwi', 'BICwi') # rows in the matrix
nodeWeightsSummed <- matrix(0, nrow = length(matrixRows), ncol = nnodes, dimnames = list(matrixRows, nodes)) # holds node weights summed; zero-filled b/c it is a sum?
More information about the Mattice-commits
mailing list