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()

Leave a Reply

Your email address will not be published. Required fields are marked *