martes, 19 de julio de 2011

Test de Mantel en R

En su sentido mas básico el Test de Mantel es un análisis estadístico de correlación entre dos matrices del mismo rango, es usado usualmente para comparar dos matrices de distancia (por ejemplo una matriz de distancia genética entre individuos y una matriz de distancia geográfica entre estos mismos individuos).
Lógicamente en R es muy sencillo realizar este Test, en principio bajo cualquier coeficiente de correlación, pero es comúnmente mayor utilizado el coeficiente de correlación producto-momento de Pearson.


Nathan Mantel (un gordito con cara de buena onda)


Aquí utilizaremos los datos de cantidad de germenes en 4 heladerías!!! (en honor a Juliette fanatica de los helados); Se midió la distancia geográfica entre las heladerías para obtener una matriz y la otra matriz de distancia se obtuvo de los valores de especies compartidas entre pares de heladerías.


los datos los tenemos en dos archivos en formato .txt (aunque podrían estar en cualquier otro formato como .csv) de la siguiente manera:

germen.dist.txt

"hel1" "hel2" "hel3" "hel4" 
"hel1" 0 16.987 107.913 139.049
"hel2" 16.987 0 104.433 129.082 
"hel3" 107.913 104.433 0 62.435
"hel4" 139.049 129.082 62.435 0 


geo.dit.txt



"hel1" "hel2" "hel3" "hel4" 
"hel1" 0 2.54 10.735 14.054
"hel2" 2.54 0 8.254 11.624 
"hel3" 10.735 8.254 0 4.325
"hel4" 14.054 11.624 4.325 0 


Teniendo estos datos podemos empezar a calcular nuestro Test de Mantel.



> gerd <- read.table("germen.dist.txt",header=T)
# llamo mis datos de la matriz de distancia de germenes presentes. Header=T o header=TRUE quiere decir que mis datos poseen los nombres de las heladerias en la primer fila y la primer columna del archivo.

> geod <- read.table("geo.dist.txt",header=T)
#llamo mis datos de la matriz de distancia geográfica entre heladerías.
Header=T o header=TRUE quiere decir que mis datos poseen los nombres de las heladerias en la primer fila y la primer columna del archivo.

> gerd


#visualizo mis datos



        hel1    hel2    hel3    hel4
hel1   0.000  16.987 107.913 139.049
hel2  16.987   0.000 104.433 129.082
hel3 107.913 104.433   0.000  62.435
hel4 139.049 129.082  62.435   0.000

> geod
#visualizo mis datos

       hel1   hel2   hel3   hel4
hel1  0.000  2.540 10.735 14.054
hel2  2.540  0.000  8.254 11.624
hel3 10.735  8.254  0.000  4.325
hel4 14.054 11.624  4.325  0.000

> install.packages ("vegan")
#instalo el paquete vegan.


> library (vegan)
#cargo el paquete vegan a R, para utilizar la función de test de mantel.

>mantel.germenes<-mantel(geod, gerd)
# Realizo el test de mantel para examinar la correlación entre las dos matrices de datos. Sí no se especifica el método o el coeficiente de correlación a utilizar, se utilizara por defecto el coeficiente de correlación producto-momento de Pearson.

>mantel.germenes
#visualizo mis resultados.

Mantel statistic based on Pearson's product-moment correlation 
Call:                                                                                                                                                                                                                            
mantel(xdis = geod, ydis = gerd)                                                                                                                                                                          
      Mantel statistic r: 0.9634                                                                                                                                                                                                       
      Significance: 0.036                                                                                                                 
Empirical upper confidence limits of r:                                                                                                                                                                                          
  90%   95% 97.5%   99%                                                                                                                                                                                                          
0.865 0.931 0.963 0.963                                                                                                                  
Based on 999 permutations

También se pueden utilizar dos coeficientes de correlacion mas (Kendall y Spearman)
mas información. 


Adjunto tambien otra forma de hallar matrices de distancia y hacer un test de mantel, a partir de datos que tengamos:




distclim.csv



estacion,SAL,TEM,O,NH4,N3,N2,P4,S04, CLA,COL,FIT,TUR,ZOT
1,5.0,22.0,2.1,1.6,1.2,2.3,2.2,16.0,1.4,100,121,9,60
2,9.0,20.0,2.2,2.3,1.4,2.5,2.3,12.0,1.4,99,125,6,75
3,11.0,18.0,2.6,0.1,2.0,0.6,1.7,13.0,1.4,89,228,9,62
4,15.0,16.0,3.2,0.0,1.3,0.5,0.7,9.0,0.8,97,245,14,121
5,20.0,18.0,3.1,0.1,0.1,0.0,0.3,6.0,0.8,97,258,14,159
6,25.0,15.0,3.5,0.1,0.5,0.2,0.4,4.0,0.6,75,287,15,186
7,24.0,18.0,4.2,0.1,0.3,0.6,0.4,3.0,1.2,72,175,16,98
8,6.0,19.0,2.7,0.1,0.6,2.5,1.4,13.0,1.6,70,103,22,421
9,10.0,19.0,3.1,0.1,1.1,1.8,1.6,9.0,1.8,85,99,21,369
10,15.0,21.0,4.9,0.1,1.3,1.6,0.6,7.0,1.0,72,145,19,251
11,18.0,18.0,4.6,0.1,2.0,1.5,0.3,5.0,0.6,56,158,16,187
12,22.0,17.0,5.2,0.1,2.0,1.2,0.3,3.0,0.4,55,167,12,111
13,24.0,16.0,4.8,0.1,2.1,1.1,0.2,1.0,0.4,42,262,10,105
14,31.0,16.0,4.2,0.1,1.7,0.2,0.2,6.0,0.6,49,112,8,192
15,30.0,16.0,4.7,0.1,1.5,0.2,0.2,8.0,0.8,23,156,8,214
16,32.0,15.0,4.5,0.5,1.3,0.2,0.3,10.0,0.6,26,128,10,227
17,33.0,16.0,5.1,0.4,1.0,0.4,0.4,6.0,0.2,15,137,12,269
18,35.0,16.0,5.2,0.3,0.7,0.8,0.3,1.0,0.2,15,139,15,301
19,35.0,15.0,4.9,0.1,0.6,0.7,0.2,3.0,0.2,13,241,18,441
20,35.0,16.0,5.1,0.1,0.9,0.3,0.2,2.0,0.3,19,258,19,379

distesp.csv


estacion,AC,OI,SP,EA,NA,PA,TY,FR,OA,CD,CO,CA
1,15,0,0,5,12,0,0,0,0,0,0,0
2,22,3,0,3,14,0,0,0,0,0,0,0
3,37,5,0,8,17,0,0,0,0,0,0,0
4,21,3,0,1,9,0,0,0,0,0,1,0
5,14,12,1,8,16,0,2,3,0,0,0,0
6,21,6,0,14,6,6,3,5,0,2,2,2
7,19,8,0,12,11,5,1,8,0,3,0,3
8,9,3,5,3,9,11,4,0,0,0,0,5
9,15,0,9,12,7,15,10,0,0,0,0,4
10,11,0,8,9,8,18,15,0,0,0,0,7
11,8,0,14,6,5,17,6,0,0,0,0,8
12,13,0,11,5,6,9,14,0,0,0,0,3
13,24,1,9,7,4,3,5,10,0,0,0,9
14,8,2,13,4,2,5,6,12,0,0,0,1
15,6,2,18,12,7,1,9,14,0,26,0,9
16,0,8,16,13,18,2,3,13,0,23,10,14
17,0,9,14,15,14,0,0,20,18,19,2,10
18,3,10,35,18,20,0,0,18,15,17,15,15
19,0,9,22,11,14,0,0,22,19,13,13,18
20,0,7,24,15,11,0,0,17,23,18,9,17



                                             
>  distspe <- read.csv("distesp.csv",header=T) 
# (llamo mis datos de las especies en las estaciones)

>  distclima <- read.csv("distclim.csv",header=T) 
# (llamo mis datos de las observaciones climaticas por estaciones)

>  mdclima <- dist(cbind(distclima$SAL, distclima$TEM, distclima$O, distclima$NH4, distclima$N3, distclima$N2, distclima$P4, distclima$SO4, distclima$CLA, distclima$CO2, distclima$FIT, distclima$TUR, distclima$ZOT)) 
# (creo mi matriz de distancias a partir de mis datos climaticos, esta formula es la via larga en donde se escojen las variables que van a ser tomadas de los datos iniciales con el comando "cbind", lo unico que se necesita hacer es escojer las variables que deso incluir para crear mi matriz)

> mdspe <- dist(cbind(distspe$AC,distspe$OI, distspe$SP, distspe$EA, distspe$NA., distspe$PA, distspe$TY, distspe$FR, distspe$OA, distspe$CD, distspe$CO, distspe$CA )) 
# (realizo otra matriz de distamcia de mis datos de especies por estaciones por la via larga escojiendo mis variables con "cbind")

>mdclima <- dist(distclima)
 # (si tengo pereza y ademas quiero incluir todas las variables de mis datos iniciales de clima, solo utilizo la via rapida para crear mi matriz de distancias)

> mdspe <- dist(distspe) 
# (lo mismo para mis datos de especie)
 
> mantelhilton<-mantel.rtest(mdspe, mdclima, nrepet = 9999) 
# (por ultimo realizo el test de mantel con 9999 replicas simuladas) 


Hasta la vista!!!