[R] Mantel's randomization test

[Takashi Mizuno]

Dear all,

Thanks for any reply to my question.
Anyway I wrote a simple program for Mantel's test in R to my purpose.
Following is a first code, dirty, nothing some comments, ...uuum.


------------------------
Takashi Mizuno
zoono at sci.osaka-cu.ac.jp
Plant Ecology Lab.
Osaka City University
------------------------

--------------BEGIN-----------------------------
## _sdmatrix()_ make the spatial distances matrix from a given count
matrix.
## This is only for grided data that an area is divided up into
## n equal-size quadrats.
## Matrix returned is a symmetric with zero diagonal terms.
sdmatrix <- function( X )
{
  n <- dim(X)[2]
  m <- dim(X)[1]
  L <- length(X)
  Z <- matrix(NA,L,L)
  for( i in 1:(L-1) ){
    if( i%%m == 0 ){
      x1 <- m
      y1 <- i%/%m
    }
    else{
      x1 <- i%%m
      y1 <- i%/%m + 1
    }
    for(j in i:L){
      if( j > i ){
        if( j%%m == 0 ){
          x2 <- m
          y2 <- j%/%m
        }
        else{
          x2 <- j%%m
          y2 <- j%/%m + 1
        }
        Z[j,i] <- Z[i,j] <- sqrt( (x2-x1)^2 + (y2-y1)^2 )
      }
    }
  }
  Z
}

## _cdiffmatrix()_ make the count difference matrix from a given matrix.
## This is only for grided data that an area is divided up into
## n equal-size quadrats.
## Matrix returned is a symmetric with zero diagonal terms.
cdiffmatrix <- function( X )
{
  n <- dim(X)[2]
  m <- dim(X)[1]
  L <- length(X)
  Z <- matrix(NA,L,L)
  for( i in 1:(L-1) ){
    if( i%%m == 0 ){
      x1 <- m
      y1 <- i%/%m
    }
    else{
      x1 <- i%%m
      y1 <- i%/%m + 1
    }
    for(j in i:L){
      if( j > i ){
        if( j%%m == 0 ){
          x2 <- m
          y2 <- j%/%m
        }
        else{
          x2 <- j%%m
          y2 <- j%/%m + 1
        }
        Z[j,i] <- Z[i,j] <- X[x2,y2] - X[x1,y1]
      }
    }
  }
  for(i in 1:L){
    Z[i,i] <- 0
  }
  Z
}

## some functions
mantelz <- function( X, Y )
{
  n <- dim(X)[1]
  z <- 0.0
  for( i in 1:(n-1) ){
    for( j in i:n ){
      if(j>i){
        z <- z + X[j,i] * Y[j,i]
      }
    }
  }
  z
}

mantelsum <- function( X )
{
  n <- dim(X)[1]
  s <- 0.0
  for( i in 1:(n-1) ){
    for( j in i:n ){
      if(j>i){
        s <- s + X[j,i]
      }
    }
  }
  s
}

mantelpow2 <- function(X)
{
  n <- dim(X)[1]
  s <- 0.0
  for( i in 1:(n-1) ){
    for( j in i:n ){
      if(j > i){
        s <- s + X[j,i]^2
      }
    }
  }
  s
}

rsdmatrix <- function( X )
{
  Z <- 1/X
  for(i in 1:dim(X)[1] ){
    Z[i,i] <- 0
  }
  Z
}

mantelr <- function(X,Y)
{
  D <- dim(X)[1]
  D <- D*(D-1)/2
  sigmaab <- mantelz(X,Y)
  suma <- mantelsum(X)
  sumb <- mantelsum(Y)
  sigmaaa <- mantelz(X,X)
  sigmabb <- mantelz(Y,Y)
  r <- ( sigmaab - (suma*sumb)/D ) /
    sqrt( (sigmaaa - (suma^2)/D) * (sigmabb - (sumb^2)/D) )
  r
}

##
## Mantel's randomization test main part
##
mantel.test <- function(X,Y,nrpt=4999)
{
  D <- dim(X)[1]
  r1 <- mantelr(X,Y)
  res <- numeric(nrpt)
  pM <- matrix(NA,dim(X)[1],dim(X)[2])
  permut <- numeric(D)
  RNGkind("Mersenne-Twister")
  for( i in 1:nrpt ){
    rs <- .Random.seed
    permut <- rank(runif(D))
    for( j in 1:(D-1) ){
      for( k in j:D ){
        if( permut[k] > permut[j]){
          pM[k,j] <- X[permut[k],permut[j]]
        }
        else{
          pM[k,j] <- X[permut[j],permut[k]]
        }
      }
    }
    res[i] <- mantelr(pM,Y)
  }
  res <- sort(res)
  over <- res[res>=r1]
  lower <- res[res<=r1]
  on <- length(over)
  ln <- length(lower)
  p.over <- (on+1)/(nrpt+1)
  p.lower <- (ln+1)/(nrpt+1)
  list(result=res,over=on,lower=ln,p.over=p.over,p.lower=p.lower,r=r1)
}
---------------------END----------------------
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._