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