14.3 Escalonamento Não-Métrico Multimensional (NMDS)

A NMDS é uma técnica de ordenação multivariada que permite visualizar graficamente distâncias entre objetos. No R há várias funções que executam isso: isoMDS(); cmdscale(); e nmds() e bestnmds() do pacote labdsv (Roberts 2019). Vamos usar a função bestnmds() nos exemplos abaixo. Para entender, veja um exemplo para distâncias geográficas entre cidades na região norte do Brasil. Se queremos representar graficamente distâncias geográficas, estaremos de certa forma reproduzindo um mapa:

Vamos utilizar o conjunto de dados contendo coordenadas geográficas de municípios do Brasil para esta prática, utilizado na seção 3.4.2.

# visualizando distancias usando NMDS
# Vamos usar o arquivo com coordenadas dos municipios brasileiros
muni <- read.table(file = "municipiosbrasil.csv", header = T, as.is = T, sep = "\t", na.strings = c("NA", "", "NULL"))

Para simplificar, vamos filtrar apenas algumas cidades da região Norte:

cids <- c("Rio Branco", "Cruzeiro do Sul", "Tabatinga", "São Gabriel da Cachoeira", "Manaus", "Santarém", "Porto Velho", "Humaitá", "Belém", "Macapá", "Marabá", "Boa Vista")

# filtrando os dados
vl <- muni$Municipio %in% cids
muni <- muni[vl, ]

# calcula a distancia geografica entre essas cidades (em graus de latitude).
# Idealmente deveríamos converter latitude e longitude em décimos de graus para UTM para obter distancias em km ou m.
mdist <- dist(muni[, c("Longitude", "Latitude")], method = "euclidean")

# calculando um nmds em dois eixos (reduzindo a variação na matriz em dois eixos)
?bestnmds

Vamos utilizar a função bestnmds do pacote labdsv Vamos agora instalar o pacote e carregar o pacote labdsv:

onmds <- bestnmds(mdist, k = 2)

##  [1] 3.557932e-13 8.167216e-03 2.228642e+01 3.791436e+01 6.603351e-03
##  [6] 6.029306e-03 2.109747e+01 9.041554e-03 4.415954e-03 4.861851e-01
## [11] 4.001459e-03 8.141675e-03 9.111469e-03 3.737914e+01 6.265940e-03
## [16] 4.124426e+01 9.818992e+00 2.702450e+01 3.316186e-03 3.712228e+01
## 
## best result = 1
## with stress = 3.557932e-13

# veja a estrutura do resultado
str(onmds)

## List of 3
##  $ points: num [1:21, 1:2] -18.2 -12.9 17.3 -14.7 -16.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:21] "3" "14" "74" "135" ...
##   .. ..$ : NULL
##  $ stress: num 3.56e-13
##  $ type  : chr "NMDS"
##  - attr(*, "class")= chr [1:2] "dsvord" "nmds"
##  - attr(*, "call")= language bestnmds(dis = mdist, k = 2)
##  - attr(*, "timestamp")= chr "Fri Oct  6 18:30:20 2023"

# o valor do stress indica o ajuste. Se o stress for 0, o ajuste é perfeito:
# a posição dos pontos é proporcional a distância entre eles
onmds$stress

## [1] 3.557932e-13

# quanto da variação foi explicada?
# essa pergunta a gente faz com PCA, não faz sentido fazer com NMDS.
# Mas, se você quer ter uma ideia de quanto um eixo capturou da variação,
# pode correlacionar a matriz de distancia original com uma matriz de distancia
# gerada pelos valores dos eixos nmds

# pega os valores dos eixos NMDS
ptsnmds <- onmds$points
# calcula a distancia
adist <- dist(ptsnmds)
# qual a correlacao entre essas matrizes de distancia?
cor(mdist, adist)

## [1] 1

# por isso o stress é baixo

# vamos comparar graficamente:
# divide o dispositivo em duas partes
par(mfrow = c(2, 1), mar = c(5, 5, 1, 1))
# adiciona limite no eixo X e y
xl <- range(ptsnmds[, 1]) + c(-1, 10)
yl <- range(ptsnmds[, 2]) + c(-2, 2)
# plota pontos
plot(ptsnmds, type = "p", pch = 21, bg = "red", xlab = "NMDS 1", ylab = "NMDS 2", xlim = xl, ylim = yl, cex = 0.5)
# adiciona o nome das cidades
text(ptsnmds, labels = muni$Municipio, cex = 0.8, pos = 4)

xl <- range(muni$Longitude) + c(-2, 5)
yl <- range(muni$Latitude) + c(-2, 2)
plot(muni$Longitude, muni$Latitude, type = "p", xlab = "Longitude", ylab = "Latitude", xlim = xl, ylim = yl, pch = 21, cex = 0.5, col = "blue")
text(muni$Longitude, muni$Latitude, labels = muni$Municipio, cex = 0.8, pos = 4)

14.3.1 Exemplo com dados morfológicos

Um exemplo de NMDS para mostra a similaridade entre indivíduos de Iris a partir de uma matriz de distância morfológica usando os dados de iris do R.

data(iris) # #carrega o conjunto de dados de iris para a area de trabalho
head(iris)

Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
5.1	3.5	1.4	0.2	setosa
4.9	3.0	1.4	0.2	setosa
4.7	3.2	1.3	0.2	setosa
4.6	3.1	1.5	0.2	setosa
5.0	3.6	1.4	0.2	setosa
5.4	3.9	1.7	0.4	setosa

# calcula a distancia morfológica
?dist

dmorfo <- dist(iris[, 1:4], method = "eucl")

# calculando um nmds em dois eixos
onmds <- bestnmds(dmorfo, k = 2)

# ops objetos 102 e 143 devem ser identificos
iris[c(102, 143), ]

	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
102	5.8	2.7	5.1	1.9	virginica
143	5.8	2.7	5.1	1.9	virginica

# entao eliminamos 1, porque senao nao funciona
iris2 <- iris[-102, ]

# calculamos novamente a distancia
dmorfo <- dist(iris2[, 1:4], method = "eucl")
# calculando um nmds em dois eixos
onmds <- bestnmds(dmorfo, k = 2)

# veja a estrutura do resultado
str(onmds)

## List of 3
##  $ points: num [1:149, 1:2] -2.67 -2.71 -2.87 -2.74 -2.72 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:149] "1" "2" "3" "4" ...
##   .. ..$ : NULL
##  $ stress: num 2.58
##  $ type  : chr "NMDS"
##  - attr(*, "class")= chr [1:2] "dsvord" "nmds"
##  - attr(*, "call")= language bestnmds(dis = dmorfo, k = 2)
##  - attr(*, "timestamp")= chr "Fri Feb  5 10:50:06 2021"

# o valor do stress indica o ajuste. Se stress for 0, o ajuste é perfeito
onmds$stress

## [1] 2.579978

# quanto da variação foi explicada?
# essa pergunta a gente faz com PCA, não faz sentido fazer com NMDS.
# Mas, se você quer ter uma ideia de quanto um eixo capturou da variação,
# pode correlacionar a matriz de distancia original
# com uma matriz de distancia gerada pelos valores dos eixos nmds

# pega os valores dos eixos NMDS
ptsnmds <- onmds$points
# cada linha nessa tabela corresponde
# à mesma linha na tabela iris2
head(ptsnmds)

-2.673344	0.2751282
-2.705573	-0.1645018
-2.866471	-0.0970515
-2.737119	-0.2591412
-2.717433	0.2980718
-2.296804	0.6687755

# essas colunas são os dois eixos NMDS

# isso deve ser verdadeiro
nrow(ptsnmds) == nrow(iris2)

## [1] TRUE

# calcula a distancia entre os pontos pelos eixos NMDS
adist <- dist(ptsnmds)
# qual a correlacao entre essa matriz e a original?
cor(dmorfo, adist)

## [1] 0.9988093

# por isso o stress é baixo, a correlação é alta

Agora vamos fazer com outro índice de distância chamado gower, por exemplo, que é um bom índice quando se tem dados semiquantitativos na matriz (não é o caso aqui). Vamos utilizar a função vegdist() do pacote vegan (Oksanen et al. 2020).

data(iris)
iris2 <- iris[-102, ]

dmorfo2 <- vegdist(iris2[, 1:4], method = "gower")
onmds2 <- bestnmds(dmorfo2, k = 2)

# pega os valores dos eixos NMDS
ptsnmds2 <- onmds2$points
# cada linha nessa tabela corresponde
# à mesma linha na tabela iris2
head(ptsnmds2)

-0.2907098	0.0456498
-0.2712303	-0.0244397
-0.3013570	-0.0090738
-0.2933819	-0.0245441
-0.3066071	0.0494793
-0.2595641	0.1165398

# essas colunas são os dois eixos NMDS
# isso deve ser verdadeiro
nrow(ptsnmds2) == nrow(iris2)

## [1] TRUE

# calcula a distancia entre os pontos pelos eixos NMDS
adist2 <- dist(ptsnmds2)
# qual a correlacao entre essa matriz e a original?
cor(dmorfo2, adist2)

## [1] 0.9950709

# por isso o stress é baixo, a correlação é alta

# vamos visualizar os dois resultados, onmds e onmds2, graficamente
# divide o dispositivo em dois
par(mfrow = c(2, 1), mar = c(5, 5, 1, 1))
cores <- c("red", "green", "blue")[as.numeric(iris2$Species)]
plot(ptsnmds, pch = 21, bg = cores, cex = 0.8, xlab = "NMDS 1", ylab = "NMDS 2", main = "Objeto: ptsnmds")
plot(ptsnmds2, pch = 21, bg = cores, cex = 0.8, xlab = "NMDS 1", ylab = "NMDS 2", main = "Objeto: ptsnmds2")

Referências

Oksanen, Jari, F. Guillaume Blanchet, Michael Friendly, Roeland Kindt, Pierre Legendre, Dan McGlinn, Peter R. Minchin, et al. 2020. vegan: Community Ecology Package. https://CRAN.R-project.org/package=vegan.

Roberts, David W. 2019. labdsv: Ordination and Multivariate Analysis for Ecology. http://ecology.msu.montana.edu/labdsv/R.