Três aspectos essenciais da estrutura das comunidades são a riqueza de espécies, suas abundâncias e a composição de espécies. Um dos grandes desafios em ecologia é avaliar o grau de similaridade entre comunidades em relação a esses três aspectos e quais os possíveis fatores que determinam os padrões observados.

Hoje em dia sabemos que as comunidades não são homogêneas e, portanto, diferentes conjuntos de espécies são encontrados em diferentes localidades. Para entender essas diferenças, o primeiro passo é saber quais localidades são mais similares entre si e quais são mais dissimilares. Uma forma de fazer isso é utilizando os chamados métodos de ordenação. Existem vários métodos com características diferentes e adequados a conjuntos de dados específicos. Alguns textos importantes estão disponíveis no final desse roteiro para quem se interessar mais pelo assunto.

Para começarmos nosso exercício prático de hoje, crie uma pasta para você na área de trabalho do seu computador. Copie para essa pasta os três arquivos abaixo:

• especies.csv
• ambientais.csv
• coordenadas.csv

Em seguida, abra o programa R, clicando no ícone que está na área de trabalho do seu computador.

Mude o diretório de trabalho do programa R para a sua pasta. Isso pode ser feito de duas formas:

• Digitando setwd(“C:/desktop/nome da sua pasta”) no R
• Abrindo o Menu “Arquivo” (ou “File”) e selecionando “Mudar dir” (ou “Change dir”). Nesse caso escolha a sua pasta na janela que abrir. [Obs. Para Mac, essa opção está no Menu “Misc” e a opção é “Change working dir”, ou usando a tecla command+D]

Para checar se você está na pasta correta, copie e cole o comando abaixo no console do R:

getwd()

e veja se aparece o nome da sua pasta. Se sim, ótimo. Se não, chame um monitor ou professor.

Para fazer as análises vamos usar as funções já prontas de um pacote muito legal para quem quer trabalhar com análises de vegetação. O pacote se chama vegan e para instalá-lo, use o comando:

install.packages("vegan")

Se o programa perguntar se você quer instalar as dependências, diga que sim.

E agora, carregue o pacote no R

library(vegan)

Se você for usar o mesmo computador outras vezes, não precisará instalar novamente o pacote vegan (a não ser quando lançarem uma nova versão), mas toda vez que for fazer essas análises usando esse pacote você precisará carregar (com a função “library” ou “require”), ok?

Antes de aplicar os métodos de ordenação, vamos entender um pouquinho algumas das informações que eles utilizam.

Imagine que amostramos uma área de 1km x 1km, dividindo-a em 100 parcelas de 100x100m. Amostramos todas as espécies arbóreas acima de 20 cm de DAP e anotamos o número de indivíduos de cada espécie (ou seja, a abundância de cada espécie) em cada parcela.

Tabulamos os dados em uma planilha cujas linhas são as amostras (parcelas ou sites) e as colunas são as abundâncias das espécies ¹⁾ Se quiser ver nossa planilha completa de espécies por parcela abra o arquivo “especies.csv” diretamente em um editor de texto ou planilha eletrônica (p. ex. Excel) ou importe para o R e veja dentro do próprio R, como explicaremos abaixo:

O primeiro passo é importar a planilha de dados para o R e criar um objeto para essa planilha:

sp<-read.csv2 ("especies.csv")

Se quiser ter certeza que a planilha foi copiada corretamente, visualize as primeiras cinco linhas da planilha usando o comando:

head(sp)

Se quiser ver um resumo das variáveis existentes na planilha, use o comando:

summary(sp)

Se tivéssemos apenas duas espécies, seria bem fácil saber quais parcelas são mais parecidas entre si. Bastaria plotarmos as abundâncias de cada espécie em um gráfico de dispersão e medir a distância em linha reta entre as duas parcelas (essa distância é chamada de Distância Euclidiana). Veja esse exemplo hipotético de duas espécies e duas parcelas:

Se existissem mais parcelas, bastaria medir as distâncias entre cada par de parcelas ²⁾.

Se tivéssemos três espécies, ainda seria possível visualizar a semelhança entre as parcelas, mas com mais dificuldade.

Porém, acima de três espécies fica praticamente impossível visualizar as relações entre as parcelas.

Uma outra forma de tentar visualizar como as espécies estão distribuídas entre as parcelas seria olhar para as correlações entre cada par de espécies, tentando entender quais espécies tendem a ocorrer com abundâncias similares.

Então, vamos fazer a matriz de correlação entre todos os pares de espécies da nossa tabela e olhar como elas se relacionam!

Antes, precisamos importar uma função para facilitar a visualização dessas correlações. Copie tudo que está dentro desse quadro e cole no R:

## Função para colocar coeficientes de correlação em uma das diagonais
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
     {
         usr <- par("usr"); on.exit(par(usr))
         par(usr = c(0, 1, 0, 1))
         r <- abs(cor(x, y))
         txt <- format(c(r, 0.123456789), digits = digits)[1]
         txt <- paste0(prefix, txt)
         if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
         text(0.5, 0.5, txt, cex = cex.cor * r)
     }

Não aconteceu nada?? Que ótimo, não precisava acontecer mesmo. A função simplesmente foi salva e armazenada no R para ser usada depois….

Agora vamos aplicar uma função (“pairs”) que vai calcular a correlação par-a-par na nossa tabela de espécies por parcelas (que se chama “sp”) e vamos utilizar a função “panel.cor” para melhorar a visualização.

pairs (sp, upper.panel=panel.cor)

Como entender essa figura?
Essa figura é uma representação da matriz de correlação entre as abundâncias das espécies nas parcelas. O nome/código das espécies está na diagonal.

Então, para entender qual espécie está represnetada no eixo x e qual está no eixo y de cada gráfico basta encontrar os nomes das espécies nas linhas e colunas correspondentes. Por exemplo, o gráfico localizado na linha 4, coluna 3, tem o gráfico de dispersão das abundâncias da espécies 3 (valores eixo x) e 4(valores no eixo y).

Os valores que estão acima e à direita da diagonal principal são os valores de correlação (entre -1,0 e +1,0) referentes aos gráficos que se encontram na posição espelhada abaixo da diagonal. Para facilitar a visualização, quanto maior a letra do número, maior a correlação. Correlações muito baixas sequer aparecem.

Inicialmente vamos aprender um dos métodos, conhecido como Análise de Componentes Principais (Principal Component Analysis - PCA).

Se conseguíssemos visualizar as correlações entre as 10 espécies em um mesmo gráfico seria lindo, mas infelizmente não temos essa capacidade… ³⁾

Quando temos mais de 3 espécies, fica difícil identificar quais parcelas são mais parecidas com as outras, levando em consideração todas as espécies e suas abundâncias. Entretanto, matematicamente isso é possível, por meio de álgebra matricial.

A Análise de Componentes Principais (PCA) procura, na nuvem de pontos formada pelas abundâncias das 10 espécies nas 100 parcelas, um eixo imaginário que melhor expressa a similaridade e as diferenças entre as parcelas. Esse eixo é definido pelos autovalores ⁴⁾ obtidos a partir de uma matriz de correlação entre as espécies (que são aqueles valores que estavam na parte direita da matriz de correlação que vimos acima).

Esse primeiro eixo é conhecido como Componente Principal 1 (ou pela sigla PC1) e ele vai usar as correlações entre as espécies para colocar de um lado do eixo as parcelas mais parecidas entre si e no lado oposto, as parcelas mais diferentes dessas. As parcelas que não são parecidas com esses extremos ficam na região intermediária do eixo.

Então, vamos ver nessa nossa comunidade hipotética, como fica esse primeiro eixo.

Vamos realizar a Análise de Componentes Principais, utilizando uma função já pronta no pacote vegan do R. Vamos colocar um nome no objeto que vai ser gerado pela análise, pois queremos ver detalhadamente os resultados encontrados.

pca.sp <-rda(sp, scale = TRUE)

A função rda está sendo aplicada na tabela “sp” e vai gerar o objeto “pca.sp”.

Mas, o que é scale=TRUE?

Isso é um requisito do pacote “vegan” para a realização de uma PCA baseada em correlações. Com esse comando os valores que serão utilizados na análise não serão os valores reais das abundâncias das espécies nas parcelas, mas sim os valores corrigidos (ou melhor, escalonados) pelo desvio padrão das abundâncias das espécies nas parcelas. É meio complicado de entender, mas não se preocupe com isso agora.

Bom, agora vamos visualizar o primeiro eixo gerado pela PCA e olhar quais parcelas são mais parecidas e mais diferentes.

plot(pca.sp, choices= c(1), display= c("sp","wa"), scaling = "species", cex=0.5)
## Para gravar um pdf com o grafico
pdf(file="pca1.pdf", height=9)
plot(pca.sp, choices= c(1), display= c("sp","wa"), scaling = "species", cex=0.5)
dev.off()

Esse gráfico parece estranho, mas não é complicado. Vamos entendê-lo:

Identifique quais são as parcelas (sit1, sit2,…,sit10,…,sit100?) mais parecidas entre si nesse eixo e quais são as espécies(sp1, sp2,…sp3?) que estão mais próximas dessas parcelas nesse eixo. Essas são as espécies que estão definindo a similaridade entre as parcelas de acordo com esse eixo.

Ok, mas uma PCA não produz apenas um eixo de explicação da variação dos dados. O primeiro eixo (esse que acabamos de visualizar) é sempre o eixo que explica a maior parte da variação dos dados. O segundo eixo produzido pela PCA será um eixo ortogonal ao primeiro eixo, pois assim, fica garantido que o segundo eixo é independente do primeiro. Ele será o eixo ortogonal que explica a maior parte da variação restante. E assim por diante com os outros eixos…Os últimos eixos geralmente explicam muito pouco sobre a ordenação das parcelas.

Uma visualização bem simplificada usando 3 variáveis apenas:
https://www.youtube.com/watch?v=BfTMmoDFXyE

Então, agora, vamos ver os resultados da PCA para a nossa comunidade hipotética:

pca.sp

Vamos entender alguns pontos importantes desse resultado:

Inércia (Inertia) é o termo usado no pacote vegan para indicar a variação nos dados. O valor total aqui nesse caso que estamos calculando é 10, pois estamos usando a matriz de correlação escalonada com valores de 0 a 1 para cada variável (espécies) e temos 10 variáveis. Então, o valor máximo que a inércia total pode assumir é 10.

Logo abaixo, a gente tem os valores de variação explicada por cada um dos eixos. Para sabermos a porcentagem de explicação de cada eixo, basta dividirmos esse valor pela inércia total.

E esse tal de “Unconstrained”??? Calma, mais abaixo a gente explica… por agora, não se preocupe com isso.

Mas, ainda temos um resultado muito simplista da PCA. Para ver um resultado mais completo e entender melhor o que são esses diferentes eixos, use:

summary (pca.sp)

Agora já podemos ver a proporção da variação total que cada eixo explica (“Proportion Explained”) e também a proporção acumulada à medida que adicinamos eixos (“Cumulative Proportion”).

Além disso, mais embaixo, você verá a tabela com os “Species scores”. Nas linhas estão as espécies e nas colunas os eixos da PCA.

É isso aí. Se tudo deu certo, esses valores devem ser relativamente parecidos.

Então, os escores das espécies (“Species scores”) representam a projeção da espécie em cada um dos eixos calculados pela PCA.

Antes de mostrarmos o gráfico com esses eixos plotados, vamos olhar agora para a segunda tabela, dos “Sites scores”

Ok, agora chega de imaginar. Vamos ver como ficam os dois primeiro eixos da PCA em um gráfico e como as parcelas estão organizadas nesse espaço bidimensional.

Vamos fazer o gráfico usando uma função especial para esse tipo:

biplot(pca.sp, scaling=2, type = c("text", "text"))

Salve essa figura em sua pasta. Faça o mesmo com as próximas figuras de ordenação que produzir.

Obs.: Note que o eixo 1 (aquele que a gente já visualizou antes), está agora na horizontal, como eixo x, ok?

Os vetores indicam em que direção aumenta a abundância de cada espécie. Os pontos neste gráfico representam as parcelas.

Precisamos olhar também a relação com o terceiro eixo?
Para esse nosso exercício prático não vamos olhar, pois os dois primeiros eixos (PC1 e PC2) já explicam cerca de 62% da variação e consideramos que isso é suficiente para os nossos propósitos. Porém, em alguns casos o terceiro eixo também pode trazer informações importantes⁵⁾.

Muito bem, depois de entendermos a lógica básica de funcionamento de uma PCA e interpretarmos a ordenação produzida pelos dois primeiros eixos, a pergunta que surge naturalmente é:

“Por que algumas espécies ocorrem predominantemente e em maior abundância em algumas parcelas e outras não?”

Tempo para pensar!

Uma das respostas mais rápidas que vêm à mente de um(a) ecólogo(a) é que parcelas com maior similaridade de espécies podem ter variáveis ambientais mais parecidas e que diferenciam do restante.

Como podemos avaliar isso?

Podemos ir a campo e medir variáveis ambientais em todas as parcelas, como características físicas do solo (granulometria, porosidade, etc), características químicas do solo (N, P, K, Ca, Mg, pH, Matéria Orgânica, Capacidade de Troca Catiônica, etc), topografia e luminosidade, etc. Com estes dados podemos contruir uma tabela em que a linhas são de novo as parcelas mas as colunas são a variávéis ambientais. Com esta tabela podemos fazer uma PCA. Nesse caso, ao invés das espécies, as variáveis que vão definir a ordenação das parcelas são as variáveis ambientais.

Vamos conhecer então as variáveis ambientais desse nosso exercício prático:

amb<-read.csv2 ("ambientais.csv")
head (amb)
summary(amb)

Nessa planilha temos as 5 variáveis ambientais hipotéticas e os valores de cada uma delas em cada parcela de amostragem.

Vamos ver como essas variáveis estão correlacionadas par a par:

pairs(amb, upper.panel=panel.cor)

Agora, vamos ver como fica a ordenação das parcelas de acordo com essas variáveis ambientais, aplicando a PCA nesse conjunto de dados.

pca.amb<-rda(amb, scale = TRUE)
pca.amb
summary (pca.amb)

biplot (pca.amb, scaling=2,type = c("text", "text"))

Será que essa ordenação é parecida com a ordenação obtida para as espécies?

Podemos colocar os gráficos lado a lado e comparar. Porém, fazendo dessa forma, teremos apenas uma comparação visual e não podemos avaliar o quanto essas variáveis ambientais explicam a distribuição das espécies nas parcelas.

Podemos fazer uma análise que junta essas duas matrizes de informações sobre as parcelas (abundâncias das espécies e variáveis ambientais), de modo que a ordenação das parcelas pelas espécies seja condicionada (“constrained”) pelas variáveis ambientais e assim podemos calcular o quanto as variáveis ambientais explicam a distribuição das espécies nas parcelas.

Uma dessas análises é chamada de Análise de Redundância (Redundancy Analysis - RDA), que é feita a partir da mesma lógica que usamos antes para fazer as PCAs, porém, reordenando as parcelas em função das variáveis ambientais.

O primeiro passo é analisar as relações lineares entre as abundâncias das espécies e os valores das variáveis ambientais em cada parcela.

Vamos olhar como ficam as relações no caso das nossas amostras. Vamos ver qual(is) das 10 espécies está(ão) relacionada(s) linearmente com qual(is) variável(is) ambiental(is). Para isso, copie e cole os comandos abaixo:

par(mfrow=c(10,5), mar=c(0.1,0.1,0.1,0.1))
for(i in 1:10){
    for(j in 1:5){
        plot(sp[,i] ~ amb[,j], axes=FALSE)
        box()
    }
}
par(mfrow=c(1,1))

Nessa figura cada gráfico apresenta a relação entre os valores de abundância de uma espécie (linhas), com os valores de uma variável ambiental (colunas) para todas as parcelas. Assim o primeiro gráfico em cima e à esquerda apresenta a relação entre a espécie 1 (sp1) e a variável 1 (x1). O gráfico logo abaixo apresenta a relação entre a espécie 2 (sp2) e a variável 1 (x1), e assim por diante.

Não vamos entrar nos detalhes técnicos da análise, mas, grosso modo, as variáveis ambientais é que condicionam a nova ordenação das parcelas. Essa nova ordenação será feita a partir dos valores das abundâncias das espécies ajustados às variáveis ambientais. É por esse motivo que esse tipo de análise é chamado de análise de ordenação restrita (ou constrained ordination analysis).

rda.sp.amb<-rda(sp, amb) 

OBS.: note que estamos usando a mesma função “rda”, mas agora temos duas matrizes dentro dos parênteses (“sp” e “amb”) e não estamos usando “scale=TRUE” como fizemos na análise de PCA)

Vamos olhar e interpretar os resultados obtidos

summary(rda.sp.amb)

A primeira informação importante desses resultados é essa tabela de Partição da Variância (Partitioning of variance).

A proporção do item “Constrained” representa a proporção da variação das abundâncias das espécies nas parcelas que é explicada pelas variáveis ambientais que adicionamos na análise de ordenação restrita (Constrained). Esse valor tem um significado “similar” ao do R^{2 6)} de uma análise de regressão linear.

A proporção do item “Unconstrained” representa a variação não explicada pelas variáveis ambientais escolhidas.

Continuando a analisar os resultados obtidos, temos os valores de proporção de explicação dos eixos de ordenação criados pela RDA. Note que são 5 eixos restritos (RDs) e 10 eixos não restritos (PCs). Se você somar a proporção explicada pelos 5 eixos RDs verá que o valor é igual ao valor total do item “Constrained” na tabela anterior. Já na soma dos eixos PCs o valor é igual ao do item “Unconstrained”

Abaixo estão as tabelas com os valores:

• dos escores das espécies
• dos escores das parcelas (sites) ponderados pelas abundâncias das espécies
• das escores das parcelas (sites) mas já com os ajustes impostos pelas variáveis ambientais
• dos escores das variáveis ambientais

A partir dos dados dessas tabelas, nós podemos fazer um gráfico que mostre a ordenação das parcelas, das espécies e das variáveis ambientais nos dois primeiros eixos (RDA1 e RDA2) para entender quais são as variáveis ambientais que determinam a distribuição das espécies, que por sua vez determinam a ordenação das parcelas.

Esse gráfico é um pouco mais complicado, então, vamos explicar passo-a-passo

Primeiro vamos montar um gráfico vazio… estranho, hein? Mas acredite, vai ser legal!

plot(rda.sp.amb, xlim=c(-5,13), type="none")

Depois, vamos adicionando os elementos aos poucos. Vá copiando e colando linha por linha e veja como o gráfico é montado.

points(rda.sp.amb, display="sites", cex=0.01) ##Adicionando os pontos referentes às parcelas

text(rda.sp.amb, display="sites", cex=0.5) ##Adicionando os nomes das parcelas

text(rda.sp.amb, display="bp", col="blue") ##Adicionando as setas das variáveis ambientais

spamb2.scores<-scores(rda.sp.amb, choices=1:2, display="sp") ##selecionando os escores das espécies

arrows(0,0,spamb2.scores[,1],spamb2.scores[,2], length=0, lty=1, col="red") ##adicionando linhas para as espécies

text(rda.sp.amb, display = "species", cex=0.8, adj=-0.3, col="red") ##Adicionando nomes às linhas das espécies

Como vocês já perceberam, as variáveis ambientais estão representadas pelas setas azuis, as espécies pelas linhas em vermelho e as parcelas pelos pontos pretos.

Ok, até aqui então, temos uma explicação bem razoável da distribuição das espécies nas parcelas, pois as variáveis ambientais parecem ser importantes na determinação da distribuição das espécies.

Porém, existe ainda uma outra explicação que os(as) ecólogos(as) sempre precisam avaliar antes de afirmar que a distribuição das espécies é definida pelas condições ambientais.

Tempo para pensar

Ok, vamos lá, você já deve estar começando a ficar cansado(a), então vamos ajudar.

As espécies podem estar distribuídas de forma agregada em algumas parcelas e isso pode ser simplesmente um reflexo da limitação da dispersão das espécies. Se isso estiver acontecendo, parcelas mais próximas devem apresentar maior similaridade de composição e abundância de espécies. Como parcelas mais próximas também podem ter condições ambientais semelhantes, a associação entre espécies e variáveis ambientais pode ser apenas uma correlação espúria, devido à agregação. Precisamos avaliar se isso acontece!

Então, é importante que tenhamos os dados de localização de cada uma das parcelas estudadas para podermos incluir na nossa próxima análise.

É exatamente isso que temos na nossa planilha coordenadas.csv Vamos importá-la para o R:

coords<-read.csv2 ("coordenadas.csv")
head (coords)
summary(coords)

Se usarmos essa planilha diretamente na nossa próxima análise, a estrutura espacial das parcelas com abundâncias similares das espécies será avaliada considerando-se apenas os ângulos retos e as diagonais de distância entre parcelas. Porém, estruturas espaciais podem ser mais complicadas que as descritas por distâncias em linhas retas e diagonais. Na imagem abaixo são mostradas diferentes estruturas espaciais possíveis.

Nessa figura, quadrados maiores representam parcelas mais similares em termos de abundância das espécies. Note como seria complicado descrever alguns desses padrões usando sistemas simples de coordenadas.

Diante disso, para podermos incorporar outras estruturas espaciais possíveis teremos que aplicar dois procedimentos (PCNM - Principal Coordinates of Neighbour Matrices e forward selection). Esses procedimentos são muito usados em trabalhos científicos cujos objetivos envolvem a análise da importância do espaço na estruturação de comunidades. Vale a pena conhecer como funcionam, mas neste exercício não vamos entrar em detalhes sobre esses métodos e vamos usar uma tabela com as variáveis espaciais já trabalhadas.

Então, vamos pedir que acredite em nós!

Se você acreditar, copie essa tabela para a sua pasta:

tab_coords_sel.csv

Agora,importe a tabela para o R:

tab.coords.sel<-read.csv2("tab_coords_sel.csv")

Essa tabela importada (“tab.coords.sel”) é a informação que será usada na última etapa do nosso exercício para representar a estrutura espacial dos dados. Nela estarão selecionadas as coordenadas que melhor explicam a distribuição espacial das abundâncias das espécies nas parcelas.

Por fim, vamos calcular a porcentagem de variação da distribuição das espécies que é explicada:

• exclusivamente pelas variáveis ambientais
• exclusivamente pela estrutura espacial
• pelos dois efeitos em conjunto

O valor do efeito conjunto das variáveis ambientais e da estrutura espacial revela a relação entre essas duas matrizes de dados. Valores altos desse efeito conjunto indicam que parcelas mais próximas têm variáveis ambientais mais similares, ou, em outras palavras, que as variáveis ambientais estão espacialmente estruturadas.

Vamos então aplicar essa análise conhecida como Partição da Variação (Variation Partitioning) e analisar as porcentagens explicadas pelas variáveis ambientais, pela estrutura espacial e por ambas.

(particao <- varpart(sp, amb , tab.coords.sel))

A tabela mostra a variação total das abundâncias das espécies decomposta em:

• $a+b+c =$ variação total
• $a+b =$ variação total devido ao ambiente ($X1$)
• $b+c =$ variação total devido ao espaço ($X2$)

Note que

• $a = $ variação exclusivamente devido ao ambiente ($X1|X2$)⁷⁾
• $c =$ variação exclusivamente devido à estrutura espacial ($X2|X1$) ⁸⁾

Um diagrama que ajuda entender esta partição:

showvarparts(2, bg = c("hotpink","skyblue") )

• Borcard et al. (2011). Numerical Ecology with R.
• Legendre & Legendre (2003). Numerical Ecology.
• Sobre os PCNMs: Borcard, D. and P. Legendre. 2002. All-scale spatial analysis of ecological data by means of principal coordinates of neighbour matrices. Ecological Modelling 153(1-2): 51-68.
• Ótima página de Michael Palmer
• Ótima página de David Zeleny
• Tutorial do vegan
• Introdução a análise multivariada do vegan