Pacotes e formatos de dados

Existem diversos pacotes no R que fazem análises de dados em comunidades. Neste roteiro, iremos utilizar dois dos principais pacotes existentes para análises de diversidade que requerem um formato específico para os dados de abundâncias de espécies (a matriz de espécies por local) e os dados de características ambientais dos locais (a matriz das variáveis por local).

install.packages("vegan")
install.package("BiodiversityR")

library(vegan)
library(BiodiversityR)

Em ambas as matrizes, os locais ficam nas linhas, enquanto nas colunas ficam as espécies para a matriz de abundâncias e as variáveis ambientais para a matriz de caraterísticas ambientais.

Vamos usar os dados de vegetação de dunas disponíveis no pacote vegan:

#carregando os dados para o seu workspace
data("dune")
dune

Vemos que estes dados tem 30 espécies (colunas), em 20 locais (linhas). Veja o help(dune) dos dados para entender melhor de onde vieram.

A matriz dos dados ambientais, com 5 variáveis, é:

data("dune.env")
dune.env
str(dune.env)

Nunca se esqueça de fazer uma inspeção nas tabelas de dados. Reveja o roteiro 2 sobre análises exploratórias de dados para inspecionar e visualizar graficamente seus dados.

Curva de acumulação de espécies

Curvas de acumulação de espécies, algumas vezes chamadas de curva do coletor, são representações gráficas que demonstram o número acumulado de espécies registradas (S) em função do esforço amostral (n). O esforço amostral pode ser o número de indivíduos coletados, ou uma medida tal como o número de amostras (e.g., quadrados) ou tempo amostral (e.g., meses). Se as curvas de acumulação de espécies atingem um ponto em que o aumento do esforço de coleta não implica num aumento no número de espécies, isto significa que aproximadamente toda a riqueza da área foi amostrada.

Dependendo do protocolo de amostragem usado para o levantamento da riqueza de espécies podemos criar uma curva baseada em indivíduo, ou baseada em amostra. Tendo como exemplo a amostragem de árvores numa floresta, na abordagem baseada em indivíduo, as árvores seriam amostrados aleatoriamente dentro de cada parcela e registrados sequencialmente. Na abordagem baseada em amostra, as parcelas teriam seu número de espécies e identidades registrados para toda a parcela e acumularíamos o número total de espécies a cada parcela amostrada. A chave da distinção das duas abordagens é a unidade de replicação: o indivíduo vs. uma amostra de indivíduos.

Para construir a curva de acumulação de espécies vamos usar o pacote BiodiversityR, que além de poder ser usado diretamente no R, tem também uma interface amigável através do pacote Rcmdr. Veja o manual do pacote aqui.

Abaixo a curva do coletor baseada em amostra para os dados de vegetação das dunas:

Accum <- accumresult(dune, method='collector')
Accum
accumplot(Accum)

Gotelli & Colwell (2001) descrevem dois tipos de curva, a curva do coletor e a curva de rarefação. A curva acima é curva do coletor, onde a riqueza de espécies foi plotada de acordo com a ordem dos locais amostrados (argumento `method=“collector”).

Colwell & Coddington (1994) sugeriram um método que consiste em montar várias curvas adicionando-se as amostras em uma ordem aleatória. Após construir várias curvas com este método, pode-se calcular uma curva do coletor média (baseada na riqueza média para cada número de amostra) e expressar a variação possível em torno dessa média. Isso pode ser feito de duas maneiras. A primeira é utilizando a riqueza média agrupada quando 1,2… todos os locais são combinados juntos. A razão para calcular a riqueza agrupada média de espécies é que diferentes combinações vão ter diferentes riquezas. Por exemplo, nestes dados o local 1 contém 5 espécies, enquanto que o local 2 e 3 contém 10 espécies. Em outras palavras, nao existe um valor único para a riqueza de espécies em UM local e por isso, o valor médio é calculado.

Accum.1 <- accumresult(dune, method='exact')
Accum.1
accumplot(Accum.1)

O resultado de Accum.1 mostra que a riqueza média para todas as possíveis combinações de 14 locais é 28.85, equanto que para todas as possíveis combinações de 17 locais é 29.51. O sd indica o desvio padrão observado, já que algumas combinações de locais contém menores número de espécies e outras combinações maiores.

Os resultados acima são baseados nos cálculos exatos da riqueza de espécies média para combinações de locais (method= "exact"), como sugerido por Colwell & Coddington (1994). A segunda forma de fazer isso é aproximar essa riqueza média de espécies em uma abordagem de aleatorização (também chamao de ‘rarefação’ por Gotelli & Colwell 2001). Por exemplo, para encontrar uma estimativa da riqueza média para 3 locais, selecionamos 3 locais aleatoriamente de todos os 30 locais. Calcumamos a riqueza de espécies e então selecionamos outros 3 locais aleatoriamente. Repetimos isso umas 1000 vezes e pegamos a média de todas estas estimativas.

Os resultados da curva de acumulação de espécies baseado nesta aleatorização para as espécies das dunas será:

# usando method="random"
Accum.2 <- accumresult(dune, method='random', permutations=1000)
Accum.2
accumplot(Accum.2)

# usando method = "rarefaction"
Accum.21 <- accumresult(dune, method='rarefaction')
Accum.21
accumplot(Accum.21)

OBS: para ver os diferentes métodos de criação das curvas de acumulação, veja o argumento method em help(accumresult).

Para alguns conjuntos de dados, a informação não é coletada por local, e sim por indivíduo observado. Frequentemente, levantamentos de vertebrados são obtidos dessa maneira, registrando as informações das espécies para cada indivíduo observado sequencialmente. Entretanto, há dados baseados em amostra que podem ser usados para fazer curvas de acumulação baseadas tanto em indivíduos quanto em amostras.

Vejamos a curva baseada em indivíduo para os dados das dunas:

# primeiro sabendo o número de indivíduos em cada site
 dune.env$site.totals <- apply(dune,1,sum)

# depois usando isso como 'scale' na curva:
Accum.5 <- accumresult(dune, y=dune.env, scale='site.totals', method='rarefaction')
Accum.5
accumplot(Accum.5, xlab='pooled individuals')

Neste gráfico percebemos que os resultados não indicam todas as 685 possíveis combinações de indivíduos (número total de indivíduos amostrados). Os resultados foram reportados apenas para o número médio de indivíduos em cada local, 685/20 = 34.25 (arredondando para 34 indivíduos). Neste exemplo, nós re-escalonamos o eixo x para número de indivíduos baseado no tamanho médio de indivíduos por amostra. Quando usamos as curvas de rarefação baseadas em amostra para comparar riqueza de espécies em níveis comparativos de esforço amostral, o número de espécies deve ser plotado como uma função do número acumulado de indivíduos (como feito no gráfico acima), porque os conjuntos de dados podem diferir quanto ao número médio de indivíduos por amostra (Gotelli & Colwell 2001).

Usando a curva de acumulação de espécies para comparar riqueza em subgrupos de dados

As curvas de acumulação de espécies (usando método de rarefação) são especialmente úteis quando comparamos a riqueza de espécies para subgrupos nos dados quando os tamanhos amostrais são diferentes destes subgrupos.

Vamos calcular a curva de acumulação de espécies para cada categoria de manejo dos dados das dunas (tabela de variáveis ambientais, coluna Management):

Accum.6 <- accumcomp(dune, y = dune.env, factor = 'Management', method = 'exact', legend = F)
Accum.6

Olhando para o gráfico acima, podemos observar uma riqueza de espécies de 16 para um tamanho amostral de 3 locais no manejo o “BF”, e uma riqueza total de 21 para um tamanho amostral de 5 locais para o manejo “HF”. Os resultados agora permitem comparar a riqueza de espécies em um mesmo tamanho amostral. Por exemplo, quando comparamos as 4 categorias de manejo no tamanho amostral 3, vemos que HF tem uma maior riqueza de espécies do que as outras 3 categorias.

Estimando o número total de espécies da área

Na maioria das situações, nós registramos o número de espécies de uma quantitade de locais que não cobre toda a área que estamos interessados. Nós não saberemos a composição de espécies das áreas que não amostramos.

Diversas fórmulas já foram desenvolvidas para estimar o número total de espécies de uma área. Em termo de curvas de acumulação de espécies, nós precisamos extrapolar a curva até o ponto do eixo x que corresponde ao tamanho total da área. Estes métodos de estimação são baseados em diferentes premissas sobre como as espécies se acumulam além da área amostrada. Provavelmente um método vai ser mais preciso e uma situação, e outro em outra.

Vejamos as propriedades de alguns estimadores (para mais detalhes ver Santos 2003 e Provete et al. 2011):

• Jackniffe 1 e 2: estimam a riqueza total somando a riqueza obsevada a um parâmetro calculado a partir do número de espécies raras e do número de amostras. As duas equações diferem basicamente em relação ao critério pelo qual se consideram uma espécie como rara.

\[ S_{jack1} = S_{obs} + Q_1\frac{m-1}{m} \]

\[ S_{jack2} = S_{obs} + Q_1\frac{2m-3} {m} - Q_2\frac{(m-2)^2} {m(m-1)} \]

\(S_{obs}\): número de espécies observadas;
\(Q_j\): número de espécies que ocorrem em j amostras;
\(m\): número de amostras;

• Chao 1 e 2: estimadores baseados no número de espécies raras dentro de uma amostra. A riqueza estimada pelo chao 1 é igual à riqueza observada, somada ao quadrado do número de espécies representadas por apenas um indivíduo (singletons), dividido pelo dobro do número de espécies com apenas dois indivíduos (doubletons). A mesma equação foi adaptada para utilizar o número de espécies que ocorrem respectivamente em uma ou em duas unidades amostrais (uniques e duplicates - chao 2). Chao 1 requer dados de abundância das espécies, enquanto chao 2 trabalha com dados de presença/ausência.

\[ S_{chao} = S_{obs} + \frac{F_1^2} {2F_2^2} \]

\(F_j\): número de espécies que tem exatamente j indivíduos em todas amostras juntas (chao1) ou que ocorrem em exatamente j amostras (chao2)

• Bootstrap: se difere dos demais por utilizar dados de todas as espécies coletadas para estimar a riqueza total, não se restringindo às espécies raras. É calculada somando-se a riqueza observada à soma do inverso da proporção de amostras em que ocorre cada espécies.

\[ S_{boot} = S_{obs} + \sum_{k=1}^{S_{obs}}(1-p_k)^m \]

\(p_k\): proporção de amostras que contém a espécie k;

Se queremos fazer esta extrapolação, então a melhor abordagem é reportar a variação nos valores obtidos com os diferentes métodos e esperar que o valor total correto esteja dentro desse intervalo.

Para os dados das dunas, fizemos a predições da riqueza total de espécies usando 4 estimadores diferentes:

Diversity.5 <- diversityresult(dune, index='jack1')
Diversity.5

Diversity.6 <- diversityresult(dune, index='jack2')
Diversity.6

Diversity.7 <- diversityresult(dune, index='chao')
Diversity.7

Diversity.8 <- diversityresult(dune, index='boot')
Diversity.8

Podemos ver que as estimativas foram similares, variando entre 31-34 espécies. Isso não nos surpreende, pois já vimos que o formato da nossa curva de acumulação de espécies praticamente atingia um platô, sugerindo que a amostragem conseguiu capturar quase toda as espécies da área de estudo.

Índice de Simpson

A medida mais simples para caracterizar a comunidade, e que leva em consideração tanto o padrão de abundância quanto a riqueza de espécies é o Índice de diversidade de Simpson. Ele é calculado obtendo-se para casa espécie a proporção de indivíduos em relação ao total da amostra, ou seja, a proporção \(P_i\) para a i-ésima espécie. Pode ser calculado de duas formas:

\[ D_1 = 1 - \sum_{i=1}^S P_i^2 \]

ou o inverso do Simpson:

\[ D_2 = \frac{1}{\sum_{i=1}^S P_i^2} \]

em que \(S\) é a riqueza de espécies na comunidade. Para uma dada riqueza, \(D\) aumenta com a equabilidade (ou equitabilidade, o contrário de dominância), e para uma dada equabilidade, \(D\) aumenta com a riqueza.

A equabilidade também pode ser quantificada (entre 0 e 1), obtendo-se a proporção entre o índice de Simpson e o valor máximo que este assumiria caso os indivíduos fossem distribuídos uniformemente entre as espécies. \(D_{max} = S\), portanto:

\[ E = \frac{D} {D_{max}} = \frac{1}{\sum_{i=1}^S P_i^2} x \frac{1}{S}\]

O índice de Simpson também é conhecido como índice de Dominância. Sua interpretação é simples e ele significa a probabilidade de 2 indivíduos sorteados de uma comunidade pertencerem à mesma espécie. É considerado um índice robustro e significartivo, captura bem a variação das distribuições de abundância e ele estabiliza com menores tamanhos amostrais. Porém é um índice que não dá muito peso à espécies raras, portanto inadequado para comunidades com muitas espécies raras (geralmente o caso dos ambientes tropicais).

No R podemos usar a função diversity do pacote vegan, ou diversityresult do pacote BiodiversityR (que é mais abrangente). Abaixo, colocando a matriz de espécies obtemos um índice para cada local:

# Simpson 1-D
diversity.simp <- diversity(dune, index='simpson')
diversity.simp

# inverso de simpson:
diversity.invsimp <- diversity(dune, index='invsimpson')
diversity.invsimp

# usando pacote BiodiversityR
Diversity.2 <- diversityresult(dune, index='Simpson' , method='s')
Diversity.2

Mas também podemos obter o índica para todos os locais agrupados:

divers.all <- diversityresult(dune, index = "Simpson", method='all')
divers.all

Índice de Shannon

Um outro índice frequentemente utilizado e que possui essencialmente as mesmas propriedades é o de Shannon, \(H\).

\[ H = -\sum_{i=1}^S P_i lnP_i \]

Cuja equabilidade é:

\[J = \frac{H} {H_{max}} = \frac{-\sum_{i=1}^S P_i lnP_i} {lnS} \]

onde \(ln\) significa o logarítimo de base natural \(e\).

Cabe lembrar que os índices descritos acima levam em consideração a riqueza e a equabilidade e se aplicam a amostras da comunidade. As premissas desses índices são de que a comunidade é infinitamente grande (amostras grandes e representativas), e que os indivíduos são amostrados aleatoriamente.

Para o índice de Shannon, as características mais atrativas são:

• Índice mais utilizado na literatura

• Geralmente valores entre 1.5 e 3.5 (raramente acima de 5)

• É sensível à espécies raras

• É sensível a variações nas abundâncias

Para calcular o índice de Shannon, usamos as mesmas funções no R, mudando apenas o argumento index:

diversity.sha <- diversity(dune, index='shannon')
diversity.sha

diversity.sha2 <- diversityresult(dune, index='Shannon', method='s')
diversity.sha2


# Equabilidade de Shannon:
diversity.shaE <- diversityresult(dune, index='Jevenness', method='s')
diversity.shaE

Índice de Dominância de Berger-Parker

Este índice é geralmente recomendado por ser intuitivamente simples e fácil de calcular e por significado biológico:

\[ d = \frac{N_{max}}{N} \]

onde \(N_{max}\) é o número de indivíduos da espécie mais abundante e \(N\) é o número total de indivíduos da comunidade.

Porém é um índice problemático apra comunidades com poucas espécies (S<15 espécies).

diversity.ber <- diversityresult(dune, index='Berger', method='s')
diversity.ber

Perfil de diversidade de Rényi

Perfis de diversidade de Rényi são curvas que provém informação sobre a riqueza e equabilidade das espécies. Esta é uma técnica de ordenação da diversidade. O formato da curva é o indicativo da equabilidade. Uma linha horizontal indica que todas as espécies tem o mesmo número de indivíduos. Quanto menos horizontal um perfil é, menos equitativa é a distribuição das espécies. O início da curva à esquerda é um indicativo da riqueza de espécies. Perfis que iniciam em um alto nível tem maior riqueza.

Uma vantagem do perfi de Rényi é que os locais podem ser facilmente ordenados da maior para a menor diversidade. Se o perfil de um local está sempre acima de outro, então este é mais diverso que o segundo. Mas se os perfis se cruzam, é provável que um local tenha mais riqueza, enquanto outro mais equabilidade.

Para obter mais detalhes de como o perfi de Rényi é calculado veja o Capítulo 5 do manual do pacote BidiversityR (esse manual é muito bom para todo o tipo de análise de diversidade).

Os índices de diversidade são um método mais compacto de comparar diversidades. Porém, um índice de diversidade frequentemente não irá prover a informação suficiente para ordenar locais da maior para menor diversidade. Apenas técnicas de ordenação da diversidade, como o perfil de Rényi, nos permitirá concluir que um local é mais diverso que o outro.

Vamos ver o perfil de Rényi para os dados das comunidades hipotéticas descritas na tabela acima

abcd <- matrix(c(1,1,4,2,1,1,2,2,1,1,1,2,0,1,0,0,0,1,0,0), nrow=4, 
               dimnames=list(c("A", "B","C","D"), c("sp1","sp2","sp3","sp4","sp5")))
abcd

reu <- renyiresult(abcd, method='s')
renyiplot(reu, legend=F)

Com este perfil podemos ver que três locais tem perfis horizontais, o que indica que as espécies estão igualmente distribuídas por local (equabilidade máxima). O local C tem a curva declinando da esquerda para a direita indicando que a distribuição das espécies nesse local é desigual (possui dominância). O perfil que começa mais alto no gráfico indica maior riqueza de espécies. Concluimos que a ordem de diversidade nos locais é \(B > A = D > C\).

Dado que a forma do perfil de diversidade de Rényi é influenciado pela equabilidade, podemos comparar a equabilidade entre os locais olhando para a forma das curvas. O Pefil de equabilidade de Rényi é um método mais direto de comparar equabilidade. Vejamos o gráfico do exemplo acima:

renyiplot(reu, legend=F, evenness = T)

A forma de interpretar o perfil de equabilidade é igual à do perfil de diversidade. Um local de maior equabilidade terá um perfil acima das áreas de menor equabilidade. Pefils que se cruzam indicam que não há uma ordenação de equabilidade possível para estas comunidades.

Renyi.1 <- renyiresult(dune)
Renyi.1
#gráfico do perfil de diversidade
renyiplot(Renyi.1,legend=F)

#gráfico do perfil de equitabilidade
renyiplot(Renyi.1, evenness=TRUE, legend=F)

# Calculando o perfil para cada local separadamente:
Renyi.2 <- renyiresult(dune, method='s')
Renyi.2
renyiplot(Renyi.2, legend=F)
renyiplot(Renyi.2, evenness=TRUE, legend=F)

Renyi.3 <- renyicomp(dune, y = dune.env, factor = 'Management',
                     permutations = 100, legend = F)

Renyi.3

Distância Ecológica

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.

# 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")

#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)

#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)

vegdist(mat, method = 'bray')

Análise da distância ecológica por Ordenação

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).

Escalonamento multidimensional não paramétrico (NMDS)

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)

Análise de componentes principais (PCA)

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).

screeplot(pca1)
points(bstick(pca1), col="red", type='o')

# também compara os eixos com o bstick
PCAsignificance(pca1) # pacote BidiversityR

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')

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

cleanplot.pca(pca1)

# 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)