R para Chibchombianos: julio 2011

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

martes, 19 de julio de 2011

Test de Mantel en R