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