[Phylobase-commits] r420 - in pkg: . R inst/doc man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Jan 1 03:39:49 CET 2009


Author: bbolker
Date: 2009-01-01 03:39:49 +0100 (Thu, 01 Jan 2009)
New Revision: 420

Added:
   pkg/R/class-phylomats.R
   pkg/man/phylomat-class.Rd
Modified:
   pkg/DESCRIPTION
   pkg/inst/doc/phylobase.Rnw
   pkg/man/as-methods.Rd
Log:
   added var-cov matrix class; extended vignette
  


Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2009-01-01 02:37:20 UTC (rev 419)
+++ pkg/DESCRIPTION	2009-01-01 02:39:49 UTC (rev 420)
@@ -9,6 +9,6 @@
 Maintainer:  Ben Bolker <bolker at ufl.edu>
 Description: Provides a base S4 class for comparative methods, incorporating one or more trees and trait data
 License: GPL
-Collate: phylo4.R checkdata.R class-multiphylo4.R class-oldclasses.R class-phylo4.R class-phylo4d.R  methods-multiphylo4.R methods-oldclasses.R methods-phylo4.R methods-phylo4d.R setAs-Methods.R pdata.R subset.R prune.R treePlot.R identify.R treestruc.R treewalk.R readNexus.R tbind.R zzz.R printphylo-deprecated.R
+Collate: phylo4.R checkdata.R class-multiphylo4.R class-oldclasses.R class-phylo4.R class-phylo4d.R  class-phylomats.R methods-multiphylo4.R methods-oldclasses.R methods-phylo4.R methods-phylo4d.R setAs-Methods.R pdata.R subset.R prune.R treePlot.R identify.R treestruc.R treewalk.R readNexus.R tbind.R zzz.R printphylo-deprecated.R
 Encoding: UTF-8
 URL: http://phylobase.R-forge.R-project.org

Added: pkg/R/class-phylomats.R
===================================================================
--- pkg/R/class-phylomats.R	                        (rev 0)
+++ pkg/R/class-phylomats.R	2009-01-01 02:39:49 UTC (rev 420)
@@ -0,0 +1,87 @@
+
+## define class for phylogenetic var-cov matrices
+setClass("phylo4vcov",
+         representation("matrix",
+                        edge.label="character",
+                        order="character"))
+
+## phylo4 -> var-cov: simply wrap ape::vcv.phylo
+##  and add other slots
+as_phylo4vcov <- function(from,...) {
+  m <- ape::vcv.phylo(as(from,"phylo"),...)
+  new("phylo4vcov",
+      m,
+      edge.label=from at edge.label,
+      order=from at order)
+}
+setAs("phylo4","phylo4vcov",
+      function(from,to) {
+        as_phylo4vcov(from)})
+
+## var-cov to phylo4
+setAs("phylo4vcov","phylo4",
+      function(from,to) {
+        matrix2tree <- function(v,reorder=TRUE) {
+          ## no polytomies allowed
+          va <- v
+          tipnames <- rownames(v)
+          ntip <- nrow(v)
+          dimnames(v) <- list(as.character(1:ntip),
+                              as.character(1:ntip))
+          diag(va) <- 0
+          edgemat <- matrix(ncol=2,nrow=0)
+          ## termlens <- diag(v)-colSums(va)
+          edgelens <- numeric(0)
+          ## maxnode <- ntip
+          curnode <- 2*ntip ## one greater than total number of nodes
+          ## can we do this in a different order?
+          while (nrow(v)>1) {
+            mva <- max(va)  ## find pair with max shared evolution
+            nextpr <- if (nrow(v)==2) c(1,2) else which(va==mva,arr.ind=TRUE)[1,]
+            ## maxnode <- maxnode+1  ## new node
+            curnode <- curnode-1
+            ## points to both of current identified nodes
+            ##   (indexed by names)
+            edgemat <- rbind(edgemat,
+                             c(curnode,as.numeric(rownames(v)[nextpr[1]])),
+                             c(curnode,as.numeric(rownames(v)[nextpr[2]])))
+            ## descending edges are amount of *unshared* evolution
+            edgelens <- c(edgelens,
+                          diag(v)[nextpr]-mva)
+            ## this clade has total evolution = shared evolution
+            diag(v)[nextpr] <- mva
+            ## assign new node name
+            rownames(v)[nextpr[1]] <- colnames(v)[nextpr[1]] <- curnode
+            ## drop rows/cols from matrix
+            v <- v[-nextpr[2],-nextpr[2],drop=FALSE]
+            va <- va[-nextpr[2],-nextpr[2],drop=FALSE]
+          }
+          ## switch order of node numbers to put root in the right place:
+          ##  much plotting code seems to assume root = node # (ntips+1)
+          ## browser()
+          reorder <- FALSE
+          if (reorder) {
+            nn <- nrow(edgemat)
+            nnode <- nn-ntip+1
+            newedge <- edgemat
+            for (i in 2:nnode) {
+              newedge[edgemat==(ntip+i)] <- nn-i+2
+            }
+            edgemat <- newedge
+          }
+          list(edgemat=edgemat,
+               edgelens=edgelens)
+        }
+        temptree <- matrix2tree(from)
+        ## browser()
+        ## add explicit root
+        rootnode <- which(tabulate(temptree$edgemat[,2])==0)
+        ## add root node to edge matrix and branch lengths
+        temptree$edgemat <- rbind(temptree$edgemat,c(NA,rootnode))
+        temptree$edgelens <- c(temptree$edgelens,NA)
+        phylo4(temptree$edgemat,edge.length=temptree$edgelens,
+               tip.label=rownames(from),
+               edge.label=from at edge.label,order=from at order)
+      })
+
+

Modified: pkg/inst/doc/phylobase.Rnw
===================================================================
--- pkg/inst/doc/phylobase.Rnw	2009-01-01 02:37:20 UTC (rev 419)
+++ pkg/inst/doc/phylobase.Rnw	2009-01-01 02:39:49 UTC (rev 420)
@@ -329,31 +329,55 @@
 tree topology and branch lengths:
 <<rtree2>>=
 set.seed(1001)
-tree <- rcoal(12)
+tree <- as(rcoal(12),"phylo4")
 @
 
+
 Next we generate the phylogenetic variance-covariance
-matrix (\code{ape::vcv.phylo}) and pick a single set
-of traits (\code{MASS:mvrnorm}).  Conveniently, the
-tip names of the original tree
-are inherited consistently by the variance-covariance
-matrix and the trait matrix:
+matrix (by coercing the tree to 
+a \code{phylo4vcov} object and pick a single set
+of traits (\code{MASS:mvrnorm}).
 <<vcvphylo>>=
-vmat <- vcv.phylo(tree,cor=TRUE)
+vmat <- as(tree,"phylo4vcov")
+vmat <- cov2cor(vmat)
 library(MASS)
 trvec <- mvrnorm(1,mu=rep(0,12),Sigma=vmat)
 @
 
-The last step (easy) is to create the \code{phylo4d}
-object and plot it:
+The last step (easy) is to 
+convert the \code{phylo4vcov} object
+back to a \code{phylo4d} object:
 <<plotvcvphylo,fig=TRUE>>=
 treed <- phylo4d(tree,tip.data=as.data.frame(trvec))
 plot(treed)
 @
 
-\subsubsection{The hard way?}
+\subsubsection{The hard way}
 
-Find root, traverse tree:
+<<keep.source=TRUE>>=
+## add node labels so we can match to data
+nodeLabels(tree) <- as.character(sort(nodeId(tree)))
+## ordering will make sure that we have ancestor value 
+## defined before descendant
+tree <- reorder(tree,"preorder")
+edgemat <- edges(tree)
+## set aside space for values
+nodevals <- numeric(nrow(edgemat))
+## label data in edge matrix order
+names(nodevals) <- labels(tree,"all")[nodeId(tree,"all")]
+## variance is proportional to edge length; drop first
+## element of edge length, which is NA
+dvals <- rnorm(nrow(edgemat)-1,sd=edgeLength(tree)[-1]^2)
+## indexing: ind[node number] gives position in edge matrix
+ind <- order(nodeId(tree,"all"))
+for (i in 2:nrow(edgemat)) {
+  ## value of ancestor node plus change
+  nodevals[i] <- nodevals[ind[edgemat[i,1]]]+dvals[i-1]
+  }
+nodevals <- data.frame(nodevals)
+treed2 <- phylo4d(tree,all.data=nodevals)
+}
+@ 
 
 
 % ========================================

Modified: pkg/man/as-methods.Rd
===================================================================
--- pkg/man/as-methods.Rd	2009-01-01 02:37:20 UTC (rev 419)
+++ pkg/man/as-methods.Rd	2009-01-01 02:39:49 UTC (rev 420)
@@ -11,6 +11,8 @@
 \alias{as,phylo4d,phylo-method}
 \alias{as,phylo4,data.frame-method}
 \alias{as,phylo4d,data.frame-method}
+\alias{as,phylo4vcov,phylo4-method}
+\alias{as,phylo4,phylo4vcov-method}
 \alias{coerce-methods}
 \alias{coerce,phylo,phylo4-method}
 \alias{coerce,phylo,phylo4d-method}
@@ -21,6 +23,8 @@
 \alias{coerce,phylo4d,phylo-method}
 \alias{coerce,phylo4,data.frame-method}
 \alias{coerce,phylo4d,data.frame-method}
+\alias{coerce,phylo4vcov,phylo4-method}
+\alias{coerce,phylo4,phylo4vcov-method}
 %\alias{coerce,phylo4d,phylo4-method}
 
 

Added: pkg/man/phylomat-class.Rd
===================================================================
--- pkg/man/phylomat-class.Rd	                        (rev 0)
+++ pkg/man/phylomat-class.Rd	2009-01-01 02:39:49 UTC (rev 420)
@@ -0,0 +1,53 @@
+\name{phylomat-class}
+\docType{class}
+\alias{phylo4vcov-class}
+\alias{as_phylo4vcov}
+\title{matrix classes for phylobase}
+\description{Classes representing phylogenies as matrices}
+\section{Objects from the Class}{
+  These are square matrices (with rows and columns corresponding
+  to tips, and internal nodes implicit) with different
+  meanings depending on the type (variance-covariance matrix,
+  distance matrix, etc.).
+}
+\section{Slots}{
+  \describe{
+    \item{\code{.Data}:}{square, numeric matrix with row and column
+      labels corresponding to the tip labels}
+    \item{\code{edge.label}:}{character vector of edge labels}
+    \item{\code{order}:}{character describing original ordering of edge matrix}
+  }
+}
+\usage{
+as_phylo4vcov(from, \dots)
+}
+\arguments{
+  \item{from}{a \code{phylo4} object}
+  \item{\dots}{optional arguments, to be passed to \code{vcov.phylo} in
+    \code{ape} (the main useful option is \code{cor}, which can be set
+    to \code{TRUE} to compute a correlation rather than a
+    variance-covariance matrix)}
+}
+\author{Ben Bolker}
+\examples{
+  example("read.tree")
+  o2 <- as(tree.owls,"phylo4")
+  ov <- as(o2,"phylo4vcov")
+  o3 <- as(ov,"phylo4")
+  ## these are not completely identical, but are
+  ## topologically identical ...
+  
+  ## edge matrices are in a different order:
+  ## cf. o2 at edge and o3 at edge
+  ## BUT the edge matrices are otherwise identical
+  identical(o2 at edge[order(o2 at edge[,2]),],
+            o3 at edge[order(o3 at edge[,2]),])
+  
+  ## There is left/right ambiguity here in the tree orders:
+  ## in o2 the 5->6->7->1 lineage
+  ## (terminating in Strix aluco)
+  ## is first, in o3 the 5->6->3 lineage 
+  ## (terminating in Athene noctua) is first.
+
+}
+\keyword{classes}



More information about the Phylobase-commits mailing list