Como comparar duas ordenações de nmds

Veremos a partir de agora como as diferenças na composição de espécies pode ser investigada através do cálculo das distâncias ecológicas (ou similaridades) entre dois locais. Os métodos podem ser aplicados à matriz inteira de espécies, resultando então em uma outra matriz de distância (ou similaridade) entre os pares de locais. Essa matriz de distância é a base para as análises multivariadas de ordenação que veremos em sequência, o PCA e o NMDS.

É importante ressaltar que as análises multivariadas podem ser divididas em dois tipos: análise de agrupamento e ordenação. Análises de agrupamento em geral tentam agrupar objetos (observações) em grupos de maneira que objetos do mesmo grupo sejam mais semelhantes entre si do que objetos de outros grupos. Por outro lado, a análise de ordenação é uma operação pela qual os objetos são posicionados num espaço que contém menos dimensões que o conjunto de dados original; a posição dos objetos em relação aos outros também podem ser usadas para agrupá-los. Não veremos as análises de agrupamento neste roteiro, mas as medidas de distância ecológica tratadas abaixo também se aplicam a estas análises.

Uma boa medida de distância ecológica entre comunidades (locais) descreve a diferença na composição de espécies de forma que locais que compartilhem muitas espécies devem ter distância ecológica pequena. Existe na literatura uma grande variedade de métodos para calcular estas distâncias, apresentamos aqui duas das mais utilizadas para dados ecológicos.

A distância euclidiana aplicada à matriz de espécies é calculada usando cada espécie como um eixo diferente em um gráfico (número de dimensões igual à riqueza total). Os locais são plotados neste gráfico multidimensional para que as distâncias entre os locais sejam calculadas.

Vejamos o seguinte exemplo:

# matriz de espécies por locais mat <- matrix(c(1,5,0,1,5,0,0,0,1), nrow=3, dimnames = list(c("A","B","C"), c("sp1","sp2","sp3"))) mat #matriz de distância entre os locais dist(mat, method = "euclidean") #usando pacote vegan que contém mais métodos de distância vegdist(mat, method="euclidean")

Vamos plotar o gráfico 3D dos dados com as espécies nos eixos para ver as distâncias euclidianas entre os locais:

#pacote usado install.packages("plot3D") library(plot3D) scatter3D(x = mat[,1], y = mat[,2], z = mat[,3], xlim = c(0, 6), ylim = c(0, 6), zlim = c(0, 4), xlab = "sp1", ylab = "sp2", zlab = "sp3", colkey = F, ticktype = "detailed", phi = 10, pch = 16, cex=2) text3D(x = mat[ , 1] + 0.2, y = mat[ , 2] + 0.2, z = mat[ , 3] + 0.3, labels = c("A", "B", "C"), add = T) lines3D(mat [ , 1], mat[ , 2], mat[ , 3], add = T, colkey = F) lines3D(mat[-2,1], mat[-2, 2], mat[-2, 3], add = T, colkey = F)

Vemos aqui que a distância euclidiana dos dados brutos não foi capaz de descrever as distâncias ecológicas muito bem. Vemos que os locais A e B tem as mesmas espécies (mas com abundâncias diferentes), enquanto que o local C tem uma espécie diferente. A distância euclidiana depende muito das abundâncias em de cada espécie e não apenas das espécies que são compartilhadas. A e B, que compartilham as mesmas espécies, estão distantes pela distância euclidiana porque o local A tem 2 indivíduos e o local B tem 10.

Embora a distância euclidiana não seja tão boa para investigar como as espécies são compartilhadas entre os locais, ela ainda é usada em algumas técnicas de ordenação e agrupamento (p.ex. PCA). Uma das soluções para dar maior peso às diferenças na composição de espécies é calcular a proporção de espécies em cada local, ou seja, dividindo as abundâncias das espécies pela abundância total do local (cada linha). Outro método seria usar a matriz de presença-ausência. Veremos na PCA uma outra transformação mais apropriada para dados de abundância (transformação de Hellinger)

#abundância total por local abund.t <- apply(mat, 1, sum) #matriz transformada para abundancias relativas mat2 <- mat /abund.t mat2 #distância euclidiana dist(mat2)

Com os dados transformados vemos que a distância entre o local A e B agora é zero, pois os locais tem as mesmas espécies em igual proporção.

A distância de Bray-Curtis é muito utilizada em dados ecológicos e seus valores estão restritos entre 0 e 1. Quando a distância é zero, os dois locais são completamente similar. Quando é 1 eles são totalmente dissimilares, significando que os locais não possuem nenhuma espécie em comum.

Vejamos as distância pra nossa matriz:

vegdist(mat, method = 'bray')

Podemos ver agora que a distância entre o local C e os locais A e B é 1, pois este local não compartilha nenhuma espécie com os demais locais.

Existem muitos outro métodos para calcular a distância ecológica entre o comunidades e também formas de transformar os dados antes de calcular as distâncias, mas este não é o foco deste roteiro. Veja mais alguns exemplos no Capítulo 8 do manual do pacote BiodiversityR.

Os métodos de ordenação arranjam as comunidades (locais) de forma que as distâncias entre elas no gráfico representem suas distâncias ecológicas. O resultado da ordenação é geralmente visualizado como um gráfico de duas dimensões. Nestes gráficos, cada local é apresentado e os locais próximos são interpetados como tendo composição de espécies similar, enquanto que locais afastados contém composição de espécies diferente. Técnicas de ordenação são particularmente bem adaptadas para analisar dados de comunidades ecológicas, que estão geralmente estruturadas em gradientes (muitas vezes não percebidos pelos pesquisadores).

As análises de ordenação sempre seguem os mesmos passos:

  1. Os dados podem ser transformados e o tipo de transformação depende da questão formulada e análise a ser realizada.

  2. É construída uma matriz de associação das distâncias entre os objetos (locais no nosso caso).

  3. O “programa” então arranja os objetos ao longo de um ou mais eixos (usualmente dois eixos ou duas dimensões) que melhor refletem os padrões encontrados nos dados. Estes eixos refletirão as variáveis ecológicas (ou gradientes) que causaram os padrões nos dados.

Vamos descrever a seguir duas das mais utilizadas análises de ordenação. Exemplificaremos o Escalonamento Multidimensional Não Paramétrico (NMDS) e a Análise de Componentes Principais, que são métodos não-restritivos (‘unconstrained’ em inglês) que usa as distâncias ecológicas baseada apenas na matriz de espécies (ou apenas na matriz de variáveis ambientais para PCA com dados ambientais).

O escalonamento multidimensional não paramétrico (NMDS) pode ser considerado uma análise indireta de gradiente, ou seja, útil para tentar descobrir se há padrões nos dados.

Se a prioridade do pesquisador é representar tão bem quanto possível a relação ordenada entre os objetos (locais) e um número pequeno de eixos, o NMDS pode ser a solução. O NMDS pode produzir a ordenação de objetos de qualquer matriz de distâncias (ou seja usando qualquer medida de distância).

Comecemos com um exemplo, temos uma tabela de dados com as abundâncias das espécies em 12 locais. Estes dados foram criados para que haja uma relação entre as espécies e a precipitação de cada local propositalmente. Entretanto, vamos fingir que não sabemos da existência desse gradiente de precipitação que está estruturando as comunidades e ver se o NMDS é capaz de capturar este padrão nos dados.

# dados fictícios tab <- matrix(c(9,0,0,0,4,0,8,6,2,0,0,10, 1,1,0,2,0,1,0,0,0,2,0,0, 0,0,0,1,0,0,0,0,2,1,4,0, 0,0,0,1,0,0,0,0,2,1,4,0, 0,0,4,0,6,0,7,3,3,0,1,1, 0,1,0,2,0,0,0,0,1,1,0,0), nrow=12) colnames(tab) <- c("spA", "spB", "spC", "spD", "spE","spF") tab # dados de precipitação prec <- c(270, 315, 200, 255, 190, 290, 150, 125, 230, 290, 240, 100)

O objetivo de nossa nálise é descrever o padrão apresentado pelas seis espécies, em menos dimensões do que a tabela de dados (lógico!). Lembremos que cada espécie representa uma dimensão. Como sabemos que existe um gradiente de chuva, vamos fazer uma ordenação dos locais (objetos) em relação às espécies (atributos).

Vamos utilizar a matriz de distância de bray-curtis como medida de associação entre os locais. Esta nos dará a informação de quão distante cada local está dos outros em relação à composição de espécies.

distmatrix <- vegdist(tab, method = 'bray')

Vamos agora usar a análise de NMDS para ordernar os locais ao longo de dois eixos, de forma tal que, tanto quanto possível, as distâncias entre os locais ao longo dos eixos sejam proporcionais às distâncias da nossa matriz de distâncias. O computador vai usar uma matemática complicada para reordenar os locais e verificar o quanto as distâncias entre eles ficaram proporcionais às da matriz de distâncias. Em seguida ele vai re-arranjar os pontos e verificar se o resultado melhorou, até que não possa chegar mais perto das proporções originais.

Vamos fazer o NMDS no R (pacote BidiversityR) escolhendo duas dimensões (podemos escolher quantas dimensões queremos), pois a princípio não sabemos que nossos dados tem apenas um padrão de gradiente (precipitação, uma dimensão) e geralmente duas dimensões é a mais recomendável.

initNMS <- NMSrandom(distmatrix, perm=100, k=2) model1 <- postMDS(initNMS, distmatrix) # adicionando os escores das espécies no modelo model1 <- add.spec.scores(model1, tab, method='wa.scores') # permite colocar a ordenação das espécies no mesmo gráfico model1$stress #esses valor de stress está em % porcentagem #plotando o NMDS plot5 <- ordiplot(model1, type='text')

Vejamos abaixo um gráfico do NMDS com os dados de precipitação categorizados (A - alta; B - baixa), para visualisarmos se o eixo MDS1 (primeira dimensão do NMDS) foi capaz de capturar o padrão existente nos dados:

# categorizando os dados de precipitação prec.cat <- c("A","A","B","A","B","A", "B","B","B","A","A","B") plot6 <- ordiplot(model1, type='none') text(x=model1$points[,1], y=model1$points[,2], labels=prec.cat) text(x=model1$cproj[,1], y=model1$cproj[,2], labels=colnames(tab), col='red' )

Parece que o MDS1 foi capaz de representar bem a precipitação, embora nenhuma observação direta da chuva tenha sido usada na ordenação. A análise fez isso usando somente as distâncias (ou dissimilaridades) entre os locais em termos de abundância de espécies.

Se colocarmos em um gráfico o gradiente verdadeiro (precipitação) contra o gradiente predito pela análise NMDS, veremos que o eixo hipotético conseguiu predizer cerca de 67% (R² da regressão) da variação da precipitação.

m<-lm(prec~model1$points[,1]) summary(m) plot(x=model1$points[,1], y=prec) abline(m, col='red')

É importante dizer que em análises de NMDS é muito raro que um padrão possa ser discernível em mais do que duas dimensões, e permanecer com muitas dimensões contraria o propósito primário da análise, que é reduzir a dimensionalidade. Frequentemente se diz que os eixos de NMDS não têm interpretação lógica, e por isso podemos girar os eixos para qualquer posição no plano dos dados, sem mudar a distância relativa entre os pontos. Usando a função postMDS, nós giramos os eixos para obter a maior correlação entre as variáveis originais e o primeiro eixo MDS1. O gradiente de precipitação nesta imagem (plot5) parece correr como uma diagonal através do espaço bidimensional determinado pelos eixos MDS.

O STRESS (Standard Residuals Sum of Squares) reportado nas análises é uma medida do quanto as posições de objetos e uma configuração tridimensional desviam-se das distâncias originais (a matriz de distâncias). O valor de stress é dado em porcentagem e dependendo da função usada pode ir de 0 - 100 ou 0 - 1. Logo, o stress pode ser usado como uma medida de quão adequada a análise é. Uma regra geral sugere que:

  • stress < 0.05 representação excelente;

  • stress <0.01 boa ordenação

  • stress < 0.2 ordenação razoável

  • stress >0.2 ordenação inviável (interpretação pode ficar comprometida)

Uma outra maneira de acessar a adequação da análise é comparar as distâncias entre objetos na ordenação com as distâncias originais, chamado de Shepard diagram. Usando a função stressplot nós obtemos o gráfico e o ajuste à ordenação medido em R² tanto para uma regressão linear quanto para não linear das distâncias do NMDS para as distâncias originais.

# NMDS usando as funções do pacote vegan spe <- metaMDS(tab, distance = 'bray',k=2) spe spe$stress # esse valor de stress está entre 0 e 1 plot6 <- ordiplot(spe, type='text') # Shepard diagram stressplot(spe)

A análise de componentes principais (PCA) é principalmente usada para reduzir a dimensionalidade dos dados, e também verificar como as amostras se relacionam, ou seja, quão semelhante são segundo as variáveis utilizadas. Diferentemente de outras análises de ordenação, só é possível usar a distância euclidiana como coeficiente de similaridade na PCA.

A principal aplicação do PCA em ecologia é a ordenação dos locais com base nas variáveis ambientais quantitativas. Embora a PCA tenha um longo histórico de método para analisar variáveis ambientais, a recente introduçao dos dados de espécies pré-transformados permitiram usar esta técnica para análise de dados da comunidade (veremos adiante como fazer).

A PCA tem como principal vantagem retirar a multicolinearidade das variáveis, pois permite transformar um conjunto de variáveis originais intercorrelacionadas em um novo conjunto de variáveis não correlacionadas (os componentes principais).

O primeiro passo é extrair a matriz de distância euclidiana dos dados, mas se os dados estiverem em escalas/unidades diferentes precisamos primeiro padronizá-los. Ou seja, o PCA deve ser computado em uma tabela de variáveis dimensionalmente homogêneas. A razão para isso é que é a soma das variâncias que são particionadas nos autovalores (veja abaixo significado). Assim, as variáveis precisam estar na mesma unidade física para produzir uma soma de variâncias com significado, ou então os dados devem estar sem dimensão, que é o caso da padronização das variáveis.

Antes de partir para as análises vejamos alguns conceitos importantes:

  • Combinações lineares: equação que agrupa as diferentes variáveis, como em uma regressão múltipla.

  • Componentes principais: são as combinações lineares das variáveis, eixos ortogonais (independentes) que resumem (explicam) a variação dos objetos, e como tal podem ser consideradas como “novas” variáveis e usadas em análises posteriores. O número de componentes principais é igual ao número de variáveis. O primeiro componente principal resume a maior variação dos dados, o segundo, a segunda direção de maior variação dos dados e asim por diante.

  • Autovalores (eigenvalues): esses valores representam a medida de importância (variância) dos componentes principais (eixos) e traz a porcentagem de explicação de cada eixo. O número de autovalores é o mesmo do número de variáveis. Os autovalores serão maiores para aquelas variáveis que forem mais importantes na formação do eixo.

  • Autovetores (eigenvectors): o mesmo que Loading, ou seja, coeficientes de combinação linear. Os autovetores são os eixos principais de dispersão da matriz e medem a importância de uma variável em cada eixo. Desse modo, representam o peso de uma variável para a construção de um eixo e variam de -1 a 1 (correlação de Pearson).

  • Escores (Z1, Z2, Zn): posição das unidades amostrais ao longo de um eixo de ordenação, pode se referir tanto à unidades amostrais quanto à variáveis. Escores são fornecidos pela substituição dos valores assumidos pelas variáveis originais nas combinações lineares. São utilizados para ordenar as unidades amostrais em um diagrama uni, bi ou tridimensional.

  • Inércia: a soma de todas as correlações das variáveis com elas mesmas, mede a quantidade de variância total que é explicada por um eixo.

  • Loadings (coeficiente de estrutura): correlação de Pearson entre os escores e as variáveis.

Vamos utilizar os dados de abundância de peixes e das variáveis ambientais em 30 pontos no rio Doubs (Europa). Esses dados foram tirados do material do livro “Numerical Ecology with R” (link para o pdf em Material Extra) A matriz de peixes possuim 27 espécies e a de variáveis ambientais 11 variáveis.

#matriz de abundâncias por local spe <- read.csv("DoubsSpe.csv", row.names=1) # matriz de variáveis ambientais por local env <- read.csv("DoubsEnv.csv", row.names=1)

Existem diversas funções que fazem PCA no R, neste exemplo vamos utiliza a função rda (e outras acessórias) do vegan. Para exemplos com outras funções, veja os links para roteiros online em Material Extra.

Vamos começar fazendo uma PCA para a matriz de variáveis ambientais:

pca1 <- rda(env, scale=TRUE) # o argumento scale=T padroniza os dados antes de fazer a análise summary(pca1)

Alguns vocabulários extra do summary:

  • Species scores: coordenadas das pontas das setas das variáveis. Por razões históricas, as variáveis resposta são sempre chamadas de ‘species’ no vegan , não importando o que representam (nesse caso as variáveis ambientais).

  • Site scores: coordenadas dos locais no diagrama de ordenação. Os objetos são sempre chamados de ‘sites’ no vegan, e nesse caso são os locais mesmo (mas poderiam ser outras coisas).

O próximo passo é selecionar quais são os eixos que foram os mais importantes, ou seja, aquele que resumem a maior quantitade de variação dos dados. A decisão pode ser arbitrária, por exemplo interpretar o número de eixos necessários para representar 75% da variação dos dados, ou seguindo algum método mais objetivo. Existem diversos métodos, porém vamos focar no chamado Broken Stick que sugere considerar apenas os eixos maiores que o valor predito pelo modelo de Broken Stick. Este modelo divide aleatoriamente um “bastão” em tantos pedaços quanto eixos do PCA, e depois ordena os pedaços em comprimente decrescente. Vejamos como fazer:

screeplot(pca1) points(bstick(pca1), col="red", type='o') # também compara os eixos com o bstick PCAsignificance(pca1) # pacote BidiversityR

Usando o método do Broken Stick para decidir quantos eixos são importantes, escolhemos os dois primeiros, que acumulam cerca de 75% da variabilidade dos dados.

Existem duas escalas usadas para plotar os gráficos da PCA:

  • escala 1 (scaling 1): ótima para observar as relações entre os objetos (locais). As distâncias entre os objetos no biplot são aproximações das distâncias euclidianas no espaço multidimensional.

  • escala 2 (scaling 2) ótima para observar as relações entre as variáveis. Os ângulos entre as variáveis no biplot reflete suas correlações.

pcplot1 <- biplot(pca1, scaling=1, type="text") pcplot2 <- biplot(pca1, scaling=2, type="text") pcplot1 <- ordiplot(pca1, scaling=1, type="text") ordiequilibriumcircle(pca1, pcplot1 , col='red')

OBS: uma outra forma de fazer os plots acima em um só comando é usando a função cleanplot.pca disponível nesse arquivo (também retirada do livro “Numerical Ecology with R”. Baixe o arquivo e coloque ele na pasta da sua área de trabalho:

source("cleanplot.pca.r") # para carregar a função no seu workspace cleanplot.pca(pca1)

Agora vamos interpretar os resultados. Primeiro a proporção de variância acumulada para os dois primeiro eixos é de 0.74 ou 74%. Este valor alto nos faz confiantes de que nossa interpretação do primeiro par de eixos extrai a informação mais relevante dos dados.

No primeiro biplot, o círculo da contribuição do equilíbrio significa um comprimento de vetor representando uma variável que contribuiria igualmente com todas as dimensões do PCA. As variáveis que tem vetores (setas) mais longos que o raio do círculo contribuem mais do que a média e podem ser interpretadas com segurança.

O primeiro biplot mostra um gradiente da esquerda para a direita começando com um grupo formado pelos locais 1-10 que apresentam os maiore valores de altitude (alt) e inclinação (pen), e os menore valores em descarga do rio (deb), e distância da fonte (das) e dureza (dur). Podemos continuar a interpretação desse plot observando a formação de grupos de locais e quais variáveis ambientais estão relaciondas com estes.

O segundo biplot mostra que as variáveis também estão relacionadas em grupos, assim como os locais analisados pelo biplot 1. Por exemplo o canto inferior esquerdo mostra que a altitude e a inclinação são altamente e positivamente correlacionadas, e que estas duas são altamente negativamente correlacionadas com os vetores das variáveis que vão em direção oposta (deb, dur, das). Cabe observar que nitrato (nit) e pH tem setas praticamente ortogonais (ângulo de 90 graus) indicando uma correlação próxima à zero.

Este exemplo mostra quão útil uma representação do biplot pode ser em sumarizar as características principais dos dados. Podemos ver agrupamentos e gradientes, assim com as correlações entre as variáveis.

A PCA, por ser um método que utiliza as distâncias euclidianas entre os locais, não é naturalmente adaptada à análise de dados de abundâncias de espécies. Entretanto, hoje em dia aplica-se a transformação dos dados para poder realizar a PCA. A pré-transormação assegura que os dados de espécies sejam tratados de acordo com sua especificidade, ou seja, sem dar importância indevida a presenças de muitos zeros ou abundância muito grandes. Um exemplo de transformação muito utilizada é a transformação de Hellinger, que nada mais é do que dividir as abundâncias de cada espécie em um local pela abundância total do local (linhas) e então fazer a raiz quadrada deste valor. Esta transformação reduz o a importância de grandes abundâncias no cálculo da distância euclidiana.

Vejamos como fazer com os dados de abundâncias de peixes.

# Transformando os dados spe.h <- decostand(spe, "hellinger") spe.h.pca <- rda(spe.h) summary(spe.h.pca) screeplot(spe.h.pca) points(bstick(spe.h.pca), col="red", type='o') cleanplot.pca(spe.h.pca)

Vemos aqui que as espécies não formam grupos claros como visto para as variáveis ambientais. Vemos que 8 espécies estão contribuindo mais fortemente nos eixos 1 e 2 (as setas que estão para fora do círculo). O biplot 1 nos revela um gradiente por trás estruturando a comunidade, os locais são ordenados ao longo dos eixos de acordo com suas posições ao longo do gradiente. O biplot2 nos mostra a relação entre as espécies.

Cabe ressaltar que a PCA pode ser aplicada também a dados binários de presença/ausência após a transformação de Hellinger.