Não gosto daquele gráfico com cara de banana do vegan!!!
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
Recentemente, durante uma disciplina aqui da UFAL (a Feramentas Computacionais em Ecologia), eu e os alunos estávamos estudando as curvas de acúmulo de espécie, e brincando com as funções do pacote vegan
e o método de rarefação. Até aí tudo bem, tudo muito divertido (juro, é bem legal!), mas rolou um problema… :P
É que o gráfico padrão do vegan pra essas curvas é, tadinho, super feioso. Olha ele aqui, ó:
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.4-1
data(dune)
resultado<-specaccum(dune)
plot(resultado)
Bleh!!!!! Não gostei…
Daí começamos a escrever um script bacaninha, que pudesse fazer uma curva legalzinha. A coisa começou a andar, o novo gráfico ficou super bonitinho, mas nesta disciplina o objetivo era puxar um pouco as coisas, então propus: vamos transformar nosso script em uma função, de forma que a gente possa automatizar a construção do gráfico! Bora? Bora!!!!
E fizemos. Agora, neste tutorialzinho, quero apresentar e disponibilizar a função para vocês usarem, curtirem, modificarem e aproveitarem como acharem melhor. ;)
Vamos lá. Primeiro, o nosso objetivo: criar uma função que leia o resultado da função specaccum
do vegan. Claro que você tem que sacar o que a função faz, o que não é o objetivo deste tutorial, beleza? Em outro momento posso explorar mais coisas sobre isso, mas neste momento pressuponho que você sabe usar a função.
Se vamos usar a função, temos que entender o que tem “dentro” do objeto que ela gera. O truque aqui, super útil em outras situações similares, é usar a função str
para olhar o objeto por dentro, assim, ó:
str(resultado)
## List of 7
## $ call : language specaccum(comm = dune)
## $ method : chr "exact"
## $ sites : int [1:20] 1 2 3 4 5 6 7 8 9 10 ...
## $ richness: num [1:20] 9.85 15.11 18.51 20.94 22.75 ...
## $ sd : num [1:20] 2.35 1.88 1.57 1.45 1.39 ...
## $ perm : NULL
## $ freq : Named num [1:30] 7 10 2 8 6 6 5 1 1 2 ...
## ..- attr(*, "names")= chr [1:30] "Achimill" "Agrostol" "Airaprae" "Alopgeni" ...
## - attr(*, "class")= chr "specaccum"
Olha que lindeza, né não?!
O que o R está nos mostrando são as “partes” do objeto que chamamos de resultado, e que foi gerado pela specaccum
. Os pedacinhos que nos interessão são os valores de riqueza e o erro padrão: richness
e sd
(sim, sd é o erro padrão neste caso; looooonga história, fica pra outra hora). Ah sim, e podemos pegar a medida de esforço amostral, guardada em samples
, pois servirá como nosso eixo horizontal.
Quer conferir? Fuce um pouco, é sempre saudável fuçar as coisas assim no R:
resultado$richness
## [1] 9.85000 15.11053 18.51053 20.93746 22.75432 24.14959 25.23956
## [8] 26.10355 26.79824 27.36496 27.83398 28.22755 28.56204 28.84964
## [15] 29.09959 29.31909 29.51404 29.68947 29.85000 30.00000
resultado$sd
## [1] 2.3510636 1.8763851 1.5722711 1.4469584 1.3901594 1.3530349 1.3164796
## [8] 1.2749034 1.2282010 1.1763410 1.1193437 1.0564537 0.9874094 0.9115998
## [15] 0.8286890 0.7380921 0.6333903 0.5139710 0.3570714 0.0000000
resultado$Samples
## NULL
Se quiser, confira o que tem dentro de tudo, é um aprendizado legal. Se você não sabia, agora sabe que a função str
é uma super mão na roda.
Agora conhecemos os pedacinhos, e podemos seguir adiante.
Bom, podemos começar salvando cada pedacinho num objeto. Isso não é verdadeiramente necessário (a gente pode continuar usando resultado$pedaço-que-queremos), mas acredito que fique um pouco mais fácil de acompanhar a coisa toda, e os nomes ficam menores. Então vamos lá:
#Pegando os resultados que queremos:
riqueza<-resultado$richness
desvio<-resultado$sd
amostras<-resultado$sites
Beleza, cada elemento é um objeto com nome fácil de lembrar e escrever. ;)
Então, o nosso gráfico vai começar com a curva dos valores de riqueza. Até aí tudo bem, mas tem um probleminha: como além da riqueza vamos incluir o erro padrão, o R pode acabar criando uma escala no eixo y que não alcança os valores necessários para as barras de erro, e pode acabar “comendo” um pedaço do gráfico.
Resulver isso é simples, basta a gente estabelecer os limites do gráfico com o bom e velho argumento ylim
. Mas aí tem uma coisa: eu quero que isso seja feito automaticamente, pra não ficarmos tendo que estabelecer isso manualmente toda santa vez…
Então vamos usar os próprios valores de riqueza como guia, e depois “espichar” os limites pra um pouco mais que a maior riqueza e um pouco menos que a menor riqueza. Sacou?
Fazendo beeeeem passo a passo, pra ficar fácil de entender, o esquema é:
#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)
Ufa, pronto. Agora teremos um objeto, o limites.y
que guarda a informação do tamanho do eixo y. Prático, né?
Bora, bora, bora!
Vamos fazer o gráfico em três etapas: curva, barras de erro e pontos. Lá vai:
#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-desvio, amostras, riqueza+desvio, angle=90, code=3, length=0.05)
## Warning in arrows(amostras, riqueza - desvio, amostras, riqueza + desvio, :
## zero-length arrow is of indeterminate angle and so skipped
#E adicionando os pontos:
points(riqueza~amostras, pch=21, bg="white")
Nóóóóó, ficou bonitinho, não foi não? :D (o nóóó foi o meu espírito mineiro falando)
Quer saber mais sobre a função arrows
? Veja este script aqui, nele eu explico.
Mas falta uma coisaaaaa!!!!11!um!!!onze!!!!!
hehe
Falta o seguinte: bora colocar essa paradinha toda dentro de uma função! Assim fica tudo mais fácil, a gente “chama” a função e pronto. Aqui como a função ficou, completinha:
curvacum<-function(resultado, nome.x, nome.y) {
riqueza<-resultado$richness
desvio<-resultado$sd
amostras<-resultado$sites
y.mínimo<-min(riqueza)
y.máximo<-max(riqueza)
y.baixo<-y.mínimo*0.80
y.alto<-y.máximo+(y.mínimo*0.20)
limites.y<-c(y.baixo, y.alto)
plot(riqueza~amostras, type="l", ylim=limites.y, las=1, xlab=nome.x, ylab=nome.y)
arrows(amostras, riqueza-desvio, amostras, riqueza+desvio, angle=90, code=3, length=0.05)
points(riqueza~amostras, pch=21, bg="white")
}
Do jeitinho que está, a nossa função curvacum
pede o seguinte: o resultado da função specaccum
do pacote vegan
, o nome do eixo x e o nome do eixo y. Legal, né?
Se você executar o código acima, a função agora está na memória do R, e aí você pode fazer assim (com outros dados agora):
data(BCI)
resultado2<-specaccum(BCI)
curvacum(resultado2, "Esforço amostral", "Riqueza de espécies")
Maneiríssimo, né? Mas pode ficar ainda mais prático!
Se você salvar o script com a função, pode “chamá-lo” a qualquer hora com a função source
. DAí basta você copiar o script, colar na sua pasta de trabalho do momento, usar a função e ser feliz (sem ter que incluir toda a função cada vez que usar). Bom, né?
Então você pode baixar o script prontinho aqui, e depois fazer assim, ó:
#Carrega a função pronta:
source("curvacum.R")
#Vamos usar ainda outro exemplo?
data(mite)
#Vegan, eu escolho você!
resultado3<-specaccum(mite)
#curvacum, vai!
curvacum(resultado3, "Esforço amostral", "Riqueza de espécies")
E é isso aí. Espero que vocês tenham se divertido. Ou pelo menos gstado um pouquinho. Ou pelo menos usem a função depois, sei lá.
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