[Adephylo-commits] r93 - in pkg: R man

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

Author: jombart
Date: 2008-12-13 14:49:45 +0100 (Sat, 13 Dec 2008)
New Revision: 93

Added:
   pkg/R/moran.R
   pkg/man/moranIdx.Rd
Log:
Added a Moran's index function and its doc. Tests are ok.


Added: pkg/R/moran.R
===================================================================
--- pkg/R/moran.R	                        (rev 0)
+++ pkg/R/moran.R	2008-12-13 13:49:45 UTC (rev 93)
@@ -0,0 +1,31 @@
+############
+# moran.idx
+############
+moran.idx <- function(x, prox){
+
+
+    ## handle arguments
+    if(any(is.na(x))) stop("NA entries in x")
+    if(!is.numeric(x)) stop("x is not numeric")
+
+    W <- as.matrix(prox)
+    if(!is.matrix(W)) stop("prox is not a matrix")
+    if(ncol(W) != nrow(W)) stop("prox is not a square matrix")
+    if(any(is.na(W))) stop("NA entries in prox")
+    diag(W) <- 0
+
+    n <- nrow(W)
+
+
+    ## main computations
+    sumW <- sum(W)
+    num <- n * sum(x * (W %*% x) )
+    denom <- sumW * sum(x*x)
+
+    if(denom < 1e-14) stop("denominator equals zero")
+
+    res <- num/denom
+
+    return(res)
+
+} # end moran.idx

Added: pkg/man/moranIdx.Rd
===================================================================
--- pkg/man/moranIdx.Rd	                        (rev 0)
+++ pkg/man/moranIdx.Rd	2008-12-13 13:49:45 UTC (rev 93)
@@ -0,0 +1,44 @@
+\name{Moran's index}
+\alias{moran.idx}
+\title{Computes Moran's index for a variable}
+\description{
+  This simple function computes Moran's index of autocorrelation given a
+  variable and a matrix of proximities among observations.
+}
+\usage{
+moran.idx(x, prox)
+}
+\arguments{
+  \item{x}{a numeric vector whose autocorrelation is computed.}
+  \item{prox}{a matrix of proximities between observations, as computed
+   by the \code{\link{proxTips}}. Off-diagonal terms must be positive or
+   null, while diagonal terms must all equal zero.}
+}
+\value{
+  The numeric value of Moran's index.
+}
+\author{ Thibaut Jombart \email{jombart at biomserv.univ-lyon1.fr} }
+\seealso{\code{\link{proxTips}} which computes phylogenetic proximities
+  between tips of a phylogeny.
+}
+\examples{
+
+## example(maples)
+W <- proxTips(tre, met="Abouheif")
+moran.idx(dom, W)
+moran.idx(bif, W)
+moran.idx(rnorm(nTips(tre)), W)
+
+## build a simple permutation test for 'bif'
+sim <- replicate(499, moran.idx(rnorm(nTips(tre)), W)) # permutations
+sim <- c(sim, moran.idx(bif, W))
+
+pval <- mean(sim>=sim[1]) # right-tail p-value
+pval
+
+hist(sim, col="grey", main="Moran's I Monte Carlo test for 'bif'") # plot
+mtext("Histogram of permutations and observation (in red)")
+abline(v=sim[1], col="red", lwd=3)
+
+}
+\keyword{manip}