[R] partial mantel

[Marcelino de la Cruz]

----- Original Message -----
From: Bouget Christophe
Sent: 9/12/2003 7:19:14 AM
To: r-help at stat.math.ethz.ch
Subject: [R] partial mantel
 > Dear all,
 > Has anyone written R code for partial Mantel Tests- as described for
 > instance in Legtendre & Legendre (1998) ?
 > In other words, in a community ecology analysis, I would like to calculate
 > the correlation between two dissimilarity matrices, controlling for a third
 > distance matrix representing geographical distances between sites.
 > Thanks!
 >
 > Christophe Bouget
 > Biodiversité et gestion des forêts de plaine - Ecosylv
 > Ecosystèmes forestiers et paysages
 > Institut de recherches pour l'ingénierie de l'agriculture et de
 > l'environnement - CEMAGREF
 > Domaine des Barres
 > F-45 290 Nogent-sur-Vernisson
 >

Some time ago I wrote the following naive function to implement the method 
2 of Legendre in  J. Statist. Comput. Simul. , 2000, Vol. 67, pp. 37 – 73 
(permutation of residuals of null model).
Acording to Legendre, it can always be used, except when highly skewed data 
are combined with small  sample size (n < 20, or n < 50 in the presence of 
outliers). I have also R code for methods 1 and 3 available.
I hope it can help.

***********************************************************************************************
PARTIAL MANTEL TEST
METHOD 2
************************************************************************************************************
#a,b,c are distance matrix (from dist())
#nsim are required simulations
************************************************************************************************************
partial.mantel2<-function(a,b,c,nsim=100){
library(mva)
library(cluster)
m<- matrix(0,(dim(as.matrix(a))[1]),dim(as.matrix(a))[1])
m[lower.tri(m)] <- residuals(lm(as.vector(a)~as.vector(c)))
resA.C<-as.dist (m)
a.sd<-((a-mean(a))/sqrt(var(a)))  #estandarización de las matrices de 
distancias
b.sd<-((b-mean(b))/sqrt(var(b)))
c.sd<-((c-mean(c))/sqrt(var(c)))
n<- length(a)-1
rmAB<-drop(crossprod(a.sd,b.sd)/n)
rmAC<-drop(crossprod(a.sd,c.sd)/n)
rmBC<-drop(crossprod(b.sd,c.sd)/n)
rmAB.C<-(rmAB-(rmAC*rmBC))/(sqrt(1-rmAC^2)*sqrt(1-rmBC^2))
distribucion<-rep(-999,nsim)
resA.C.m<-as.matrix(resA.C)
resA.C.m.sim<-as.matrix(resA.C)

for (h in 1:nsim){
     x<-sample(1:dim(as.matrix(a.sd))[1]) # permutación de la matriz a 
“nsim” veces
     for (i in 1:length(x)){
        for (j in 1:length (x)){
         resA.C.m.sim[i,j]<-resA.C.m[x[i],x[j]]
         }
     }
    rmRB<-drop(crossprod(as.dist(resA.C.m.sim),b.sd)/n)
    rmRC<- drop(crossprod(as.dist(resA.C.m.sim),c.sd)/n)
    distribucion[h]<-(rmRB-(rmRC*rmBC))/(sqrt(1-rmRC^2)*sqrt(1-rmBC^2))
}
distribucion<-c(rmAB.C,distribucion)

if(rmAB.C>0) 
p<-length(which(sort(distribucion)>=rmAB.C))/length(distribucion) else p<- 
length(which(sort(distribucion)<=rmAB.C))/length(distribucion)


print(c(“rMab.c:”,rmAB.C))
if(rmAB.C>0) print(c("prop rM < r*M:",p)) else print(c("prop rM > r*M:", p))
}


*****************************************************************************************************************************








***************************************
Marcelino de la Cruz Rot
Depto. Biología Vegetal
ETS Ingenieros Agrónomos
Universidad Politécnica de Madrid
28040-Madrid

Tel: 91.336.56.63.