Forrageamento ótimo: Como predadores ativos minimizam o tempo de busca?

Bem, vamos pensar um pouco sobre outra teoria bem legal em ecologia, a de como ser um bom predador.
Nem só de garras e presas vive um predador, a maioria dos predadores terrestres é limitada pela disponibilidade de alimentos a que está submetido. Isso significa que aqueles que conseguem mais que os outros da mesma espécie deixam mais descendentes (filhotinhos, prole, como quiser chamar) que aqueles indivíduos menos aptos a encontrar, capturar e matar suas presas.

Ao longo das gerações, o comportamento de forrageamento de uma população tende a ficar mais e mais eficiente. Varias adaptações podem contribuir para diminuir:

(1) o tempo gasto caçando
(2) a energia gasta caçando

Ou ambos, essa área de estudo é conhecida como Teoria do forrageamento ótimo (optimal foraging theory).

De forma simplificada, a gente pode definir duas categorias gerais de predadores: Os que fazem busca ativa (ficam perambulando procurando o que comer) e os que sentam e esperam (sit-and-wait).

Um predador do tipo busca ativa, se move através de seu habitat e captura as presas que ele quer quando precisa conseguir comida, como uma aranha que fica andando pela casa, pulando (Salticidae é conhecida como a aranha saltadora) atrás das presas.

Um predador que senta e espera, bem ele senta e espera, se posiciona em algum ponto de vantagem e espera por uma presa inocente passar ali por perto, como uma aranha parada numa flor (Uma aranha da família das Thomizidae por exemplo).

Aranhas são legais porque apresentam os dois tipos de estrategias entre suas várias famílias, mas vamos pensar primeiro num predador de busca ativa, cuja estratégia de forrageamento ideal é adquirir mais presas no menor espaço de tempo.

Considere uma aranha saltadora então, que se move através da serrapilheira a procura de insetos, larvas e outros pequenos invertebrados para comer. Sua atividade de forrageamento a torna mais vulnerável ​​a seus próprios predadores, como algumas aves por exemplo. Então quanto mais tempo uma aranha gasta forrageamento, maior a chance de morrer durante esse processo. Assim, uma aranha de maior sucesso (aquela que terá a barriga cheia para deixar mais filhotinhos) toma decisões que lhe dão mais presas na menor quantidade de tempo, já que cada segundo desperdiçado pode significar a morte.
Ela faz isso continuamente, decidir se deve parar e comer uma presa especial ou se vai seguir em frente em busca de algo melhor. A decisão é baseada em tempo:

Quanto tempo vai demorar para encontrar outra presa? Será uma presa diferente?
Quanto tempo vai demorar para capturar essa presa em comparação com o que ele já encontrou?

Estes dois componentes que formam o tempo total de forrageamento e são chamados de tempo de busca e manuseio.

Primeiro vamos pensar no tempo de busca, bem não é difícil visualizar que ele é inversamente proporcional à abundância das presas. Quanto mais presas, menos tempo se leva para encontrar uma presa.

A nossa amiga aranha, de certa forma, controla a abundância de suas presas, controlando o número de espécies de presas mantém no cardápio, um cardápio mais restrito pode ser mais raro enquanto um cardápio mais variado pode facilitar encontrar comida pela serrapilheira onde a dona aranha vive.

Se a dona aranha come apenas larvas, sua presa é bem menos abundantes e seu tempo de acaba sendo maior. Agora se no cardápio está incluso larvas e besouros por exemplo, a busca é mais rápida, dependendo da abundância de cada uma dessas presas na serapilheira.

Porém existe um tempo de manipulação da presa, para comer, e esse tempo depende do comportamento da presa. É mais rápido capturar uma larva que um besouro, podemos imaginar isso, porque um besouro está mais consciente da presença da aranha e pode se mover mais rápido, ou ter adaptações para sua defesa.

Mas a partir da discussão que começamos podemos fazer um modelo que represente o tempo médio que leva para encontrar uma presa, e esse será mais ou menos o inverso da abundância dessa presa, então:

{TeP} = {1\over{AP}}

Onde
TeP significa tempo para encontrar uma presa,
AP é abundância dessa presa.

Temos que usar umas siglas senão a formula fica muito grande e não cabe na tela. Mas está bem simples até aqui.

Certo, mas como gostamos mais de gráficos que de formulas, vamos visualizar isso num gráfico.

01

O que acontece, estamos simulando aqui o tempo de captura em relação a a abundância da presa, onde a abundância é expressa como ocorrência por minuto. Por exemplo, uma abundância média de larvas, 100 larvas por metro quadrado significa que a aranha leva 2 minutos mais ou menos para encontrar uma larva, isso da uma taxa de meia larva por minuto. Se tivermos 200 larvas por metro quadrado, bem abundânte, podemos ter uma taxa de quase 1 larva encontrada por minuto. O Contrario, sei la, 10 lavas por metro quadrado, pode dar uma taxa de 0.1 larvas por minuto, a aranha leva 10 minutos para achar uma larva. Esse valores eu estou inventando agora, veja que eu nem coloquei valores no eixo x do gráfico.
A ideia é que estamos modelando a vida da aranha, e estamos pensando sempre numa unidade para fazer comparação, e essa unidade é presas capturadas pelo tempo gasto. Como a gente sempre ta pensando numa unidade, podemos mais para frente fazer comparações, baseadas nessa unidade. Então essa unidade é como dinheiro, dinheiro sabemos que mais dinheiro é melhor, aqui, menos tempo é melhor. Temos que lembrar disso sempre, menos tempo é melhor.

Mas como dissemos, existe ainda outro componente, o tempo de manipulação da presa, e o tempo total de forrageamento torna-se então:

{TEC} = {{1\over{AP}} + TM}

Onde:
TEC é o tempo de encontrar e capturar a presa
AP é a abundância da presa
TM é o tempo que tem que ser gasto manipulando essa presa.

Que podemos olhar num gráfico assim:

01

Basicamente, somar aquele número ali, o tempo de manipulação aumenta o tempo geral para encontrar e capturar a presa, que no final das contas só pioras as coisas para a aranha, então um menor tempo de manipulação é melhor certo? So lembrar que nossa unidade é presa por minuto, e em qualquer abundância de presa ali, um tempo de manipulação faz gastarmos mais tempo total com aquela presa.

Mas estamos falando de um cardápio simples aqui, vamos supor dois itens de cardápio então, larvas e besouros.

Sabemos que, pelo que vimos até agora, temos que ter algo assim:

{TEC_l} = {{1\over{AP_l}} + TM_l}

Colocamos ali um elezinho subescrito, para lembrar que estamos falando de larva e

{TEC_b} = {{1\over{AP_b}} + TM_b}

um bezinho subescrito para lembrar que estamos falando do besouro.

Quando uma aranha decide comer apenas larvas, apenas besouros ou os dois?

Podemos expandir nosso modelo para considerar as duas presas juntas. Considere agora o tempo de forrageamento de uma aranha que decide deixar no cardápio larvas e besouros. A primeira coisa a considerar é que a abundância total vai ser a abundância de besouros mais a abundância de larva e consequentemente o tempo médio de busca vai ser menor.

{TEC} = {1\over{AP_l + AP_b}}

Certo, neste modelo temos o tempo que uma aranha leva para encontrar e capturar uma presa, mas qual a chance dessa presa ser uma larva ou um besouro?
Estamos perguntando mais especificamente: quando uma aranha encontra uma presa aceitável (larva ou besouro), qual é a probabilidade de que ela seja uma larva?

(A probabilidade, neste caso de duas presas, de que a presa seja um besouro vai ser apenas um menos a probabilidade de que ele é uma larva)

Estamos querendo saber a abundância de larvas em relação à abundância de todas as presas possíveis, larvas e besouros:

{p_l} = {AP_l\over{AP_l + AP_b}}

Achou que ia ser algo complicado né ^^. Não, a resposta é apenas a abundância de larvas divido pela abundância total. Para o besouro vai ser:

{p_b} = {AP_b\over{AP_l + AP_b}}

Ou , como dissemos la em cima:

{p_b} = 1 - p_l

Mas para que queremos essa informação?

Porque imaginamos que o tempo de manipulação seja diferente entre besouros e larvas, é mais fácil comer uma larva que um besouro, logo o tempo de manipulação para as larvas da dieta é:

 {TM_l} = { {AP_l\over{AP_l + AP_b}} \cdot{M_l}}

E nos podemos usar a mesma ideia para calcular o tempo para o besouro

 {TM_b} = { {AP_b\over{AP_l + AP_b}} \cdot{M_b}}

Agora é so juntar tudo isso, e temos o tempo médio nossa aranha leva para encontrar e capturar uma presa qualquer do seu cardápio.

{TEC} = {  {1\over{AP_l + AP_b}} + TM_l + TM_b }

E substituímos para cada caso, e ficamos com:

{TEC} = {  {1\over{AP_l + AP_b}} + { {AP_l\over{AP_l + AP_b}} \cdot{M_l}} + { {AP_b\over{AP_l + AP_b}} \cdot{M_b}} }

Pode parecer uma formula grande, mas continuamos ainda com a mesma coisa, podemos falar que o Tempo de encontrar e capturar uma presa depende da abundância os espécies do cardápio, aqui larvas e besouros, mais o tempo de manipular elas, mas ponderamos a manipulação pela abundância da presa em questão (larva ou besouro) em relação a todos as presa que a aranha pode encontrar.

Agora vamos fazer uma tentativa aqui, uma pequena simulação.
Vamos atribuir abundâncias aos vermes e besouros, sempre pensando a abundância de um grupo em relação ao outro, em vários cenários. Ai vamos procurar em qual cenário a aranha se daria melhor. Qual a melhor escolha de cardápio?

01

Bem, como o besouro tem um tempo de manipulação muito maior, a primeira coisa que vemos é que comer somente larvas é melhor que comer somente besouros. Mas algo talvez não tão obvio, é que em baixas abundâncias, manter larvas e besouros no cardápio é melhor que somente larvas, ja que a abundância de ambos é baixo, a aranha demora para encontrar algo para comer, então mesmo demorando para capturar o besouro, vale a pena parar e comer o besouro.

Mas veja que a linha vermelha, que representa somente larvas, se cruza com a linha preta, que inclui larvas e besouros no cardápio. Isso acontece porque nem sempre incluir o besouro será uma boa estratégia. Como assim? Se o tempo que você leva capturando um besouro for maior que o tempo para encontrar e manipular uma larva, a partir dai não vale mais a pena se dar o trabalho de comer besouros, porque no tempo de capturar um a gente consegue comer uma larva. E a partir dai, conforme aumentamos a densidade das larvas, passa sempre a valer a pena comer só larvas, veja que no final do eixo x, a linha vermelha está sempre mais perto de zero, que é o que queremos, mais presas por segundo.

Agora vamos manipular um pouco a abundância dos besouros, que está como a metade da abundância das larvas.

01

Agora fica mais fácil de ver o que estávamos comentando no gráfico anterior, veja que fomos aumentando a abundância do besouro para ver o que acontecia. Colocamos a abundância igual a abundância das larvas, depois o dobro, quatro vezes mais e ainda assim as linhas se cruzam no mesmo ponto. Que é o tempo de manipulação do besouro. Se a larva tem uma abundância muito alta, sempre vai acabar valendo a pena tirar o besouro do cardápio, agora veja que quanto maior a abundância do besouro, olhe a linha azul, passa a valer muito a pena incluir ele se a abundância de larva é media ou baixa. Quando a larva tem uma abundancia muito baixo, um besouro no prato garante a refeição de todos os dias.

Bem interessante, principalmente como não existe uma escolha simples, haverá sempre um balanço entre tempo que se passa manipulando uma presa e tempo para encontrar ela. E lembre que tempo de manipulação, a gente está falando de modo geral, mas em manipulação, inclua o tempo que você demora para digerir uma a presa, certos tecidos são bem mais difíceis de digerir que outros, e se seu estomago está cheio, não da para comer mais nada.

Ja a abundância pode mudar com a temperatura ambiental, besouros e larvas são animais ectotérmicos (de sangue frio), assim todo mundo vai ficar menos ativo no frio, mais ativo no calor, e podemos imaginar que a atividade dessas presas vai influenciar nas chances delas serem encontradas, então podemos incluir isso também como parte da “Abundância”.

E além de tudo isso, ainda temos o ganho energético proporcionado por cada tipo de presa, que não foi considerado aqui. Até aqui estamos dando pesos iguais para larvas e besouros, mas talvez um besouro seja mais gordinho e suculento, e isso compensa todo os problemas para capturar ele. Mas isso vai ficar para o próximo post.

Bem ficamos por aqui, o script está no repositório recologia, e até o próximo post.

Referência:
Bernstein, R. 2003 – Population Ecology – An Introduction to Computer Simulations. John Wiley & Sons, Ltd 159pp.

#Tempo para encontrar uma presa.
tep<-function(ap) {
    return (1/ap)
}
 
sequencia<-seq(0.2,1,0.01)
 
#
plot(sequencia,tep(sequencia),pch=19,type="b",frame=F,xaxt="n",
     ylab="Tempo para encontrar uma presa",xlab="Abundância da presa no ambiente")
axis(1,at=c(0.2,0.4,0.6,0.8,1),labels=F)
mtext("Muito raro", side = 1, line = 1, at = 0.2)
mtext("Raro", side = 1, line = 1, at = 0.4)
mtext("Mais ou menos", side = 1, line = 1, at = 0.6)
mtext("Abundante", side = 1, line = 1, at = 0.8)
mtext("Muito abundante", side = 1, line = 1, at = 1)
 
#
tepm<-function(ap,tm) {
    return ((1/ap)+tm)
}
 
#
curve(tep(x),0.1,1,frame=F,ylab="Tempo para encontrar uma presa",xlab="Abundância da presa no ambiente",xaxt="n")
axis(1,at=c(0.1,0.4,0.6,0.8,1),labels=F)
mtext("Muito raro", side = 1, line = 1, at = 0.1)
mtext("Raro", side = 1, line = 1, at = 0.4)
mtext("Mais ou menos", side = 1, line = 1, at = 0.6)
mtext("Abundante", side = 1, line = 1, at = 0.8)
mtext("Muito abundante", side = 1, line = 1, at = 1)
tm<-seq(1,5,1)
for(i in tm){
    curve(tepm(x,i),col=i,add=T,lwd=2,lty=2)
}
legend("topright",col=c(1,tm),lwd=c(1,rep(2,5)),lty=c(1,rep(2,5)),legend=paste("Tempo =",c(0,tm)),bty="n")
 
#
tiw<-function(aw,hw) {
    return(1/aw + hw)
    }
#Larva
hw<-1
aw<-seq(0.000,0.03,0.0005)
#Besouro
ab<-0.5*aw
hb<-60
 
plot(aw,tiw(aw,hw),pch=19,cex=0.5,type="b",col=2,frame=F,xlab="Abundância",ylab="Tempo total manipulação",
     ylim=c(0,200),xlim=c(0,0.03))
points(aw,tiw(ab,hb),pch=19,cex=0.5,type="b",col=3)
points(aw,tiwb(aw,ab,hw,hb),type="l",col=1,lwd=1,lty=3)
 
legend("topright",col=c(2,3,1),cex=1,pch=19,lty=c(1,1,3),bty="n",legend=c("Larva","Besouro","Ambos"))
 
tiwb<- function(aw,ab,hw,hb) {
    return(1/(aw+ab)+hw*aw/(aw+ab)+hb*ab/(aw+ab))
}
#
plot(aw,tiwb(aw,ab,hw,hb),type="l",col=1,lwd=2,lty=2,frame=F,xlab="Abundância",ylab="Tempo total manipulação",
     ylim=c(0,200),xlim=c(0,0.03),main="Tempo de manipulação:\n Larvas = 1\n Besouros=60")
#
ab<-aw;
points(aw,tiwb(aw,ab,hw,hb),type="l",col=2,lwd=2,lty=2)
#
ab<-2*aw;
points(aw,tiwb(aw,ab,hw,hb),type="l",col=3,lwd=2,lty=2)
#
ab<-4*aw;
points(aw,tiwb(aw,ab,hw,hb),type="l",col=4,lwd=2,lty=2)
#
points(aw,tiw(aw,hw),type="l",col=1,lwd=1,lty=3)
legend("topright",col=c(1,2,3,4,1),lwd=c(2,2,2,2,1),lty=c(2,2,2,2,3),bty="n",title="Relações de abundância",
       legend=c("larvas=0.5*besouros","larvas=besouros","larvas=2*besouros","larvas=4*besouros","Somente larvas"))

Usando recursão para desenhar a árvore H

Esquemas simples para desenhos recursivos podem gerar figuras incrivelmente complexas.
A recursão, que ja vimos a muito tempo atrás aqui, resolvendo problemas do Rosalind é muito interessante.
O conceito de recursão é bem simples até. Mas as vezes é duro de conseguir uma intuição. Mas quando algo é dificil a gente desenha, e pra mim desenhar sempre funciona para ajudar a entender algo.

Então vamos considerar o seguinte desenho, um H.

1

Como esse desenho é feito?

Primeiro a gente inicia algumas variáveis e um plot

x<-5
y<-5
tamanho<-2
plot(x,y,type="n",xlim=c(x-1.5*tamanho,x+1.5*tamanho),ylim=c(y-.5*tamanho,y+1.5*tamanho),frame=F)

Aqui, o x e y é onde queremos desenhar, um ponto de centro para o H, e o tamanho é o tamanho que damos as partes dele. Veja como os lados do H medem 2, vão do 4 ao 6 no eixo y por exemplo.
Ok, mas o desenho foi feito com essa função.

desenheH<-function(x,y,tamanho) {
    #Calcule as coordenadas das 4 pontas do H
    x0 <- x - tamanho/2
    x1 <- x + tamanho/2
    y0 <- y - tamanho/2
    y1 <- y + tamanho/2
 
    #Desenhe 3 linhas que formam o H
    points(c(x0, x0),c(y0, y1),type="l",lty=1)    # Linha Esquerda
    points(c(x1, x1),c(y0, y1),type="l",lty=1)    # Linha Direita
    points(c(x0,  x1),c(y,  y),type="l",lty=1)    # Linha do meio
}

Veja que no fundo ela é bem simples, ela recebe o ponto central da H que a gente vai desenhar, os valores de X e Y, e o tamanho, calcula os pontos que precisamos para desenhar o H, o que a gente faz com a função points do R, usando ela como tipo linha (type=”l”).

Agora e se a gente desenhar um H em cada ponta do H inicial. Tipo isso:

2

Haha legal né?
Agora o que a gente fez aqui? Veja que assim como calculamos as pontas do H para desenhar o H, a gente pode usar esses pontos como centro de um H , chamando a função desenhaH com metade do tamanho anterior e produzir essa figura.

3

Mais legal ainda, é fazer mais uma rodada.
E como fazer isso? Usando Recursão!!!
Podemos fazer o seguinte.

desenhe<-function(n,x,y,tamanho) {
 
    if (n == 0) {
        return("Fim da Recursão")
    } else {
 
        desenheH(x, y, tamanho)
        #Sys.sleep(0.5)
 
        #Calcule as coordenadas das 4 pontas do H
        x0 = x - tamanho/2
        x1 = x + tamanho/2
        y0 = y - tamanho/2
        y1 = y + tamanho/2
 
        #Desenhe, recursivamente, 4 arvores H com metade do tamanho e de ordem n-1
        desenhe(n-1, x0, y0, tamanho/2)    #// Árvore da esquerda de baixo
        desenhe(n-1, x0, y1, tamanho/2)    #// Árvore da esquerda de cima
        desenhe(n-1, x1, y0, tamanho/2)    #// Árvore da direita de baixo
        desenhe(n-1, x1, y1, tamanho/2)    #// Árvore da direita de cima
    }
}

O que essa função vai fazer? Agora além do x e y e tamanho do H vamos entrar com mais um argumento, um n
que vai dizer quantas vezes queremos fazer o H dentro do H

4

Veja que o lance, é que a função desenhe chama a própria função desenhe. E chama ela 4 vezes para cada uma das pontas do H. Mas a gente sempre está chamando ela com n-1. Ou seja, o n vai diminuir a cada chamada, e quando chegar a zero, ele vai entrar no if (n == 0) e isso vai fazer a recursão acabar (return(“Fim da Recursão”)).

5

Como o desenho acaba sendo feito muito rápido pelo computador. A gente pode usar o comando Sys.sleep para ver o desenho sendo feito. Eu deixei o comando marcado com um #, mas se a gente tirar, haverá uma pausa de meio segundo entre cada chamado do desenhoH. Assim da para ver o desenho sendo feito. É bem legal, mas tome cuidado. Serão desenhados \sum\limits_{i=1}^n 4^{n-1} letras H na figura. Então quando a gente chama a função usando n=1, teremos 4^{1-1} e 4^{0}=1. Veja como a primeira figura tem um H.
Na segunda figura, temos n=2, então teríamos 4^1+4^0=5, cinco H, que é o que vemos la.
Assim da para saber quantos desenhos teremos na figura, o que cresce bem rápido, então não use um n muito grande, ja que vai demorar um bocado para terminar. Com n=6 ja da 5460.

6

Bem é isso ae. Ficamos por aqui, o script está aqui e no repositório recologia.
Esse exemplo eu tirei do livro Introduction to Programming in Java enquanto tento aprender java. Mas é legal replicar o código em R.
Uma coisa que ainda dava para fazer, é criar uma outra função, que usa a função plot dentro dela e depois chamar a função desenhe dentro dela. Assim para usar a função não precisaria usar o comando plot separado. Outra coisa legal é deixar a chamada da função desenheH no final da função, e ver o desenho final ser feito de traz para frente.

desenheH<-function(x,y,tamanho) {
    #Calcule as coordenadas das 4 pontas do H
    x0 <- x - tamanho/2
    x1 <- x + tamanho/2
    y0 <- y - tamanho/2
    y1 <- y + tamanho/2
 
    #Desenhe 3 linhas que formam o H
    points(c(x0, x0),c(y0, y1),type="l",lty=1)    # Linha Esquerda
    points(c(x1, x1),c(y0, y1),type="l",lty=1)    # Linha Direita
    points(c(x0,  x1),c(y,  y),type="l",lty=1)    # Linha do meio
}
 
desenhe<-function(n,x,y,tamanho) {
 
    if (n == 0) {
        return("Fim da Recursão")
    } else {
 
        desenheH(x, y, tamanho)
        #Sys.sleep(0.5)
 
        #Calcule as coordenadas das 4 pontas do H
        x0 = x - tamanho/2
        x1 = x + tamanho/2
        y0 = y - tamanho/2
        y1 = y + tamanho/2
 
        #Desenhe, recursivamente, 4 arvores H com metade do tamanho e de ordem n-1
        desenhe(n-1, x0, y0, tamanho/2)    #// Árvore da esquerda de baixo
        desenhe(n-1, x0, y1, tamanho/2)    #// Árvore da esquerda de cima
        desenhe(n-1, x1, y0, tamanho/2)    #// Árvore da direita de baixo
        desenhe(n-1, x1, y1, tamanho/2)    #// Árvore da direita de cima
    }
}
 
 
x<-5
y<-5
tamanho<-2
 
plot(x,y,type="n",xlim=c(x-1.5*tamanho,x+1.5*tamanho),ylim=c(y-.5*tamanho,y+1.5*tamanho),frame=F)
desenhe(4, x, y, tamanho)

Overdispersion em um modelo para contagens (Poisson)

Ola, após quase um mês enrolado, o blog está de volta. O ultimo host nos baniu rapidinho hehehe e ainda no meio das provas.
Mas agora vamos ser no que da, achei um novo host sugerido no Reddit, o biz.nf. Espero que pelo menos no fim de ano o blog fique ativo.

Mas chega de blablabla. Vamos continuar onde paramos no último post quando falamos de GLM usando distribuição de Poisson.
Um dos problemas, não é bem um problema, mais um característica da distribuição de poisson é que ele so tem um parâmetro, o \lambda (lambda), ele determina onde está a média e a variação é igual também a esse valor. Daqui a alguns post espero escrever sobre como é o processo de poisson, mas essas características, como falamos muitas vezes, fazem todo o sentido se pensar que quanto contamos poucas coisas, erramos por pouquinho, agora quando contamos muito erramos por muito. Porém as vezes , quando estamos pegando dados reais, muitas vezes a variação pode ser maior ou menor que isso, por exemplo se existe algum fator que não estamos levando em conta, isso faz com que a variação seja diferente de \lambda, apesar de média estar certa. Isso vai levar a gente a fazer previsões erradas, ja que se parar para pensar um momento, se a variação a mais ou a menos não for levada em conta, a gente pode ter confiança demais num valor, ou de menos.

Vamos tentar ilustrar com um exemplo.
Imagine que visitamos várias fazendas e reservas, e sempre parávamos em alguns locais e contávamos quantos coelhos víamos e ficamos com a suspeita de que coelhos são mais comuns onde a terra está arada, foi trabalhada do que em locais nativos. Então aqui esta um box-plot de 20 contagens para cada um níveis que estamos considerando (2 níveis, arado e nativo), 20 para cada, da um total de 40 amostras.

Uma tabelinha assim:

Local Contagem 1 Nativo 0 2 Nativo 2 3 Nativo 3 . . . 39 Arado 2 40 Arado 2

O que nos da um gráfico assim:

Coelhos

Certo. A primeira coisa que vem a cabeça é um glm com distribuição de Poisson, mas agora podemos ou não pensar que há uma maior variação nos dados.

Bem vamos fazer os dois usando o comando GLM do R e ver o que acontece.

> summary(glm.fit.no.OD) Call: glm(formula = C.OD ~ x, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -2.00000 -0.78339 -0.02871 0.65787 2.27670 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.6931 0.1581 4.384 1.17e-05 *** xArado 0.4220 0.2035 2.074 0.0381 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 54.580 on 39 degrees of freedom Residual deviance: 50.182 on 38 degrees of freedom AIC: 154.28 Number of Fisher Scoring iterations: 5

Esse primeiro caso ta considerando o mesmo que fizemos no post teste t de poisson. Vemos as estimativas de \lambda e um valor p significativo, legal.

> summary(glm.fit.with.OD) Call: glm(formula = C.OD ~ x, family = quasipoisson) Deviance Residuals: Min 1Q Median 3Q Max -2.00000 -0.78339 -0.02871 0.65787 2.27670 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.6931 0.1757 3.945 0.000332 *** xArado 0.4220 0.2261 1.867 0.069681 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 1.234685) Null deviance: 54.580 on 39 degrees of freedom Residual deviance: 50.182 on 38 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5

E agora, olhando o número vermelhinho ali destacado, o valor p não é mais significativo.
Mas o que mudou? Se a gente olhar as estimativas de \lambda aqui, são iguais ao modelo inicial. Mas veja que os erros(“Std. Error”) do valor de \lambda são maiores, agora estamos levando em conta uma gama maior de valores, por isso a menor certeza da diferença, e a diferença nos valores p. Se pensar em valor p como a chance de observamos novamente um valor tão extremo, claro que uma variação maior diminui nossas certezas.
Mas veja também que temos um outro valor diferente ainda, o “Dispersion parameter”. Agora é maior que 1. Não é difícil chutar que ele ta maior porque estamos levando em conta a maior variação nos dados.

Bem podemos repetir essa análise usando uma análise bayesiana com o Openbugs para melhorar nossa intuição sobre o assunto.

Então vamos olhar como fica o modelinho na linguagem bugs.

sink("overdispersion.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 10)
 tau <- 1 / (sigma * sigma)
 
# Likelihood
 for (i in 1:n) {
    C.OD[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i] + eps[i]
    eps[i] ~ dnorm(0, tau)
 }
}
",fill=TRUE)
sink()

Veja que o modelo é igual a todos que temos feito, a não ser por uma adição, que é o parâmetro de dispersão eps, que esta contabilizando a “overdispersion”.

Muito cuido na leitura a partir daqui, que eu posso ter entendido as coisas erradas, mas vamos indo.

Bem, para diferenciar os lugares temos:

C.OD[i] ~ dpois(lambda[i])
log(lambda[i]) <- alpha + beta *x[i]

Ou seja, cada local depende dos seus parâmetros locais, que é o \alpha(alpha), e o \beta(beta), lembrando que o x vai ser 0 para locais arado e 1 para nativo. Com isso construímos o \lambda.
Mas ai adicionamos a maior variação no \lambda

log(lambda[i]) <- alpha + beta *x[i] + eps[i]

Além disso tudo, temos mais alguma coisa. Que coisa? A coisa que influência a quantidade de coelho e que não conseguimos controlar, ou seja a quantidade de coelhos observada não é um processo de Poisson perfeito, temos outras coisinhas, ou coisonas acontecendo durantes as observações que não controlamos. Talvez a quantidade de alimento disponível, presença de predadores, a ação desses predadores. O mundo é mais complexo que estamos descrevendo nesse modelinho. Muitas coisas acontecendo ao mesmo tempo. So que não temos como contabilizar tudo isso, não há tempo viável para ficar vasculhando o local para disponibilidade de alimento, nem ficar vários dias la olhando o que aparece de predador, como eles agem, como são todas as interações. E tem outra, os dados ja estão coletados, não vamos pular do barco e não analisar os dados so porque não existem outras 5 milhões de variáveis disponíveis, tempo que sempre tentar fazer o melhor possível com o que temos a mão, e não reclamar do que não temos.

Mas não é por causa disso que precisamos “mentir” para o nosso modelo, podemos dizer, modelo esta acontecendo alguma coisa (essa coisa é o eps ali), que sempre acontece mas eu não tehho controle, leve isso em consideração ok? E ficamos assim:

log(lambda[i]) <- alpha + beta *x[i] + eps[i]
eps[i] ~ dnorm(0, tau)

Veja que temos um erro de distribuição normal, com média zero e um desvio padrão sigma, e a variância que chamamos de tau. É isso que tamos fazendo. Se for um processo perfeito, esperamos que esse valor seja próximo de zero, se não esperamos que esse valor contabilize essa maior, ou menor dispersão dos dados.

Mas depois disso separamos os dados numa lista no formato certo para mandar para o bugs, Fazemos uma função para iniciar valores para os parâmetros de interesse. configuramos como queremos a cadeia e pimba, temos o resultado.

> print(out, dig = 3) Inference for Bugs model at "overdispersion.txt", Current: 3 chains, each with 3000 iterations (first 1000 discarded), n.thin = 5 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 0.625 0.185 0.224 0.506 0.635 0.755 0.956 1.001 6000 beta 0.431 0.236 -0.022 0.274 0.427 0.586 0.916 1.001 6000 sigma 0.308 0.168 0.029 0.181 0.300 0.420 0.655 1.001 5300 deviance 144.011 7.084 129.600 139.000 144.600 149.600 155.700 1.001 6000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 10.3 and DIC = 154.3 DIC is an estimate of expected predictive error (lower deviance is better).

E aqui estamos, as estimativas do \lambda são bem próximas daquelas usando o comando GLM. A conclusão também, assim como o valor p não era significativo, se a gente olhar o \beta e sua distribuição, vemos que o intervalo de 95% vai de -0.022 a 0.916. Ou seja o zero, ou zero de diferença está dentro daquilo que consideramos como verdadeiro, então estamos inclinados a que tanto o locais arados como nativos tem a mesma quantidade de coelhinhos, tem um \lambda igual.

Mas agora como pensar nesse desvio, nesse “Overdispersion”?

Eu não tenho muita certeza desses gráficos, por causa do kernel que usei, mas é algo mais ou menos assim, vamos colocar esses valores num gráfico e olhar para locais nativos e arados.

Arado

Nativo

Bem veja a diferença entre como é o ajuste com e sem overdispersion, veja que o nos dois casos, apesar de não dar para ver direito para locais arado, a distribuição azul, que é quando consideramos o “overdispersion” é mais abertinha, que a linha vermelha, considerando somente a variação esperada no processo de Poisson, pode parecer pouca diferença, mas veja que fez uma grande diferença nas nossas conclusões, sendo em um caso com valor p menor que 0.05 e no outro não, ao abrir as perninhas, a distribuição fica mais sobrepostas, uma na outra, entre arado e nativo, fazendo com que começamos a levar mais em consideração a afirmação que elas são iaguais, e que a diferença que estamos vendo é ao acaso. Mas a idéia é essa, estamos assumindo que existe algo a mais na variação, correlações que não levamos em conta, algo a mais, e estamos deixando isso aparecer no modelo, e somente então tomando nossa conclusão.

Para confirmar que estamos no caminho certo, podemos gerar um conjunto de dados sem Overdispersion e ver o que acontece. Usaremos os mesmo parâmetros alpha e beta, so vamos tirar a maior variação no lambda e usar o mesmo modelo que usamos no caso acima.
Primeiro geramos os dados e vemos como fica.

exemplo2

Ai fazemos tudo denovo, juntar dados e tal e rodamos o openbugs denovo para esse novo conjunto de dados.

> print(out2, dig = 3) Inference for Bugs model at "overdispersion.txt", Current: 3 chains, each with 3000 iterations (first 1000 discarded), n.thin = 5 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 0.335 0.113 0.111 0.261 0.335 0.413 0.554 1.001 3500 beta 0.604 0.139 0.325 0.513 0.604 0.696 0.881 1.001 6000 sigma 0.171 0.101 0.009 0.095 0.161 0.238 0.387 1.120 49 deviance 397.125 8.296 377.100 392.600 399.000 403.100 408.900 1.006 410 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 10.4 and DIC = 407.5 DIC is an estimate of expected predictive error (lower deviance is better).

Agora vamos olhar o que aconteceu com o valor de sigma. Veja agora que ele é baixinho, comparado com o caso anterior. E ele vai de 0.009 a 0.387. Bem temos praticamente um zero dentro da amostra, ja que 0.009 está bem próximo de zero. Esse valor na verdade é 0.00852687, mas ta arredondado para 3 casas depois do zero, a gente sempre imprime o resultado print(out2, dig = 3), 3 dígitos de precisão.

Mas parece que essa idéa funciona, e é bem interessante para produzir modelos melhores. E se o sigma for algo perto de zero, a gente pode excluir ele do modelo e fazer um modelo mais simples. E assim começamos a tentar ver um mundo um pouco maior além dos glm convencionais. Tudo isso tem pronto no R como vimos ao usar family=quasipoison no R, então uma vez entendendo, vale a pena considerar a possibilidade de uso, dependendo do que estamos considerando.

Num próximo post, vamos estender essa ideia para os modelos com zero-inflados, que ja fizemos um post aqui a milhões de anos atras, mas agora vamos voltar nele com estatística bayesiana.

Bem ficamos por aqui, legal ter o blog devolta, espero que esse host não chute minha bunda também como os anteriores, deixei até uma propagandinha ali em cima na lateral do blog pra ver se faço algum dinheiro para pagar um servidor hehe.

O script vai estar la no repositório recologia do github além de aqui embaixo, e é isso ai, até a próxima.

### Gerando os dados
set.seed(123)
n.locais <- 20
x <- gl(n = 2, k = n.locais, labels = c("Nativo", "Arado"))
eps <- rnorm(2*n.locais, mean = 0, sd = 0.5)
lambda.OD <- exp(0.60 +(0.35*(as.numeric(x)-1) + eps) )
 
C.OD <- rpois(n = 2*n.locais, lambda = lambda.OD)
 
data.frame(Local=x,Contagem=C.OD)
 
#GRafico 1
boxplot(C.OD ~ x, col = "grey", xlab = "Uso da terra",frame=F,
ylab = "Contagem de coelhos", las = 1, ylim = c(0, max(C.OD)))
 
 
### Analise usando o R
glm.fit.no.OD <- glm(C.OD ~ x, family = poisson)
glm.fit.with.OD <- glm(C.OD ~ x, family = quasipoisson)
summary(glm.fit.no.OD)
summary(glm.fit.with.OD)
 
anova(glm.fit.no.OD, test = "Chisq")
anova(glm.fit.with.OD, test = "F")
 
### Analise com OpenBugs
 
sink("overdispersion.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 10)
 tau <- 1 / (sigma * sigma)
 
# Likelihood
 for (i in 1:n) {
    C.OD[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i] + eps[i]
    eps[i] ~ dnorm(0, tau)
 }
}
",fill=TRUE)
sink()
 
# Juntando os dados
bugs.data <- list(C.OD = C.OD, x = as.numeric(x)-1, n = length(x))
 
# Função de inicialização
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1), sigma = rlnorm(1))}
 
# Parametros a serem guardados na saida.
params <- c("alpha", "beta", "sigma")
 
# Configurando o MCMC
nc <- 3		# Número de cadeias
ni <- 3000	# Tamanho da cadeia
nb <- 1000	# burn-in, quanto do inicio da cadeia vai ser jogado fora
nt <- 5		# Thinning rate (A cada uma amostra que pegamos, jogamos 5 fora, para evitar auto-correlação)
 
# Iniciando o Gibbs sampling
library(R2OpenBUGS)
out <- bugs(data=bugs.data, inits=inits, parameters.to.save=params,model.file="overdispersion.txt", n.thin=nt,
            n.chains=nc,n.burnin=nb, n.iter=ni)
print(out, dig = 3)
 
 
#O que realmente estamos estimando.
 
#Grafico 2
hist((C.OD[which(x=='Arado')]),freq=F,main="Arado",xlab="Histograma de Contagens",ylab="Frequência da contagem")
lines(density(rpois(1000,out$summary[1,1]),adjust = 4),lwd=2,col="red",lty=2)
valores<-rep(out$summary[1,1],1000)+rnorm(1000,0,out$summary[3,1])
valores[valores<0]<-0
lines(density(rpois(1000,valores),adjust = 3),lwd=2,col="blue",lty=2)
legend("topright",col=c("red","blue"),lty=2,lwd=2,legend=c("Sem overdispersion","Com overdispersion"),bty="n")
 
#Grafico 3
hist((C.OD[which(x=='Nativo')]),freq=F,main="Nativo",xlab="Histograma de Contagens",ylab="Frequêcia da contagem")
lines(density(rpois(1000,out$summary[1,1]+out$summary[2,1]),adjust = 3),lwd=2,col="red",lty=2)
valores<-rep(out$summary[1,1]+out$summary[2,1],1000)+rnorm(1000,0,out$summary[3,1])
valores[valores<0]<-0
lines(density(rpois(1000,valores),adjust = 3),lwd=2,col="blue",lty=2)
legend("topright",col=c("red","blue"),lty=2,lwd=2,legend=c("Sem overdispersion","Com overdispersion"),bty="n")
 
 
#Exemplo 2
n.locais <- 60
x <- gl(n = 2, k = n.locais, labels = c("Nativo", "Arado"))
lambda.OD <- exp(0.60 +(0.35*(as.numeric(x)-1)) )
 
C.OD <- rpois(n = 2*n.locais, lambda = lambda.OD)
 
boxplot(C.OD ~ x, col = "grey", xlab = "Uso da terra",frame=F,
ylab = "Contagem de coelhos", las = 1, ylim = c(0, max(C.OD)))
 
bugs.data2 <- list(C.OD = C.OD, x = as.numeric(x)-1, n = length(x))
 
params <- c("alpha", "beta", "sigma")
 
out2 <- bugs(data=bugs.data2, inits=inits, parameters.to.save=params,model.file="overdispersion.txt", n.thin=nt,
            n.chains=nc,n.burnin=nb, n.iter=ni)
print(out2, dig = 3)