Cantinho do R
Trazendo os resultados do EstimateS pro R
Porque o vegan nem sempre faz o que eu quero…
Por Marcos Vital, do LEQ-UFAL
Este material foi construído com a ajuda de muitas pessoas que acreditam no LEQ e em Ciência Livre. Muito obrigado!
Para mais material, visite o Cantinho do R
Então, se você trabalhar com biodiversidade e já teve que construir uma curva de acúmulo de espécies por rarefação ou calcular um estimador de riqueza, aposto que conhece o bom e velho EStimateS, programa gratuito muito bacana mantido pelo Robert Colwell.
O programa do Colwell é excelente. Dá um trabalhinho às vezes para se aprender como funciona a entrada de dados, mas a documentação dele, com um manual lá na página principal, é muito bem feita e dá todos os detalhes (além de arquivos de exemplo e tudo mais). Mas uma coisa que chateia muita gente é que o EstimateS apenas calcula as coisas, mas não gera, sozinho, gráficos. Isso você tem que fazer sozinho, pegando os resultados do programa e levando para o seu ambiente de gráficos favorito.
E aí, por que não usar o R pra isso? ;)
A proposta deste script é criar uma função que leia os resultados gerados pelo EstimateS, leve para o R e crie um gráfico para você. Prático, né?
Mas peraí, por que não usar o R pra tudo, incluindo fazer os cálculos?
Pois é, por que?!
Na verdade eu diria que este é o melhor caminho. Mas a questão é que nem sempre ele é trivial. Por exemplo, se usamos o pacote vegan para isso, teremos algumas diferenças, e algumas coisas que o EstimateS faz, ele não faz. Eu recomendo fortemente conhecer o pacote iNEXT, mantido pela Anne Chao, pois ele vai fazer várias coisas além. Mas de vez em quando, ainda assim, você pode estar naquela situação de já ter uma análise pronta, tudo feito, e quer simplesmente fazer o gráfico. E como o EStimateS é um programinha gratuito e super tradicional na área, refazer tudo no R até seria interessante, mas não seria algo absolutamente necessário.
Então bora lá?
Ah sim, antes que eu me esqueça: estou escrevendo este script partindo do pressuposto que você já conhece o EstimateS, sabe usá-lo e leu com atenção as instruções do Guia do Usuário. Tendo isso em mente, podemos seguir adiante.
Lendo o arquivo do EstimateS
O output do EstimateS (quando o programa é usado para gerar estimadores, índices e a curva de acúmulo) segue um formatinho padrão, o que é a coisa mais legal que poderíamos esperar. Afinal, se o formato é sempre igual, fica mais fácil padronizar a leitura.
E mais: ele gera o resultado em formato txt, o que é ótimo pra gente ler diretão no R. Então vamos lá.
Pra começar, eu vou usar a versão mais atual do EStimateS (no momento em que vos escrevo, é a 9.1, de Agosto de 2016) para processar um dos dados de exemplos do próprio programa, de forma que qualquer pessoa pode reproduzir em casa (os dados de exemplo vem junto com o EstimateS quando você o instala no seu computador).
Eu fiz o seguinte: carreguei o Single-Sample-Based-Example.txt no EStimateS (basta mandar ler o arquivo de dados, e usar todas as opções iniciais, pois ele carrega diretão) e rodei o programa sem mexer em qualquer configuração do Diversity Settings. Foi tipo uma rodada de exemplo básica, sem mexer em nada, certo?
O arquivo com o resultado foi este aqui, pode baixar se quiser. Mas se você preferir, pode rodar o EstimateS e gerar o seu próprio resultado, beleza?
Bom, o resultado do EstimateS vem sempre em um arquivo com um cabeçalho que ocupa duas linhas; depois dele segue uma linha em branco e aí, finalmente, os resultados. E aí tem um detalhe importante: os nomes das colunas contém espaços.
Sabendo disso, vamos lá!
dados.estimates<-read.table("resultado_EstimateS.txt", skip=3, header=T, sep="\t")Aqui usamos a boa e velha
read.table, já que o arquivo é um txt;O argumento mágico da vez é
skip, que usamos para pedir para o R ignorar as três primeiras linhas do arquivo;Uma coisa importante é usar o
sep, que indicou para o R que a separação das colunas é dada apenas por tabulações. Em consequência disso, o R substituí os espaços por pontos nos nomes das variáveis.
Vamos conferir?
str(dados.estimates)## 'data.frame': 11 obs. of 46 variables:
## $ Samples : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Individuals..computed. : num 29.9 59.8 89.7 119.6 149.6 ...
## $ S.est. : num 8.36 12.73 15.87 18.24 20.12 ...
## $ S.est..95..CI.Lower.Bound : num 5.66 9.75 12.82 15.24 17.2 ...
## $ S.est..95..CI.Upper.Bound : num 11.1 15.7 18.9 21.2 23 ...
## $ S.est..SD : num 1.38 1.52 1.56 1.53 1.49 1.45 1.41 1.38 1.36 1.36 ...
## $ S.Mean..runs. : num 8.64 12.5 16.13 18.33 20.5 ...
## $ Singletons.Mean : num 2.99 4.12 5.02 5.22 5.52 5.53 5.24 5.08 4.92 4.59 ...
## $ Singletons.SD..runs. : num 1.51 2.21 2.15 1.93 1.89 1.96 1.71 1.54 1.39 1.02 ...
## $ Doubletons.Mean : num 2.12 2.98 3.72 4.15 4.51 4.54 4.54 4.55 4.62 4.6 ...
## $ Doubletons.SD..runs. : num 1.42 1.48 1.71 1.75 1.62 1.63 1.67 1.43 1.32 1.08 ...
## $ Uniques.Mean : num 8.64 8.6 9.69 9.48 9.6 9.22 8.57 8.03 7.52 6.88 ...
## $ Uniques.SD..runs. : num 2.71 3.17 2.63 2.37 2.1 1.9 1.84 1.81 1.67 1.21 ...
## $ Duplicates.Mean : num 0 3.9 3.82 5 5.29 5.56 5.95 6.23 6.61 7.17 ...
## $ Duplicates.SD..runs. : num 0 1.39 1.76 2.07 2.13 1.92 1.71 1.54 1.46 1.25 ...
## $ ACE.Mean : num 12.3 18.5 22.4 23.9 26.1 ...
## $ ACE..SD..runs. : num 4.26 8.38 6.62 5.41 5.36 5.11 4.07 2.82 2.17 1.22 ...
## $ ICE.Mean : num 8.64 41.48 36.72 33.24 33.76 ...
## $ ICE.SD..runs. : num 2.71 23.52 15.81 10.41 7.34 ...
## $ Chao.1.Mean : num 10.1 15.2 19.3 21.3 23.6 ...
## $ Chao.1.95..CI.Lower.Bound : num 8.84 12.98 16.72 18.91 21.07 ...
## $ Chao.1.95..CI.Upper.Bound : num 20.8 29 34.8 36 38.2 ...
## $ Chao.1.SD..analytical. : num 2.23 3.09 3.52 3.34 3.38 3.45 3.23 2.95 2.73 2.38 ...
## $ Chao.2.Mean : num 8.64 16.72 23.84 24.92 27.36 ...
## $ Chao.2.95..CI.Lower.Bound : num 13.7 13.4 18 19.9 22.1 ...
## $ Chao.2.95..CI.Upper.Bound : num 13.7 32 49.2 46.8 49.9 ...
## $ Chao.2.SD..analytical. : num 0 3.84 6.58 5.67 5.86 5.34 4.43 3.84 3.34 2.68 ...
## $ Jack.1.Mean : num 8.64 16.8 22.59 25.44 28.18 ...
## $ Jack.1.SD..analytical. : num 0 1.73 2.83 3.09 3.25 3.15 2.96 2.75 2.51 2.34 ...
## $ Jack.2.Mean : num 0 16.8 25.2 28.5 31.6 ...
## $ Jack.2.SD..runs. : num 0 4.74 5.03 5.24 4.53 3.96 4.06 4.21 3.92 3.1 ...
## $ Bootstrap.Mean : num 8.64 14.65 19.14 21.65 24.09 ...
## $ Bootstrap..SD..runs. : num 2.71 3.99 3.42 3.44 2.59 2.02 1.95 1.75 1.32 0.83 ...
## $ MMRuns.Mean : num 0 25.6 35 68.3 51.2 ...
## $ MMMeans..1.run. : num 0 26.6 28.4 29.7 30.6 ...
## $ Cole.Rarefaction : num 11.1 15.3 18.1 20 21.6 ...
## $ Cole.SD..analytical. : num 1.98 2.01 1.93 1.8 1.65 1.48 1.3 1.1 0.88 0.61 ...
## $ Alpha.Mean : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Alpha.SD..analytical. : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Shannon.Mean : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Shannon.SD..runs. : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Shannon.Exponential.Mean : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Shannon.Exponential.SD..runs.: int 0 0 0 0 0 0 0 0 0 0 ...
## $ Simpson.Inv.Mean : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Simpson.Inv.SD..runs. : int 0 0 0 0 0 0 0 0 0 0 ...
## $ X : logi NA NA NA NA NA NA ...Perfeito! :D
No final apareceu uma variável X, que eu confeço não saber direito de onde veio. Talvez seja algum detalhe do arquivo de saída do EstimateS que eu não tenha percebido. Mas note que ele não altera em nada o restante dos dados, e nós temos agora todos os resultados do EstimateS dentro do R!
Os nomes das variáveis ficam com aqueles pontinhos todos porque o EstimateS gera vários espaços nos nomes das colunas, e o read.table os substituiu por pontos. Você pode eliminá-los se quiser (veja aquele recente script dos pokemons para um exemplo, usando a função gsub) mas é mais uma questão estética, não vai nos afetar na prática.
Pegando as variáveis que nos interessam
Aqui a coisa é simples: vamos criar objetos com as variáveis que nos interessam, e deixá-las prontinhas para criarmos nossos gráficos. Vou aproveitar o arcabouço criado naquele nosso script de curvas de acúmulo do vegan, e apenas modificar um pouquinho as coisas para este novo fim.
Neste exemplo vou usar os dados da curva de acúmulo por rarefação mesmo. Mas você pode, se quiser, trabalhar com os estimadores, índices ou o que quiser. Os dados estão aí, basta se divertir com eles!
amostras<-dados.estimates$Samples
riqueza<-dados.estimates$S.est.
erro<-dados.estimates$S.est..SDOlha que moleza, foi fácil, não foi?
Agora vou aproveitar a função que criamos para o vegan, vai ficar assim:
#Achando o valor máximo e mínimo:
y.mínimo<-min(riqueza)
y.máximo<-max(riqueza)
#Estabelecendo valores um pouco abaixo do mínimo e um pouco acima do máximo:
y.baixo<-y.mínimo*0.80
y.alto<-y.máximo+(y.mínimo*0.20)
#Salvando tudo num objeto que vai guardar estes limites:
limites.y<-c(y.baixo, y.alto)
###
#Criando a curva:
plot(riqueza~amostras, type="l", ylim=limites.y, las=1, xlab="Esforço amostral", ylab="Riqueza de espécies")
#Adicionando as barras de erro:
arrows(amostras, riqueza-erro, amostras, riqueza+erro, angle=90, code=3, length=0.05)
#E adicionando os pontos:
points(riqueza~amostras, pch=21, bg="white")
Yahoooooooo, deu certo!!!! :D
Adicionando um estimador
Legal, conseguimos fazer tudo. Mas podemos aproveitar que estamos usando os resultados do EstimateS e dar aquela espiada nos estimadores de riqueza, não é? Então, como um exemplo prático, eu vou fazer o seguinte: pegar o resultado do estimador Jacknife de primeira ordem e adicionar no gráfico como um ponto de comparação. Bora?
Achando o resultado que queremos
Na tabela do EstimateS é toda uma série de dados para tudo que ele calcula, inclusive os estimadores. Mas no caso deles, dificilmente precisamos da curva completa, pois o mais provável é que vamos querer comparar a riqueza observada com a estimada, usando apenas a estimativa final (com todos os dados) gerada pelo estimador. Então vamos lá. Se quiser, rode de novo o str lá de cima pra relembrar, ok?
O estimador Jacknife está na coluna Jack.1.Mean, e seu erro padrão logo ao lado, em Jack.1.SD..analytical.. Se ficar em dúvida, consulte o manual do EstimateS, belezinha? ;)
Nós queremos apenas o último valor, isto é, que está na última linha. No nosso caso isso significa a linha 11. Mas nós podemos usar uma forma mais genérica, de forma que você pode usar o mesmo comando paraum resultado com mais ou menos linhas. Vamos lá:
#Podemos ver o número de linhas e colunas assim:
dim(dados.estimates)## [1] 11 46#E se quisermos apenas o número delinhas, tá aqui:
dim(dados.estimates)[1]## [1] 11#Então basta dizer que:
linha.final<-dim(dados.estimates)[1]
#De forma que o valor de Jacknife é:
estimador<-dados.estimates$Jack.1.Mean[linha.final]
#E o erro padrão é:
erro.estimador<-dados.estimates$Jack.1.SD..analytical.[linha.final]
#Olhando os valores:
estimador## [1] 31.45erro.estimador## [1] 2.07Maravilha!
Adicionando no gráfico
Vamos ter que rodar o gráfico de novo, pois desta vez queremos dar uma leve esticada no tamanho dos eixos, vocês já vão sacar por que. E aí, com ele aberto, podemos adicionar o estimador. A coisa toda vai ficar assim:
#Alterando o tamanho do eixo y, para caber o estimador:
limites.y<-c(y.baixo, estimador+erro.estimador)
#Criando a curva (agora aumentando os eixos):
plot(riqueza~amostras, type="l", ylim=limites.y, las=1, xlab="Esforço amostral", ylab="Riqueza de espécies", xlim=c(0, max(amostras)+2))
#Adicionando as barras de erro:
arrows(amostras, riqueza-erro, amostras, riqueza+erro, angle=90, code=3, length=0.05)
#E adicionando os pontos:
points(riqueza~amostras, pch=21, bg="white")
##
#Agora adicionando o ponto do estimador:
points(max(amostras)+1, estimador, pch=16)
#E adicionando a barrinha de erro do estimador:
arrows(max(amostras)+1, estimador-erro.estimador, max(amostras)+1, estimador+erro.estimador, angle=90, code=3, length=0.05)
Aê, olha só que boniteza! :D
Empacotando tudo em uma função, pode ser?
Bom, chegamos no final, e achei que podia ser legal automatizar a parada toda em uma funçãozinha esperta. Nada sofisticado e cheio de opções: só uma função pra ler o resultado do EstimateS e fazer o gráfico que fizemos ali em cima. Se você quiser alterar o que ela mostra (na curva e no ponto com o estimador) vai precisar abrir e mexer nela por dentro, beleza? Mas tenho certeza que você consegue. ;)
Eu vou salvar tudo numa função e deixar prontinha pra vocês. Basta baixar daqui e usar. Qualquer coisa abra pra dar uma espiada, mas não tem nada de muito novo, os procedimentos da função são mais ou menos o que fizemos aqui por etapas.
A função vai usar apenas um argumento, o nome do arquivo do EstimateS (entre aspas e com a extensão, belezinha?). Ela vai ler os dados, mostrar os valores de riqueza e do estimador, e fazer o nosso gráfico. Para testar, vou pegar um arquivo novo do EstimateS, o resultado do arquivo de exemplo Seedbank.txt
O arquivo vocês podem pegar aqui, basta baixar; mas se preferirem, podem seguir adiante e rodar o EstimateS de novo, que é bem mais legal. ;)
Lá vai:
source("curvaestimates.R")
curva.estimates("resultado estimates seedbank.txt")## [1] "Valores de riqueza"
## [1] 3.81 6.62 8.80 10.57 12.06 13.34 14.48 15.50 16.43 17.29 18.08
## [12] 18.82 19.51 20.15 20.76 21.32 21.86 22.36 22.84 23.29 23.72 24.12
## [23] 24.51 24.87 25.22 25.55 25.87 26.17 26.45 26.72 26.98 27.23 27.47
## [34] 27.70 27.91 28.12 28.32 28.51 28.70 28.87 29.04 29.20 29.36 29.51
## [45] 29.66 29.80 29.93 30.06 30.18 30.30 30.42 30.53 30.64 30.74 30.85
## [56] 30.94 31.04 31.13 31.22 31.31 31.39 31.47 31.55 31.63 31.70 31.77
## [67] 31.84 31.91 31.98 32.04 32.10 32.16 32.22 32.28 32.34 32.39 32.45
## [78] 32.50 32.55 32.60 32.65 32.70 32.74 32.79 32.83 32.88 32.92 32.96
## [89] 33.00 33.04 33.08 33.12 33.16 33.20 33.23 33.27 33.30 33.34 33.37
## [100] 33.41 33.44 33.47 33.50 33.53 33.56 33.59 33.62 33.66 33.68 33.71
## [111] 33.74 33.77 33.79 33.82 33.85 33.87 33.90 33.92 33.95 33.98 34.00
## [1] "Estimador e seu erro padrão"
## [1] 36.98 1.70
Yeah, deu certin!
Bom, pessoal, é isso aí. Neste exemplo eu não entrei nas minúncias das coisas, mas isso acontece porque ele tem elementos de vários tutoriais passados. Qualquer coisa volte lá no Cantinho do R e veja um pouco o material, beleza? Pode ajudar a entender tudo direitinho.
Então tá bão! Um abração, e até a próxima!
Este documento faz parte do material que disponibilizo no meu blog, o “Cantinho do R”, e foi foi criado utilizando o software R, o programa RStudio e a linguagem Markdown. Para saber mais sobre eles, acesse:
Quer comentar ou perguntar? Faça isso no blog, aqui
