Problemas de convergencia en el escalamiento multidimensional no métrico (NMDS) ¿Qué hacer?

El escalamiento multidimensional no métrico (NMS, MDS, NMDS o NMMDS) es una técnica multivariante de interdependencia que trata de representar en un espacio geométrico de pocas dimensiones las proximidades existentes entre un conjunto de objetos. El NMDS es un método de ordenación adecuado para datos que no son normales o que están en una escala discontinua o arbitraria. Una ventaja del NMDS frente a otras técnicas de ordenación es que, al estar basada en rangos de distancias, tiende a linealizar la relación entre las distancias ambientales y las distancias biológicas (esto es, calculadas a partir de una matriz de sitios x especies). Una de las desventajas de esta técnica es la di cultad para alcanzar una solución estable única. A pesar de ello, el NMDS es una técnica ampliamente utilizada en ecología para detectar gradientes en comunidades biológicas.

El NMDS se implementa de la siguiente forma:

1. Se calcula la matriz de disimilaridad X a partir de la matriz de datos de sitios x especies. Esta matriz nos indica cómo de iguales son cada par de sitios utilizando para ello la similaridad entre sus especies. Supongamos que tenemos tres especies (sp1, sp2, sp3) y tres sitios (A, B, C). El sitio A tiene sp1 = 3, sp2 = 0 y sp3 = 8. El sitio B tiene sp1 = 3, sp2 = 0 y sp3 = 6. El sitio C tiene sp1 = 0, sp2 = 5 y sp3 = 1. Por tanto, podemos calcular una matriz de disimilaridad que nos indique con números que los sitios A y B son muy iguales, mientras que los sitios A y C y B y C son muy distintos entre sí. Cuando se trata de datos biológicos la distancia más usada es la distancia de Sorensen (Bray-Curtis) en vez de la distancia Euclídea.
2. Se asignan los sitios (unidades muestrales) a una con guración inicial aleatoria en un espacio k-dimensional (dónde k es el número de especies), aunque en realidad, la ordenación se va a realizar principalmente sobre unas pocas dimensiones (2 o 3).
3. Se calculan las distancias sobre este nuevo espacio geométrico y se calcula una matriz de distancia Y .
4. Se comparan las matrices de distancia X e Y y se mide cómo son de parecidas entre ellas (stress).
5. A partir de la con guración inicial, se reasignan los sitios (unidades muestrales) para reducir las distancias con la matriz X.
6. Se repite este proceso de manera iterativa hasta que se consigue una solución óptima en dónde la matriz de distancias Y es muy parecida a la matriz de distancias X. Esto es, se minimiza el stress.

En R tenemos una implementación de esta función (metaMDS) en el paquete vegan.

Problemas de convergencia

Uno de los problemas que pueden surgir a la hora de realizar un NMDS en R es el de la convergencia a la hora de encontrar una solución óptima. Esto ocurre porque existen dos sitios (o más) que tienen exactamente la misma composición de especies. Es decir, todos los sitios tienen que tener una composición de especies única (puede ser muy parecida pero no exactamente iguales). Todavía no se muy bien si esto es una premisa del NMDS como tal o simplemente un resultado de su implementación en R. Lo que está claro es que el NMDS está basado en distancias y puede darse el caso de que la distancia entre dos sitios sea 0 (independientemente del método que utilicemos para calcular dichas distancias), lo que implicaría que su composición de especies sea exactamente igual. Por lo tanto no es un problema de cálculo de distancias sino del algoritmo que realiza el NMDS.

Aunque desconozco todavía el motivo de este error, existe una forma de solucionar el problema. La idea es sencilla. Consiste en incluir un pequeño ruido aleatorio en los números que forman la matriz de datos. Este error debe de ser lo suficientemente pequeño como para no alterar las distancias multidimensionales entre pares de sitios, pero conseguirán que los valores de composición no sean exactamente iguales. Y aunque al tratarse de abundancias de especies los datos deberían de ser enteros, no supone un problema análitico el meter un ruido aleatorio que sea de tipo decimal. Veámoslo con un ejemplo.

Creamos una matriz de datos de abundancia con 10 sitios y 10 especies:

m <- matrix(rpois(100, 1), nrow=10, ncol=10)
colnames(m)<- paste("sp", c(1:10), sep="")
rownames(m) <- paste("sitio", c(1:10), sep="_")
m

sp1 sp2 sp3 sp4 sp5 sp6 sp7 sp8 sp9 sp10
sitio_1 1 1 0 4 4 2 1 1 1 0
sitio_2 2 1 2 1 0 3 1 0 1 0
sitio_3 0 0 2 1 0 1 0 0 0 1
sitio_4 2 1 2 0 2 1 1 0 0 0
sitio_5 0 1 0 1 0 0 0 1 0 0
sitio_6 1 1 0 0 0 1 2 1 1 0
sitio_7 2 1 1 1 0 2 0 2 1 0
sitio_8 1 0 3 1 0 2 2 1 1 0
sitio_9 0 1 1 2 1 3 0 1 0 0
sitio_10 0 2 6 1 1 2 0 0 1 1

La analizamos en R con la función metaMDS del paquete vegan.

library(vegan)
set.seed(0)
nmds1 <- metaMDS(m)
plot(nmds1)

Como vemos, todo va bien. Ahora vamos a crear una nueva fila (sitio) en la matriz m que tendrá la misma abundancia que el sitio_1.

m2 <- rbind(m, m[1,])
rownames(m2) <- paste("sitio", c(1:11), sep="_")
m2

sp1 sp2 sp3 sp4 sp5 sp6 sp7 sp8 sp9 sp10
sitio_1 1 1 0 4 4 2 1 1 1 0
sitio_2 2 1 2 1 0 3 1 0 1 0
sitio_3 0 0 2 1 0 1 0 0 0 1
sitio_4 2 1 2 0 2 1 1 0 0 0
sitio_5 0 1 0 1 0 0 0 1 0 0
sitio_6 1 1 0 0 0 1 2 1 1 0
sitio_7 2 1 1 1 0 2 0 2 1 0
sitio_8 1 0 3 1 0 2 2 1 1 0
sitio_9 0 1 1 2 1 3 0 1 0 0
sitio_10 0 2 6 1 1 2 0 0 1 1
sitio_11 1 1 0 4 4 2 1 1 1 0

Y analizamos esta nueva matriz de datos con la función metaMDS.

set.seed(0)
nmds2 <- metaMDS(m2)

Error in metaMDSdist(comm, distance = distance, autotransform = autotransform, :
Zero dissimilarities are not allowed

Resolvemos el problema. Vamos a meter ruido únicamente en los valores > 0.

m3 <- m2
for(i in 1:dim(m2)[2]){
m3[, i] <- ifelse(m2[,i] > 0, jitter(m2[,i],
factor=1e-05), m2[,i])
}

Y ahora analizamos la matriz m3 con la función metaMDS().

set.seed(0)
nmds3 <- metaMDS(m3)
plot(nmds3)

El único problema a este enfoque es la replicabilidad exacta de los resultados. Al utilizar la función set.seed(0) hacemos que la configuración inicial del NMDS sea siempre la misma. Sin embargo, al meter un ruido aleatorio con la función jitter() cada vez que corremos el código vamos a crear matrices ligeramente distintas y el gráfico de ordenación que obtengamos será distinto cada vez, aunque la posición relativa de los puntos va a ser muy parecida siempre. Si queremos evitar este problema podríamos salvar la matriz m3 como un archivo de texto y leer esos mismos datos cada vez que vayamos a repetir el análisis. Esto solucionaría el problema.

Si alguien averigua por qué ocurre este problema de convergencia en R que por favor me lo diga. Se puede encontrar más información sobre análisis multivariante aquí.