# [R] Mantel's randomization test

Takashi Mizuno zoono at sci.osaka-cu.ac.jp
Tue May 1 04:25:59 CEST 2001

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

```