* Pós-graduação em Saúde Pública pela Faculdade de Saúde Pública da Universidade de São Paulo (FSP/USP) * aralphms@usp.br

Ao se fazer um experimento ou tomar um conjunto de observações na natureza, busca-se inferir sobre os padrões que podem ter gerado os dados obtidos. A inferência estatística permite que o pesquisador opte por refutar ou aceitar hipóteses a partir das evidências geradas pelos dados (Hilborn & Mangel 1997, Lewin-Koh et al. 2004). Ao longo do último século, a confrontação de hipóteses, que constrói o conhecimento científico, tem se pautado principalmente na abordagem frequentista. Nesta, os dados obtidos são confrontados contra uma hipótese nula que tenta explicar da forma mais simples os padrões observados nos dados, isto significa atribuir a variação nos dados à aleatoriedade ou ao erro de medidas (Gotelli & Elison, 2011). Assim, pergunta-se com que frequência seriam encontrados valores iguais ou mais extremos que os observados se tomássemos infinitas amostras ou repetíssemos infinitos experimentos, tendo como base uma distribuição de probabilidades conhecida (como a Gaussiana, por exemplo). Caso a probabilidade seja baixa (em geral menor que 0,05) rejeita-se a hipótese nula (Hilborn & Mangel 1997).

A inferência estatística por verossimilhança é uma abordagem que tem sido usada como alternativa à abordagem frequentista (Hobbs & Hilborn, 2006). Em determinadas situações é possível gerar duas ou mais hipóteses que possam explicar o comportamento de uma variável aleatória X qualquer. Ao se tomar uma observação desta variável pergunta-se: qual a probabilidade desta observação ter sido obtida dada a hipótese A? Fazendo o mesmo em relação a hipótese B ou para as demais hipóteses geradas. A Lei da Verossimilhança afirma que a observação (X = x) é uma evidência que favorece a hipótese A sobre a hipótese B se e somente se “PA(x)>PB(x)”. Ao se tomar diversas observações de uma variável aleatória a probabilidade de se obter a amostra é igual ao produto (multiplicação) das probabilidades de cada observação dada a hipótese de interesse (A, B,…). Portanto, na inferência por verossimilhança confrontam-se hipóteses e define-se qual é a mais plausível ou verossímil dado os dados observados, sempre de forma relativa, ou seja, uma em relação à outra, não sendo a hipótese mais plausível necessariamente a melhor aproximação da realidade, mas sim a hipótese cujo os dados geraram a maior evidência a favor (Batista, 2009).

Para um melhor entendimento, a probabilidade de uma observação qualquer está condicionada aos parâmetros (θ) de uma função que descreve uma distribuição de probabilidades: f (x)=P(X=x|θ). Em outras palavras, a uma observação será atribuído um valor de probabilidade P(X=x) que é dependente do valor dos parâmetros (θ) da função probabilística f(x), a exemplo de μ e σ que são parâmetros da distribuição Gaussiana e atribuem probabilidades a uma determinada observação. Na função de verossimilhança inverte-se esta lógica, agora os parâmetros da função são condicionados aos dados:

L (θ|X=x)

ou

L (θ|X=x1) ∙L (θ|X=x2) ∙ … ∙L (θ|X=xn) para múltiplas observações.

A Estimativa de Máxima Verossimilhança (MLE – Maximum Likelihood Estimate) é obtida pelos valores dos parâmetros da distribuição que melhor se ajustam aos dados observados. Para facilitar os cálculos utiliza-se a log-verossimilhança negativa, esta é obtida pelo negativo do logaritmo natural da função de verossimilhança: -ln [L (θ|X=x)] . Deste modo, evita-se trabalhar com valores muito pequenos (obtidos pela multiplicação das probabilidades de muitas observações) facilitando o tratamento matemático dos dados (Bolker, 2008; Batista, 2009).

Um pesquisador pode construir diferentes modelos que tentem explicar os padrões observados em seus dados, sendo que cada modelo é construído a partir de uma hipótese que tente elucidar a questão de interesse. Para cada modelo é possível obter o MLE para os parâmetros e posteriormente selecionar o que (ou quais) melhor explique seus dados (Burnham and Anderson, 2002). Ao selecionar um modelo o que possuímos é a informação mais próxima do modelo ou processo real que gerou nossos dados. No mundo real desconhecemos o modelo verdadeiro e os parâmetros de qualquer modelo que melhor aproximem o verdadeiro, o que possuímos são apenas modelos propostos com base em nosso conhecimento sobre o sistema. O critério de informação de Akaike (AIC – Akaike Information Criteria) leva em conta esta perda de informação que há entre um modelo verdadeiro teórico e o modelo proposto (Burnham and Anderson, 2002). Espera-se que quanto maior o número de parâmetros melhor tenda a ser o ajuste no modelo em relação àqueles com menos parâmetros. Uma vez que isto não é considerado na razão de verossimilhança o calculo do AIC tende a penalizar os modelos com maior número de parâmetros, assim:

,

sendo K o número de parâmetros do modelo.

A partir do AIC (ou seu ajuste para pequenas amostras - AICc) é possível ranquear os modelos candidatos e selecionar aquele ou aqueles que melhor expliquem seus dados.

Em minha pesquisa busco avaliar as variações nos padrões de diversidade de culicídeos em diferentes parques urbanos da cidade de São Paulo. A família Culicidae compreende insetos dípteros popularmente conhecidos como mosquitos, pernilongos, carapanãs ou muriçocas. São insetos holometábolos, com formas imaturas aquáticas e os adultos terrestres. O ciclo de vida compreende ovo, quatro estádios larvais, pupa e adultos alados. As fêmeas se alimentam do sangue de vertebrados, o que contribui para o desenvolvimento de seus ovos, mas permite que, em algumas ocasiões, veiculem patógenos entre seus hospedeiros, como os agentes causadores da dengue, da malária e da filariose, por exemplo (Forattini, 2002).

Com a intensificação do processo de urbanização e as altas taxas de crescimento populacional, a cidade de São Paulo sofreu ao longo do século XX um forte processo de descentralização, com o deslocamento e assentamento da população, principalmente dos mais pobres, em áreas periféricas, promovendo gradualmente a supressão da cobertura vegetal existente (ISA, 2008). Atualmente, uma grande parte das áreas verdes remanescentes no município situa-se dentro de parques públicos e unidades de conservação.

Os parques urbanos podem ser bastante distintos em suas características, indo desde pequenas áreas contendo jardins e bosques até grandes unidades que preservam importantes fragmentos de mata natural. Esta heterogeneidade de ambientes pode refletir-se em diferenças consideráveis na riqueza, abundância e composição de mosquitos. Uma das questões que norteiam minha pesquisa é saber como os padrões de diversidade de culicídeos variam entre os parques da cidade e quais fatores podem influenciar tais variações?

A diversidade beta (β) é a extensão com que as composições de espécies em assembleias de duas ou mais unidades espaciais diferem entre si. Existe uma série de métodos para medir a diversidade β, entre os mais fáceis e intuitivos há os coeficiente de similaridade/dissimilaridade, como os conhecidos índices de Jaccard ou de Sorensen (Magurran, 2004). Estas medidas de diversidade podem refletir dois diferentes fenômenos, a rotatividade (turnover) e o aninhamento (nestedness), que resultam da substituição de espécies e da perda de espécies, respectivamente (Baselga, 2010).

O índice de dissimilaridade de Sorensen (Sørensen, 1948) entre duas assembleias é obtido por:

βsor = b+c/2a+b+c

Onde, a é o número de espécies comum aos dois locais, b é o número de espécies que ocorre no primeiro mas não no segundo local e c é o número de espécies que ocorre no segundo mas não no primeiro local. Este índice incorpora tanto a rotatividade quanto diferenças na riqueza de espécies (Koleff et al., 2003). Por sua vez, o índice de Simpson (Simpson, 1943) é capaz de descrever a rotatividade sem a influência dos gradientes de riqueza. Este é obtido por:

βsimp = min(b,c)/a+min(b,c)

Onde, a,b e c são as mesmas variáveis definidas para o índice de dissimilaridade de Sorensen.

Baselga (2010) mostrou que uma simples subtração do primeiro índice pelo segundo gera o índice de aninhamento (β nestedness), que reflete o quanto assembleias menos ricas são compostas por subconjunto de espécies das assembleias mais ricas, logo:

βnes = (b+c/2a+b+c)-( min(b,c)/a+min(b,c))

Como exemplo, vamos criar um cenário hipotético no R mostrando como os padrões de rotatividade e aninhamento podem estar associados a diferentes variáveis preditoras. Para tanto, utilizaremos a função de verossimilhança e o critério AIC para seleção dos melhores modelos dentre os propostos.

Primeiro simularemos a localização de cinco parques urbanos. Supondo que desconfiemos que a distância e a diferença entre os tamanhos das áreas sejam boas preditoras para a variação na composição de espécies, vamos mensura-las.

px=c(6,8.5,4,6.5,8)
py=c(7,10,12,10,6.5)
rotulos=c("p1","p2","p3","p4","p5")
plot(px,py,ylim=c(6,13),xlim=c(3,10))
text(y=py-0.3,x=px+0.0,labels=rotulos)

dist.geo=dist(cbind(px,py),method="euclidean") ###distância hipotética entre as áreas de estudo
dist.geo.v=dist.geo[1:10]

tam.area=c(50,45,60,90,75)###tamanho hipotético da área (mil m²) de cada parque####
d.area=dist(tam.area,method="euclidian")
dif.area.v=(d.area[1:10])###diferença entre os tamanhos

Agora criaremos uma matriz de presença ou ausência de espécies para cada uma das áreas de estudo.

parq1=c(1,0,0,0,0,0,1,1,1,1,1,0)
parq2=c(0,0,0,1,0,0,1,1,1,1,0,0)
parq3=c(1,1,1,1,1,0,0,1,0,0,0,0)
parq4=c(1,0,1,1,1,1,1,1,1,1,0,1)
parq5=c(1,0,0,1,0,1,1,1,1,1,1,0)
tot.parq=rbind(parq1,parq2,parq3,parq4,parq5)
cont.sp=rowSums(tot.parq) ###contando o número de espécies coletadas em cada área

Vamos chamar o pacote Vegan e utilizar a função betadiver para calcular, par a par, os índices de de Sorensen (dissimilaridade) e Simpson, e após iremos subtrair o primeiro pelo segundo para obter o beta nestedness.

library(vegan)###chamando o pacote "vegan" e utilizando a função betadiver
beta.sor=1-(betadiver(tot.parq,method="sor"))##índice de dissimilaridade de Sorensen
beta.simp=betadiver(tot.parq,method="sim")###índice de Simpson
beta.nes=beta.sor-beta.simp ###"beta nestedness"(aninhamento)###

beta.sor.v=beta.sor[1:10] ###transformando em objeto da classe "numeric"
beta.simp.v=beta.simp[1:10]
beta.nes.v=beta.nes[1:10]

Com estes dados ajustaremos quatro modelos lineares: variação da composição de espécies como função da distância entre as áreas (mod2), da diferença entre os tamanhos das áreas (mod3) ou de ambas as variáveis (mod4), também um modelo de ausência de efeito (mod1), Utilizaremos a princípio o índice de dissimilaridade de Sorensen como variável resposta. Para isto, chamaremos o pacote bbmle e aplicaremos a função mle2 para obter o melhor ajuste para cada modelo e ICtab para compararmos os AICc’s de cada modelo.

library(bbmle)

mod1.sor=function(m=0.5,s=0.18)-sum(dnorm(x=beta.sor.v, mean=m, sd=s,log=T))
LL.mod1.sor=mle2(mod1.sor)

mod2.sor=function(b0=-1, b1=-2,s=10)
{ m <- b0+b1*(dist.geo.v)
-sum(dnorm(x=beta.sor.v, mean=m, sd=s,log=T))}
LL.mod2.sor=mle2(mod2.sor)

mod3.sor=function(b0=0, b1=1,s=0.2)
{ m <- b0+b1*(dif.area.v)
-sum(dnorm(x=beta.sor.v, mean=m, sd=s,log=T))}
LL.mod3.sor=mle2(mod3.sor)

mod4.sor=function(b0=-1, b1=-2,b2=-1,s=10)
{ m <- b0+b1*(dif.area.v)+b2*(dist.geo.v)
-sum(dnorm(x=beta.sor.v, mean=m, sd=s,log=T))}
LL.mod4.sor=mle2(mod4.sor)

ICtab(LL.mod1.sor,LL.mod2.sor,LL.mod3.sor,LL.mod4.sor,
    nobs=10,weights=T,type="AICc",base=T,logLik=T)

O “modelo 2” é apontado como o mais plausível. No entanto é possível se testar os mesmos modelos para os outros dois índices que correspondem à rotatividade e aninhamento.

mod1.simp=function(m=0.5,s=0.2)-sum(dnorm(x=beta.simp.v, mean=m, sd=s,log=T))
LL.mod1.simp=mle2(mod1.simp)
mod1.nes=function(m=0.5,s=0.2)-sum(dnorm(x=beta.nes.v, mean=m, sd=s,log=T))
LL.mod1.nes=mle2(mod1.nes)

mod2.simp=function(b0=-1, b1=-2,s=10)
{ m <- b0+b1*(dist.geo.v)
-sum(dnorm(x=beta.simp.v, mean=m, sd=s,log=T))}
LL.mod2.simp=mle2(mod2.simp)
mod2.nes=function(b0=-1, b1=-2,s=10)
{ m <- b0+b1*(dist.geo.v)
-sum(dnorm(x=beta.nes.v, mean=m, sd=s,log=T))}
LL.mod2.nes=mle2(mod2.nes)

mod3.simp=function(b0=0, b1=1,s=0.2)
{ m <- b0+b1*(dif.area.v)
-sum(dnorm(x=beta.simp.v, mean=m, sd=s,log=T))}
LL.mod3.simp=mle2(mod3.simp)
mod3.nes=function(b0=2, b1=3,s=10)
{ m <- b0+b1*(dif.area.v)
-sum(dnorm(x=beta.nes.v, mean=m, sd=s,log=T))}
LL.mod3.nes=mle2(mod3.nes)

mod4.simp=function(b0=0, b1=1,b2=1,s=0.2)
{ m <- b0+b1*(dif.area.v)+b2*(dist.geo.v)
-sum(dnorm(x=beta.simp.v, mean=m, sd=s,log=T))}
LL.mod4.simp=mle2(mod4.simp)
mod4.nes=function(b0=0, b1=1,b2=1,s=0.2)
{ m <- b0+b1*(dif.area.v)+b2*(dist.geo.v)
-sum(dnorm(x=beta.nes.v, mean=m, sd=s,log=T))}
LL.mod4.nes=mle2(mod4.nes)

ICtab(LL.mod1.simp,LL.mod2.simp,LL.mod3.simp,LL.mod4.simp,nobs=28,weights=T,type="AICc",base=T)
ICtab(LL.mod1.nes,LL.mod2.nes,LL.mod3.nes,LL.mod4.nes,nobs=28,weights=T,type="AICc",base=T)

Para a rotatividade o modelo dois continua sendo o mais plausível, já para o aninhamento os modelos 3 e 4 são os mais (e igualmente) plausíveis. Mostrando, segundo nossos dados hipotéticos, a relação entre substituição de espécies e distância geográfica entre as áreas e uma maior relação entre o tamanho da área e perda de espécies.

Um simples plot, poderá deixar mais evidente a relação que existe entre as variáveis.

par(mfrow=c(2,2))
plot(dist.geo.v,beta.sor.v,ylab="beta Sorensen (dissimilaridade)",xlab="distância geografica(km)")
plot(dist.geo.v,beta.simp.v,ylab="beta Simpson",xlab="distância geografica (km)")
plot(dif.area.v,beta.nes.v,ylab="beta nestedness",xlab="área (diferença em m²)")
plot(px,py,ylim=c(6,13),xlim=c(3,10))
text(y=py-0.3,x=px+0.0,labels=rotulos)

Em geral, ao se estudar padrões e diversidade de culicídeos busca-se observar como estes se comportam em diferentes gradientes ambientais, de naturais até urbanos. As respostas a estas alterações podem ser traduzidas em termos de um maior risco de veiculação de patógenos a humanos. Conhecer quais os fatores ambientais que podem promover variações na riqueza, abundância e composição destes insetos pode melhorar a nossa compreensão e predição sobre o risco de veiculação de doenças a população. A abordagem estatística da função de verossimilhança e seleção de modelos ajudará a buscar entre nossas hipóteses aquela que pode ser a melhor explicação para os padrões observados na diversidade de culicídeos, dado os dados.

Baselga, A. 2010. Partitioning the turnover and nestedness components of beta diversity. Global Ecology and Biogeography. 19(1):134-143.

Batista, J.L.F. 2009. Verossimilhança e Máxima Verossimilhança.

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

Burnham, K.P. Anderson, D.R. 2002. Model Selection and Multimodel Inference: A Practical-Theoretic Approach, 2nd ed. New York, Springer-Verlag.

Forattini O.P. 2002. Culicidologia Médica, vol 2. São Paulo: EDUSP.

Hilborn, R., Mangel, M., 1997. The ecological detective: confronting models with data, Monographs in population biology. Princeton University Press, Princeton, NJ.

Hobbs, N.T. & Hilborn, R. 2006. Alternatives to statistical hypothesis testing in ecology: A guide to self-teaching. Ecological Applications: 16(1): 5-19.

ISA - Instituto Socioambiental. 2008. Parques urbanos municipais de São Paulo: subsídios para a gestão. São Paulo: Instituto Socioambiental.

Koleff, P., Gaston, K. J., & Lennon, J. J. (2003). Measuring beta diversity for presence–absence data. Journal of Animal Ecology: 72(3): 367-382.

Lewin-Koh, N., Taper, M.L., Lele, S.R. 2004. A brief tour of statistical concepts, in: Taper, M.L., Lele, S.R. (Eds.), The Nature of Scientific Evidence: Statistical, Philosophical, and Empirical Considerations. University of Chicago Press, Chicago, pp. 119–152.

Magurran AE. 2004. Measuring biological diversity. Oxford, UK: Blackwell Publishing.

Simpson, G. G. (1943). Mammals and the nature of continents. American Journal of Science, 241(1): 1-31.

Sørensen T. 1948. A method of establishing groups of equal amplitude in plant sociology based on similarity of species and its application to analyses of the vegetation on Danish commons. Biol. skr. 5:1-34.

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

Medeiros-Sousa, A.R. 2014. Função de verossimilhança e seleção de modelos no estudo dos padrões de diversidade de culicídeos. Universidade de São Paulo. url: http://cmq.esalq.usp.br/BIE5781.