[Mattice-commits] r167 - pkg/misc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 19 16:42:41 CET 2009
Author: andrew_hipp
Date: 2009-02-19 16:42:41 +0100 (Thu, 19 Feb 2009)
New Revision: 167
Added:
pkg/misc/multiModels-dataFormat.csv
Modified:
pkg/misc/multiModels.R
Log:
updating multiModels... cleaning up parameter slotting
Added: pkg/misc/multiModels-dataFormat.csv
===================================================================
--- pkg/misc/multiModels-dataFormat.csv (rev 0)
+++ pkg/misc/multiModels-dataFormat.csv 2009-02-19 15:42:41 UTC (rev 167)
@@ -0,0 +1,12 @@
+,loglik,dof,sigma.squared,alpha,theta,optimum,optimum.uptree,optimum.downtree
+whole.ou2,x,x,x,x,,,x,x
+whole.ou1,x,x,x,x,,x,,
+whole.brown,x,x,x,,x,,,
+part.ou.uptree,x,x,x,x,,x,,
+part.ou.downtree,x,x,x,x,,x,,
+part.brown.uptree,x,x,x,,x,,,
+part.brown.downtree,x,x,x,,x,,,
+
+"paramHeader <- c('loglik', 'dof', 'sigma.squared', 'alpha', 'theta', 'optimum', 'optimum.uptree', 'optimum.downtree')",,,,,,,,
+"paramSets <- list(all = c('loglik', 'dof', 'sigma.squared'), brown = c('theta'), ou1 = c('alpha', 'optimum'), ou2 = c('alpha', 'optimum.uptree', 'optimum.downtree'))",,,,,,,,
+"modelAll = c('whole.ou2', 'whole.ou1', 'whole.brown', 'part.ou.uptree', 'part.ou.downtree', 'part.brown.uptree', 'part.brown.downtree')",,,,,,,,
Modified: pkg/misc/multiModels.R
===================================================================
--- pkg/misc/multiModels.R 2009-02-11 16:14:03 UTC (rev 166)
+++ pkg/misc/multiModels.R 2009-02-19 15:42:41 UTC (rev 167)
@@ -1,28 +1,40 @@
multiModels <- function(phy, dat, node, models = c('whole.brown', 'whole.ou1', 'whole.ou2','part.brown', 'part.ou'), ic = "BIC") {
# test the support for alternative models on simple and partitioned trees
-# currently only works on one tree; fix so it runs on a set of trees
- parameters <- c('loglik', 'dof', 'alpha', 'sigma.squared', 'theta', 'optima.nodeAncestor', 'optima.nodeDescendent')
+# currently only works on one tree; fix so it runs on a set of trees, conditioned on those trees that have the node of interest;
+# return percent of trees possessing that node as an additional value
+ paramHeader <- c('AICc', 'BIC', 'AICc.weight', 'BIC.weight', 'loglik', 'dof', 'sigma.squared', 'alpha', 'theta', 'optimum', 'optimum.uptree', 'optimum.downtree')
+ paramsAll <- c('loglik', 'dof', 'sigma.squared')
+ paramSets <- list(brown = c(paramsAll, 'theta'),
+ ou1 = c(paramsAll, 'alpha', 'optimum'),
+ ou2 = c(paramsAll, 'alpha', 'optimum.uptree', 'optimum.downtree')
+ )
+ modelsAll = c('whole.ou2', 'whole.ou1', 'whole.brown', 'part.ou.uptree', 'part.ou.downtree', 'part.ou.summed', 'part.brown.uptree', 'part.brown.downtree', 'part.brown.summed')
+ # parameters <- c('loglik', 'dof', 'alpha', 'sigma.squared', 'theta', 'optimum', 'optima.nodeAncestor', 'optima.nodeDescendent')
pSum <- c('loglik', 'dof') # parameters to sum for evaluating partitioned trees
+
aboveNodeTaxa <- phy$tip.label[-which(phy$tip.label %in% node)]
wholeTree <- ape2ouch(phy)
+
upTree <- ape2ouch(drop.tip(phy, node))
downTree <- ape2ouch(drop.tip(phy, aboveNodeTaxa))
partialTree <- list(upTree = upTree, downTree = downTree)
- outMatrix <- matrix(NA, nrow = length(models), ncol = length(parameters), dimnames = list(models, parameters))
+
+ outMatrix <- matrix(NA, nrow = length(modelsAll), ncol = length(paramHeader), dimnames = list(modelsAll, paramHeader))
modelsSplit <- strsplit(models, ".", fixed = T)
+
for(i in modelsSplit) {
model <= paste(i[1], ".", i[2], sep = "")
# as written, the following won't correctly chunk parameters into the matrix... mismatch between partitioned- and whole-tree tests -- NOW MAY WORK... DOUBLE CHECK RELATIVE TO NAMES SPIT OUT IN WHOLE MODEL, AND CHEKC PARTIALMODEL NAMES
- outMatrix[model, ] <- ifelse(i[1] == 'whole',
- wholeModel(wholeTree, dat, i[2], node, parameters)$params,
- partialModel(partialTree, dat, i[2], c('uptree', 'downtree'), parameters)$params
- )
+ if (i[1] == 'whole') outMatrix[model, ] <- wholeModel(wholeTree, dat, i[2], node, paramSets[[i[2]]], paramHeader)$params
+ if (i[1] == 'part') outMatrix[c(paste(model, '.uptree', sep = ""), paste(model, '.downtree', sep = ""), paste(model, '.summed', sep = "")), ]
+ <- partialModel(partialTree, dat, i[2], c('uptree', 'downtree'), paramSets[[i[2]]], pSum, paramHeader)$params
+
}
MAKE A WEIGHTS COLUMN using informationCriterion
return(outMatrix)
}
-wholeModel <- function(phy, dat, model, node, parameterVector) {
+wholeModel <- function(phy, dat, model, node, parameterVector, paramHeader) {
if(model == "brown") analysis <- brown(dat, phy)
if(model == "ou1") analysis <- hansen(dat, phy,
regimes = structure(rep(1, phy at nnodes),
@@ -35,14 +47,18 @@
descNum <- as.character(unique(regime))[unique(regime) != tree at root]
analysis <- hansen(dat, phy, regime, alpha = 1, sigma = 1)
} # close if
+ ## CHANGE SO THAT PARAMS IS FORMATTED USING paramHeader
params <- unlist(summary(analysis)[parameterVector])[parameterVector]
- params['optima.nodeAncestor'] <- summary(analysis)$optima[[1]][ancNum]
- params['optima.nodeDescendent'] <- summary(analysis)$optima[[1]][descNum]
+ if(model == "ou2") {
+ params['optima.nodeAncestor'] <- summary(analysis)$optima[[1]][ancNum]
+ params['optima.nodeDescendent'] <- summary(analysis)$optima[[1]][descNum]
+ }
+ if(model == "ou1") params['optimum'] <- summary(analysis)$optima[[1]]
out <- list(analysis = analysis, params = params)
return(out)
}
-partialModel <- function(phyList, dat, model, treeNames, parameterVector = NULL, pSum = NULL) {
+partialModel <- function(phyList, dat, model, treeNames, parameterVector = NULL, pSum = NULL, paramHeader) {
# sums a subset of parameters (indicated by pSum) and leaves the others separate
allParams <- pSum
for(i in parameterVector) {
@@ -58,7 +74,9 @@
levels = 1, class = 'factor'),
sigma = 1, alpha = 1)
} # close else
+ # change params so that it returns a matrix that fits -- the two partial tree models + the sum
params <- lapply(analysis, function(x) {unlist(summary(x)[pV])[pV]}, pV = allParams)
+ params['optimum'] <- summary
rawMat <- matrix(unlist(params), nrow = length(params), ncol = length(allParams), byrow = T, dimnames = list(treeNames, allParams))
params <- colSums(rawMat, na.rm = T)
out <- list(analysis = analysis, rawMat = rawMat, params = params)
More information about the Mattice-commits
mailing list