Poisson Point Process

A gente ja falou aqui sobre a distribuição dos pontos nos espaço, mas existe mais coisa a se pensar.

Eu comecei a ler o novo livro do Marc Kery, Applied hierarchical modeling in Ecology, eu ja gostava muito das coisas que ele escrevia desde os dois primeiros livros dele, Introduction to WinBUGS for Ecologists e Bayesian population analysis using WinBUGS, bem se o livro do Ben Bolker Ecological Models and Data in R abria os olhos para um monte de coisa, o Marc Kery detalhava bem modelos bayesianos desde de um começo simples, e me parece que com esse novo livro, vamos um pouquinho mais para a frente, ao mesmo estilo.

Bem uma coisa que comum no livro dele, que fala principalmente de modelos populacionais relacionados a abundância e ocupação, é o processo de pontos de poisson, ou Poisson Point Process. A gente ja falou da distribuição de poisson, que aliás é comum a muitas disciplinas, não so a ecologia. Mas vamos pensar um pouco mais sobre ela.

O Poisson point process (me perdoem usar o termo em inglês, mas as vezes é mais fácil ja ver o termo em inglês, tanto porque é mais fácil para procurar no google, achar na literatura, como muitas vezes eu traduzo as coisas erradas com meu inglês nórdico e depois vira uma confusão) tem a propriedade que cada ponto é estocasticamente independente de todos os outros pontos no processo. Apesar disso, ele é amplamente usado para modelar fenomemos representados por pontos, como fazemos direto ao usar o glm com distribuição de poisson, mas isso implica que ele nao descreve adequadamente fenomenos que tem interação suficientemente forte entre os pontos, o que é bem comum em biologia, e que explica porque a gente ve coisas como modelos quasi-poisson no glm e porque nos preocupamos com overdispersion.

O processo tem seu nome do matemático Siméon Denis Poisson e vem do fato que se uma coleção de pontos aleatórios no espaço formam um processo de poisson, então o número de pontos em uma região finita é uma variável aleatória com uma distribuição de poisson. O processo depende de um único parâmetro, que, dependendo do contexto, pode ser constante ou uma função localmente integrável (na maioria dos nossos casos a+bx ou uma função com várias médias, como numa anova). Esse parâmetro representa o valor esperado da densidade de pontos em um processo de Poisson, também conhecido como intensidade e normalmente representado pela letra grega lambda (\lambda).

A função massa de probabilidade é dada por:

p(n)=\frac{\lambda^n e^{-\lambda}}{n!}

Bem, a integral de p(n) da 1, como toda função massa de probabilidade (isso parece obvio hoje, mas eu demorei so me liguei disso depois de uma disciplina no coursera).

Para um \lambda grande, teremos o equivalente a uma distribuição normal.

Lembrando, a distribuição normal é:

p(x)=\frac{1}{\sqrt{2\sigma^2\pi}}e^{\frac{(x-\mu)^2}{2\sigma^2}}

Tudo isso basta dar ?rnorm e ?rpois, que você vai ver que ambas as distribuições são implementadas assim no R. Mas agora olha so que mágico.

Se você olhar uma parada chamada aproximação de stirling, que é uma formula fechada que mais ou menos da o resultado de um número fatorial para valores grandes, nos temos que

x! \approx \sqrt{2\pi x}e^{-x}x^x \ quando \ x \rightarrow \infty

Então, poisson é uma distribuição discreta e a normal continua, mas seja x=n=\lambda(1+\delta) onde \lambda \geq 10 (essa relação so vale para números grandes) e \delta \leq 1 (poisson é discreto, normal é contínuo). Dai substituindo temos…

p(x)=\frac{\lambda^{\lambda(1+\delta)} e^{-\lambda}}{\sqrt{2\pi \lambda(1+\delta)}e^{-\lambda(1+\delta)}(\lambda(1+\delta))^{\lambda(1+\delta)}}

Primeiro trazemos a potência do e do divisor para cima

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda(1+\delta)}(\lambda(1+\delta))^{\lambda(1+\delta) - \lambda(1+\delta)}}

Depois descemos a potência do lambda la de cima

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda(1+\delta)}(\lambda(1+\delta))}

simplificamos

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda}(\lambda(1+\delta)) (1+\delta)^{1/2}}

ai tiramos o delta de dentro da raiz

p(x)=\frac{ e^{\lambda \delta} \lambda(1+\delta)^{-1/2}}{\sqrt{2\pi \lambda}(1+\delta)}

subimos ele

p(x)=\frac{ e^{\lambda \delta} \lambda^{-1/2}}{\sqrt{2\pi \lambda}}

cortamos o um mais delta

p(x)=\frac{ \lambda e^{-(\lambda \delta)/2}}{\sqrt{2\pi \lambda}}

e arrumamos as potências

E agora tem uma parte que eu não entendo perfeitamente, mas lembre-se que a média e a variância da distribuição de poisson são iguais, então \lambda=\mu=\sigma^2.
Assim, a parte debaixo já está pronta, e a parte de cima da equação, já é praticamente a pdf da distribuição normal, só não sei exatamente o que fazer agora, mas como vemos a relação dela com a distribuição normal nesse diagrama do blog do John Cook, é algo por ae, se alguém quiser comentar o que fazer daqui para frente, eu ficaria feliz 🙂

distribution_chart

Tudo bem que eu não sou muito bom em matemática, mas se a gente fizer um gráfico das pdfs, assim:

1
2
3
plot(1:30,dpois(1:30,10),type="l",lty=3,lwd=3,col="red",frame=F,xlab="x",ylab="Probabilidade")
points(1:30,dnorm(1:30,10,sqrt(10)),type="l",lty=3,lwd=3,col="blue")
legend("topright",legend=c("Poisson","Normal"),lty=3,lwd=3,col=c("red","blue"),bty="n")

A gente vai ver o seguinte:

01

Bem é isso ai, sem script hoje, mas olhe o repositório recologia de qualquer forma, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail, que para esse post em especial, eu tenho certeza que devo ter falado algo errado!!! Então correções são bem vindas.

Referência:

John D. Cook – Diagram of distribution relationships
Marc Kéry & Andy Royle 2015 Applied hierarchical modeling in Ecology Elsevier 183pp

Método de Monte Carlo

De forma geral, o método de Monte Carlo é um método estatístico, experimental de calcular integrais usando posições aleatórias, chamadas amostras, cuja distribuição tem que ser escolhidas com muito cuidado. Aleatória, ou “Random” vem de uma palavra francesa antiga randon (to run around, andar em volta) e “sample”, amostras é derivada do latin exemplum, exemplo, um pouco de cultura inútil para alegrar o dia.

Eu as vezes estava dando uma olhada no livro Statistical Mechanics. Algorithms and Computation, mais para olhar o primeiro capitulo, que tem uma introdução muito legal sobre MCMC.

Primeiro ele fala do exemplo que já vimos a alguma tempo atrás aqui.
Mas ele explica isso como um jogo de crianças.

pebbles

Então o joguinho tem um circulo desenhado dentro de um quadrado, e os guris vão jogando pedras, e vão anotando quantas pedras caem dentro do circulo. Ai basicamente ele traz o seguinte algorítimo.

mcalg.

Podemos implementar ele no R, para testar.

1
2
3
4
5
6
7
8
9
10
11
direct_pi<- function(n=4000) {
    n_acertos<-0
    for(i in 1:n){
        x<-runif(n=1,min=-1,max=1)
        y<-runif(n=1,min=-1,max=1)
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}

E fazer alguns experimentos.

> direct_pi() [1] 3152

Se usarmos nossa função, temos o número de acertos, relativos a 4 mil tentativas, 4 mil jogadas de pedras no chão. E se repetirmos vezes o suficiente, que podemos fazer com a função replicate do R, que replica coisas.

> repeticoes<-replicate(100,direct_pi()) > mean(4*(repeticoes/4000)) [1] 3.14076 > pi pi>min(4*(repeticoes/4000)) [1] TRUE

Vemos que fazendo isso 100 vezes, a média é muito próxima do valor de pi real, também que pi está contido dentro do menor e maior valores encontrados.

Massa, mas até aqui estamos falando de Monte Carlo, agora entra o Markov Chain, que é a tal da cadeia. Aqui, ele propõe outro jogo, pensando que os adultos foram jogar, só que eles foram jogar num heliporto. A primeira grande diferença aqui é que agora não da para jogar uma pedra, e cobrir toda a área, porque agora o espaço é muito grande. Mas o que da para fazer é ir andando pela área. Então a gente joga uma pedra, vai onde ela caiu, e joga outra pedra, de olhos fechados e assim vai, devagarinho, cobrindo toda a área do heliporto.

mcmcheli

Agora, o acerto é estar dentro do circulo, então todo mundo vai com um saquinho cheio de pedras, e toda lugar onde joga sua pedra, deixa uma pedrinha onde ela caiu para contar depois.
Mas existe um problema, se estamos jogando a pedra, indo onde ela está e jogando denovo, pode acontecer de eventualmente, chegarmos na borda do heliporto e jogar a pedra pra fora. E agora, a gente vai buscar a pedra? Continua jogando ela la fora até voltar pro heliporto? Bom, uma solução, é contar denovo, deixar 2 pedrinhas no mesmo lugar para contar depois, e jogar denovo do mesmo lugar, a sua pedra de arremessar.

Nisso a gente vai ter o seguinte

mcmcheli2

Agora a gente tem uma pilhas de pedras nas beiradas, e só pedras espalhadas pelo meio. Meio louco, mas esse joguinho funciona (alias, isso é o metropolis-hasting).
E temos um algorítimo para isso também.

markov_pi

E podemos implementar ele no R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
markov_pi<-function(n=4000,passo=0.1){
    n_acertos<-0
    x<-1
    y<-1
    for(i in 1:n){
        delta_x<-runif(n=1,min=-passo,max=passo)
        delta_y<-runif(n=1,min=-passo,max=passo)
        if(abs(x+delta_x)<1 & abs(y+delta_y)<1){
            x<-x+delta_x
            y<-y+delta_y
        }
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}

Veja que agora a gente começa de um lugar (x=1,y=1), e vai andando no heliporto, andando na área, então, o antes, a gente jogava as pedrinhas, simplesmente sorteava números, agora a gente sorteia uma distância, para adicionar ao passo, que são os deltas.
Outra coisa, veja que se a gente olhar agora, antes de adicionar um acerto, se a pedra está dentro do heliporto, se estiver fora, a gente não soma os deltas, ou seja o local onde a pedra está continua igual, o que é adicionado aos acertos, ou não.

Podemos testar essa nova função, e vemos o seguinte.

> markov_pi() [1] 3397

Temos a contagem de passos.

> repeticoes<-replicate(100,markov_pi()) > mean(4*(repeticoes/4000)) [1] 3.11109

E o na média, se realizamos esse procedimento várias vezes, o valor é próximo de pi, sendo que o valor real está contido entre o maior e menor caso.

> pi pi>min(4*(repeticoes/4000)) [1] TRUE

Mas diferente do Monte Carlo, veja que aqui, nos não pegamos as amostras diretamente, como fizemos antes, aqui uma amostra depende da anterior, depende de onde você está no espaço amostral, e isso pode gerar uma certa dependência a curto prazo (que é que a gente tenta resolver quanto usa o thin no openbugs), mas ainda sim funciona.
Outra coisa é que veja que dependendo do tamanho do passo, a gente pode não explorar totalmente o espaço.

Se aumentarmos, ainda sim conseguimos andar bem em tudo

> repeticoes<-replicate(100,markov_pi(passo=0.5)) > mean(4*(repeticoes/4000)) [1] 3.13863

Mas se os passos forem muito pequenos, simplesmente podemos nunca chegar na borda, simplesmente demorar muito mais para percorrer todo o espaço.

> repeticoes<-replicate(100,markov_pi(passo=0.01)) > mean(4*(repeticoes/4000)) [1] 0.83486

Isso implica que para o MCMC, quando não retiramos nossas amostras diretamente da distribuição de interesse, temos que nos preocupar com esse ajuste do algorítimo, uma regra consensual, é que a rejeição, ou aceitação de novas amostras deve ficar em torno de 50%, o que é possível ser avaliado, é só contar o quanto estamos rejeitando amostras, para dar um parecer.

Bem é isso ai, so queria compartilhar essa explicação sobre MCMC, que achei muito legal, eu vi o cara do livro num curso do Coursera na verdade, aqui e depois eu consegui uma cópia do livro para dar uma olhada, o que é bem legal também. Alias o livro tem uma excelente página wiki, de acompanhamento, com os algorítimos em tudo que é linguagem, menos no R hehe, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Werner Krauth 2006 Statistical Mechanics – Algorithms and Computations

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#Monte Carlo
direct_pi<- function(n=4000) {
    n_acertos<-0
    for(i in 1:n){
        x<-runif(n=1,min=-1,max=1)
        y<-runif(n=1,min=-1,max=1)
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}
 
#Testes
direct_pi()
 
repeticoes<-replicate(100,direct_pi())
mean(4*(repeticoes/4000))
pi<max(4*(repeticoes/4000))
pi>min(4*(repeticoes/4000))
 
#MCMC
markov_pi<-function(n=4000,passo=0.1){
    n_acertos<-0
    x<-1
    y<-1
    for(i in 1:n){
        delta_x<-runif(n=1,min=-passo,max=passo)
        delta_y<-runif(n=1,min=-passo,max=passo)
        if(abs(x+delta_x)<1 & abs(y+delta_y)<1){
            x<-x+delta_x
            y<-y+delta_y
        }
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}
 
 
#Teste
markov_pi()
 
repeticoes<-replicate(100,markov_pi())
mean(4*(repeticoes/4000))
pi<max(4*(repeticoes/4000))
pi>min(4*(repeticoes/4000))
 
 
repeticoes<-replicate(100,markov_pi(passo=0.5))
mean(4*(repeticoes/4000))
 
repeticoes<-replicate(100,markov_pi(passo=0.01))
mean(4*(repeticoes/4000))
pi<max(4*(repeticoes/4000))
pi>min(4*(repeticoes/4000))

Data Augmentation para dados de Captura e Recaptura com efeito do tempo como uma variável aleatória

A algum tempo atrás, a gente falou de data augmentation aqui e depois usando data augmentation num modelo com efeito do tempo aqui.

Muito legal, mas particularmente no modelo com efeito do tempo, a gente cogitou (no livro do Marc Kery é cogitado na verdade) a possibilidade de tratar esse efeito como uma variável aleatória, então vamos tentar fazer isso e uma pequena modificação depois para ver o que acontece.

Bem os dados serão os mesmo que para o post de data augumentation com efeito do tempo. Mas vamos mudar o modelo. Na verdade vamos mudar um pouco a geração dos dados, vamos pegar o efeito de tempo de uma distribuição normal, para condizer com nossa analise.

Mas no mais é tudo igual, agora o modelo na linguagem Bugs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()

Basicamente ele continua o mesmo, só mudamos a parte do prior la em cima.

1
2
3
4
5
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo

Veja que agora temos ainda a geração dos efeitos do tempo, de cada tempo de amostragem, mas agora pegamos esses valores de uma distribuição normal, e não de uma distribuição uniforme como da outra vez.

Mas tentamos determinar agora, qual a média e o desvio padrão dessa distribuição.
Para a média partimos de uma distribuição normal com média zero e desvio 0.001, ou seja, um pico no zero, e para o desvio padrão saímos de uma distribuição uniforme de 0 a 5, lembre-se que o desvio padrão tem que ser um número maior que zero, e não pode ser negativo, senão o bugs vai reclamar.

Uma outra coisa da linguagem bugs é que a distribuição normal usa precisão, e não o desvio diretamente na função. Isso não tem nenhum impacto na estatística, tem haver com precisão numérica apenas eu acho. Até hoje não entendi esse conceito perfeitamente, mas isso não muda nada.

1
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)

No fim, tau vai ser 1 dividido pela variância, que é o desvio padrão ao quadrado, veja que o desvio só vai aumentar, e um sobre um valor entre zero e alguma coisa, vai ser cada vez mais um valor pequenininho, e bem é isso.

Fazemos umas cadeias maiores, adicionamos na função de geradores de valores iniciais nossos novos priors

1
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),sigma.tempo = runif(1,0, 5))

E é isso, agora é so rodar o modelo.

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 5000 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 13500 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.825 0.995 99.000 99.000 100.000 100.000 102.000 1.001 7100 p[1] 0.656 0.047 0.562 0.625 0.658 0.689 0.743 1.001 14000 p[2] 0.140 0.034 0.080 0.115 0.138 0.162 0.212 1.002 3000 p[3] 0.070 0.025 0.029 0.051 0.067 0.085 0.124 1.002 2400 p[4] 0.559 0.049 0.464 0.525 0.559 0.593 0.653 1.001 14000 p[5] 0.931 0.027 0.869 0.916 0.935 0.951 0.974 1.002 2000 med.tempo 0.470 0.306 -0.123 0.325 0.469 0.612 1.073 1.001 14000 sigma.tempo 0.579 0.381 0.236 0.359 0.471 0.663 1.608 1.001 14000 omega 0.402 0.031 0.341 0.381 0.402 0.423 0.464 1.001 5100 deviance 436.569 10.169 424.600 428.000 435.400 441.800 461.100 1.001 4200 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 = var(deviance)/2) pD = 51.7 and DIC = 488.3 DIC is an estimate of expected predictive error (lower deviance is better). >

E apesar da estimativa do tamanho populacional estar ok, as estimativas do tempo ficam bem ruins, confira os valores que geraram os dados com os que foram estimados.

> variacao_Tempo [1] 1.517343 -1.058690 -1.568697 1.357698 3.792000

Veja que a estimativa do desvio padrão ficou bem abaixo do real usado, que era 2. Ela ficou algo entre zero e 1.6, enquanto o desvio real era 2. A média ficou correta quanto ao usado para gerar os dados, englobando o zero, mas o desvio deixou a desejar.

Se a gente olhar a figura de das estimativas do tamanho da população.

01

Veja que a gente perdeu bastante em variabilidade, já que o tempo variou menos, porque ficou “preso” dentro do nosso desvio padrão subestimado.

Mas ainda existe uma possibilidade. Lembre-se que aqui temos poucas medidas de de tempo, apenas 5. Lembra de alguma distribuição que tentamos tratar do número de amostras para controlar o desvio da média? Se pensou na distribuição t, pensou certo.

Podemos modelar a variável aleatória que gera o tempo como uma distribuição t, isso vai adicionar um parâmetro a mais no modelo, mas no entanto, teremos como ter a possibilidade de ter os valores mais longe da média como valores mais comuns para o valor do tempo. Ja que se você ler o post sobre a distribuição t, vera que isso que o novo parâmetro no t faz, o grau de liberdade, quanto mais graus de liberdade, mais parecido com uma distribuição normal, mas quanto menos, deixando a distribuição mais “uniforme”, com uma maior chance de valor mais longe da média, e como temos apenas 5 valores, isso pode ser uma propriedade interessante. É meio maluco pensar que vamos estimar esse parâmetro, já que na distribuição t ele é simplesmente o valor de amostras, mas ok, podemos estimar um valor para ele, para tentar ajustar um modelo mais realista.

De todo aquele modelo, vamos mudar apenar a parte dos prior, mais especificamente onde está a variável aleatória que é o tempo.

1
2
3
4
5
6
7
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
        k.tempo ~ dunif(2, 30) # Hyperprior da variacao do tempo
 
        for (i in 1:tempo){
            p[i] ~ dt(med.tempo, tau.tempo, k.tempo)
        }

Veja que usamos o dt ao invés de dnorm dentro do loop, e temos um parâmetro chamado de k.tempo agora, que vamos tirar de uma distribuição uniforme de 2 a 30, quando ele for 30, temos uma distribuição normal, quanto menor, deixamos tudo mais uniforme. Começamos no 2 porque o valor de k tem que ser no mínimo 2 para a linguagem bugs, veja a definição aqui.
E é isso, agora é só colocar ele no inicializador.

1
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),sigma.tempo = runif(1,0, 5),k.tempo= runif(1,2, 30))

E pronto, podemos rodar. Não esquecendo que temos que colocar ele entre os parâmetros que temos que acompanhar.

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 5000 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 13500 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.866 1.010 99.000 99.000 100.000 100.000 102.000 1.001 11000 med.tempo 0.470 0.294 -0.100 0.322 0.469 0.613 1.065 1.001 14000 sigma.tempo 0.558 0.361 0.217 0.345 0.462 0.648 1.519 1.001 5900 k.tempo 16.565 7.884 2.979 9.950 16.660 23.360 29.310 1.001 5100 omega 0.402 0.031 0.343 0.380 0.402 0.423 0.465 1.001 14000 deviance 437.026 10.301 424.700 428.200 435.800 442.900 461.900 1.001 14000 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 = var(deviance)/2) pD = 53.1 and DIC = 490.1 DIC is an estimate of expected predictive error (lower deviance is better).

E olhando para a figura temos a mesma coisa

02

E no final das contas, o valor de k ficou entre 2 e 30, ou seja, não conseguimos estimar nada para ele, ja que de 2 a 30 era as possibilidades que abrimos. As estimativas de média para a distribuição de parâmetros continua ok, mas o desvio ainda está subestimado. Então apesar dessa ideia de usar a distribuição t seja legal, ela não contribuiu em nada praticamente, aumentamos um parâmetro e não melhoramos muito o modelo. Veja que o deviance continua a mesma coisa. Bom mas valeu a tentativa.

Eu imagino que dado o pequeno número de tempos, essa ideia de usar ele como uma variável aleatória não seja uma boa ideia, ja que estamos trocando 5 parâmetros por 3, ou 2 no primeiro exemplo. Além disso, precisamos lembrar que essa é uma afirmativa forte, pensar que que o tempo é uma variável aleatória, talvez num span pequeno de tempo, pode estar tudo ok, mas num span de tempo maior, dependendo da espécie que estamos trabalhando, podemos pensar o tanto que falamos em aquecimento global, em degradação do ambiente, em tantas mudanças que falar que nada acontece no tempo é algo bem “forte”. Mas apesar de tudo, talvez para mais medidas e um span pequeno de tempo, essa seja uma possibilidade bem legal.

Bem é isso ai, a gente tem que tentar fazer as coisas para aprender, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Kéry e Shaub 2011 – Bayesian population analysis using WinBUGS. A hierarchical perspective – Academic Press 554pp

John K. Kruschke 2014 Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan.2nd Edition. Academic Press / Elsevier.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
library(R2OpenBUGS)
 
##################################
## Gerando dados
##################################
set.seed(51)
n = 100
probabilidade_detecao = 0.25
tempo = 5
variacao_Tempo = rnorm(tempo,0,2)
 
yreal<- matrix(NA, nrow=n,ncol=tempo)
 
for (j in 1:tempo){
    p <- plogis(log(probabilidade_detecao / (1-probabilidade_detecao)) + variacao_Tempo[j])
    yreal[,j] <- rbinom(n = n, size = 1, prob = p)
}
 
detectado_pelo_menos_uma_vez <- apply(yreal, 1, max)
C <- sum(detectado_pelo_menos_uma_vez)
 
yobservado <- matrix(NA, nrow=n,ncol=tempo)
yobservado <- yreal[detectado_pelo_menos_uma_vez == 1,]
 
##################################
## Modelo com tempo
##################################
 
# Data Augumentation
nz <- 150
yaug <- rbind(yobservado, array(0, dim = c(nz, tempo)))
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()
 
# Juntando os dados
bugs.data <- list(yaug = yaug, amostras = nrow(yaug), tempo = ncol(yaug))
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),
                  sigma.tempo = runif(1,0, 5))
 
# Parametros para saida
params <- c("N", "p","med.tempo","sigma.tempo","omega")
 
# Configurações do MCMC
ni <- 5000
nt <- 2
nb <- 500
nc <- 3
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura1
jpeg("01.jpg")
hist(out_tempo$sims.list$N,breaks=seq(90,110,by=1), col = "gray", main = "", xlab = "Tamanho Populacional", las = 1,xlim=c(95,105))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
dev.off()
#####################################
## Efeito aleatorio com distribuiçao t
#####################################
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
 
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
        k.tempo ~ dunif(2, 30) # Hyperprior da variacao do tempo
 
        for (i in 1:tempo){
            p[i] ~ dt(med.tempo, tau.tempo, k.tempo)
        }
 
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),
                         sigma.tempo = runif(1,0, 5),k.tempo= runif(1,2, 30))
inits()
# Parametros para saida
params <- c("N","med.tempo","sigma.tempo","k.tempo","omega")
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura2
jpeg("02.jpg")
hist(out_tempo$sims.list$N,breaks=seq(90,110,by=1), col = "gray", main = "", xlab = "Tamanho Populacional", las = 1,xlim=c(95,105))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
dev.off()

Data Augmentation para dados de Captura e Recaptura com efeito do tempo

Ano passado, a gente olhou aqui a ideia de data augumentation. Que é muito interessante. Ai eu queria fazer posts sobre todos esses modelos, que são muito legais, mas fiquei enrolando, mas agora sai mais um e sabe-se la quando o resto.

Mas vamos la, esse modelo assume que a probabilidade de detecção p varia por ocasião de coleta, se estivermos coletando passarinhos, talvez a chuva faça com que eles não saiam para voar e cair na rede, ou para pequenos mamíferos, as condições do clima estejam ruins, o dia para ficar na cama abaixo das cobertas, ou ainda diferentes armadilhas ou mecanismos de detecção foram usados.

Então vamos gerar dados, pensando em passarinhos caindo numa rede, que colocamos no campo, e conseguimos ir pegando os passarinhos e anilhando eles, conforme caem na rede, então identificamos eles perfeitamente.

Então imagine:

1
2
3
4
n = 100
probabilidade_detecao = 0.25
tempo = 5
variacao_Tempo = runif(tempo,-2,2)

Cem passarinhos (n), coletados em 5 campanha diferentes (tempo), mas além desses parâmetros, sabemos que… , bem nos não sabemos, mas a mãe natureza estabelece uma probabilidade desse passarinho cair na rede (probabilidade_detecao), mas essa probabilidade não é constante ao longo das campanhas, as vezes ela é mais um pouquinho, quando o clima ta bom, sem vento, e as vezes menos um pouquinho (esses valores estão no vetor variacao_Tempo).

Ok, agora vamos fazer a matriz de dados que a mãe natureza vê, ela vê todos dos 100 passarinhos que la residem, e ela observa as 5 campanhas também, ou seja uma matriz de n pelo tempo.

1
yreal<- matrix(NA, nrow=n,ncol=tempo)

Para simular os dados fazemos o seguinte

1
2
3
4
5
6
7
#Num loop para todas as campanha
for (j in 1:tempo){
    #Qual a probabilidade de detecção nessa campanha j?
    p <- plogis(log(probabilidade_detecao / (1-probabilidade_detecao)) + variacao_Tempo[j])
    #Ai simulamos todos os 100 individuos de uma vez
    yreal[,j] <- rbinom(n = n, size = 1, prob = p)
}

E pronto, a simulação é isso, gerar números a partir de uma distribuição binomial, na verdade.
Só que quem a gente nunca viu em nenhuma campanha, quem nunca foi anilhado não existe pra gente, só foi anilhado quem foi pego pelo menos uma vez

1
2
detectado_pelo_menos_uma_vez <- apply(yreal, 1, max)
C <- sum(detectado_pelo_menos_uma_vez)

E assim esses são os cara que observamos

1
2
yobservado <- matrix(NA, nrow=n,ncol=tempo)
yobservado <- yreal[detectado_pelo_menos_uma_vez == 1,]

Agora basta a gente estabelecer nosso modelo para a linguagem bugs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dunif(0, 1)
        }
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        p.medio <- mean(p[])
    }

E a principio o likelihood pode parecer idêntico ao do post anterior, mas veja que a diferença está em usar um vetor para as probabilidade de detecção da campanha

1
p.eff[i,j] <- z[i] * p[j]

Veja que no prior, inicializamos um vetor agora, e não um único valor, para contabilizar essa variabilidade entre campanha

1
2
3
for (i in 1:tempo){
     p[i] ~ dunif(0, 1)
}

Ai agora é só estimar os parâmetros, temos que colocar os dados no formado do openbugs que estamos usando aqui, lembrando das outras opções aqui, estabelecer valores iniciais para as cadeia do MCMC, seu parâmetros como tamanho, valores iniciais a serem queimados, número de cadeia e rodar.

Aqui o resultado fica o seguinte:

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.604 4.646 92.000 96.000 99.000 102.000 110.000 1.001 6000 p[1] 0.542 0.055 0.436 0.505 0.542 0.579 0.649 1.001 6000 p[2] 0.108 0.031 0.055 0.086 0.105 0.127 0.179 1.001 6000 p[3] 0.138 0.034 0.078 0.113 0.135 0.160 0.211 1.001 6000 p[4] 0.672 0.056 0.560 0.635 0.673 0.711 0.776 1.001 6000 p[5] 0.099 0.030 0.049 0.078 0.096 0.117 0.164 1.001 6000 p.medio 0.312 0.022 0.268 0.297 0.312 0.327 0.354 1.001 6000 omega 0.418 0.037 0.349 0.392 0.417 0.443 0.495 1.002 2700 deviance 470.030 20.741 434.300 455.500 468.300 483.100 514.302 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 = var(deviance)/2) pD = 215.1 and DIC = 685.2 DIC is an estimate of expected predictive error (lower deviance is better).

Bem existem mais ou menos entre 92 a 110 indivíduos por ali, o intervalo de confiança, o que cobre a verdade (já que sabemos os valores usados pela mãe natureza foi 100), e temos um valor de detecção por campanha e duas quantidades derivadas, o p médio e o omega, que falamos depois.

Podemos tentar olhar esse resultado numa figura das nossas estimativas.

01

Beleza, mas será que valeu a pena o esforço, o que aconteceria se a gente ajusta-se o modelo considerando um único valor de detecção. Bem, podemos ajustar e ver o que acontece, por curiosidade.

E a saída é o seguinte

> print(out, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 112.268 7.386 100.000 107.000 112.000 117.000 129.000 1.002 1700 p 0.274 0.026 0.224 0.256 0.274 0.292 0.326 1.002 1700 omega 0.471 0.044 0.390 0.440 0.469 0.499 0.561 1.002 1400 deviance 657.828 23.208 616.100 641.300 656.900 672.425 707.800 1.002 1700 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 = var(deviance)/2) pD = 269.1 and DIC = 926.9 DIC is an estimate of expected predictive error (lower deviance is better).

Aqui o intervalo de confiança ainda pegou o valor real, que era de 100, o tamanho populacional teve uma estimativa de 100 a 129 indivíduos, o que está ok, com um p de 274, a probabilidade média de captura. Agora o que isso quer dizer?
Vamos olhar os valores da mãe natureza.

Veja que os valores de detecção foram

> plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo) [1] 0.50132242 0.09139408 0.12948952 0.69882099 0.09683113

Da média original, na primeira campanha foi bem mais, na segunda e terceira bem menos, na quarta a mais, e na quinta a menos. se fizermos a média disso da

> mean(plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo)) [1] 0.3035716

O valor que aquele p realmente queria chegar era 0.304, mas ele ficou em 0.274, isso quer dizer, que na segunda, terceira e quinta campanha, a gente deveria ter pensado em achar menos passarinhos, mas o modelo fez simulação com valores altos, sempre achando mais, por isso a gente super-estima muito, veja que comparando as duas distribuições, nesse exemplo o modelo sem levar o tempo foi la para frente.

figura2

Veja que isso é a quantidade derivada la do modelo com o tempo, que ficou próximo desse valor real.

p.medio 0.312 0.022 0.268 0.297 0.312 0.327 0.354 1.001 6000

Um pouquinho maior, mas próximo.

Veja que isso vai depender claro do número de campanhas, pensando que a não ser que haja uma tendência nesses valores de detecção, por exemplo, a não ser que eles sempre melhorem, ou piorem, se eles simplesmente forem um processo estocástico com média zero, com muitas campanhas, os modelo sem o tempo conseguiria um valor bom de estimativa populacional, mas basta pensar um pouco para saber que isso não vai acontecer, porque?

Porque campanhas de coletas são caras, dificilmente o dinheiro vai dar para 30 ou 90 campanhas, e com poucas campanhas, a chance de você considerar um valor fixo de detecção e dar errado é maior, chance no sentido de que o modelo sempre vai ser pior, então talvez começar com um modelo assim seja mais seguro, e é fácil você simplificar se for possível.

Veja o deviance do modelo sem o tempo é muito pior:

deviance 657.828 23.208 616.100 641.300 656.900 672.425 707.800 1.002 1700

Compare como o deviance do modelo com o tempo é muito melhor

deviance 470.030 20.741 434.300 455.500 468.300 483.100 514.302 1.001 6000

Veja o quanto, quando comparamos as estimativas, o desvio padrão das estimativas para o modelo com tempo é menor.

figura3

A importante contribuição da incerteza sobre o tamanho da população vem da incerteza sobre a probabilidade p. Só que a gente paga em parâmetros a mais que precisamos para estimar p. Mas ignorar esses parâmetros levam a erros, que podem ser muito ruins. Talvez se a variação fosse pequena, esse modelo nem vale-se a pena, mas na dúvida, compensa começar do modelo mais complexo, e ir simplificando, vendo se é possível simplificar, afinal, se todas as probabilidades de detecção fossem iguais, daria para ver nos parâmetros estimados e não precisaríamos de um vetor de probabilidade de detecção, mas pero que si, pero que no, compensa tentar todas as possibilidades.

Nesse modelo, nos assumimos que os efeitos do tempo são fixos, mas poderíamos ainda ajustar um modelo estocástico também, pensando que ao invés de valores fixos de p, eles viriam de uma distribuição como a normal ou uniforme, que é o caso dos dados que geramos aqui.

Bem é isso ai, o script vai estar la no repositório recologia, falta uns dois modelos a mais sobre esse tema no livro muito legal do Kery, mas um dia chegamos la, isso tem no capítulo 6. Se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e até mais.

Referência:
Kéry e Shaub 2011 – Bayesian population analysis using WinBUGS. A hierarchical perspective – Academic Press 554pp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
#
library(R2OpenBUGS)
 
##################################
## Gerando dados
##################################
set.seed(51)
n = 100
probabilidade_detecao = 0.25
tempo = 5
variacao_Tempo = runif(tempo,-2,2)
 
yreal<- matrix(NA, nrow=n,ncol=tempo)
 
for (j in 1:tempo){
    p <- plogis(log(probabilidade_detecao / (1-probabilidade_detecao)) + variacao_Tempo[j])
    yreal[,j] <- rbinom(n = n, size = 1, prob = p)
}
 
detectado_pelo_menos_uma_vez <- apply(yreal, 1, max)
C <- sum(detectado_pelo_menos_uma_vez)
 
yobservado <- matrix(NA, nrow=n,ncol=tempo)
yobservado <- yreal[detectado_pelo_menos_uma_vez == 1,]
 
##################################
## Modelo com tempo
##################################
 
# Data Augumentation
nz <- 150
yaug <- rbind(yobservado, array(0, dim = c(nz, tempo)))
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dunif(0, 1)
        }
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        p.medio <- mean(p[])
    }
    ",fill = TRUE)
sink()
 
# Juntando os dados
bugs.data <- list(yaug = yaug, amostras = nrow(yaug), tempo = ncol(yaug))
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1))
 
# Parametros para saida
params <- c("N", "p","p.medio", "omega")
 
# Configurações do MCMC
ni <- 2500
nt <- 2
nb <- 500
nc <- 3
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura1
hist(out_tempo$sims.list$N,nclass=30, col = "gray", main = "", xlab = "Tamanho Populacional", las = 1, xlim = c(80, 120))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
 
##################################
## Modelo sem variação de tempo
##################################
 
# Modelo na linguagem bugs
sink("modelo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        p ~ dunif(0, 1)
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
    }
    ",fill = TRUE)
sink()
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(1, 0, 1))
 
#
params <- c("N", "p", "omega")
 
#Rodando o modelo
out <- bugs(bugs.data, inits, params, "modelo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out, dig = 3)
print(out_tempo, dig = 3)
 
#Figura 2
hist(out_tempo$sims.list$N, nclass = 40, col = "darkgray", main = "", xlab = "Tamanho Populacional", las = 1, xlim = c(70, 150))
hist(out$sims.list$N, nclass = 40, col = "lightgray",add=T)
abline(v = mean(out_tempo$sims.list$N), col = "blue", lwd = 3)
abline(v = mean(out$sims.list$N), col = "red", lwd = 3)
legend("topright",c("Modelo com tempo","Modelo sem tempo"),lty=1,lwd=3,col=c("blue","red"),bty="n")
 
#
plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo)
mean(plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo))
 
#Figura 3
par(mfrow=c(2,1))
hist(out_tempo$sims.list$N, nclass = 40, col = "darkgray", main = "Modelo com variação no tempo",
     xlab = "Tamanho Populacional", las = 1, xlim = c(70, 150),ylim=c(0,600))
text(130,400,paste("Desvio padrão =",round(sd(out_tempo$sims.list$N),3)))
hist(out$sims.list$N, nclass = 40, col = "darkgray", main = "Modelo sem variação no tempo",
     xlab = "Tamanho Populacional", las = 1, xlim = c(70, 150),ylim=c(0,600))
text(130,400,paste("Desvio padrão =",round(sd(out$sims.list$N),3)))

Data Augmentation para dados de Captura e Recaptura

Acho que um tipo de dados muito comum em ecologia são os de Captura e Recaptura. Desde aqueles feitos com armadilhas, onde elas são vistoriadas periodicamente num grid, até as armadilhas de câmeras, que tiram fotos por sensor de movimento, redes capturando passarinhos, peixes, borboletas, talvez mais algumas formas de conseguir dados assim.

A ideia é que capturando indivíduos de alguma espécie, e recapturando depois (basta dar um jeito de identificar o um indivíduo) a gente pode calcular uma serie parâmetros populacionais que temos muito interesse, taxas de crescimento, sobrevivência, migração e assim vai.

Mas esses modelo sempre foram meio misticos pra mim, um monte de formula fechada que vem em um programa, a gente usa sem saber o que é, e usa da forma fechadinha. Mas se a gente conseguir ver esse tipo de modelo como um glm, as portas se abrem para um monte de possibilidades. Nesse post aqui a gente viu como estimar a colonização e extinção, usando o pacote unmarked e um modelo Bugs equivalente.

No modelo final, a gente vê que coisas como colonização, extinção são na verdade medidas derivadas, e que a resposta propriamente dita é encontrar ou não um individuo em uma amostra, observar ele ou não, a captura. Sendo que indivíduos podem ser recapturados. Resumindo, tentamos fazer um modelo assim.

logit(p_{i,j})=\alpha_i+\beta_j+\delta \cdot F_{i,j}+ \gamma \cdot X_{i,j}

E o \alpha_i \sim Normal(\mu_{\alpha},\sigma^2_{\alpha})

Vamos entender isso tudo, primeiro aqui, a gente ta lidando com dois índices, o i e o j, aqui o i diz respeito aos indivíduos e j ao evento de captura. Por exemplo pegando um exemplo em que a gente arma uma rede para pegar passarinhos, o i é referente a identidade do passarinho, o número da anilha dele, e o j é o dia que a rede for armada, no final isso resulta numa matriz de dados onde os indivíduos são linhas e os dias de rede armada as colunas, e dentro dessa matriz a gente marca as captura com um e zero para quando não houve captura, sendo que uma linha pode ter mais de uma captura, a gente pode recapturar um mesmo passarinho.

Agora seguindo esse exemplo, indivíduos, passarinhos podem diferir em capturabilidade, a gente pode tratar isso, tratar o fato de alguns passarinhos serem mais bobinhos e outros não assumindo que eles vem de uma distribuição normal para essa variabilidade, teremos então alguns bobinhos, alguns espertalhões mas a maioria vai ter uma determinada capturabilidade, isso é o \alpha_i da formula, veja que ele usa o índice i de individuo, e dentro desse \alpha podemos modelar outras coisas dos indivíduos, sei la, sexo por exemplo, se machos são mais curiosos, derrepente sejam mais capturados.

Agora a chance de uma captura não depende só de um individuo, por exemplo um dia ensolarado pode ter mais passarinho voando que um dia de chuva, ou dia depois da chuva, ou seja, o dia da amostra pode influenciar na capturabilidade. Isso ai a gente vei colocar dentro do \beta_j, lembrando que j é o índice dos eventos de captura. Então a gente pode modelar as coisas desse tipo aqui dentro. Lembrando que tanto no \alpha quanto o \beta, a gente ta colocando uma letrinha, mas cabe uma formulinha aqui dentro, com intercepto e inclinação de por exemplo, dados que você coleta dos passarinhos, tipo a morfometria, se você acha que ela tem a ver com a chance do passarinho ser capturado, podemos tentar observar isso.

Agora as outras duas coisas, o \delta \cdot F_{i,j} tem haver com o famoso ditado, gato escaldado tem medo de água, ou o contrario. Apos primeira comida, a chance de captura pode ser alterada, tanto para menos quanto pra mais, por exemplo, se ser capturado é muito estressante, o individuo pode começar a ficar mais esperto e difícil de capturar, ou se a gente ta lidando com uma gaiola com alguma isca, é fácil pensar que passa pela cabeça do individuo a frase “boca livre”, e ele pode começar a ser capturado com mais facilidade por procurar uma boca livre, ou sei la, veja que esse parâmetro depende tanto do individuo i como da captura j, já que esse efeito a partir da segunda captura, se alguém é capturado apenas uma vez, não tem como apresentar esse efeito, por isso leva os dois índices.

O último parâmetro que tem ali é o \gamma \cdot X_{i,j}, que ai são as questões ambientais, como é o ambiente, ou a disponibilidade de algo especial. Isso a gente enquadra ai nesse \gamma, veja que ele também depende de i e j.

Certo então com o modelo ali em cima a gente conseguiria estimar uma probabilidade de captura, mas qual o interesse nisso?
A gente se interessa em quantidades mais palpáveis, o tamanho da população, sobrevivência, imigração, mas essas quantidade pode ser derivada dessa chance de captura, sendo que com o modelo acima a gente tem a chance de considerar praticamente tudo que é do nosso interesse.

Agora vamos simular uns dados aqui para ver melhor como isso acontece. Vamos supor uma população de 100 indivíduos, ai vamos supor também uma chance de captura de 50% e que coletamos por três vezes, em três eventos de coleta.

Beleza, se tudo que colocamos no modelo la em cima for constante, teríamos o seguinte.

> N = 100 > p = 0.5 > T = 3 > Real<-array(NA, dim = c(N, T)) > for (j in 1:T){ + Real[,j] <- rbinom(n = N, size = 1, prob = p) + } > Real [,1] [,2] [,3] [1,] 0 1 1 [2,] 1 0 1 [3,] 0 1 0 [4,] 0 0 1 . . . [95,] 0 1 1 [96,] 1 0 1 [97,] 0 0 0 [98,] 1 0 1 [99,] 1 1 1 [100,] 0 1 1

Certo, colocamos nossos parâmetros e meramente sorteamos de uma distribuição binomial a com a chance de captura que definimos se capturamos ou não. Temos uma matriz de 100 linhas, que são 100 indivíduos e 3 colunas, que são os três eventos de coleta, agora da uma olhada o que aconteceu na linha 97, a gente nunca capturou esse individuo nos três eventos.

Se a gente somar as linhas dessa matriz, a gente vai ver que

> detectado <- apply(Real, 1, max) > sum(detectado) [1] 86

vimos somente 86 indivíduos dos 100 que existe, existem 14 sortudos que nunca caíram na rede, simples assim.
Agora a pergunta é a seguinte, existem 100 indivíduos, mas a gente não ve uma matriz como a acima, esse valor de 100 esta oculto pela natureza, aqui a gente sabe porque inventamos ele, mas com dados reais a gente não sabe.
Os dados reais que a gente ve vão ser assim:

> dados<-array(NA, dim = c(N, T)) > dados <- Real[detectado == 1,] > dados [,1] [,2] [,3] [1,] 0 1 1 [2,] 1 0 1 [3,] 0 1 0 . . . [83,] 1 0 1 [84,] 1 0 1 [85,] 1 1 1 [86,] 0 1 1

So vão até a linha 86, porque as 14 linhas que so tem zero a gente não viu, não ve e nunca verá. Nossas observações são sempre imperfeitas. Agora como faz para chegar aquele valor de 100 indivíduos, se a gente so tem nossas observações acima?

Existiam outros estimadores, presentas em programas como o Capture e o Mark, eles estimavam tudo e depois faziam uma conta para estimar o tamanho da população, aquele valor de 100. Agora quando a gente ta usando estatística bayesiana, com MCMC e tal, existe um problema que a gente vai de iteração a iteração, e o vetor de N pode mudar de iteração a iteração, ou seja não cabe nos dados que a gente manda para o BUGS. Mas como problemas existem para serem solucionados, dois caras chamados Tanner e Wong inventaram uma técnica chamada de data augmentation, que o Royle e companhia trouxeram pros modelos de captura e recaptura usando MCMC.

Basicamente duas coisas são feitas, adicionamos um numero arbitrário de zeros a matriz de dados analisamos uma versão reparametrizada dos dados originais. Essencialmente convertendo um modelo de população fechada para um modelo de ocupação.

Então a gente está adicionando um grande numero de observações em potencial que não foram observadas. Os dados aumentados terão um número de linhas M e o mesmo numero de colunas T das outras observações, se o tamanho da população deveria ser o número de linhas, que deveria ser 100, e só temos 86 linhas, ao aumentarmos os dados teremos uma quantidade linhas agora muito muito maior que o tamanho da população, que é o que nos interessa. Nesses dados a gente ajusta um tipo de modelo de zero inflado se o tamanho da população fosse conhecido. Para isso a gente adiciona um indicador binario que indica se uma linha da matriz de dados aumentada representa um individuo real ou um que não existe na pratica. Esse indicador é dado por um prior da distribuição de Bernoulli, ou Binomial, e esse parâmetro, que podemos chamar de \omega (omega) é estimável dos dados. Uma probabilidade de inclusão. Com essa probabilidade de inclusão a gente tem o tamanho da população saindo dos dados.

Então criamos nosso modelo na linguagem bugs

# Modelo na linguagem Bugs
sink("modelo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        p ~ dunif(0, 1)
        # Likelihood
        for (i in 1:M){
            z[i] ~ dbern(omega)
            # Indicador de inclusão
            for (j in 1:T){
                dadosaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p # So pode ser detectado se z=1
                } #j
            } #i
        # Quantidades derivadas
        N <- sum(z[])
        }
    ",fill = TRUE)
sink()

Vamos observar alguma coisas, primeiro que temos dois priores

omega ~ dunif(0, 1)
p ~ dunif(0, 1)

Que são as chances daquela linha ser um indivíduo e a chance do indivíduo ter sido capturado.

Veja então que no modelo de verossimilhança, a gente tem uma chance do individuo ser existir

for (i in 1:M){
    z[i] ~ dbern(omega)

E a chance dele ter sido capturado num dado evento pra esse individuo

for (j in 1:T){
    dadosaug[i,j] ~ dbern(p.eff[i,j])
    p.eff[i,j] <- z[i] * p # So pode ser detectado se z=1

Mas veza que ele é multiplicado pelo valor de z, que vai ser zero ou um, com uma chance omega, que é a chance dele ser um individuo real.

N <- sum(z[])

Ou seja a soma do vetor de z vai ser exatamente a estimativa do tamanho populacional.

Agora é só juntar os dados, gerar valores iniciais, configurar as cadeias e rodar o esse modelo. E temos o seguinte resultado.

> print(out, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 95.451 4.114 89.000 92.000 95.000 98.000 105.000 1.001 2700 p 0.542 0.037 0.466 0.517 0.542 0.567 0.613 1.001 4100 omega 0.405 0.036 0.337 0.380 0.404 0.429 0.478 1.001 4000 deviance 395.385 18.993 363.200 383.300 393.500 407.200 436.600 1.001 2800 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 = var(deviance)/2) pD = 180.3 and DIC = 575.7 DIC is an estimate of expected predictive error (lower deviance is better).

Sempre temos que fazer aquelas validações, que tivemos uma convergência dos parâmetros, mas podemos ver que funciona bem.

Primeiro vamos olhar o N que é o tamanho populacional

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 95.451 4.114 89.000 92.000 95.000 98.000 105.000 1.001 2700

Veja que temos uma média de 95 para o parâmetro, que é perto de 100, o valor real do tamanho populacional, com um desvio padrão de 4.
Olhando o intervalo de confiança vemos que ele vai de 89 a 105, o que inclui o valor real de 100, ou seja, deu certo, esse negocio funciona mesmo para estimar o tamanho populacional.

E estimamos bem também a chance de captura, que na realidade era de 0.5 (50%), que esta dentro do intervalo de confiança.

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff p 0.542 0.037 0.466 0.517 0.542 0.567 0.613 1.001 4100

Então vamos tentar fazer um gráfico para ver melhor essa saída.

01

Certo, nossa estimativa cobre o valor real, dentro do intervalo de confiança de 95%, não sei nem se posso falar intervalo de confiança, mas é o 95% ali das saídas, as vezes esses termos tem pequenos pressupostos que a gente não sabe, até escrever errado.
Mas agora podemos ver como o \omega (omega) faz, vamos olhar o valor dele.

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff omega 0.405 0.036 0.337 0.380 0.404 0.429 0.478 1.001 4000

Média de 0.405 com desvio de 0.036, esse 40.5% é o número de linhas da matriz de dados aumentada que mandamos que são provavelmente realmente indivíduos. Como a gente tinha 86 observações, e aumentamos a matriz em 150 linhas, temos um total de 236 linhas. E dessas 40.5% consideramos como indivíduos.

> 0.405*(86+150) [1] 95.58

As outras não devem ser observações, apenas linhas. Mas o aumento deu a possibilidade de considerar esses individuos que não observamos, mas deviam estar la.
Mas uma pergunta razoavel é o quanto devemos aumentar a matriz. Porque a gente aumentou 150 aqui e deu certo, mas será que qualquer número funciona?
Bem sabemos que o tamanho da população vai ser igual ao número das observações, ou maior. Se aumentamos pouco a matriz, por exemplo, somente em 5 linhas, podemos ver o que acontece.

> print(out5, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 90.004 1.093 87.000 89.000 90.000 91.000 91.000 1.001 6000 p 0.574 0.031 0.512 0.553 0.574 0.595 0.636 1.001 5300 omega 0.978 0.019 0.930 0.970 0.984 0.992 0.999 1.002 6000 deviance 369.337 5.821 354.500 365.200 369.700 373.700 376.900 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 = var(deviance)/2) pD = 16.9 and DIC = 386.3 DIC is an estimate of expected predictive error (lower deviance is better).

Veja que o valor de omega ficou espremido, grudado no 100%, olhando o grafico isso fica mais facil de visualizar.

02

Veja que a distribuição é bem truncada, como se tivesse uma barreira ali segurando a distribuição de ir mais para direita, ou seja a gente precisava de mais “espaço”. Então se a gente ver uma distribuição truncada assim, sabemos que temos que aumentar o tamanho da matriz mais um pouco. Com 150 deu legal, podemos olhar denovo isso.

> print(out150, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 95.551 4.222 89.000 93.000 95.000 98.000 105.000 1.001 6000 p 0.542 0.038 0.467 0.517 0.543 0.568 0.615 1.002 1700 omega 0.406 0.037 0.339 0.381 0.405 0.430 0.481 1.001 5800 deviance 395.845 19.407 363.300 383.300 393.500 407.400 438.710 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 = var(deviance)/2) pD = 188.3 and DIC = 584.2 DIC is an estimate of expected predictive error (lower deviance is better).

03

Aqui a distribuição ficou legal, veja que o omega não ficou espremido, truncado nem para a esquerda nem para a direita. Então está tudo ok. Mas e se a gente aumentar demais, tem problema? Será que o precavido aqui vai se ferrar?

> print(out1500, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 95.501 4.175 89.000 92.000 95.000 98.000 105.000 1.001 6000 p 0.542 0.037 0.468 0.517 0.542 0.568 0.614 1.001 2800 omega 0.061 0.006 0.049 0.056 0.061 0.065 0.074 1.001 3400 deviance 395.610 19.188 363.300 383.300 393.600 407.100 437.600 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 = var(deviance)/2) pD = 184.1 and DIC = 579.7 DIC is an estimate of expected predictive error (lower deviance is better).

04

E a resposta é não, não tem problema nenhum aumentar muito. O que acontece é que o omega vai ficar baixinho, mas não ficando truncado ta tudo certo. Veja que aumentando em 1500 linhas a matriz, temos exatamente a mesma estimativa de tamanho populacional, unica coisa é que aqui tivemos uma matriz gigante para processar, e aqui o 6.1% de 1586 linhas que é o total de indivíduos.

Eu achei essa solução muito mirabolante. Primeiro porque a gente cai em um glm, o que torna um modelo de captura e recaptura algo bem próximo dos glm que comumente olhamos e essa sacada de aumentar a matriz é muito foda. Melhor sobrar que faltar. Agora aqui conseguimos ter uma noção de que data augumentation funciona, mas em um modelo simples.

No inicio do post colocamos o seguinte modelo.
logit(p_{i,j})=\alpha_i+\beta_j+\delta \cdot F_{i,j}+ \gamma \cdot X_{i,j}
E o \alpha_i \sim Normal(\mu_{\alpha},\sigma^2_{\alpha})

Mas o que ajustamos foi algo como

logit(p_{i,j})=\alpha

Somente um parâmetro, so lembrar do modelo, que so tem uma linha aproximando nossas observações.

dadosaug[i,j] ~ dbern(p.eff[i,j])

Agora temos que começar a expandir esse parametro para o modelo completo. Sendo que não somos obrigados a tentar ajustar tudo, podemos usar somente aquilo que nos interessa, ou que temos dados coletados.

Bem é isso ai, hora de ir dormir. O script vai estar la no repositório recologia, se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, ou mande um e-mail e boa noite.

Referencia:
Kéry e Shaub 2011 – Bayesian population analysis using WinBUGS. A hierarchical perspective – Academic Press 554pp

#Gerando dados
N = 100
p = 0.5
T = 3
Real<-array(NA, dim = c(N, T))
for (j in 1:T){
    Real[,j] <- rbinom(n = N, size = 1, prob = p)
    }
#A visão da natureza
Real
 
detectado <- apply(Real, 1, max)
sum(detectado)
 
dados<-array(NA, dim = c(N, T))
dados <- Real[detectado == 1,]
#A nossa visão
dados
 
# Aumentando os dados com 150 individuos em potencial
nz <- 150
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
 
# Modelo na linguagem Bugs
sink("modelo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        p ~ dunif(0, 1)
        # Likelihood
        for (i in 1:M){
            z[i] ~ dbern(omega)
            # Indicador de inclusão
            for (j in 1:T){
                dadosaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p # So pode ser detectado se z=1
                } #j
            } #i
        # Quantidades derivadas
        N <- sum(z[])
        }
    ",fill = TRUE)
sink()
 
# Juntando os dados para mandar para o bugs
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
# Função geradora de valores iniciais
inits <- function() list(z = rep(1, nrow(dadosaug)), p = runif(1, 0, 1))
# Parametros a serem monitorados
params <- c("N", "p", "omega")
 
# Configurando o MCMC
ni <- 2500
nt <- 2
nb <- 500
nc <- 3
# Chamando o openbugs do R
library(R2OpenBUGS)
out <- bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
 
 
# Distribuições posteriores
print(out, dig = 3)
 
str(out)
out$summary
 
#Figura 1
hist(out$sims.list$N,nclass=50,col="gray",main="",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))
 
 
nz <- 5
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
 
out5 <- bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
print(out5, dig = 3)
 
#Figura 2
hist(out5$sims.list$N,nclass=50,col="gray",main="Dados aumentados em 5",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out5$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))
 
 
nz <- 150
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
 
 
out150 <-bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
print(out150, dig = 3)
 
#Figura 3
hist(out150$sims.list$N,nclass=50,col="gray",main="Dados aumentados em 150",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out150$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))
 
 
nz <- 1500
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
 
out1500 <- bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
print(out1500, dig = 3)
 
#Figura 4
hist(out1500$sims.list$N,nclass=50,col="gray",main="Dados aumentados em 1500",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out1500$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))