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

Leave a Reply

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