[Mattice-commits] r20 - pkg/R

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Nov 13 21:50:48 CET 2008 

Author: andrew_hipp
Date: 2008-11-13 21:50:48 +0100 (Thu, 13 Nov 2008)
New Revision: 20

Modified:
   pkg/R/revisedBatchHansen.R
Log:
Cleaned up batchHansen... now should call hansen correctly

Modified: pkg/R/revisedBatchHansen.R
===================================================================
--- pkg/R/revisedBatchHansen.R	2008-11-13 20:50:16 UTC (rev 19)
+++ pkg/R/revisedBatchHansen.R	2008-11-13 20:50:48 UTC (rev 20)
@@ -33,9 +33,11 @@
 #  "ouchTrees" = list of OUCH-style trees; if a single tree, send as 'list(TREE)'
 #  "characterStates" = vector of character states, either extracted from an ouch-style tree data.frame or a named vector
 #  REMOVED: "SEM"= standard error of the mean, vector extracted from an ouch-style tree data.frame
-#  REMOVED: "scalingFactor" = factor to multiply against (times / max(times)) -- choose based on trial analyses
+#  "rescale" = factor to multiply against (times / max(times)) -- choose based on trial analyses; set at <= 0 if you don't want to rescale trees
 #  "cladeMembersList" = list of vectors containing names of the members of each clade (except for the root of the tree)
-function(ouchTrees, characterStates, cladeMembersList, regimeNames = NULL, logData = FALSE) {
+#  "brown" = whether to analyse the data under a Brownian motion model
+#  "..." = additional arguments to pass along to hansen
+function(ouchTrees, characterStates, cladeMembersList, maxNodes = NULL, regimeTitles = NULL, logData = FALSE, brown = F, rescale = 1, ...) {
   ## do all the objects in ouchTrees inherit ouchtree?
   treeCheck <- unlist(lapply(ouchTrees, function(x) is(x,'ouchtree')))
   if(F %in% treeCheck) 
@@ -65,14 +67,19 @@
   stopFlag <- F  
 
   if(logData) characterStates <- log(characterStates)
-  if(identical(regimeNames, NULL)) regimeNames <- c(as.character(1:length(regimesList)), "brown")
-  
+  if(identical(regimeTitles, NULL)) {
+    regimeTitles <- as.character(1:length(regimesList))
+    if(brown) regimeTitles <- c(regimeTitles, 'brown')
+  }
   hansenBatch = vector("list", length(ouchTrees))
   treeCounter = 0
   for (i in 1:length(ouchTrees)) {
     tree <- ouchTrees[[i]]
-    regimesList = regimeVectors(tree, cladeMembersList)
+    regimesList = regimeVectors(tree, cladeMembersList, maxNodes)
     
+    ## rescale tree if requested
+    if(rescale>0) tree at times <- rescale * tree at times / max(tree at times) 
+    
     ## need to revise regimeVectors so that it only returns regimes for nodes that are supported in the tree
     ## for now, assume nodes of interest are present in all trees
 
@@ -85,10 +92,11 @@
     
     ## send it off to batchHansen and just stick the results in hansenBatch... this won't work as the number of regimes gets large, 
     ##   so there should be some option here to just hang onto the coefficients for each run (i.e., hang onto 'coef(hansen(...))' rather than 'hansen(...)')
-    hansenBatch[[i]] = batchHansen(tree, dataIn, regimesList, regimeNames, numberOfTermini = tree at nterm)
+    ##   there could also be an option to save the entire object as a series of files in addition to hanging onto 
+    hansenBatch[[i]] = batchHansen(tree, dataIn, regimesList, regimeTitles, brown, ...)
     message(paste("Tree",i,"of",length(ouchTrees),"complete"))}
     
-    ## right now no summary is returned; it should be
+    ## right now no summary is returned; one is needed, summarizing over trees what is summarized for each tree in batchHansen
   
   return(list(hansens = hansenBatch, regimes = regimesList)) }
 
@@ -98,37 +106,18 @@
 #  "node" "ancestor" "times" "data" = the standard tree specification vectors of the OUCH-style tree
 #  "regimesList" = list of regime-paintings as output from regimeVectors
 #  "scalingFactor" = factor to multiply against (times / max(times)) -- choose based on trial analyses
-# Value: a matrix with nrow = regimes + 1 and columns for u, d.f., all estimated parameters, LRvsBM, AIC, and AIC weight
-# ADDED "error" on 2 march 07 to accomodate the "me" in Pienaar's ouch
-function(tree, data, regimesList, regimeTitles, numberOfTermini) {
-  variables = c("loglik", "df", "alpha", "sigma", "theta0", "theta1", "theta2", "theta3", "theta4", "theta5","theta6","theta7","theta8","theta9", "LRvsBM", "AIC", "AICweight", "deltaAIC", "exp(-0.5*deltaAIC)", "AICc", "AICcWeight", "deltaAICc", "exp(-0.5*deltaAICc)")
-  treeData = matrix(data = NA, nrow = (length(regimesList) + 1), ncol = length(variables),
-             dimnames = list(regimeTitles,variables))
-  brownianResults = brown.fit(data, error, node, ancestor, (times/max(times))*scalingFactor)
-  treeData[length(regimesList) + 1, "AIC"] = brownianResults$aic
-  treeData[length(regimesList) + 1, "loglik"] = brownianResults$loglik
-  treeData[length(regimesList) + 1, "sigma"] = brownianResults$sigma
-  treeData[length(regimesList) + 1, "theta0"] = brownianResults$theta
-  treeData[length(regimesList) + 1, "df"] = brownianResults$df
-  treeData[length(regimesList) + 1, "AICc"] = brownianResults$aic + ((2 * brownianResults$df * (brownianResults$df + 1)) / (numberOfTermini - brownianResults$df - 1))
-  for (i in 1:length(regimesList)) {
+# Value: a matrix with nrow = regimes (+ 1 if brownian model is included) and columns for u, d.f., all estimated parameters, LRvsBM, AIC, and AIC weight
+function(tree, data, regimesList, regimeTitles, brown, ...) {
+  n <- tree at nterm
+  ## set up a matrix that returns lnL, K, sigmasq, theta0, and alpha for every model; thetas will go along into a list that is indexed by model
+  variables <- c("loglik", "dof", "sigma.squared", "theta", "alpha") # it's important that 'alpha' go last so that the matrix fills up right when the brownian motion model is used
+  treeData <- matrix(data = NA, nrow = length(regimeTitles), ncol = length(variables), dimnames = list(regimeTitles,variables))
+  if(brown) br <- brown(data, tree)
+  treeData["brown", ] <- unlist(summary(br)[variables])
+  for (i in seq(regimesList)) {
     message(paste("Running regime",i))
-    tempHansen = hansen.fit(data, error, node, ancestor, (times/max(times))*scalingFactor, regimesList[[i]])
-    treeData[i,"loglik"] = tempHansen$loglik
-    treeData[i,"df"] = tempHansen$df
-    treeData[i,"alpha"] = tempHansen$alpha
-    treeData[i,"sigma"] = tempHansen$sigma
-    treeData[i,"LRvsBM"] = 2 * ((tempHansen$loglik / 2) - (treeData[length(regimesList) + 1, "loglik"] / 2))
-    treeData[i,"AIC"] = tempHansen$aic
-    treeData[i,"AICc"] = tempHansen$aic + ((2 * tempHansen$df * (tempHansen$df + 1)) / (numberOfTermini - tempHansen$df - 1))
-    for (j in 0:(length(names(tempHansen$theta))-1)) {
-      treeData[i, paste("theta",j,sep = "")] = tempHansen$theta[[j+1]]}}
-  for (i in 1:(length(regimesList) + 1)) {
-    treeData[i, "deltaAIC"] = treeData[i, "AIC"] - min(treeData[1:(length(regimesList) + 1), "AIC"])
-    treeData[i, "exp(-0.5*deltaAIC)"] = exp(-0.5 * treeData[i,"deltaAIC"])
-    treeData[i, "deltaAICc"] = treeData[i, "AICc"] - min(treeData[1:(length(regimesList) + 1), "AICc"])
-    treeData[i, "exp(-0.5*deltaAICc)"] = exp(-0.5 * treeData[i,"deltaAICc"]) }
-  for (i in 1:(length(regimesList) + 1)) {
-    treeData[i, "AICweight"] = treeData[i, "exp(-0.5*deltaAIC)"] / sum(treeData[1:(length(regimesList) + 1), "exp(-0.5*deltaAIC)"])
-    treeData[i, "AICcWeight"] = treeData[i, "exp(-0.5*deltaAICc)"] / sum(treeData[1:(length(regimesList) + 1), "exp(-0.5*deltaAICc)"]) }
+    ## at this point, the user has to give an initial alpha and sigma for hansen to search on... this should be relaxed
+    ha = hansen(data, tree, regimesList[[i]], ...)
+    treeData[i, ] <- unlist(summary(ha)[variables])
+  }
   return(treeData) }

More information about the Mattice-commits mailing list