[Mattice-commits] r173 - in pkg: . R man misc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 24 07:19:14 CET 2009
Author: andrew_hipp
Date: 2009-02-24 07:19:14 +0100 (Tue, 24 Feb 2009)
New Revision: 173
Added:
pkg/R/multiModel.R
pkg/man/multiModel.Rd
Removed:
pkg/misc/multiModel.R
pkg/misc/multiModels.R
Modified:
pkg/DESCRIPTION
pkg/man/carex.Rd
Log:
fixing documentation; moving multiModel.R into the live folder
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2009-02-24 05:02:09 UTC (rev 172)
+++ pkg/DESCRIPTION 2009-02-24 06:19:14 UTC (rev 173)
@@ -3,13 +3,12 @@
Title: Mapping Transitions in Continuous Character Evolution
Version: 0.4
Date: 2008-12-18
-Depends: ouch
-Suggests: ape
+Depends: ouch, ape
Enhances: ouch
-Author: Andrew Hipp <ahipp at mortonarb.org>, with contributions from Marcial Escudero.
+Author: Andrew Hipp <ahipp at mortonarb.org>, with contributions from Marcial Escudero
Maintainer: Andrew Hipp <ahipp at mortonarb.org>
Description: Tools for an information-theoretic approach to estimating
- the probability of continuous character transitions on phylogenetic trees.
+ the probability of continuous character transitions on phylogenetic trees
License: GPL
URL: http://mattice.R-forge.R-project.org
Revision: 1
\ No newline at end of file
Added: pkg/R/multiModel.R
===================================================================
--- pkg/R/multiModel.R (rev 0)
+++ pkg/R/multiModel.R 2009-02-24 06:19:14 UTC (rev 173)
@@ -0,0 +1,95 @@
+multiModel <- function(phy, dat, node, models = c('whole.brown', 'whole.ou1', 'whole.ou2','part.brown', 'part.ou')) {
+# 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, conditioned on those trees that have the node of interest;
+# return percent of trees possessing that node as an additional value
+ paramHeader <- c('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')
+ 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(modelsAll), ncol = length(paramHeader), dimnames = list(modelsAll, paramHeader))
+ compareK <- rep(NA, length(models)); names(compareK) <- models
+ comparelnL <- rep(NA, length(models)); names(comparelnL) <- models
+ modelsSplit <- strsplit(models, ".", fixed = T)
+
+ for(i in modelsSplit) {
+ model <- paste(i[1], ".", i[2], sep = "")
+ if (i[1] == 'whole') {
+ outMatrix[model, ] <- wholeModel(wholeTree, dat, i[2], node, paramSets[[i[2]]], paramHeader)$params
+ comparelnL[model] <- outMatrix[model, 'loglik']
+ compareK[model] <- outMatrix[model, 'dof']
+ }
+ 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
+ compareK[model] <- outMatrix[paste(model, '.summed', sep = ""), 'dof']
+ comparelnL[model] <- outMatrix[paste(model, '.summed', sep = ""), 'loglik']
+ }
+ }
+ IC <- informationCriterion(lnL = comparelnL, K = compareK, n = length(phy$tip.label), names = models)
+ outdata <- list(IC = IC, modelMatrix <- outMatrix)
+ return(outdata)
+ }
+
+wholeModel <- function(phy, dat, model, node, parameterVector, paramHeader) {
+ dat <- dat[phy at nodelabels]; names(dat) <- phy at nodes
+ if(model == "brown") analysis <- brown(dat, phy)
+ if(model == "ou1") analysis <- hansen(dat,
+ phy,
+ regimes = structure(rep(phy at root, phy at nnodes), names = phy at nodes, levels = 1, class = 'factor'),
+ alpha = 1,
+ sigma = 1
+ )
+ if(model == "ou2") {
+ regime <- paintBranches(list(node), phy)
+ uptreeNum <- as.character(phy at root)
+ downtreeNum <- as.character(unique(regime))[unique(regime) != phy at root]
+ analysis <- hansen(dat, phy, regime, alpha = 1, sigma = 1)
+ }
+ params <- unlist(summary(analysis)[parameterVector])[paramHeader]
+ names(params) <- paramHeader
+ if(model == "ou2") {
+ params['optimum.uptree'] <- summary(analysis)$optima[[1]][uptreeNum]
+ params['optimum.downtree'] <- summary(analysis)$optima[[1]][downtreeNum]
+ }
+ 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, paramHeader) {
+# sums a subset of parameters (indicated by pSum) and leaves the others separate
+ allParams <- pSum
+ analysis <- vector('list', length(phyList))
+ if(model == "brown") {
+ for (i in seq(length(phyList))) analysis[[i]] <- brown(data = structure(dat[phyList[[i]]@nodelabels], .Names = phyList[[i]]@nodes), tree = phyList[[i]])
+ }
+ else {
+ for (i in seq(length(phyList))) analysis[[i]] <- hansen(data = structure(dat[phyList[[i]]@nodelabels], .Names = phyList[[i]]@nodes), tree = phyList[[i]],
+ regimes = structure(rep(1, phyList[[i]]@nnodes),
+ names = phyList[[i]]@nodes,
+ levels = 1, class = 'factor'),
+ sigma = 1, alpha = 1)
+ }
+ names(analysis) <- treeNames
+ params <- matrix(NA, nrow = length(c(treeNames, 'summed')), ncol = length(paramHeader), dimnames = list(c(treeNames, 'summed'), paramHeader))
+ for (i in treeNames) {
+ temp <- summary(analysis[[i]])
+ params[i, ] <- unlist(temp[paramHeader])[paramHeader]
+ if(model != "brown") params[i, 'optimum'] <- temp$optima[[1]]
+ }
+ params['summed', pSum] <- colSums(params[treeNames, pSum])
+ out <- list(analysis = analysis, params = params)
+ return(out)
+ }
+
\ No newline at end of file
Modified: pkg/man/carex.Rd
===================================================================
--- pkg/man/carex.Rd 2009-02-24 05:02:09 UTC (rev 172)
+++ pkg/man/carex.Rd 2009-02-24 06:19:14 UTC (rev 173)
@@ -28,7 +28,7 @@
\source{
Hipp, A.L. (2007)
Non-uniform processes of chromosome evolution in sedges (Carex: Cyperaceae).
- Evolution 61:2175-2194.
+ \emph{Evolution} 61:2175-2194.
}
\references{
Hipp, A.L., P.E. Rothrock, A.A. Reznicek, and J.A. Weber (2006)
@@ -47,6 +47,7 @@
ovales.tree <- ape2ouch(ovales.tree) # tree comes in in phylo format, but we need an ouchtree object
trial <- runBatchHansen(ovales.tree, ovales.data, ovales.nodes[1:3], maxNodes = 2) # for expedience, only tests for changes at up to 2 of the first 4 nodes
summary(trial) # summarizes results using summary.hansenBatch and displays the results using print.hansenSummary
+ multiModel(carex$ovales.tree, ovales.data, ovales.nodes[[2]]) # compares five different models of character change at node 2
trialSim <- ouSim(trial, ovales.tree) # simulates the evolution of the chromosome number under the model-averaged values
plot(trialSim) # plots the character simulation, with all branches black
plot(trialSim, colors = paintBranches(ovales.nodes[1:4], ovales.tree)) # plots the character simulation, with branch colors changing at all 8 nodes
Added: pkg/man/multiModel.Rd
===================================================================
--- pkg/man/multiModel.Rd (rev 0)
+++ pkg/man/multiModel.Rd 2009-02-24 06:19:14 UTC (rev 173)
@@ -0,0 +1,42 @@
+\name{multiModel}
+\alias{multiModel}
+\title{Test multiple models of character transition at a single node}
+\description{
+ Compares up to five different models of character evolution: two models of no transition (Brownian motion and a single-
+ optimum Ornstein-Uhlenbeck model); one model allowing transition at a single node in a whole-tree context; and two models
+ that split the tree into two subtrees and treat character evolution separately in the two trees (the 'censored' approach -- see
+ discussion below).
+}
+\usage{
+multiModel(phy, dat, node, models = c("whole.brown", "whole.ou1", "whole.ou2", "part.brown", "part.ou"))
+}
+\arguments{
+ \item{phy}{ An \pkg{ape}-style tree }
+ \item{dat}{ A vector of data, with names corresponding to tips in \code{phy} }
+ \item{node}{ A single node, defined by the tips that are descendent from it }
+ \item{models}{ The vector of models, where \dQuote{whole} indicates a model in which the tree is treated as a single entity, and
+ \dQuote{part} indicates a model in which the tree is subdivided and treated as two entities }
+}
+\details{
+ The 'censored' approach was presented in O'Meara et al. 2006
+}
+\value{
+ A list with two items:
+ \item{IC}{A matrix of information criterion statistics, generated by \code{informationCriterion}}
+ \item{modelMatrix}{A matrix of model parameters for each of the whole tree, partial tree, and summed-partial-tree models}
+}
+\references{
+ Hipp, A.L. (2007)
+ Non-uniform processes of chromosome evolution in sedges (Carex: Cyperaceae).
+ \emph{Evolution} 61:2175--2194.
+
+ O'Meara, B. C., C. Ané, M. J. Sanderson, and P. C. Wainwright (2006)
+ Testing for different rates of continuous trait evolution using likelihood.
+ \emph{Evolution} 60:922--933.
+}
+\author{Andrew L. Hipp \email{ahipp at mortonarb.org}}
+\note{
+ In \pkg{ouch} version 2.x, the \code{ouchtree} function rescales trees to depth = 1. As a consequence, the alpha and sigma parameters
+ are not directly comparable between subtrees, as they were in previous versions of \pkg{ouch}.
+}
+\seealso{ \code{\link{carex}} for an example}
\ No newline at end of file
Deleted: pkg/misc/multiModel.R
===================================================================
--- pkg/misc/multiModel.R 2009-02-24 05:02:09 UTC (rev 172)
+++ pkg/misc/multiModel.R 2009-02-24 06:19:14 UTC (rev 173)
@@ -1,85 +0,0 @@
-multiModel <- 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, 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')
- 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(modelsAll), ncol = length(paramHeader), dimnames = list(modelsAll, paramHeader))
- modelsSplit <- strsplit(models, ".", fixed = T)
-
- for(i in modelsSplit) {
- model <- paste(i[1], ".", i[2], sep = "")
- 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, paramHeader) {
- if(model == "brown") analysis <- brown(dat, phy)
- if(model == "ou1") analysis <- hansen(dat,
- phy,
- regimes = structure(rep(phy at root, phy at nnodes), names = phy at nodes, levels = 1, class = 'factor'),
- alpha = 1,
- sigma = 1
- )
- if(model == "ou2") {
- regime <- paintBranches(list(node), phy)
- uptreeNum <- as.character(phy at root)
- downtreeNum <- as.character(unique(regime))[unique(regime) != phy at root]
- analysis <- hansen(dat, phy, regime, alpha = 1, sigma = 1)
- }
- params <- unlist(summary(analysis)[parameterVector])[paramHeader]
- names(params) <- paramHeader
- if(model == "ou2") {
- params['optimum.uptree'] <- summary(analysis)$optima[[1]][uptreeNum]
- params['optimum.downtree'] <- summary(analysis)$optima[[1]][downtreeNum]
- }
- 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, paramHeader) {
-# sums a subset of parameters (indicated by pSum) and leaves the others separate
- allParams <- pSum
-# for(i in parameterVector) {
-# if(i %in% pSum) next
-# else allParams <- c(allParams, paste(treeNames, i, sep = ""))
-# }
- if(model == "brown") analysis <- lapply(phyList, brown, data = dat)
- else {
- analysis <- vector('list', length(phyList))
- for (i in seq(length(phyList))) analysis[[i]] <- hansen(data = dat, tree = phyList[[i]],
- regimes = structure(rep(1, phyList[[i]]@nnodes),
- names = phyList[[i]]@nodes,
- levels = 1, class = 'factor'),
- sigma = 1, alpha = 1)
- } # close else
- names(analysis) <- treeNames
- params <- matrix(NA, nrow = length(c(treeNames, 'summed')), ncol = length(paramHeader), dimnames = list(c(treeNames, 'summed'), paramHeader))
- for (i in treeNames) {
- temp <- summary(analysis)
- params[i, ] <- unlist(temp[paramHeader])[paramHeader]
- if(model != "brown") params[i, 'optimum'] <- summary(temp)$optima[[1]]
- }
- params['summed', pSum] <- colSums(params[treeNames, pSum])
- out <- list(analysis = analysis, params = params)
- return(out)
- }
-
\ No newline at end of file
Deleted: pkg/misc/multiModels.R
===================================================================
--- pkg/misc/multiModels.R 2009-02-24 05:02:09 UTC (rev 172)
+++ pkg/misc/multiModels.R 2009-02-24 06:19:14 UTC (rev 173)
@@ -1,85 +0,0 @@
-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, 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(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
- 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, paramHeader) {
- if(model == "brown") analysis <- brown(dat, phy)
- if(model == "ou1") analysis <- hansen(dat, phy,
- regimes = structure(rep(1, phy at nnodes),
- names = phy at nodes, levels = 1,
- class = 'factor'),
- alpha = 1, sigma = 1)
- if(model = "ou2") {
- regime <- paintBranches(list(node), phy)
- ancNum <- as.character(tree at root)
- 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]
- 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, paramHeader) {
-# sums a subset of parameters (indicated by pSum) and leaves the others separate
- allParams <- pSum
- for(i in parameterVector) {
- if(i %in% pSum) next
- else allParams <- c(allParams, paste(treeNames, i, sep = ""))
- }
- if(model == "brown") analysis <- lapply(phyList, brown, data = dat),
- else {
- analysis <- list(length(phyList))
- for (i in seq(length(phyList))) analysis[[i]] <- hansen(data = dat, tree = phyList[[i]],
- regimes = structure(rep(1, phyList[[i]]@nnodes),
- names = phyList[[i]]@nodes,
- 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)
- return(out)
- }
-
\ No newline at end of file
More information about the Mattice-commits
mailing list