[Adephylo-commits] r48 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 27 15:52:40 CET 2008
Author: jombart
Date: 2008-11-27 15:52:39 +0100 (Thu, 27 Nov 2008)
New Revision: 48
Added:
pkg/R/orthobasis.R
pkg/man/orthobasis.Rd
Modified:
pkg/DESCRIPTION
pkg/R/proximities.R
Log:
orthobasis function and doc. Needs somes fixes yet.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2008-11-27 12:11:39 UTC (rev 47)
+++ pkg/DESCRIPTION 2008-11-27 14:52:39 UTC (rev 48)
@@ -9,4 +9,4 @@
Description: Multivariate tools to analyze comparative data, i.e. a phylogeny and some traits measured for each taxa.
License: GPL (>=2)
LazyLoad: yes
-Collate: utils.R partition.R s.phylo4d.R distances.R proximities.R
\ No newline at end of file
+Collate: utils.R partition.R s.phylo4d.R distances.R proximities.R orthobasis.R
\ No newline at end of file
Added: pkg/R/orthobasis.R
===================================================================
--- pkg/R/orthobasis.R (rev 0)
+++ pkg/R/orthobasis.R 2008-11-27 14:52:39 UTC (rev 48)
@@ -0,0 +1,63 @@
+###################
+# orthobasis.phylo
+###################
+orthobasis.phylo <- function(x=NULL, prox=NULL,
+ method=c("brlength","nNodes","Abouheif","sumDD"), a=1){
+ if(!require(phylobase)) stop("phylobase package is not installed")
+ if(!require(ade4)) stop("ade4 package is not installed")
+
+ ## handle arguments
+ method <- match.arg(method)
+
+ if(is.null(prox)){ # have to compute prox
+ x <- as(x, "phylo4")
+ if (is.character(checkval <- check_phylo4(x))) stop(checkval)
+ W <- proxTips(x, tips="all", method=method, a=a, normalize="row", symmetric=TRUE)
+ } else { # prox is provided
+ W <- as.matrix(prox)
+ if(!is.matrix(W)) stop("W is not a matrix")
+ if(ncol(W) != nrow(W)) stop("W is not a square matrix")
+ diag(W) <- 0
+ W <- 0.5 * (t(W) + W) # re-symmetrization
+ }
+
+ n <- nrow(W)
+
+
+ ## main computation
+ W0 <- bicenter.wt(W)
+ decomp <- eigen(W0, sym=TRUE)
+ E <- decomp$vectors # Moran eigenvectors (ME)
+ ## must be re-orthogonalized
+ temp <- cbind(rep(1,n) , E)
+ temp <- qr.Q(qr(temp))
+ E <- as.data.frame(temp[,-1])*sqrt(n)
+ row.names(E) <- rownames(W)
+ names(E) <- paste("ME", 1:ncol(E), sep=".")
+
+ ## retrieve Moran' I for each ME
+ f1 <- function(vec){
+ res <- mean( vec * (W0 %*% vec) )
+ return(res)
+ }
+
+ Ival <- apply(E, 2, f1)
+
+ ## build output
+ res <- E
+ attr(res,"values") <- Ival
+ attr(res,"weights") <- rep(1/n,n)
+ attr(res,"call") <- match.call()
+ attr(res,"class") <- c("orthobasis","data.frame")
+
+ return(res)
+} # end orthobasis.phylo
+
+
+
+
+
+###########
+# me.phylo
+###########
+me.phylo <- orthobasis.phylo
Modified: pkg/R/proximities.R
===================================================================
--- pkg/R/proximities.R 2008-11-27 12:11:39 UTC (rev 47)
+++ pkg/R/proximities.R 2008-11-27 14:52:39 UTC (rev 48)
@@ -14,7 +14,6 @@
N <- nTips(x)
if(tips[1]=="all") { tips <- 1:N }
tips <- getnodes(x, tips)
- tips.names <- names(tips)
## some checks
if (is.character(checkval <- check_phylo4(x))) stop(checkval)
Added: pkg/man/orthobasis.Rd
===================================================================
--- pkg/man/orthobasis.Rd (rev 0)
+++ pkg/man/orthobasis.Rd 2008-11-27 14:52:39 UTC (rev 48)
@@ -0,0 +1,66 @@
+\name{Moran's eigenvectors}
+\alias{orthobasis.phylo}
+\alias{me.phylo}
+\title{Computes Moran's eigenvectors from a tree or a phylogenetic
+ proximity matrix}
+\description{
+ The function \code{orthobasis.phylo} (also nicknamed \code{me.phylo})
+ computes Moran's eigenvectors (ME) associated to a tree. If the tree has
+ 'n' tips, (n-1) vectors will be produced. These vectors form an
+ orthonormal basis: they are centred to mean zero, have unit variance,
+ and are uncorrelated. Each vector models a different pattern of
+ phylogenetic autocorrelation. The first vectors are those with maximum
+ positive autocorrelation, while the last vectors are those with
+ maximum negative autocorrelation. ME can be used, for instance, as
+ regressors to remove phylogenetic autocorrelation from data (see
+ references).\cr
+
+ ME can be obtained from a tree, specifying the
+ phylogenetic proximity to be used. Alternatively, they can be obtained
+ directly from a matrix of phylogenetic proximities as constructed by
+ \code{\link{proxTips}}.
+}
+\usage{
+me.phylo(x=NULL, prox=NULL, method=c("brlength","nNodes","Abouheif","sumDD"), a=1)
+orthobasis.phylo(x=NULL, prox=NULL, method=c("brlength","nNodes","Abouheif","sumDD"), a=1)
+}
+\arguments{
+ \item{x}{A tree of class \code{\link[pkg:ape]{phylo}},
+ \linkS4class{phylo4} or \linkS4class{phylo4d}.}
+ \item{prox}{a matrix of phylogenetic proximities as returned by \code{\link{proxTips}}.}
+ \item{method}{a character string (full or abbreviated without
+ ambiguity) specifying the method used to compute proximities;
+ possible values are:\cr
+ - \code{brlength}: (inversed sum of) branch length \cr
+ - \code{nNodes}: (inversed) number of nodes on the path between the nodes \cr
+ - \code{Abouheif}: Abouheif's proximity (see details) \cr
+ - \code{sumDD}: (inversed) sum of direct descendants of all nodes on the path
+ (see details in \code{\link{proxTips}}).}
+ \item{a}{the exponent used to compute the proximity (see ?\code{\link{proxTips}}).}
+}
+\value{
+ An object of class \code{orthobasis}. This is a data.frame with
+ Moran's eigenvectors in column, with special attributes:\cr
+ - attr(...,"values"): Moran's index for each vector
+ - attr(...,"weights"): weights of tips; current implementation uses
+ only uniform weights
+}
+\author{ Thibaut Jombart \email{jombart at biomserv.univ-lyon1.fr} }
+\seealso{\code{\link{proxTips}} which computes phylogenetic proximities
+ between tips.
+}
+\examples{
+if(require(ape) & require(ade4)){
+## make a tree
+x <- rtree(50)
+
+## compute Moran's eigenvectors
+ME <- me.phylo(x)
+ME
+
+## plot the 5 first vectors
+obj <- phylo4d(x, ME[,1:5])
+s.phylo4d(obj, cent=FALSE, scale=FALSE)
+}
+}
+\keyword{manip}
More information about the Adephylo-commits
mailing list