[Vegan-commits] r892 - in pkg/vegan: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 15 00:33:33 CEST 2009
Author: gblanchet
Date: 2009-07-15 00:33:32 +0200 (Wed, 15 Jul 2009)
New Revision: 892
Added:
pkg/vegan/R/mantel.correlog.R
pkg/vegan/R/plot.mantel.correlog.R
pkg/vegan/R/print.mantel.correlog.R
pkg/vegan/man/mantel.correlog.Rd
Log:
Add functions mantel.correlog.R, print.mantel.correlog.R, and plot.mantel.correlog.R and documentation
Added: pkg/vegan/R/mantel.correlog.R
===================================================================
--- pkg/vegan/R/mantel.correlog.R (rev 0)
+++ pkg/vegan/R/mantel.correlog.R 2009-07-14 22:33:32 UTC (rev 892)
@@ -0,0 +1,106 @@
+'mantel.correlog' <-
+ function(D.eco, D.geo=NULL, XY=NULL, n.class=0, break.pts=NULL, cutoff=TRUE, r.type="pearson", nperm=999)
+{
+require(vegan)
+
+epsilon <- .Machine$double.eps
+D.eco = as.matrix(D.eco)
+
+# Geographic distance matrix
+if(!is.null(D.geo)) {
+ if(!is.null(XY)) stop("You provided both a geographic distance matrix and a list of site coordinates. Which one should the function use?")
+ D.geo = as.matrix(D.geo)
+ } else {
+ if(is.null(XY)) {
+ stop("You did not provided a geographic distance matrix nor a list of site coordinates.")
+ } else {
+ D.geo = as.matrix(dist(XY))
+ }
+ }
+
+n = nrow(D.geo)
+if(n != nrow(D.eco)) stop("Numbers of objects in D.eco and D.geo are not equal.")
+n.dist = n*(n-1)/2
+vec.D = as.vector(as.dist(D.geo))
+vec.DD = as.vector(D.geo)
+
+# Number of classes
+if(n.class > 0) {
+ if(!is.null(break.pts)) stop("You provided both a number of classes and a list of break points. Which one should the function use?")
+ } else {
+ if(is.null(break.pts)) {
+ # Use Sturge's rule to determine the number of classes
+ n.class = round(1 + log(n.dist, base=2))
+ start.pt = min(vec.D)
+ end.pt = max(vec.D)
+ width = (end.pt - start.pt)/n.class
+ break.pts = vector(length=n.class+1)
+ break.pts[n.class+1] = end.pt
+ for(i in 1:n.class) { break.pts[i] = start.pt + width*(i-1) }
+ } else {
+ # Use the list of break points
+ n.class = length(break.pts) - 1
+ }
+ }
+half.cl = n.class %/% 2
+
+# Move the first breakpoint a little bit to the leftepsilon <- .Machine$double.eps
+break.pts[1] = break.pts[1] - epsilon
+
+# Find the break points and the class indices
+class.ind = break.pts[1:n.class] + (0.5*(break.pts[2:(n.class+1)]-break.pts[1:n.class]))
+
+# Create the matrix of distance classes
+vec2 = vector(length=n^2)
+for(i in 1:n^2) { vec2[i] = min( which(break.pts >= vec.DD[i]) ) - 1 }
+
+# Start assembling the vectors of results
+class.index = NA
+n.dist = NA
+mantel.r = NA
+mantel.p = NA
+# check.sums = matrix(NA,n.class,1)
+
+# Create a model-matrix for each distance class, then compute a Mantel test
+for(k in 1:n.class) {
+ class.index = c(class.index, class.ind[k])
+ vec3 = rep(0, n^2)
+ sel = which(vec2 == k)
+ vec3[sel] = 1
+ mat.D2 = matrix(vec3,n,n)
+ diag(mat.D2) = 0
+ n.dis = sum(mat.D2)
+ n.dist = c(n.dist, n.dis)
+ if(n.dis==0) {
+ mantel.r = c(mantel.r, NA)
+ mantel.p = c(mantel.p, NA)
+ } else {
+ row.sums = apply(mat.D2,1,sum)
+ # check.sums[k,1] = length(which(row.sums == 0))
+ if((cutoff==FALSE) | !(cutoff==TRUE & k>half.cl & any(row.sums == 0))) {
+ temp = mantel(mat.D2, D.eco, method=r.type, permutations=nperm)
+ mantel.r = c(mantel.r, temp$statistic)
+ temp.p = temp$signif
+ if(temp$statistic >= 0) {
+ temp.p = ((temp.p*nperm)+1)/(nperm+1)
+ } else {
+ temp.p = (((1-temp.p)*nperm)+1)/(nperm+1)
+ }
+ mantel.p = c(mantel.p, temp.p)
+ } else {
+ mantel.r = c(mantel.r, NA)
+ mantel.p = c(mantel.p, NA)
+ }
+ }
+ }
+
+mantel.res = cbind(class.index, n.dist, mantel.r, mantel.p)
+mantel.res=mantel.res[-1,]
+rownames(mantel.res) = rownames(mantel.res,do.NULL = FALSE, prefix = "D.class.")
+colnames(mantel.res) = c("class.index", "n.dist", "Mantel.r", "Pr(Mantel)")
+
+# Output the results
+res = list(mantel.res=mantel.res, n.class=n.class, break.pts=break.pts, call=match.call() )
+class(res) = "mantel.correlog"
+res
+}
Added: pkg/vegan/R/plot.mantel.correlog.R
===================================================================
--- pkg/vegan/R/plot.mantel.correlog.R (rev 0)
+++ pkg/vegan/R/plot.mantel.correlog.R 2009-07-14 22:33:32 UTC (rev 892)
@@ -0,0 +1,10 @@
+'plot.mantel.correlog' <- function(x, ...)
+{
+how.many.classes = which(x$mantel.res[,3] != "NA")
+lim = max(how.many.classes)
+plot(x$mantel.res[1:lim,1],x$mantel.res[1:lim,3], xlab="Distance class index", ylab="Mantel r", pch=22)
+not.signif = which((x$mantel.res[1:lim,4] <= 0.05))
+points(x$mantel.res[not.signif,1], x$mantel.res[not.signif,3], pch=22, bg="black")
+abline(a=0, b=0, col="red")
+invisible()
+}
Added: pkg/vegan/R/print.mantel.correlog.R
===================================================================
--- pkg/vegan/R/print.mantel.correlog.R (rev 0)
+++ pkg/vegan/R/print.mantel.correlog.R 2009-07-14 22:33:32 UTC (rev 892)
@@ -0,0 +1,9 @@
+'print.mantel.correlog' <- function(x, ...)
+{
+ cat('\nMantel Correlogram Analysis\n')
+ cat('\nCall:\n','\n')
+ cat(deparse(x$call),'\n')
+ cat('\n')
+ printCoefmat(x$mantel.res, P.values=TRUE, signif.stars=TRUE, Pvalues = TRUE)
+ invisible(x)
+}
Added: pkg/vegan/man/mantel.correlog.Rd
===================================================================
--- pkg/vegan/man/mantel.correlog.Rd (rev 0)
+++ pkg/vegan/man/mantel.correlog.Rd 2009-07-14 22:33:32 UTC (rev 892)
@@ -0,0 +1,73 @@
+\name{mantel.correlog}
+\alias{mantel.correlog}
+\alias{print.mantel.correlog}
+\alias{plot.mantel.correlog}
+\title{ Mantel correlogram }
+\description{
+Function \code{mantel.correlog} computes a multivariate Mantel correlogram.
+}
+\usage{
+mantel.correlog(D.eco, D.geo=NULL, XY=NULL, n.class=0, break.pts=NULL,
+cutoff=TRUE, r.type="pearson", nperm=999)
+
+\method{plot}{mantel.correlog}(x,...)
+}
+
+\arguments{
+ \item{D.eco}{ An ecological distance matrix, with class either \code{dist} or \code{matrix}. }
+ \item{D.geo}{ A geographic distance matrix, with class either \code{dist} or \code{matrix}. Provide either \code{D.geo} or \code{XY}. Default: \code{D.geo=NULL}. }
+ \item{XY}{ A file of Cartesian geographic coordinates of the points. Default: \code{XY=NULL}. }
+ \item{n.class}{ Number of classes. If \code{n.class=0}, the Sturge equation will be used unless break points are provided. }
+ \item{break.pts}{ Vector containing the break points of the distance distribution. Default: \code{break.pts=NULL}. }
+ \item{cutoff}{ For the second half of the distance classes, \code{cutoff = TRUE} limits the correlogram to the distance classes that include all points. If \code{cutoff = FALSE}, the correlogram includes all distance classes. }
+ \item{r.type}{ Type of correlation in calculation of the Mantel statistic. Default: \code{r.type="pearson"}. Other choices are \code{r.type="spearman"} and \code{r.type="kendall"}, as in functions \code{\link{cor}} and \code{\link{mantel}}. }
+ \item{nperm}{ Number of permutations for the tests of significance. Default: \code{nperm=999}. For large data files, permutation tests are rather slow. }
+ \item{x}{ Output of \code{mantel.correlog}. }
+ \item{...}{ Parameters passed to other functions. }
+}
+
+\details{
+The Mantel tests are done by \code{vegan}'s function \code{\link{mantel}}. Function \code{\link[ape]{mantel.test}} of the \code{ape} package could also have been used.
+
+The \code{print.mantel.correlog} function prints out the correlogram. See example.
+}
+
+\value{
+
+ \item{mantel.res }{A table with the distance classes as rows and the class indices, number of distances per class, Mantel statistics, and p-values as columns. }
+ \item{n.class }{The number of distance classes. }
+ \item{break.pts }{The break points provided by the user or computed by the program. }
+ \item{call }{The function call. }
+}
+
+\author{ Pierre Legendre, Universite de Montreal }
+
+\references{
+Legendre, P. & L. Legendre. 1998. Numerical ecology, 2nd English edition. Elsevier Science BV, Amsterdam.
+
+Oden, N. L. & R. R. Sokal. 1986. Directional autocorrelation: an extension of spatial correlograms to two dimensions. Syst. Zool. 35: 608-617.
+
+Sokal, R. R. 1986. Spatial data analysis and historical processes. 29-43 in: E. Diday et al. [eds.] Data analysis and informatics, IV. North-Holland, Amsterdam.
+}
+
+\examples{
+# Mite data from "vegan"
+data(mite)
+data(mite.xy)
+mite.hel = decostand(mite, "hellinger")
+mite.hel.D = dist(mite.hel)
+
+mite.correlog = mantel.correlog(mite.hel.D, XY=mite.xy, nperm=99)
+summary(mite.correlog)
+mite.correlog
+plot(mite.correlog)
+
+mite.correlog2 = mantel.correlog(mite.hel.D, XY=mite.xy, cutoff=FALSE,
+r.type="spearman", nperm=99)
+summary(mite.correlog2)
+mite.correlog2
+plot(mite.correlog2)
+
+}
+
+\keyword{ multivariate }
More information about the Vegan-commits
mailing list