[Mattice-commits] r63 - pkg/R

[noreply at r-forge.r-project.org]

Author: andrew_hipp
Date: 2008-12-02 09:46:00 +0100 (Tue, 02 Dec 2008)
New Revision: 63

Modified:
   pkg/R/batchHansen.R
Log:
batchHansen now creates a model by branch theta matrix for each tree, to be used for model-averaging theta on branches for each tree

Modified: pkg/R/batchHansen.R
===================================================================
--- pkg/R/batchHansen.R	2008-12-02 07:57:12 UTC (rev 62)
+++ pkg/R/batchHansen.R	2008-12-02 08:46:00 UTC (rev 63)
@@ -50,7 +50,9 @@
   if(!identical(di, NULL)) dir.create(di)
   ar = regimeVectors(ouchTrees, cladeMembersList, maxNodes)
   hansenBatch <- list(length(ouchTrees))
+  thetas <- list(length(ouchTrees))
   for (i in 1:length(ouchTrees)) {
+    fP <- NULL
     if(!identical(filePrefix, NULL)) fP <- paste(filePrefix, ".t", i, ".", sep = "")
     if(!identical(di, NULL)) fP <- paste(di, "/", fP, sep = "")
     tree <- ouchTrees[[i]]
@@ -73,47 +75,57 @@
     ## 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(...)')
     ##   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, ar$regList[[i]], regimeTitles, brown, fP, alpha, sigma, ...)
+    hb <- batchHansen(tree, dataIn, ar$regList[[i]], regimeTitles, brown, fP, alpha, sigma, ...)
+    hansenBatch[[i]] <- hb$treeData
+    thetas[[i]] <- hb$thetas
     message(paste("Tree",i,"of",length(ouchTrees),"complete", "\n-----------------------------"))
   }
     
     ## right now no summary is returned; perhaps one is needed, summarizing over trees what is summarized for each tree in batchHansen
-  outdata <- list(hansens = hansenBatch, regList = ar$regList, regMatrix = ar$regMatrix, nodeMatrix = ar$nodeMatrix, brown = brown, N = ouchTrees[[i]]@nterm, analysisDate = date(), call = match.call())
+  outdata <- list(hansens = hansenBatch, thetas = thetas, regList = ar$regList, regMatrix = ar$regMatrix, nodeMatrix = ar$nodeMatrix, brown = brown, N = ouchTrees[[i]]@nterm, analysisDate = date(), call = match.call())
   class(outdata) <- 'hansenBatch'
   return(outdata)}
 
 batchHansen <-
 # Runs hansen and brown on a tree over a batch of selective regimes
 # Arguments:
-#  "node" "ancestor" "times" "data" = the standard tree specification vectors of the OUCH-style tree
+#  "tree" = the standard OUCH-style (S4) 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 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, filePrefix, alpha, sigma, ...) {
+function(tree, data, regimesList, regimeTitles, brown, filePrefix = NULL, alpha, sigma, ...) {
   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
+  ## set up a matrix that returns lnL, K, sigmasq, theta0, and alpha for every model
+  ## thetas go into a models-by-branch matrix
   hansenOptima <- list(length(regimeTitles))
   variables <- c("loglik", "dof", "sigma.squared", "theta / alpha") # only display variables... set the selecting variables in the next two lines
   brVars <- c("loglik", "dof", "sigma.squared", "theta")
   haVars <- c("loglik", "dof", "sigma.squared", "alpha")
+  if(brown) thetaModels <- regimeTitles[1: (length(regimeTitles) - 1)]
+  else thetaModels <- regimeTitles
+  thetas <- matrix(NA, 
+                   nrow = length(thetaModels), 
+                   ncol = tree at nnodes, 
+                   dimnames = list(thetaModels, tree at nodes))
   treeData <- matrix(data = NA, nrow = length(regimeTitles), ncol = length(variables), dimnames = list(regimeTitles,variables))
   if(brown) {
     br <- brown(data, tree)
+    if(!identical(filePrefix, NULL)) save(br, file = paste(filePrefix, 'b.Rdata', sep = ""))
     treeData["brown", ] <- unlist(summary(br)[brVars])
     }
   for (i in seq(regimesList)) {
     if(any(is.na(regimesList[[i]]))) {
       message(paste("skipping regime", i))
       treeData[i, ] <- rep(NA, dim(treeData)[[2]])
-      hansenOptima[[i]] <- NA
       }
     else {
       message(paste("Running regime",i))
       ## 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 = data, tree = tree, regimes = regimesList[[i]], alpha = alpha, sigma = sigma, ...)
-      if(!identical(filePrefix, NULL)) save(ha, paste(filePrefix, 'r', i, '.Rdata', sep = ""))
       treeData[i, ] <- unlist(summary(ha)[haVars])
-      hansenOptima[[i]] <- summary(ha)$optima[[1]]
+      thetas[i, ] <- ha at theta$data[ha at regimes[[1]]]
+      if(!identical(filePrefix, NULL)) save(ha, file = paste(filePrefix, 'r', i, '.Rdata', sep = ""))
       }
   }
-  return(treeData) }
\ No newline at end of file
+  outdata <- list(treeData = treeData, thetas = thetas)
+  return(outdata) }
\ No newline at end of file