R/mantel.test.R

   1 ## mantel.test.R (2012-03-16)
   2
   3 ##   Mantel Test for Similarity of Two Matrices
   4
   5 ## Copyright 2002-2011 Ben Bolker and Julien Claude
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 perm.rowscols <- function(m1, n)
  11 {
  12     s <- sample(1:n)
  13     m1[s, s]
  14 }
  15
  16 ### calculate the Mantel z-statistic for two square matrices m1 and m2
  17 mant.zstat <- function(m1, m2) sum(lower.triang(m1 * m2))
  18
  19 lower.triang <- function(m)
  20 {
  21     if (diff(dim(m))) print("Warning: non-square matrix")
  22     m[col(m) <= row(m)]
  23 }
  24
  25 mantel.test <- function(m1, m2, nperm = 999, graph = FALSE,
  26                         alternative = "two.sided", ...)
  27 {
  28     alternative <- match.arg(alternative, c("two.sided", "less", "greater"))
  29     n <- nrow(m1)
  30     realz <- mant.zstat(m1, m2)
  31     nullstats <- replicate(nperm, mant.zstat(m1, perm.rowscols(m2, n)))
  32     pval <- switch(alternative,
  33                    "two.sided" =  2 * min(sum(nullstats >= realz), sum(nullstats <= realz)),
  34                    "less" = sum(nullstats <= realz),
  35                    "greater" = sum(nullstats >= realz))
  36     pval <- (pval + 1) / (nperm + 1) # 'realz' is included in 'nullstats'
  37     if (graph) {
  38         plot(density(nullstats), type = "l", ...)
  39         abline(v = realz)
  40     }
  41     list(z.stat = realz, p = pval, alternative = alternative)
  42 }