Regressão segmentada ou piecewise.

Depois de um tempo sumido, espero voltar finalmente a postar aqui. Estou devagar para estudar ultimamente, é de dezembro meu último post aqui, mas agora que acabou minha faculdade, vou voltar a estudar novas coisas e postar regularmente, existe muito coisa que pensei em postar e parei.

Mas para retomar as atividades, vamos falar sobre a piecewise regression, ou regressão segmentada. A gente ta acostumado a modelos de diferenças de médias, tipo anova, e modelos de regressão linear, ou a combinação desses dois. Eles resolvem bem quase todos os casos, mas as vezes a gente já sabe um pouco mais sobre o sistema que estamos estudando, e podemos fazer melhor, utilizando esse conhecimento.

Por exemplo, sabemos que plantas muitas vezes são limitadas, no seu crescimento, pelo mínimo de alguns nutrientes, a lei de Liebig. Podemos expandir esse conceito e pensar assim, enquanto não há um mínimo de nutrientes no solo, o crescimento é estagnado, mas a partir de algum ponto começamos a ter um incremento no crescimento de uma espécie dado um acréscimo nos nutrientes.

Podemos simular alguns dados para essa situação, então imaginando que temos uma medida dos nutrientes num vaso e do tamanho da planta após um determinado tempo fixo, podemos visualizar esses dados da seguinte forma.

Veja então que no eixo x temos a concentração de nutrientes, e no eixo y temos o tamanho das plantulas, então enquanto não temos uma certa concentração parece que todo mundo tem o mesmo tamanho, mas a partir de uma concentração, quanto maior a concentração, maior ficam as plantulas. Se a gente pensar, isso pode ocorrer em um bocado de situações, o que interessa é que, a partir de um ponto, as coisas mudam, algo começa a ocorrer, ou para, aqui temos um ponto crítico, algo entre o 4 e o 6 em que abaixo dele, a concentração de nutrientes é irrelevante.

Podemos representar isso matematicamente como uma equação do tipo piecewise, nesse caso, o que usamos para simular os dados foi:

  f(x) = \left\{          \begin{array}{ll}              8 & \quad x\ \textless \  5 \\              1+2x & \quad x \geq 5          \end{array}      \right.

Sendo x o nosso preditor, concentração de nutriente e f(x) o tamanho, para a parte determinística do modelo.
Esses são os valores que usamos para gerar essas dados de exemplo, que no caso é partido no 5, antes do 5 temos uma reta f(x)=8+0x(eu omiti o 0x, porque tudo que é multiplicado por zero da zero, mas temos uma reta, só que sem inclinação, ou melhor, com inclinação zero) e depois do 5 temos uma reta f(x)=1+2x.

Nesse caso, o partição, o valor 5 tem um significa biológico importante, ele é o limite mínimo de nutrientes, para começarmos a observar alguma coisa, antes dele nada acontece, mas como podemos recuperar esse valor dos dados? A primeira opção é usar um tipo de força bruta, testar o modelo acima.

Para aquela equação podemos ajustar ela usando uma variável tipo dummy usando a função lm, ou usar a função nls e ajustar diretamente a equação piecewise usando um ifelse, aqui vamos usar diretamente a equação com nls que é mais direto a interpretação dos parâmetros(uma reta é a+bx, então vamos ter dois a, o a1 e a2 e dois b, o b1 e b2 para as duas retas), ajustamos duas retas e um valor de quebra, no caso ajustamos o modelo para muitas valores de quebra e vemos qual da o melhor ajuste usando logLikelihood

Então criamos um vetor de possíveis valores de quebra.

1
2
quebras <- seq(1,9,0.5)
minimo_quadrado <- vector(length = length(quebras))

E testamos todos eles

1
2
3
4
for(i in 1:length(quebras)){
 piecewise <- nls(tamanho ~ ifelse(nutriente < quebras[i],a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1))
 minimo_quadrado[i] <- logLik(piecewise)
}

E podemos avaliar graficamente o resultado.

Certo, como esperávamos, o valor da simulação é recuperado, veja que o pico de ajusta está bem no 5.

Podemos então ajustar o modelo usando 5 como ponto de quebra e ver se recuperamos os coeficientes corretamente

> modelo<-nls(tamanho ~ ifelse(nutriente < 5,a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1)) > summary(modelo) Formula: tamanho ~ ifelse(nutriente < 5, a1 + b1 * nutriente, a2 + b2 * nutriente) Parameters: Estimate Std. Error t value Pr(>|t|) a1 7.7307 0.3410 22.668 <2e-16 *** a2 3.2449 1.1431 2.839 0.0063 ** b1 0.2193 0.1249 1.756 0.0846 . b2 1.7430 0.1369 12.734 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9221 on 56 degrees of freedom Number of iterations to convergence: 1 Achieved convergence tolerance: 1.57e-08

E temos as inclinações, zero para a primeira equação esta no intervalo do b1, que é o nome que demos para o parâmetro, e assim como o 2 está no intervalo do b2, os interceptos também são recuperados. Vamos ver ele graficamente, e além disso, vamos comparar com um modelo linear.

Veja que ele fica bem mais certinho, mas o modelo linear também detectaria uma tendência aqui, porém não nos forneceria essa informação, de a partir de onde, de qual concentração é relevante, seriamos enganados pensando que a concentração de nutrientes é relevante em qualquer intensidade quando ela so importa em uma parte do preditor, isso mesmo avaliando os resíduos dos modelos.

Temos um ajuste menor (lembre-se que estamos falando de números negativos, então -78 é maior que -111), mas não temos uma tendência clara nos resíduos para esse caso, rejeitar o modelo linear, somente observando os resíduos, então aqui interessa muito o que pensamos a priori da construção do modelo.

Como podemos ter casos onde a busca precisa ser mais precisa para o ponto de quebra, temos funções prontas que vão cuidar dessas atividades, por exemplo o pacote segmented, que ajusta esse tipo de modelo diretamente, fazendo um update do modelo linear.

Basicamente temos o seguinte resultado:

> modelo_pieciwise<- segmented(modelo_linear, seg.Z = ~nutriente, psi=1) > summary(modelo_pieciwise) ***Regression Model with Segmented Relationship(s)*** Call: segmented.lm(obj = modelo_linear, seg.Z = ~nutriente, psi = 1) Estimated Break-Point(s): Est. St.Err 3.799 0.304 Meaningful coefficients of the linear terms: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.8733 0.4533 17.369 <2e-16 *** nutriente 0.1274 0.1967 0.648 0.52 U1.nutriente 1.9158 0.2222 8.621 NA --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.045 on 56 degrees of freedom Multiple R-Squared: 0.9576, Adjusted R-squared: 0.9554 Convergence attained in 7 iterations with relative change 3.485019e-16

Veja que temos o Estimated Break-Point, que da 3.799, um pouco inferior ao que realmente usamos, e temos apenas um intercepto, mas a primeira inclinação fica em zero, e a segunda em 0.12+1.91, diferença de inclinação para depois do ponto de quebra, que é 2, então ele estima bem os parâmetros, mas é um modelo um pouco diferente do que usamos (só um intercepto), eu não li muito sobre esse pacote, para saber ajustar exatamente como fizemos ali em cima, mas é um pacote a ser explorado, para quem quer utilizar esse tipo de modelo e o sistema de plot dele é bem simples.

Bem é isso ai, depois de um tempo enroscado, de volta a atividade, o script vai estar lá no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Michael J. Crawley 2012 The R Book, 2nd ed. Wiley 1076pp

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
set.seed(31415)
 
nutriente <- runif(60,0,10)
 
fn_bs <- function(x){
    saida<-vector(length=length(x))
    for(i in 1:length(x)){       
        if(x[i]<5){
            saida[i] <- 8+0*x[i]
        }else{
            saida[i] <- 1+2*x[i]
        }
    }
    saida <- saida+rnorm(length(x),0,1)
    return(saida)
}
 
tamanho <- fn_bs(nutriente)
 
plot(tamanho~nutriente,xlab="Concentração de nutrientes",ylab="Tamanho da plantula (cm)",frame=F,ylim=c(0,30),xlim=c(0,12),pch=19)
 
#####################################
## Força bruta
#####################################
quebras <- seq(1,9,0.5)
minimo_quadrado <- vector(length = length(quebras))
 
for(i in 1:length(quebras)){
 piecewise <- nls(tamanho ~ ifelse(nutriente < quebras[i],a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1))
 minimo_quadrado[i] <- logLik(piecewise)
}
 
minimo_quadrado <- unlist(minimo_quadrado)
 
plot(quebras,minimo_quadrado,type="b",pch=19,xlab="Valor de quebra",ylab="LogLikelihood",frame=F,xlim=c(0,10))
 
modelo<-nls(tamanho ~ ifelse(nutriente < 5,a1+b1*nutriente,a2+b2*nutriente),start = c(a1=1,a2=1,b1=1,b2=1))
summary(modelo)
 
modelo_linear<- lm(tamanho~nutriente)
summary(modelo_linear)
 
plot(tamanho~nutriente,xlab="Concentração de nutrientes",ylab="Tamanho da plantula (cm)",frame=F,ylim=c(0,30),xlim=c(0,12),pch=19)
 
abline(modelo_linear,lwd=3,lty=3,col="blue")
 
curve(coef(modelo)[1]+coef(modelo)[3]*x,0,5,add=T,lwd=3,lty=3,col="red")
curve(coef(modelo)[2]+coef(modelo)[4]*x,5,12,add=T,lwd=3,lty=3,col="red")
abline(v=5,lwd=1,lty=2,col="red")
legend("topleft",lwd=3,lty=3,col=c("blue","red"),legend=c("Modelo Linear","Modelo piecewise"),bty="n")
 
par(mfrow=c(2,1))
plot(resid(modelo_linear),main=paste("Loglikelihood do modelo liner =",round(logLik(modelo_linear),3)),frame=F,xlab="Amostras",ylab="Resíduo",ylim=c(-5,5))
abline(h=0,lty=2,lwd=2)
plot(resid(modelo),main=paste("Loglikelihood do modelo piecewise =",round(logLik(modelo),3)),frame=F,xlab="Amostras",ylab="Resíduo",ylim=c(-5,5))
abline(h=0,lty=2,lwd=2)
 
#####################################
## Pacote Segmented
#####################################
##install.packages("segmented")
library(segmented)
 
modelo_pieciwise<- segmented(modelo_linear, seg.Z = ~nutriente, psi=1)
summary(modelo_pieciwise)
 
plot(tamanho~nutriente,xlab="Concentração de nutrientes",ylab="Tamanho da plantula (cm)",frame=F,ylim=c(0,30),xlim=c(0,12),pch=19)
plot(modelo_pieciwise,add=T,lwd=2,lty=2,col="red")

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