----- Original Message ----- From: Bouget Christophe Sent: 9/12/2003 7:19:14 AM To: [EMAIL PROTECTED] 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.
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
