Análisis de Comunidades: ANOSIM y NMDS
junio 24, 2022
1 ANÁLISIS DE COMUNIDADES: ANOSIM y NMDS
Autor: Álvaro Alonso Fernández
Departamento de Ciencias de la Vida
Universidad de Alcalá (España)
Citation:
Alonso A (2022) Análisis de comunidades: ANOSIM y NMDS. Rpubs, accesed FECHA, https://rpubs.com/aafernandez1976/ANOSIMyNMDS
2 Nuestro objetivo
Nuestro objetivo es analizar la posible existencia de diferencias en la estructura de las comunidades de peces entre diferentes puntos de muestreo a lo largo de un río (tramo alto, medio y bajo). Para ello analizaremos la abundancia de 11 especies de peces en 12 puntos de muestreo, cuatro en tramo alto, cuatro en medio y cuatro en bajo. El primer paso consistirá en un análisis ANOSIM (análisis de similitudes), que es una prueba estadística no paramétrica muy empleada para el análisis de comunidades. Al ser no paramétrica nuestros datos no tiene el requisito de la normalidad (podemos tener muchos ceros en nuestra matriz sin preocuparnos por ello). El estadístico resultante (R) nos indica con valores cercanos a 1 la mayor similitud entre comunidades, cuanto más cercano a cero más parecidas son las comunidades entre ellas.
3 Nuestros datos
Cargamos nuestra matriz de datos, que tendrá el siguiente aspecto:
setwd(dir = "F:/R/MARKDOWN/ANOSIMyNMDS/")
datos<-read.table("ANOSIMyNMDS.csv", sep = ";", header = TRUE, dec = ".")
datos## sample tramo sp1 sp2 sp3 sp4 sp5 sp6 sp7 sp8 sp9 sp10 sp11
## 1 1 upper 3 8 9 3 8 10 2 5 10 2 15
## 2 2 upper 6 6 7 6 5 9 10 6 12 6 17
## 3 3 upper 3 9 10 1 7 8 9 4 10 6 20
## 4 4 upper 3 7 6 9 4 9 10 4 14 5 18
## 5 5 medium 10 9 9 8 12 11 14 8 9 16 11
## 6 6 medium 12 14 11 16 15 13 8 15 10 8 11
## 7 7 medium 10 14 13 11 16 16 9 8 8 13 10
## 8 8 medium 13 15 14 14 17 10 11 8 5 13 10
## 9 9 lower 10 20 13 21 11 13 10 10 5 15 1
## 10 10 lower 11 16 17 20 19 19 20 21 8 14 6
## 11 11 lower 21 22 22 15 10 20 21 21 7 21 3
## 12 12 lower 11 14 21 17 20 22 19 15 9 11 4En cuanto a su estructura, tenemos un identificador de la muestra sampley del tramo. el resto de columnas son cada una de las 11 especies estudiadas:
str(datos)## 'data.frame': 12 obs. of 13 variables:
## $ sample: int 1 2 3 4 5 6 7 8 9 10 ...
## $ tramo : chr "upper" "upper" "upper" "upper" ...
## $ sp1 : int 3 6 3 3 10 12 10 13 10 11 ...
## $ sp2 : int 8 6 9 7 9 14 14 15 20 16 ...
## $ sp3 : int 9 7 10 6 9 11 13 14 13 17 ...
## $ sp4 : int 3 6 1 9 8 16 11 14 21 20 ...
## $ sp5 : int 8 5 7 4 12 15 16 17 11 19 ...
## $ sp6 : int 10 9 8 9 11 13 16 10 13 19 ...
## $ sp7 : int 2 10 9 10 14 8 9 11 10 20 ...
## $ sp8 : int 5 6 4 4 8 15 8 8 10 21 ...
## $ sp9 : int 10 12 10 14 9 10 8 5 5 8 ...
## $ sp10 : int 2 6 6 5 16 8 13 13 15 14 ...
## $ sp11 : int 15 17 20 18 11 11 10 10 1 6 ...Nuestros datos muestran 11 especies de peces (las columnas) y hemos muestreado en 12 sitios (las filas). Parece que algunos peces son más abundantes en algunos tramos que en otros, pero de momento no sabemos muy bien lo que ocurre en nuestra comunidad.
4 Representación gráfica preliminar
Hagamos una figura sencilla que muestre como es la abundancia de cada especie en cada uno de los tramos:
library("tidyr")#Agrupamos las columnas con las especies en una única columna para poder representar la figura
datos2 <- gather(datos,
key = "species",#nombre para la variable resultante de transponer todas las columnas de los diferentes tiempos (t1, t2, t3, etc.)
value = "abundancia",#variable que estamos agrupando y el nombre que tendrá la nueva columna
sp1:sp11)#desde que columna hasta que columna queremos transponer
datos2## sample tramo species abundancia
## 1 1 upper sp1 3
## 2 2 upper sp1 6
## 3 3 upper sp1 3
## 4 4 upper sp1 3
## 5 5 medium sp1 10
## 6 6 medium sp1 12
## 7 7 medium sp1 10
## 8 8 medium sp1 13
## 9 9 lower sp1 10
## 10 10 lower sp1 11
## 11 11 lower sp1 21
## 12 12 lower sp1 11
## 13 1 upper sp2 8
## 14 2 upper sp2 6
## 15 3 upper sp2 9
## 16 4 upper sp2 7
## 17 5 medium sp2 9
## 18 6 medium sp2 14
## 19 7 medium sp2 14
## 20 8 medium sp2 15
## 21 9 lower sp2 20
## 22 10 lower sp2 16
## 23 11 lower sp2 22
## 24 12 lower sp2 14
## 25 1 upper sp3 9
## 26 2 upper sp3 7
## 27 3 upper sp3 10
## 28 4 upper sp3 6
## 29 5 medium sp3 9
## 30 6 medium sp3 11
## 31 7 medium sp3 13
## 32 8 medium sp3 14
## 33 9 lower sp3 13
## 34 10 lower sp3 17
## 35 11 lower sp3 22
## 36 12 lower sp3 21
## 37 1 upper sp4 3
## 38 2 upper sp4 6
## 39 3 upper sp4 1
## 40 4 upper sp4 9
## 41 5 medium sp4 8
## 42 6 medium sp4 16
## 43 7 medium sp4 11
## 44 8 medium sp4 14
## 45 9 lower sp4 21
## 46 10 lower sp4 20
## 47 11 lower sp4 15
## 48 12 lower sp4 17
## 49 1 upper sp5 8
## 50 2 upper sp5 5
## 51 3 upper sp5 7
## 52 4 upper sp5 4
## 53 5 medium sp5 12
## 54 6 medium sp5 15
## 55 7 medium sp5 16
## 56 8 medium sp5 17
## 57 9 lower sp5 11
## 58 10 lower sp5 19
## 59 11 lower sp5 10
## 60 12 lower sp5 20
## 61 1 upper sp6 10
## 62 2 upper sp6 9
## 63 3 upper sp6 8
## 64 4 upper sp6 9
## 65 5 medium sp6 11
## 66 6 medium sp6 13
## 67 7 medium sp6 16
## 68 8 medium sp6 10
## 69 9 lower sp6 13
## 70 10 lower sp6 19
## 71 11 lower sp6 20
## 72 12 lower sp6 22
## 73 1 upper sp7 2
## 74 2 upper sp7 10
## 75 3 upper sp7 9
## 76 4 upper sp7 10
## 77 5 medium sp7 14
## 78 6 medium sp7 8
## 79 7 medium sp7 9
## 80 8 medium sp7 11
## 81 9 lower sp7 10
## 82 10 lower sp7 20
## 83 11 lower sp7 21
## 84 12 lower sp7 19
## 85 1 upper sp8 5
## 86 2 upper sp8 6
## 87 3 upper sp8 4
## 88 4 upper sp8 4
## 89 5 medium sp8 8
## 90 6 medium sp8 15
## 91 7 medium sp8 8
## 92 8 medium sp8 8
## 93 9 lower sp8 10
## 94 10 lower sp8 21
## 95 11 lower sp8 21
## 96 12 lower sp8 15
## 97 1 upper sp9 10
## 98 2 upper sp9 12
## 99 3 upper sp9 10
## 100 4 upper sp9 14
## 101 5 medium sp9 9
## 102 6 medium sp9 10
## 103 7 medium sp9 8
## 104 8 medium sp9 5
## 105 9 lower sp9 5
## 106 10 lower sp9 8
## 107 11 lower sp9 7
## 108 12 lower sp9 9
## 109 1 upper sp10 2
## 110 2 upper sp10 6
## 111 3 upper sp10 6
## 112 4 upper sp10 5
## 113 5 medium sp10 16
## 114 6 medium sp10 8
## 115 7 medium sp10 13
## 116 8 medium sp10 13
## 117 9 lower sp10 15
## 118 10 lower sp10 14
## 119 11 lower sp10 21
## 120 12 lower sp10 11
## 121 1 upper sp11 15
## 122 2 upper sp11 17
## 123 3 upper sp11 20
## 124 4 upper sp11 18
## 125 5 medium sp11 11
## 126 6 medium sp11 11
## 127 7 medium sp11 10
## 128 8 medium sp11 10
## 129 9 lower sp11 1
## 130 10 lower sp11 6
## 131 11 lower sp11 3
## 132 12 lower sp11 4library(ggplot2)
gr1<-ggplot(datos2,aes(species,abundancia))+geom_boxplot(aes(fill=tramo))+
geom_vline(xintercept = c(1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5),linetype = "dashed")+
scale_x_discrete(drop = FALSE)
gr1
Parece que algunas especies son más abundantes en tramos altos (por ejemplo la especie 9) y otras en tramos bajos (por ejemplo la especie 7).
5 ANOSIM
Comparemos el conjuntos de las comunidades en función del tramo estudiado. Sería similar a un ANOVA de un factor, pero en este caso tenemos varias variables (las especies) y las queremos comparar entre tramos (un factor con tres niveles: alto, medio y bajo).
5.1 Generamos una matriz solo con las especies
Para este análisis necesitamos generar una matriz con las abundancias de nuestras 11 especies. Para ello cargamos el paquete vegan y aplicamos la función as.matrix de la siguiente forma:
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-2#Ahora solo necesitamos las abundancias que son de la tercera columna en adelante. Pasamos de un dataframe a una matriz de datos:
matriz <- datos[,3:ncol(datos)]
matriz <- as.matrix(matriz)El aspecto de nuestra matriz es el siguiente:
matriz## sp1 sp2 sp3 sp4 sp5 sp6 sp7 sp8 sp9 sp10 sp11
## [1,] 3 8 9 3 8 10 2 5 10 2 15
## [2,] 6 6 7 6 5 9 10 6 12 6 17
## [3,] 3 9 10 1 7 8 9 4 10 6 20
## [4,] 3 7 6 9 4 9 10 4 14 5 18
## [5,] 10 9 9 8 12 11 14 8 9 16 11
## [6,] 12 14 11 16 15 13 8 15 10 8 11
## [7,] 10 14 13 11 16 16 9 8 8 13 10
## [8,] 13 15 14 14 17 10 11 8 5 13 10
## [9,] 10 20 13 21 11 13 10 10 5 15 1
## [10,] 11 16 17 20 19 19 20 21 8 14 6
## [11,] 21 22 22 15 10 20 21 21 7 21 3
## [12,] 11 14 21 17 20 22 19 15 9 11 45.2 Aplicamos ANOSIM
Aplicamos la función anosim a nuestra matriz de datos en función del factor ‘tramo’ de la siguiente forma:
analisis<-anosim(matriz, datos$tramo, distance="bray",permutations=9999)#hay que indicarle el factor que es el tramo
analisis##
## Call:
## anosim(x = matriz, grouping = datos$tramo, permutations = 9999, distance = "bray")
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.8796
## Significance: 3e-04
##
## Permutation: free
## Number of permutations: 9999Nuestro estadístico R muestra una gran disimilitud (0.8796) (cercano a 1) y además hay diferencias significativas entre las comunidades de los tramos estudiados.
5.3 Influencia de las especies
Para analizar que especies difieren entre los tramos aplicamos el ‘Indicator species analysis’. Se requiere un vector indicando el factor (en este caso tramo) y se aplica la función multipatt de la siguiente forma:
library(indicspecies)
tramo<-datos$tramo#vector con el tramo
tramo## [1] "upper" "upper" "upper" "upper" "medium" "medium" "medium" "medium"
## [9] "lower" "lower" "lower" "lower"#analizamos las especies:
speciesdiferencias<-multipatt(matriz,tramo,func="r.g",control = how(nperm=9999))
speciesdiferencias## $call
## multipatt(x = matriz, cluster = tramo, func = "r.g", control = how(nperm = 9999))
##
## $func
## [1] "r.g"
##
## $cluster
## [1] "upper" "upper" "upper" "upper" "medium" "medium" "medium" "medium"
## [9] "lower" "lower" "lower" "lower"
##
## $comb
## lower medium upper lower+medium lower+upper medium+upper
## [1,] 0 0 1 0 1 1
## [2,] 0 0 1 0 1 1
## [3,] 0 0 1 0 1 1
## [4,] 0 0 1 0 1 1
## [5,] 0 1 0 1 0 1
## [6,] 0 1 0 1 0 1
## [7,] 0 1 0 1 0 1
## [8,] 0 1 0 1 0 1
## [9,] 1 0 0 1 1 0
## [10,] 1 0 0 1 1 0
## [11,] 1 0 0 1 1 0
## [12,] 1 0 0 1 1 0
##
## $str
## lower medium upper lower+medium lower+upper medium+upper
## sp1 0.5467673 0.26149742 -0.8082647 0.8082647 -0.26149742 -0.5467673
## sp2 0.7461760 0.02407019 -0.7702462 0.7702462 -0.02407019 -0.7461760
## sp3 0.8021795 -0.13170111 -0.6704784 0.6704784 0.13170111 -0.8021795
## sp4 0.7405434 0.05696488 -0.7975083 0.7975083 -0.05696488 -0.7405434
## sp5 0.4095142 0.40951418 -0.8190284 0.8190284 -0.40951418 -0.4095142
## sp6 0.7961815 -0.12841637 -0.6677651 0.6677651 0.12841637 -0.7961815
## sp7 0.7363971 -0.18684703 -0.5495501 0.5495501 0.18684703 -0.7363971
## sp8 0.7616010 -0.08016853 -0.6814325 0.6814325 0.08016853 -0.7616010
## sp9 -0.4716666 -0.25941665 0.7310833 -0.7310833 0.25941665 0.4716666
## sp10 0.5922620 0.22349508 -0.8157570 0.8157570 -0.22349508 -0.5922620
## sp11 -0.8376578 0.00000000 0.8376578 -0.8376578 0.00000000 0.8376578
##
## $A
## NULL
##
## $B
## NULL
##
## $sign
## s.lower s.medium s.upper index stat p.value
## sp1 1 1 0 4 0.8082647 0.0056
## sp2 1 1 0 4 0.7702462 0.0167
## sp3 1 0 0 1 0.8021795 0.0176
## sp4 1 1 0 4 0.7975083 0.0120
## sp5 1 1 0 4 0.8190284 0.0104
## sp6 1 0 0 1 0.7961815 0.0165
## sp7 1 0 0 1 0.7363971 0.0316
## sp8 1 0 0 1 0.7616010 0.0168
## sp9 0 0 1 3 0.7310833 0.0275
## sp10 1 1 0 4 0.8157570 0.0056
## sp11 0 0 1 3 0.8376578 0.0108
##
## attr(,"class")
## [1] "multipatt"summary(speciesdiferencias)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: r.g
## Significance level (alpha): 0.05
##
## Total number of species: 11
## Selected number of species: 11
## Number of species associated to 1 group: 6
## Number of species associated to 2 groups: 5
##
## List of species associated to each combination:
##
## Group lower #sps. 4
## stat p.value
## sp3 0.802 0.0176 *
## sp6 0.796 0.0165 *
## sp8 0.762 0.0168 *
## sp7 0.736 0.0316 *
##
## Group upper #sps. 2
## stat p.value
## sp11 0.838 0.0108 *
## sp9 0.731 0.0275 *
##
## Group lower+medium #sps. 5
## stat p.value
## sp5 0.819 0.0104 *
## sp10 0.816 0.0056 **
## sp1 0.808 0.0056 **
## sp4 0.798 0.0120 *
## sp2 0.770 0.0167 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Si nos fijamos en el último resultado (Multilevel pattern analysis) vemos que hay una serie de especies asociadas a los tramos. Por ejmeplo, las species 3, 6, 8 y 7) se encuentran asociadas al tramo bajo. El orden en que se muestra en en su grado de afinidad por ese tramo, vemos que se ordena de mayor a menor valor del estadístico. Recordemos nuestra figura del principio para ver si esto es correcto:
Parece que efectivamente esas especies son más abundantes en el tramo bajo. Mirar el resto de especies y sus afinidades por los tramos. Hay algunas que tienen afinidad por dos tramos.
6 NMDS Escalamiento multidimensional no métrico
Se trata de un análisis multivariante que nos permite representar gráficamente una series de puntos de muestreo en función de la similitud o disimilitud de la composición de las comunidades de esos puntos de muestreo. Cuanto más similares sean las comunidades más cercanos aparecen los puntos en el espacio, cuanto más diferentes más alejados. Es una forma relativamente sencilla de comprobar qué tramos de nuestro estudio se parecen más entre si en función de la abundancia de las diferentes especies de peces de nuestra comunidad.
6.1 Realizamos el análisis NMDS
Aplicamos la función metaMDS del paquete simED
library(simEd)## Loading required package: rstream##
## Attaching package: 'simEd'## The following objects are masked from 'package:base':
##
## sample, set.seedmatriz## sp1 sp2 sp3 sp4 sp5 sp6 sp7 sp8 sp9 sp10 sp11
## [1,] 3 8 9 3 8 10 2 5 10 2 15
## [2,] 6 6 7 6 5 9 10 6 12 6 17
## [3,] 3 9 10 1 7 8 9 4 10 6 20
## [4,] 3 7 6 9 4 9 10 4 14 5 18
## [5,] 10 9 9 8 12 11 14 8 9 16 11
## [6,] 12 14 11 16 15 13 8 15 10 8 11
## [7,] 10 14 13 11 16 16 9 8 8 13 10
## [8,] 13 15 14 14 17 10 11 8 5 13 10
## [9,] 10 20 13 21 11 13 10 10 5 15 1
## [10,] 11 16 17 20 19 19 20 21 8 14 6
## [11,] 21 22 22 15 10 20 21 21 7 21 3
## [12,] 11 14 21 17 20 22 19 15 9 11 4set.seed(123)#saca siempre el mismo resultado del NMDS para estos datos
nmds<-metaMDS(matriz,distance="bray")## Wisconsin double standardization
## Run 0 stress 0.04589633
## Run 1 stress 0.04949178
## Run 2 stress 0.04589633
## ... Procrustes: rmse 5.783155e-05 max resid 0.0001091742
## ... Similar to previous best
## Run 3 stress 0.04589634
## ... Procrustes: rmse 4.768201e-05 max resid 8.352725e-05
## ... Similar to previous best
## Run 4 stress 0.04589632
## ... New best solution
## ... Procrustes: rmse 1.963139e-05 max resid 3.48667e-05
## ... Similar to previous best
## Run 5 stress 0.04513076
## ... New best solution
## ... Procrustes: rmse 0.1063576 max resid 0.2214685
## Run 6 stress 0.04513068
## ... New best solution
## ... Procrustes: rmse 0.0001012067 max resid 0.000290964
## ... Similar to previous best
## Run 7 stress 0.04513072
## ... Procrustes: rmse 5.774534e-05 max resid 0.0001623201
## ... Similar to previous best
## Run 8 stress 0.04817096
## Run 9 stress 0.05741258
## Run 10 stress 0.04949171
## Run 11 stress 0.04589632
## Run 12 stress 0.04366891
## ... New best solution
## ... Procrustes: rmse 0.02278309 max resid 0.06090961
## Run 13 stress 0.04589633
## Run 14 stress 0.04589632
## Run 15 stress 0.05741267
## Run 16 stress 0.04817111
## Run 17 stress 0.04949185
## Run 18 stress 0.05741257
## Run 19 stress 0.04817102
## Run 20 stress 0.05741279
## *** No convergence -- monoMDS stopping criteria:
## 20: stress ratio > sratmaxnmds##
## Call:
## metaMDS(comm = matriz, distance = "bray")
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(matriz)
## Distance: bray
##
## Dimensions: 2
## Stress: 0.04366891
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(matriz)'plot(nmds) 
Las cruces rojas son las especies y los puntos negros los sitios. La figura no nos dice gran cosa. Vamos a mejorarla para poder entender lo que ocurre con nuestras comunidades.
6.2 Mejoramos la figura
Ahora sobre las coordenadas de nuestra figura anterior añadimos el tramo y la muestra:
coordenadas <- as.data.frame(scores(nmds)$sites)
coordenadas## NMDS1 NMDS2
## 1 -0.30513956 0.13099014
## 2 -0.22617762 -0.06176535
## 3 -0.31684082 0.02027648
## 4 -0.32920636 -0.08700135
## 5 0.01579458 -0.07815066
## 6 0.06277971 -0.02197830
## 7 0.07610254 0.04309527
## 8 0.15617844 -0.04578167
## 9 0.27144682 -0.07093430
## 10 0.17265815 0.04286058
## 11 0.27863944 0.02717179
## 12 0.14376468 0.10121737#Le añadimos a las coordenadas una columna con los tramos:
coordenadas$tramo = datos$tramo
coordenadas$sample = datos$sample
head(coordenadas)## NMDS1 NMDS2 tramo sample
## 1 -0.30513956 0.13099014 upper 1
## 2 -0.22617762 -0.06176535 upper 2
## 3 -0.31684082 0.02027648 upper 3
## 4 -0.32920636 -0.08700135 upper 4
## 5 0.01579458 -0.07815066 medium 5
## 6 0.06277971 -0.02197830 medium 6str(coordenadas)## 'data.frame': 12 obs. of 4 variables:
## $ NMDS1 : num -0.3051 -0.2262 -0.3168 -0.3292 0.0158 ...
## $ NMDS2 : num 0.131 -0.0618 0.0203 -0.087 -0.0782 ...
## $ tramo : chr "upper" "upper" "upper" "upper" ...
## $ sample: int 1 2 3 4 5 6 7 8 9 10 ...Y ahora hacemos el gráfico:
gr2<- ggplot(coordenadas, aes(x = NMDS1, y = NMDS2))+
geom_point(size = 4, aes( shape = tramo, colour = tramo))+
geom_text(hjust=0.5, vjust=1.5, label=datos$sample)
gr2
Ahora cada tramo se representa por un color y forma. En cada punto de muestreo aparece el número de ese sitio. Se puede ver como los tramos altos aparecen separados del resto y las diferencias entre el tramo medio y bajo son menores que respecto al tramo alto que hemos estudiado. Con esta figura y el análisis ANOSIM ya podemos sacar bastantes conclusiones sobre cómo son nuestras comunidades de diferentes (o de similares).
7 CRÉDITOS
Departamento de Ciencias de la Vida
Universidad de Alcalá (España)
Citation:
Alonso A (2022) Análisis de comunidades: ANOSIM y NMDS. Rpubs, accesed FECHA, https://rpubs.com/aafernandez1976/ANOSIMyNMDS
Truchas en un tramo de cabecera de un río