[Mattice-commits] r255 - in pkg: R inst/doc

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

Author: andrew_hipp
Date: 2012-05-03 05:17:48 +0200 (Thu, 03 May 2012)
New Revision: 255

Modified:
   pkg/R/batchHansen.R
   pkg/inst/doc/maticce.Rnw
Log:
summary.batchHansen is misbehaving when the Brownian motion model is included; this option is now blocked and the vignette rewritten so there is no longer an error message.

Modified: pkg/R/batchHansen.R
===================================================================
--- pkg/R/batchHansen.R	2012-05-03 02:18:31 UTC (rev 254)
+++ pkg/R/batchHansen.R	2012-05-03 03:17:48 UTC (rev 255)
@@ -95,6 +95,7 @@
 #  "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 = NULL, sqrt.alpha, sigma, ...) {
+  if(brown) stop("Including the Brownian motion model has been discontinued in batchHansen")
   n <- tree at nterm
   ## set up a matrix that returns lnL, K, sigmasq, theta0, and sqrt.alpha for every model
   ## thetas go into a models-by-branch matrix

Modified: pkg/inst/doc/maticce.Rnw
===================================================================
--- pkg/inst/doc/maticce.Rnw	2012-05-03 02:18:31 UTC (rev 254)
+++ pkg/inst/doc/maticce.Rnw	2012-05-03 03:17:48 UTC (rev 255)
@@ -39,7 +39,7 @@
 In case you aren't familiar with \code{R}, the following commands will get you started.
 
 <<startAnalysis, fig=FALSE>>=
-library(maticce) # load maticce and required packages
+# library(maticce) # load maticce and required packages
 data(carex) # load dataset
 attach(carex) # attach dataset to search path
 # convert the Bayes consensus to an ouchtree object...
@@ -123,16 +123,16 @@
 <<runBatch, fig=FALSE, echo=TRUE>>=
 # First, analyze with maxNodes set to 2
 ha.4.2 <- runBatchHansen(ovales.tree, ovales.data, 
-          ovales.nodes[1:4], maxNodes = 2, brown = T)
+          ovales.nodes[1:4], maxNodes = 2)
 print(summary(ha.4.2))
 # Then, analyze with maxNodes set to 4
 ha.4.4 <- runBatchHansen(ovales.tree, ovales.data, 
-          ovales.nodes[1:4], maxNodes = 4, brown = T)
+          ovales.nodes[1:4], maxNodes = 4)
 print(summary(ha.4.4))
 # Then, assess the effects of phylogenetic uncertainty by 
 #   analyzing over a set of trees
 ha.4.2.multi <- runBatchHansen(trees, ovales.data, 
-                ovales.nodes[1:4], maxNodes = 2, brown = T)
+                ovales.nodes[1:4], maxNodes = 2)
 print(summary(ha.4.2.multi))
 @
 
@@ -156,13 +156,15 @@
 Then we can find the likelihood and parameter estimates of these models on a given tree:
 
 <<haLnl, fig = FALSE, echo = TRUE>>=
-ha.4.2[['hansens']][[1]][c(7, 11, 'brown'), ]
+# ha.4.2[['hansens']][[1]][c(7, 11, 'brown'), ]
+ha.4.2[['hansens']][[1]][c(7, 11), ]
 @
 
 or the information criterion weights:
 
 <<haWeights, fig=FALSE, echo=TRUE>>=
-summary(ha.4.2)[['modelsMatrix']][[1]][c(7, 11, 'brown'), ]
+# summary(ha.4.2)[['modelsMatrix']][[1]][c(7, 11, 'brown'), ]
+summary(ha.4.2)[['modelsMatrix']][[1]][c(7, 11), ]
 @
 
 Considering just these models, model 7 is not overwhelmingly supported (BIC weight = 0.530, AICc weight = 0.369), but it is much more strongly supported than the Brownian motion model or the OU model with no change. This points to the utility of model-averaging as a means of localizing character transitions on a phylogenetic tree. Moreover, the fact that a character transition is strongly supported only for node 2 tells us little about whether each node, analyzed on its own, would support a character transition model over a no-transition model. In fact, in the sample data, nodes 1, 2, 3, 4, and 7 all support a transition over no-transition model. You can investigate this node-by-node using the \code{multiModel} function.