[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