[R] Mantel's randomization test

[ben@zoo.ufl.edu]

  Years ago I coded this in S-PLUS ... I think it's been posted to the
S-news list before.  It seems to work just fine in R.  I have not tested
it at all, so it comes with the usual warranty (i.e., none).

  Just cut on the dotted lines and source() the file.

  One example (someday I'll package this all up and write a help file
etc.).

q1 _ matrix(runif(36),nrow=6)
q2 _ matrix(runif(36),nrow=6)
mantel(q1,q2,graph=TRUE,rpt=200)

  [Note to R web page maintainers: perhaps there should be a
general-purpose "How do I find out if R can ... ?" entry, which tells
people to try

  (a) searching the documentation
  (b) searching the list of contributed packages
  (c) searching the mailing lists
  (d) searching the S-news mailing lists, which sometimes have usable
code

  Yes, I've never gotten around to it (I don't run my own web server), but
someone could do a service to the community by settin up an ht://dig
server that indexed all of those sources ... ]


--------------- cut here -----------------------
# Mantel test for similarity of two matrices via permutation
# and calculation of Mantel "z statistic"
#  (and associated utility functions)
#
#  from Bryan F.J. Manly, <-Multivariate statistical methods: a primer<-
#         (Chapman & Hall, 1986)
####
#  TO DO:
#      (mostly bells and whistles)
#
#   add normalized Z-stat (or option for it)
#   add asymptotic conf. limits/tests (Sokal and Rolf (Biometry, 3d ed.) cite
#       Mantel (1967), "The detection of disease clustering ..." Cancer
#       Res. 27:209-220, and suggest the asymptotic results should be OK
#        for matrices as small as 25x25)
#   add 3- or multi-way extensions of Mantel (S&R cite Smouse, Long, and
#        Sokal (1986) "Multiple regression and correlation extensions of the
#        Mantel test of matrix correspondence" Syst. Zool. 35:627-632 and
#        Oden and Sokal (1992) "An investigation of three-matrix permutation
#        tests" J. Classif. 9:275-290.
#   documentation

### utility functions

factorial <- function(n)gamma(n+1)

perm.rowscols <- function(m1,n,p=NULL)
{
# permute rows and columns of a matrix
# m1 is a (square) nxn matrix
# if p is numeric then the pth permutation will be taken, otherwise
#   a random permutation will be taken
   if (is.numeric(p))
     s <- gen.perm(p,n)
   else
     s <- sample(1:n)
   m1[s,s]
}

mant.zstat <- function(m1,m2)
{
#  calculate the Mantel z-statistic for two square matrices m1 and m2
   sum(lower.triang(m1*m2))
}

gen.perm <- function(k,n,vec=1:n)
{
#  generate a specified permutation of a vector
#  k is an integer between 1 and n!
#  n is the length of the sample
#  vec is the vector to permute, 1:n by default
   p <- as.list(vec)
   perm <- numeric(n)
   for (i in (n:1)) {
      v <- (k%%i)+1
      perm[i] <- p[[v]]
      p[[v]] <- NULL
      k <- k %/% i
   }
   perm
}

all.perm <- function(n,vec=1:n)
{
# generate all permutations of a list of length n
#  (1:n by default)
   k <- matrix(NA,nrow=factorial(n),ncol=n)
   for (i in (1:factorial(n)))
       k[i,] <- gen.perm(i,n,vec)
   k
}

lower.triang <- function(m)
{
   d <- dim(m)
   if (d[1] != d[2]) print("Warning: non-square matrix")
   m[col(m)<row(m)]
}

mantel <- function(m1,m2,nperm=250,graph=F,rpt=F)
{
   n <- dim(m1)[1]
   realz <- mant.zstat(m1,m2)
   enum<-F
   if (nperm>factorial(n)) {
      cat("Warning: nperm>n!, enumerating permutations =",factorial(n),"\n")
      nperm <- factorial(n)
      enum<-T
   }
   nullstats <- numeric(nperm)
   for (i in (1:nperm)) {
      if (enum)
         nullstats[i] <- mant.zstat(m1,perm.rowscols(m2,n,i))
      else
         nullstats[i] <- mant.zstat(m1,perm.rowscols(m2,n))
      if (is.numeric(rpt))
         if (i %% rpt == 0)
            cat(i,"...\n")
   }
   nullstats <- sort(nullstats)
   pval <- 1-((1:nperm)[nullstats>realz])[1]/nperm
   if (graph) {
      plot(density(nullstats),type="l",
       main="Distribution of Mantel z-statistic",
       xlab="Z statistic",ylab="# of permutations",
       sub=paste("Actual z-stat =",round(realz,3),
         ": p<",round(pval,4),",",nperm," permutations"))
      abline(v=realz)
   }
   list(z.stat=realz,p=pval)
}

-------------------- cut here ------------------------

On Sun, 29 Apr 2001, Takashi Mizuno wrote:

> Dear all,
>
> Dose anyone know whether there is a good R packege or
> program for Mantel's randomization test?
>
>
> Thanks in advance.
>
> ------------------------
> Takashi Mizuno
> zoono at sci.osaka-cu.ac.jp
> Plant Ecology Lab.
> Osaka City University
> ------------------------
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._