[Vegan-commits] r943 - pkg/vegan/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Aug 23 19:25:05 CEST 2009
Author: jarioksa
Date: 2009-08-23 19:25:04 +0200 (Sun, 23 Aug 2009)
New Revision: 943
Modified:
pkg/vegan/R/mantel.correlog.R
pkg/vegan/R/plot.mantel.correlog.R
Log:
formatting in standard style
Modified: pkg/vegan/R/mantel.correlog.R
===================================================================
--- pkg/vegan/R/mantel.correlog.R 2009-08-23 13:14:07 UTC (rev 942)
+++ pkg/vegan/R/mantel.correlog.R 2009-08-23 17:25:04 UTC (rev 943)
@@ -1,140 +1,153 @@
-'mantel.correlog' <-
- function(D.eco, D.geo=NULL, XY=NULL, n.class=0, break.pts=NULL, cutoff=TRUE, r.type="pearson", nperm=999, mult="holm", progressive=TRUE)
+`mantel.correlog` <-
+ function(D.eco, D.geo=NULL, XY=NULL, n.class=0, break.pts=NULL,
+ cutoff=TRUE, r.type="pearson", nperm=999, mult="holm",
+ progressive=TRUE)
{
-require(vegan)
+ r.type <- match.arg(r.type, c("pearson", "spearman", "kendall"))
+ mult <- match.arg(mult, c("sidak", p.adjust.methods))
-r.type <- match.arg(r.type, c("pearson", "spearman", "kendall"))
-mult <- match.arg(mult, c("sidak", p.adjust.methods))
+ epsilon <- .Machine$double.eps
+ D.eco = as.matrix(D.eco)
-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 {
+ ## 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))
- }
- }
+ 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)
-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 {
+ ## 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
+ ## 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 left
-epsilon <- .Machine$double.eps
-break.pts[1] = break.pts[1] - epsilon
+ ## Move the first breakpoint a little bit to the left
+ 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]))
+ ## 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 }
+ ## 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)
+ ## 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)
- }
- }
- }
+ ## 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*n)
+ 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,]
-mantel.res = cbind(class.index, n.dist, mantel.r, mantel.p)
-mantel.res = mantel.res[-1,]
-
-## Note: vector 'mantel.p' starts with a NA value
-mantel.p = mantel.p[-1]
-n.tests = length(which(mantel.p != "NA"))
-
-if(mult=="none") {
- colnames(mantel.res) = c("class.index", "n.dist", "Mantel.cor", "Pr(Mantel)")
- } else {
+ ## Note: vector 'mantel.p' starts with a NA value
+ mantel.p <- mantel.p[-1]
+ n.tests <- length(which(!is.na(mantel.p)))
+
+ if(mult == "none") {
+ colnames(mantel.res) <-
+ c("class.index", "n.dist", "Mantel.cor", "Pr(Mantel)")
+ } else {
## Correct P-values for multiple testing
- if(progressive) {
- p.corr = mantel.p[1]
- if(mult == "sidak") {
- for(j in 2:n.tests) { p.corr = c(p.corr, 1-(1-mantel.p[j])^j) }
- } else {
- for(j in 2:n.tests) {
- temp = p.adjust(mantel.p[1:j], method=mult)
- p.corr = c(p.corr, temp[j])
- }
- }
- } else {
- ## Correct all p-values for 'n.tests' simultaneous tests
- if(mult == "sidak") {
- p.corr = 1-(1-mantel.p[1:n.tests])^n.tests
- } else {
- p.corr = p.adjust(mantel.p[1:n.tests], method=mult)
- }
- }
- temp = c(p.corr, rep(NA,(n.class-n.tests)))
- mantel.res = cbind(mantel.res, temp)
- colnames(mantel.res) = c("class.index", "n.dist", "Mantel.cor", "Pr(Mantel)", "Pr(corrected)")
- }
-rownames(mantel.res) = rownames(mantel.res,do.NULL = FALSE, prefix = "D.cl.")
+ if(progressive) {
+ p.corr <- mantel.p[1]
+ if(mult == "sidak") {
+ for(j in 2:n.tests) {
+ p.corr <- c(p.corr, 1-(1-mantel.p[j])^j) }
+ } else {
+ for(j in 2:n.tests) {
+ temp <- p.adjust(mantel.p[1:j], method=mult)
+ p.corr <- c(p.corr, temp[j])
+ }
+ }
+ } else {
+ ## Correct all p-values for 'n.tests' simultaneous tests
+ if(mult == "sidak") {
+ p.corr <- 1 - (1 - mantel.p[1:n.tests])^n.tests
+ } else {
+ p.corr <- p.adjust(mantel.p[1:n.tests], method=mult)
+ }
+ }
+ temp <- c(p.corr, rep(NA,(n.class-n.tests)))
+ mantel.res <- cbind(mantel.res, temp)
+ colnames(mantel.res) <-
+ c("class.index", "n.dist", "Mantel.cor", "Pr(Mantel)", "Pr(corrected)")
+ }
+ rownames(mantel.res) <-
+ rownames(mantel.res,do.NULL = FALSE, prefix = "D.cl.")
+
+ ## Output the results
+ res <- list(mantel.res=mantel.res, n.class=n.class, break.pts=break.pts,
+ mult=mult, n.tests=n.tests, call=match.call() )
+ class(res) <- "mantel.correlog"
+ res
+}
-# Output the results
-res = list(mantel.res=mantel.res, n.class=n.class, break.pts=break.pts, mult=mult, n.tests=n.tests, call=match.call() )
-class(res) = "mantel.correlog"
-res
-}
\ No newline at end of file
Modified: pkg/vegan/R/plot.mantel.correlog.R
===================================================================
--- pkg/vegan/R/plot.mantel.correlog.R 2009-08-23 13:14:07 UTC (rev 942)
+++ pkg/vegan/R/plot.mantel.correlog.R 2009-08-23 17:25:04 UTC (rev 943)
@@ -1,15 +1,17 @@
-'plot.mantel.correlog' <- function(x, alpha=0.05, ...)
+`plot.mantel.correlog` <-
+ function(x, alpha=0.05, ...)
{
-lim = max(x$n.tests)
-plot(x$mantel.res[1:lim,1],x$mantel.res[1:lim,3], xlab="Distance class index", ylab="Mantel correlation", pch=22)
-if(x$mult=="none") {
- signif = which((x$mantel.res[1:lim,4] <= alpha))
- } else {
- signif = which((x$mantel.res[1:lim,5] <= alpha))
- }
-lines(x$mantel.res[1:lim,1], x$mantel.res[1:lim,3])
-points(x$mantel.res[1:lim,1], x$mantel.res[1:lim,3], pch=22, bg="white")
-points(x$mantel.res[signif,1], x$mantel.res[signif,3], pch=22, bg="black")
-abline(a=0, b=0, col="red")
-invisible()
-}
\ No newline at end of file
+ lim <- max(x$n.tests)
+ plot(x$mantel.res[1:lim,1],x$mantel.res[1:lim,3],
+ xlab="Distance class index", ylab="Mantel correlation", pch=22)
+ if(x$mult == "none") {
+ signif <- which((x$mantel.res[1:lim,4] <= alpha))
+ } else {
+ signif <- which((x$mantel.res[1:lim,5] <= alpha))
+ }
+ lines(x$mantel.res[1:lim,1], x$mantel.res[1:lim,3])
+ points(x$mantel.res[1:lim,1], x$mantel.res[1:lim,3], pch=22, bg="white")
+ points(x$mantel.res[signif,1], x$mantel.res[signif,3], pch=22, bg="black")
+ abline(a=0, b=0, col="red")
+ invisible()
+}
More information about the Vegan-commits
mailing list