BIE5781 - Modelagem Estatística para Ecologia e Recursos Naturais

• Pós-Graduação em Recursos Florestais, ESALQ-USP
• Centro de Métodos Quantitativos (CMQ)
• rafaelapnaves@gmail.com

Contextualização do projeto

As florestas amazônicas têm um potencial madeireiro enorme, estimado em 60 bi m³ (VERISSIMO et al., 2012). A exploração antigamente era feita de forma convencional (EC), onde não existia planejamento acurado das operações. A partir da década de 80 surgiu a exploração de impacto reduzido (EIR) para essas áreas, com o objetivo de tornar a exploração mais sustentável, minimizando impactos ecológicos e sociais das operações (SOGOBAL et al., 2000).

Modelos para o volume dessas áreas apontam que com a EIR, tem-se a recuperação dessa variável mais cedo que com a EC, embora seja argumentado que mesmo com EIR o ciclo de corte deva ser aumentado (atualmente o ciclo é de 30 anos) (CASTRO, 2012). Neste sentido estudos que comparam a evolução das florestas remanescentes submetidas aos diferentes tratamentos trazem a real eficácia dos métodos. O desenvolvimento de modelos de crescimento para florestas com grande diversidade deve estar voltado para as mudanças na estrutura e composição da floresta e refletir as respostas das operações silviculturais. Ou seja, esses estudos devem ser respaldados nas mudanças ocorridas na comunidade.

O objetivo do projeto será avaliar a trajetória da diversidade de comunidades arbóreas ao longo de 20 anos. Essas áreas foram submetidas a três tipos de tratamentos: sem exploração, exploração convencional e exploração de impacto reduzido. O estudo das áreas foi feito através de parcelas permanentes que foram instaladas antes de exploração e remedidas periodicamente.

O que é e como comparar comunidades?

Para estudo em comunidade é necessário delimitação temporal, espacial e taxonômica desta (TOKESHI,1993). Nesse caso serão as comunidades arbóreas (com diâmetro mínimo de inclusão) medidas em um dado ano para um dado tratamento.

Como comparar essas áreas?

Três abordagens foram definidas na comparação de comunidades (MCGILL et al., 2007). Essas abordagens variam no grau de complexidade de suas metodologias e no poder de explicar a variação dos dados observados.

A primeira abordagem é a mais difundida pela facilidade de aplicação, consiste em comparar riqueza, índices de diversidade, tabelas de fitossociologia entre outros. Apesar de fácil execução, muita informação é perdida quando se resume comunidades em poucos números, além de ser dependente da composição de espécies, ou seja, o resultado vai variar conforme a localização das parcelas em campo.

A segunda abordagem faz uso de métodos estatísticos multivariados, como análise de correspondência canônica, sendo possível incorporar variáveis preditoras (ambientais, por exemplo). Apesar de explicar melhor as comunidades, essa abordagem é de difícil execução e interpretação, além de ainda ser dependente da composição.

Uma terceira abordagem que tem nível intermediário de exigência metodológica, mas que explica melhor a variação e não depende da composição de espécie é a modelagem da distribuição das abundâncias das espécies, com sigla SAD (do inglês Species Abundance Distribution) (MCGILL et al., 2007).

Um dos poucos padrões que se mantém em quase todas as comunidades é a dominância numérica de poucas espécies, enquanto a maioria é representada por poucos indivíduos, padrão conhecido como J invertido (MCGILL et al., 2007). Logo se este é um padrão que se mantém, é possível encontrar funções que atribuam tais probabilidades aos valores de abundância (PRADO,2010).

A função de densidade de probabilidade

Um modelo estatístico para SAD é probabilístico, baseado na abundância das espécies na comunidade toda e no processo de amostragem (FISHER et al.,1943, PRESTON 1948).

A abundância na comunidade toda pode então ser descrita com uma distribuição de probabilidade f(n) que permita essa forma assimétrica de J’ invertido. Embora dados de contagem sejam do tipo discreto, uma aproximação contínua é aceitável, pois o número de indivíduos é elevado e assim facilita os cálculos. Fisher et al. (1943) Preston (1948) foram os primeiros a darem tratamento matemático a essa abordagem, sendo que o primeiro usou distribuição Gamma e o segundo Log-normal. Este modelo nos dará a probabilidade de uma dada espécie ter abundância entre dois valores, que é a área de f(n) neste intervalo.

Entretanto, dificilmente vamos medir a comunidade toda, faremos amostragens, e provavelmente o número de indivíduos de uma espécie na amostra pode não ser exatamente a proporção esperada, além de muitas espécies de baixa abundância não estarem presentes em nossa amostra. Para isso é adicionada outra função de probabilidade, que diz a chance da espécie ocorrer com abundância x na amostra, dado que no todo, ela tinha n. Para efeito de amostragem Fisher et al. (1943) e Preston (1948) adicionaram a série Poisson, que reflete a ocorrência de um evento em um intervalo.

A distribuição da abundância de espécies na amostra será então a integral do produto desses dois eventos independentes. A solução dessa integral dá a distribuição das abundâncias das espécies na amostra, dadas as distribuições na comunidade e o processo amostral (PRADO, 2010)

Conhecendo as propriedades de diferentes distribuições é possível combiná-las para melhor se ajustar aos dados.

Abaixo script de um conjunto real de dados de abundância de espécies e o gráfico da densidade empírica.

dado=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,4,4,5,5,5,5,6,6,6,7,7,7,7,8,8,10,11,11,11,11,11,12,14,18,19,21,27,28,33,33,46,51,66,79,83,88,123,135,194,226)
plot(density(dado,bw=5),xlab="Abundância",ylab="Densidade ",main="Densidade empírica SAD")

Vamos supor que seja de nosso interesse ajustar a distribuição Poisson exponencial truncada, que possui dois parâmetros (“lambda” e “a”) (PRADO,2010). A distribuição Exponencial (1) tem lambda como parâmetro. A distribuição Poisson (2) tem theta como parâmetro que representa a fração da comunidade amostrada (a) pelo número total (n). Fazendo a integral do produto dessas duas funções ficamos com h(x) (3):

Façamos a distribuição truncada no R:

fun.densidade.pexp=function(lambda,x,a)
{(lambda*a^(x-1))/(lambda+a)^x}

Para o parâmetro “a” ficaremos com 0.1, que é a fração amostrada da comunidade. Vemos na figura abaixo que o parâmetro “lambda” da distribuição muda o formato da curva da poisson exponencial truncada. Com qual lambda devemos ficar para ter melhor ajuste aos dados?

plot(density((dado-1),bw=5),xlab="Abundância",ylab="Densidade ",main="Distr_Poisson_Exponencial_Truncada")
curve(fun.densidade.pexp(x,lambda=0.001,a=0.1),add=T,col="blue",lwd=2)
curve(fun.densidade.pexp(x,lambda=0.01,a=0.1),add=T,col="green",lwd=2)
curve(fun.densidade.pexp(x,lambda=0.2,a=0.1),add=T,col="red",lwd=2)
legend(150,0.05, c("lambda=0.001", "lambda=0.01", "lambda=0.2"),cex=1, bty="n", fill=c("blue", "red", "green"))

A verossimilhança

Um cenário estocástico gerou esses dados e nós nunca teremos acesso ao modelo real, mas a partir do padrão observado podemos construir modelos que possuem parâmetros, sendo possível construir infinitos modelos (BOLKER,2008). Um modelo pode ser melhor ou pior quando comparado a outro. É o próprio dado que vai indicar, num contexto de comparação, qual modelo é mais plausível, ou seja, toda a informação necessária está contida nos dados e essa comparação é sempre vis-à-vis (BATISTA,2009).

Cada modelo atribui uma probabilidade a cada observação. A função de verossimilhança é o produto dessas probabilidades, ou seja, agora os dados são fixos e variamos os parâmetros (BOLKER,2008).

Cada modelo (ou explicação) atribui uma probabilidade diferente aos dados. A explicação mais plausível será aquela que atribuir maior probabilidade aos dados. A lei da verossimilhança diz que a razão de verossimilhança é a medida da força de evidência a favor de uma hipótese contra outra. Como o produto de probabilidades gera números muito pequenos, por conveniência matemática trabalhamos com a log verossimilhança, que transforma esse produto em soma, essa função deveria então ser maximizada para encontrar seu maior valor. Em geral, procedimentos de otimização minimizam funções, então, trabalha-se com a log-verossimilhança-negativa (-LL), para que a função possa ser minimizada (BOLKER,2008).

Vamos construir a função de -LL no R e através do gráfico veremos em qual lambda a função tem seu valor mínimo.

#função_de_verossimilhança
fun.vero=function(lambdas){ -sum(log(fun.densidade.pexp(lambda=lambdas,x=dado,a=0.1)))} 
fun.vero.vetorizada=Vectorize(fun.vero,vectorize.args=c("lambdas")) #vetorizando_a_função
lam=seq(0.005,0.01,length=100)
plot(fun.vero.vetorizada (lam)~lam,type="l",col="pink",lwd=3,xlab="lambdas",ylab="log-veros.negativa")
abline(v=lam[ fun.vero.vetorizada(lam) == min(fun.vero.vetorizada(lam) )],col="green",lwd=2)

Podemos obter esses valores da -LL minimizada através da função mle2 do pacote bbmle http://cran.r-project.org/web/packages/bbmle/index.html (BOLKER,2014):

library(bbmle)
(modelo_pexp=mle2(fun.vero,start=list(lambdas=(1/mean(dado)) )))

E como esses conceitos mudaram minha maneira de comparar as comunidades

A morte de indivíduos pelas atividades da operação florestal pode levar a redução da abundância de algumas espécies, o que pode mudar a forma da curva da SAD. A abertura de grandes clareiras pela exploração cria novos nichos que podem vir a ser ocupados pelas espécies oportunistas/pioneiras, aumentando sua abundância, que geralmente é intermediária (DORNELLAS et al., 2009), podendo alterar o parâmetro da distribuição novamente. A exploração de impacto reduzido visa tornar a operação mais sustentável, portanto, é de se esperar que cause menor morte de indivíduos que a exploração convencional.

Será possível então, estabelecer alguma relação entre os impactos causados pela exploração florestal na distribuição da abundância das espécies? Como o (s) parâmetro (s) da distribuição varia (m) ao longo do tempo? Para isso serão criados modelos onde os parâmetros não são fixos, eles variam em função de variáveis preditoras relacionadas à dinâmica da floresta e aos impactos da exploração. Por exemplo, se a distribuição de densidade probabilística escolhida for a Poisson exponencial truncada teremos modelos descritos da seguinte forma:

X=Abundância ~ Poisson-Exponencial-truncada(xi |λ=f(mortos, tempo, incremento…),a=constante)

Através da seleção de modelos veremos quais variáveis explicam a SAD. Para a seleção será usado o Critério de Informação de Akaike (AIC). O AIC mede a distância de cada modelo ao modelo ideal, sendo o melhor modelo àquele que apresentar o menor valor da AIC. Essa métrica é derivada da verossimilhança máxima e também penaliza modelos com muitos parâmetros. Desde que a diferença entre o AIC de dois modelos seja menor que dois, ambos são considerados plausíveis (BATISTA, 2009). Com essa abordagem saímos da simples comparação de índices e composição e passamos a explorar relações.

Neste ensaio é apresentada a distribuição exponencial para fins didáticos, mas grande parte dos estudos empíricos se ajustou melhor aos modelos Poisson-gama ou Poisson-log-normal (HUBELL, 2001), e parte do procedimento também pode ser feito com o uso da função fitsad do pacote sads http://cran.r-project.org/web/packages/sads/sads.pdf (PRADO e MIRANDA,2014).

BATISTA, J.L.F. 2009. Verossimilhança e Máxima Verossimilhança. Apostila 27 p.

BOLKER, B.M. 2008 Ecological Models and Data in R Princeton: Princeton University Press.

BOLKER, B. and R Development Core Team (2014). bbmle: Tools for general maximum likelihood estimation. R package version 1.0.17. http://CRAN.R-project.org/package=bbmle

CASTRO, T, N. 2012. Comparação de diferentes formas de colheita florestal na Amazônia brasileira através da modelagem da produção e do crescimento. 2012. Dissertação (Mestrado em Recursos Florestais) - Escola Superior de Agricultura Luiz de Queiroz, Fundação de Amparo à Pesquisa do Estado de São Paulo. 72 p.

DORNELAS, M.et al. 2009. Species abundance distributions reveal environmental heterogeneity in modified landscapes. Journal of Applied Ecology. 46, 666–672.

FISHER, R. A.et al 1943. The relation between the number of the species and the number of individuals in a random sample from animal population. Journal of Animal Ecology 12:42-58.

HUBBEL, S. P. 2001. The Unified Neutral Theory os Biodiversity and Biogeography. Princeton, Princeton University Press.

MCGILL, B. et al. 2007. Species abundance distributions: moving beyond single prediction theories to integration within an ecological framework. Ecology Letters, 10, 995-1015.

PRADO, P. I. 2010. Distribuições de Abundâncias de Espécies: avanços analíticos para entender um padrão básico em ecologia. Revista Ciência e Ambiente. 20 p.

PRADO,P. I, MIRANDA, M. D.(2014). sads: Maximum Likelihood Models for Species Abundance Distributions. R package version 0.1.10. http://CRAN.R-project.org/package=sads

PRESTON, F. W., 1948. The commonness and rarity of species. Ecology 29:254-283.

SABOGAL, C., SILVA, J., ZWEEDE, J., PEREIRA JÚNIOR, R., BARRETO, P., & GUERREIRO, C. 2000. Diretrizes técnicas para a exploração de impacto reduzido em operações florestais de terra firme na Amazônia brasileira. Documentos 64. Embrapa: Belém. 57p.

TOKESHI, M. 1993. Species Abundance Patterns and Community. Advances in ecological research. 24: 111.

Este ensaio é um produto de disciplina da pós-graduação da Universidade de São Paulo. Para citá-lo

Naves, R. P. Diversidade em áreas de manejo florestal na Amazônia. In: Prado , P.I & Batista, J.L.F. Modelagem Estatística para Ecologia e Recursos Naturais. Universidade de São Paulo. url: http://cmq.esalq.usp.br/BIE5781.