[R] Mantel's randomization test

[ben@zoo.ufl.edu]

  This looks OK (I really can't comment on style points a whole lot).  The
distance matrix is obviously specific to your case ... the only thing I
can point out is that you can often be more efficient in R with vector
operations than using for-loops ... (it's my understanding that implicitly
vectorized operations and matrix operators are faster than apply() and
friends, which are in turn faster than for-loops).

  Ben Bolker

On Tue, 1 May 2001, Takashi Mizuno wrote:

> 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
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>

-- 
318 Carr Hall                                bolker at zoo.ufl.edu
Zoology Department, University of Florida    http://www.zoo.ufl.edu/bolker
Box 118525                                   (ph)  352-392-5697
Gainesville, FL 32611-8525                   (fax) 352-392-3704

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._