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







7 comentarios:

  1. hola me gusto tu blog. pero me quedo una duda con este test, los datos deben ser transformados bajo alguna condición o son utilizados directamente. me quedo esa duda por los valores que utilizas de temperatura, oxigeno y amonio.
    .

    Pablo Muñoz Ortiz
    tesista Biología Marina
    Universidad de Concepción.
    pabmunoz@udec.cl

    ResponderEliminar
  2. Hola, definitivamente la transformación de los datos es algo inherente a sus necesidades y tipos de datos; los datos a los que se refiere son un ejemplo de cualquier matriz de datos y lo que buscaba era en realidad explicar la implementación del test en R, yo los use directamente en el test, pero obviamente pueden ser transformados, siguiendo las reglas básicas de transformación, para que la transformación no suene como un truco sucio para acomodar los resultados :)

    ResponderEliminar
  3. me ayudaste mucho en mi tesis.. estoy muy agradecido

    ResponderEliminar
  4. Hola: ¿Es posible hacer una prueba de permutaciones en R para medir el grado de correlación? Tengo 10 columnas en mi base de datos de excel 4 de ellas variables independientes y el resto dependientes. Quisiera ver cómo se introducen los datos y la interpretación en R. Gracias desde ya.

    ResponderEliminar
  5. Este comentario ha sido eliminado por el autor.

    ResponderEliminar
  6. Hola, si es posible hacer este tipo de pruebas en R para medir correlación entre variables. Para introducir datos en R le recomiendo otro post de "R para chibchombianos" http://rchibchombia.blogspot.com/2011/12/entrada-y-manejo-de-datos-en-r-guia.html
    con este, todas sus dudas de como introducir su matriz de datos seran resuelta, sí no entiende algo, me escribe de nuevo y yo le ayudare con eso; ahora para hacer test de permutaciones para correlación en R, hay un manual muy sencillo, que esta en ingles pero es muy fácil de utilizar; http://www.psych.umn.edu/faculty/waller/downloads/mpt/mptcorr.pdf.
    Cualquier cosa que depronto le genere confusión o no este funcionando, acá estaré.
    mucha suerte.

    ResponderEliminar
  7. estimado, me puedes ayudar a interpretar mis resultados..
    Mantel statistic r: 0.09782
    Significance: 0.676

    Si no es significativo (0.676) entonces las variables no cambian según su posición geográfica?
    Y el r que indica?
    Muchas gracias!

    ResponderEliminar