Duas espécies em metapopulações de Levins

Ja vimos metapopulações para uma espécie aqui, aqui e aqui, mas será que poderíamos adicionar mais uma espécie ao modelo de Levins?

Primeiro, vamos lembrar como é a metapopulação de levins.

\frac{dO}{dt}=cO(1-O)-eO

onde:

  • O – Ocupação
  • c – taxa de colonização
  • e – taxa de extinção

Então basicamente temos dois termos,

cN(1-N)

que faz a população crescer a uma taxa c, mas quanto maior a ocupação, mais propagulos são emitidos, mas veja que quanto maior a ocupação, menor o crescimento, que é a multiplicação por (1-N), conforme o N cresce, o valor tende a zero, porque cada vez teremos 1 menos alguma coisa, que da menos que um, estaremos multiplicando por uma fração menor.

eN

Cada mancha ocupada tem uma chance de extinção e, logo quanto mais manchas ocupadas, mais manchas veremos sendo extintas.

Agora vamos colocar outra espécie, ela vai ocupar as manchas vagas pela primeira espécie, ou seja, vagou, ela tem chance de entrar ocupar, mas a se a espécie um ocupar novamente, ela perde a competição.

Como são duas espécies, precisamos agora acompanhar a mudança na ocupação de cada uma, então precisamos de duas equações

E espécie 1 continua como a equação anterior
\frac{dO_1}{dt}=c_1O_1(1-O_1)-eO_1
mas a segunda fica
\frac{dO_2}{dt}=c_2O_2(1-O_1-O_2)-c_1*O_1*O_2-eO_2

Ou seja, a metapopulação 2 expande sua ocupação, diminuindo a expansão de acordo com o o quanto ja está ocupado pela ocupação da espécie 1 e dela mesmo, e tiramos as manchas que ela perde na competição, além do que ela perde com a extinção.

Agora, quem lembra do pacote desolve, vamos usar ele para ver como isso fica ao longo do tempo, então temos que transformar numa função para usar com a função ode.

1
2
3
4
5
6
7
8
9
10
11
metpop2sp <- function(t, y, p) {
    {
        o1 <-y[1] 
        o2 <-y[2]
    }
    with(as.list(p), {
        do1.dt <- c1*o1*(1-o1)-e*o1
        do2.dt <- c2*o2*(1-o1-o2)-c1*o1*o2-e*o2
        return(list(c(do1.dt, do2.dt)))
    })
}

Agora temos que determinar por quanto tempo queremos observar a ocupação e os paramtros definir os parametros de ocupação e extinção

1
2
3
4
5
6
tempo <- seq(0, 200, by =5)
c1 = 0.4
c2 = 0.5 
o1 = 0.1
o2 = 0.4 
e = 0.25

E vemos o seguinte

01

A espécie um rapidamente domina o ambiente, retirando a espécie 2 da jogada.

Mesmo quando mudamos os parametros um pouco, deixando ambas com a mesma taxa de colonização, maior a taxa de extinção, lembrando que para uma espécie, quando a extinção ficava maior que a colonização, ela tendia a extinção, sendo a diferença o quão rápido isso ia acontecer.

1
2
3
4
5
c1 = 0.1 
c2 = 0.1
o1 = 0.05 
o2 = 0.05 
e = 0.05

02

Agora se mudarmos o jogo para a espécie 1 extinguir, a espécie 2 consegue dominar, salvo que ela também não tenha a colonização menor que a extinção

1
2
3
4
5
c1 = 0.05
c2 = 0.07 
o1 = 0.05 
o2 = 0.05
e = 0.06

03

E existem também os casos onde podemos ter uma coexistencia

1
2
3
4
5
c1 = 0.1
c2 = 0.5
o1 = 0.05 
o2 = 0.05
e = 0.05

04

Agora, as possibilidades de combinações de parâmetros são grandes, tornando mais complexo definir quando uma ou outra espécie ganha, ou quando temos a cooexistencia, mas podemos ter um feeling da coisa já. além disso, nossas equações garantem a vitoria da espécie 1, mas isso não precisa ser assim.

Bem é isso ai, lembrando um pouco de equações diferencias, que não falamos a um bom tempo aqui e vendo a parte essa expansão para a metapopulação, 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:
A. A. de Oliveira e P. I. Prado – Coexistência em Metapopulações – Roteiro no EcoVirtual
Robert May & Angela McLean 2007 – Theoretical Ecology: Principles and Applications. Oxford University 272 pp

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
##
##install.packages("deSolve")
library(deSolve)
rm(list=ls())
 
##Função
metpop2sp <- function(t, y, p) {
    {
        o1 <-y[1] 
        o2 <-y[2]
    }
    with(as.list(p), {
        do1.dt <- c1*o1*(1-o1)-e*o1
        do2.dt <- c2*o2*(1-o1-o2)-c1*o1*o2-e*o2
        return(list(c(do1.dt, do2.dt)))
    })
}
 
##Tempo
tempo <- seq(0, 200, by =5)
 
##Caso 1
c1 = 0.4
c2 = 0.5 
o1 = 0.1
o2 = 0.4 
e = 0.25
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
 
##jpeg("01.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 2
c1 = 0.1 
c2 = 0.1
o1 = 0.05 
o2 = 0.05 
e = 0.05
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("02.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 3
c1 = 0.05
c2 = 0.07 
o1 = 0.05 
o2 = 0.05
e = 0.06
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("03.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 4
c1 = 0.1
c2 = 0.5
o1 = 0.05 
o2 = 0.05
e = 0.05
 
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("04.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()

Estimando colonização e extinção de uma espécie usando OpenBugs e Unmarked

Bem, nesse post aqui, a muito tempo atrás, vimos como estimar parâmetros do modelo de metapopulações. Mais especificamente dentro do nosso interesse, colonização e extinção.
Para relembrar um pouco sobre metapopulações, podemos dar uma olhada nesse post aqui, para ter uma ideia do modelo.

A ideia aqui é so começar a olhar, como estimar, os mesmos parâmetros, como no unmarked, para colonização e extinção. O mesmo que faz a função colext() desse pacote.Mas levando em conta também a detecção.

Pra falar a verdade, vamos primeiro olhar aqui como fica o modelo no na linguagem bugs, e depois em outros posts olhamos melhor pedacinho por pedacinho, mas vamos fazer ele e deixar funcionando para começar a realizar alguns testes.

Bem primeiro vamos olhar como são os dados, então temos a ocupação ao longo do tempo. Vamos simular alguns dados e obtemos o seguinte.

01

Veja que a ocupação que observamos é menor que a ocupação real do ambiente, a realizada realmente. Isso porque a detecção não é 100%, ou seja, se pensarmos em aranhas em arbustos, não é porque eu não vi a aranha que é 100% garantido que o arbusto esta vazio, e pior que isso, pode ser um processo estocástico, quero dizer, não é garantido que sempre eu “não consigo ver igual” a aranha, numa coleta num dia de vento pode ser pior que uma coleta num dia sem vento, no dia de vento a aranha se escondeu e ai eu posso achar que existem menos aranhas no ambiente que realmente tem. E assim vai.
Então ali a linha azul é como se fosse a ocupação (pode chamar de prevalência também) das aranhas se eu visse todas. Mas temos a linha vermelha pontilhada que é a chance de eu ver essas aranhas, ou seja quanto menor ela esta, menos aranha eu vejo, e no final o que eu realmente vejo no campo é a linha preta. Veja bem que em dados reais, tudo que a gente ve é a linha preta, as outras você nunca vê, a gente vê aqui porque simulou os dados, mas elas são ocultas, e é isso que tentamos estimar. Todas essas probabilidades.

E por curiosidade como são esses dados?

> data$y[,,1] [,1] [,2] [,3] [1,] 0 1 0 [2,] 1 0 0 [3,] 0 0 0 [4,] 0 0 0 . . . [27,] 0 1 1 [28,] 0 1 0 [29,] 0 1 0 [30,] 0 0 0

É tipo uma matriz de 3 dimensões, onde as linhas são os meus arbustos, eu tenho 30 nesse caso, eu ia três vezes ao campo todo ano, olha se minha arainha estava la ou não, e por vários anos, anos ficam como a terceira dimensão dessa matriz, ou array seria um termo melhor, da para organizar de outras formas os dados, mas assim fica mais simples de trabalhar com os índices, que ficam [amostras, repetições , anos]

Certo então vamos definir nosso modelo na linguagem bugs.

model {
 
        # Priors
        psi1 ~ dunif(0, 1)
        for (k in 1:(nyear-1)){
            phi[k] ~ dunif(0, 1)
            gamma[k] ~ dunif(0, 1)
            p[k] ~ dunif(0, 1)
        }
        p[nyear] ~ dunif(0, 1)
 
        # Modelo ecologico, definindo estado real
        for (i in 1:nsite){
            z[i,1] ~ dbern(psi1)
            for (k in 2:nyear){
                muZ[i,k]<- z[i,k-1]*phi[k-1] + (1-z[i,k-1])*gamma[k-1]
                z[i,k] ~ dbern(muZ[i,k])
            } #k
        } #i
 
        # Modelo referente as observações
        for (i in 1:nsite){
            for (j in 1:nrep){
                for (k in 1:nyear){
                    muy[i,j,k] <- z[i,k]*p[k]
                    y[i,j,k] ~ dbern(muy[i,j,k])
                } #k
            } #j
        } #i
 
# Parametros derivados
        psi[1] <- psi1
        n.occ[1]<-sum(z[1:nsite,1])
        for (k in 2:nyear){
            psi[k] <- psi[k-1]*phi[k-1] + (1-psi[k-1])*gamma[k-1]
            n.occ[k] <- sum(z[1:nsite,k])
            growthr[k-1] <- psi[k]/psi[k-1]
            turnover[k-1] <- (1 - psi[k-1]) * gamma[k-1]/psi[k]
        }
    }

E a parte mais importante é aqui, onde realmente as coisas estão acontecendo

 
        # Modelo ecologico, definindo estado real
        for (i in 1:nsite){
            z[i,1] ~ dbern(psi1)
            for (k in 2:nyear){
                muZ[i,k]<- z[i,k-1]*phi[k-1] + (1-z[i,k-1])*gamma[k-1]
                z[i,k] ~ dbern(muZ[i,k])
            } #k
        } #i

Veja que aqui a gente ta trabalho com o locais e as repetições dentro de anos.
Usamos distribuição binomial, ou seja nada mais é que um tipo de regressão logístico.
E se a gente lembrar do modelo do Levins:

{dN\over{dt}}={c \dot N \dot (1-N)-e\dot N}

Temos alguns comentários de muito tempo atras aqui.

Veja que isso é o que estamos usando nessa linha aqui:

muZ[i,k]<- z[i,k-1]*phi[k-1] + (1-z[i,k-1])*gamma[k-1]

So trocar o gamma por c, o z por N e o phi por e.

Mas essa quantidade Z pode estar subestimada, pode estar contata a menos, e isso a gente corrige aqui.

        # Modelo referente as observações
        for (i in 1:nsite){
            for (j in 1:nrep){
                for (k in 1:nyear){
                    muy[i,j,k] <- z[i,k]*p[k]
                    y[i,j,k] ~ dbern(muy[i,j,k])
                } #k
            } #j
        } #i

Juntando todos os dados no Bugs, temos o resultado:

> print(out1, dig = 2) Inference for Bugs model at "Dynocc.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 4 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff psi[1] 0.89 0.08 0.69 0.84 0.90 0.95 0.99 1.00 6000 psi[2] 0.83 0.07 0.68 0.78 0.83 0.87 0.94 1.00 2600 ... psi[9] 0.58 0.12 0.35 0.50 0.58 0.67 0.83 1.00 720 psi[10] 0.55 0.11 0.33 0.47 0.54 0.62 0.78 1.00 2400 phi[1] 0.87 0.07 0.71 0.82 0.88 0.93 0.99 1.00 6000 ... phi[9] 0.68 0.15 0.37 0.57 0.68 0.78 0.96 1.00 1200 gamma[1] 0.45 0.27 0.02 0.22 0.43 0.66 0.95 1.00 730 ... gamma[9] 0.35 0.19 0.03 0.21 0.34 0.47 0.73 1.00 1900 p[1] 0.30 0.05 0.20 0.26 0.29 0.33 0.41 1.00 6000 ... p[9] 0.40 0.09 0.24 0.34 0.40 0.47 0.60 1.00 6000 p[10] 0.44 0.09 0.27 0.37 0.43 0.50 0.62 1.00 5700 n.occ[1] 27.45 1.92 23.00 26.00 28.00 29.00 30.00 1.00 4200 ... n.occ[10] 16.58 2.71 13.00 14.00 16.00 18.00 23.00 1.00 1200 growthr[1] 0.94 0.11 0.75 0.87 0.93 0.99 1.17 1.00 2800 ... growthr[9] 0.97 0.29 0.53 0.77 0.93 1.13 1.67 1.00 680 turnover[1] 0.06 0.07 0.00 0.01 0.04 0.09 0.25 1.00 2200 ... turnover[9] 0.28 0.16 0.02 0.15 0.27 0.39 0.61 1.00 1100 deviance 771.25 20.97 732.50 756.60 770.60 785.02 813.20 1.00 1100 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 = 219.5 and DIC = 990.7 DIC is an estimate of expected predictive error (lower deviance is better).

Cada parâmetro por ano, lembrando que quando vemos por exemplo 9 parâmetros gamma, é porque ele é usado na transição de um tempo para o próximo, por isso temos 10 anos de dados simulados, mas 9 parâmetros estimados, não temos como estimar o gamma para a primeira observação.

Podemos testar a ocupação estimada com as dos dados originais e ver se obtivemos algum sucesso.

02

Bem o vermelhinho é o que estimamos, o pretinho são os dados que tinhamos, e o azul são os valores que fomos tirar la da simulação, o real, utilizado na simulação, com exceção do primeiro ano, o intervalo de confiança cobriu os valores reais.

É, parece que esse tipo de modelagem funciona, mas é preciso uma quantidade grande de dados, e aqui melhorar os dados é determinado por mais replicas, mais locais e mais tempo. Somente um desses três não melhora muito a parametrização.
Mas agora a gente ja sabe mais ou menos o que a função colext esta tentando fazer. Isso ja ajuda a pensar um pouco.

Agora é continuar estudando. Esse script, esta presente aqui, junto com outros varios exemplos legais. Existem também esse site referente ao pacote unmarked e um google group de discussão aqui. Esses são os pontos primeiros para buscar informação e exemplos para entender melhor tudo isso.

#install.packages("unmarked")
library(unmarked)
 
data.fn<-function(R=30,J=3,K=10,psi1=0.4,range.p=c(0.2, 0.4),range.phi=c(0.6,0.8),range.gamma=c(0, 0.1)) {
 
    # Iniciando os vetores
    site <- 1:R # Locais
    year <- 1:K # Anos
    psi <- rep(NA, K) # Probabilidade de ocupação
    muZ <- z <- array(dim = c(R, K)) # Ocorrencia esperada e realizada
    y <- array(NA, dim = c(R, J, K)) # Historico de deteção
 
    # Determinando a ocupação inicial e seus parametros
    psi[1] <- psi1 # Probabilidade de ocupação inicial
    p <- runif(n = K, min = range.p[1], max = range.p[2])
    phi <- runif(n = K-1, min = range.phi[1], max = range.phi[2])
    gamma <- runif(n = K-1, min = range.gamma[1], max = range.gamma[2])
 
    # Estados latentes de ocupação
    # Primeiro ano
    z[,1] <- rbinom(R, 1, psi[1])
    # Estado de ocupação inicial
    # Proximos anos
    for(i in 1:R){
        # Loop sobre locais
        for(k in 2:K){
            # Loop sobre anos
            muZ[k] <- z[i, k-1]*phi[k-1] + (1-z[i, k-1])*gamma[k-1] # Prob for occ.
            z[i,k] <- rbinom(1, 1, muZ[k])
            }
        }
    # Grafico da ocupação realizada
    plot(year,apply(z,2,mean),type="l",xlab="Ano",ylab="Ocupação",col="blue",xlim=c(0,K+1),
         ylim=c(0,1),lwd=2,lty=1,frame.plot=FALSE,las = 1)
    lines(year,p,type="l",col="red",lwd=2,lty=2)
    # Gerando dados de detecção/não detecção
    for(i in 1:R){
        for(k in 1:K){
            prob <- z[i,k] * p[k]
            for(j in 1:J){
                y[i,j,k] <- rbinom(1, 1, prob)
                }
            }
        }
    # Computando ocupação anual
    for (k in 2:K){
        psi[k] <- psi[k-1]*phi[k-1] + (1-psi[k-1])*gamma[k-1]
        }
    # Plotando ocupação aparente
    psi.app <- apply(apply(y,c(1,3), max), 2, mean)
    lines(year, psi.app, type= "l", col = "black", lwd = 2)
    legend("topright",lty=c(1,2,1),lwd=2,col=c("blue","red","black"),
           legend = c("Ocupação Verdadeira","Detecçao"," Ocupação Observada"))
 
 
# Retornando dados
    return(list(R=R,J=J,K=K,psi=psi,psi.app=psi.app,z=z,phi=phi,gamma=gamma,p=p,y=y))
    }
 
set.seed(123)
jpeg("01.jpg")
data<-data.fn(R=30,J=3,K=10,psi1=0.6,range.p=c(0.1,0.9),range.phi=c(0.7,0.9),range.gamma=c(0.1,0.5))
str(data)
dev.off()
 
yy <- matrix(data$y, data$R, data$J* data$K)
str(yy)
head(yy)
 
year <- matrix(c('01','02','03','04','05','06','07','08','09','10'),
nrow(yy), data$K, byrow=TRUE)
head(year)
 
simUMF <- unmarkedMultFrame(y = yy,yearlySiteCovs = list(year = year),numPrimary=data$K)
summary(simUMF)
 
fm1 <- colext(psiformula = ~1,      #Ocupação no primeiro ano
              gammaformula = ~ 1,   #Colonização
              epsilonformula = ~ 1, #Extinção
              pformula = ~ 1,       #Detecção
              data = simUMF)
 
summary(fm1)
system.time(summary(fm2 <- colext(psiformula=~1, gammaformula=~year-1,epsilonformula=~1, pformula =~1, data = simUMF,
                                  control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm3 <- colext(psiformula=~1, gammaformula=~1,epsilonformula=~year-1, pformula =~1, data = simUMF,
                                  control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm4 <- colext(psiformula=~1, gammaformula=~1,epsilonformula=~1, pformula =~ year-1, data = simUMF,
                                  control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm5 <- colext(psiformula=~1, gammaformula=~year-1,epsilonformula=~year-1, pformula =~1,
                                  data = simUMF,control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm6 <- colext(psiformula=~1, gammaformula=~year-1,epsilonformula=~1, pformula =~year-1,
                                  data = simUMF,control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm7 <- colext(psiformula=~1, gammaformula=~1,epsilonformula=~year-1, pformula =~year-1,
                                  data = simUMF,control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
system.time(summary(fm8 <- colext(psiformula=~1, gammaformula=~year-1,epsilonformula=~year-1, pformula =~year-1,
                                  data = simUMF,control=list(trace=TRUE, REPORT=1, maxit = 500), se = TRUE)))
 
 
modelos <- fitList(
'psi(.)gam(.)eps(.)p(.)' = fm1,
'psi(.)gam(Y)eps(.)p(.)' = fm2,
'psi(.)gam(.)eps(Y)p(.)' = fm3,
'psi(.)gam(.)eps(.)p(Y)' = fm4,
'psi(.)gam(Y)eps(Y)p(.)' = fm5,
'psi(.)gam(Y)eps(.)p(Y)' = fm6,
'psi(.)gam(.)eps(Y)p(Y)' = fm7,
'psi(.)gam(Y)eps(Y)p(Y)' = fm8)
modSel(modelos)
 
# Bundle data
bugs.data <- list(y = data$y, nsite = dim(data$y)[1], nrep = dim(data$y)[2],nyear = dim(data$y)[3])
 
data$y[,,1]
 
# Specify model in BUGS language
sink("Dynocc.txt")
cat("
    model {
 
        # Priors
        psi1 ~ dunif(0, 1)
        for (k in 1:(nyear-1)){
            phi[k] ~ dunif(0, 1)
            gamma[k] ~ dunif(0, 1)
            p[k] ~ dunif(0, 1)
        }
        p[nyear] ~ dunif(0, 1)
 
        # Modelo ecologico, definindo estado real
        for (i in 1:nsite){
            z[i,1] ~ dbern(psi1)
            for (k in 2:nyear){
                muZ[i,k]<- z[i,k-1]*phi[k-1] + (1-z[i,k-1])*gamma[k-1]
                z[i,k] ~ dbern(muZ[i,k])
            } #k
        } #i
 
        # Modelo referente as observações
        for (i in 1:nsite){
            for (j in 1:nrep){
                for (k in 1:nyear){
                    muy[i,j,k] <- z[i,k]*p[k]
                    y[i,j,k] ~ dbern(muy[i,j,k])
                } #k
            } #j
        } #i
 
# Parametros derivados
        psi[1] <- psi1
        n.occ[1]<-sum(z[1:nsite,1])
        for (k in 2:nyear){
            psi[k] <- psi[k-1]*phi[k-1] + (1-psi[k-1])*gamma[k-1]
            n.occ[k] <- sum(z[1:nsite,k])
            growthr[k-1] <- psi[k]/psi[k-1]
            turnover[k-1] <- (1 - psi[k-1]) * gamma[k-1]/psi[k]
        }
    }
    ",fill = TRUE)
sink()
 
# Valores iniciais da cadeia
zst <- apply(data$y, c(1, 3), max) # Observações ocorridas para iniciar Z
inits <- function(){ list(z = zst)}
# Parametros que queremos salvar
params <- c("psi", "phi", "gamma", "p", "n.occ", "growthr", "turnover")
# Definições do MCMC
ni <- 2500
nt <- 4
nb <- 500
nc <- 3
# Chamando Openbus (Pode demorar um pouco)
library(R2OpenBUGS)
out1 <- bugs(bugs.data,inits,params,"Dynocc.txt",n.chains=nc,n.thin=nt,n.iter=ni,n.burnin=nb)
 
# Sumario dos resultados
print(out1, dig = 2)
 
str(data)
jpeg("02.jpg")
plot(1:data$K,data$psi,type="l",xlab="Ano",ylab="Probabilidade de ocupação",col="blue",xlim=c(0,data$K+1),ylim=c(0,1),
     lwd=2,lty=1, frame.plot=FALSE,las = 1)
lines(1:data$K, data$psi.app, type = "l", col = "black", lwd = 2)
points(1:data$K, out1$mean$psi, type = "l", col = "red", lwd = 2)
segments(1:data$K, out1$summary[1:data$K,3], 1:data$K,out1$summary[1:data$K,7], col = "red", lwd= 1)
dev.off()

Metapopulações de Levins e o equilíbrio estável

Aproveitando que falamos que na análise bayesiana que quase tudo converge para um ponto, vamos olhar mais uma vez a equação de Levins e a dinâmica de metapopulações.

Desta vez façamos o seguinte, vamos começar a população com diferentes níveis de ocupação e ver o que acontece. Vamos usar 20%, 40%, 60% e 80% de ocupação, representando cada dinâmica com uma cor.

grafico

Bem, parece que independente de onde a população esta começando, ela vai para o mesma taxa de ocupação, ela caminha para uma estabilidade, como se tivesse sendo atraída para la.

Alias o ponto de equilíbrio é um tópico interessante para conservação e uso sustentável de recursos, e não é a toa que até hoje aprendemos equação de lotka-volterra, já que ela desdobra em muitas outras coisas, como aqui a teoria de metapopulações, só que simplificado já que não olhamos a abundância em cada população, apenas se ela existe ou não, 0 ou 1.

Mas aqui pensamos como o consumo sensato de recursos naturais é plenamente possível, sabendo como as populações do recurso, seja uma planta ou animal se comportam, podemos prever o ponto de equilíbrio, e quanto mais perto do equilíbrio menor o crescimento, então podemos estabelecer um meio termo entre manter a população abaixo do equilíbrio o suficiente para ela ter sempre um crescimento alto, já que quanto mais perto do equilíbrio, menor o crescimento. Claro que isso só vai funcionar dependendo dos parâmetros de colonização e extinção da espécie, como ja falamos que se a colonização é menor ou igual a extinção a espécies esta é caminhando para a extinção.

Ainda bem que sso já é bem estabelecido e conhecido como “Maximum sustainable yield”, e talvez a única forma de escapar da tragédia dos comuns, que ao que tudo indica hoje em dia, estamos falhando fortemente em sérios casos, como com a pesca.

Estimando a colonização e extinção a partir de dados de censo

Ja vimos em vários posts do que se trata a teoria de metapopulações e como podemos aprender muitas coisas sobre uma espécie modelando ela como uma metapopulação. No entanto eu mesmo já me frustrei muito ao tentar fazer a ligação entre dados e teoria.
Por exemplo como estimar as taxas de colonização e extinção de uma espécie a partir de dados que podemos coletar? Que tipo de dados precisamos e principalmente onde nesses dados esta essa informação que pode nos gerar tantas conclusões.
A resposta é que se a gente coleta várias vezes em um lugar a informação esta la e é mais simples do que pelo menos o que eu esperava.
Olhando esse diagrama fica mais simples, vamos pensar em aranhas em plantinhas como estávamos fazendo nas simulações anteriores.
Então temos:

Na primeira coleta eu olho várias plantas, mas somente duas coisas podem acontecer, ou eu vejo a aranha la ou eu não vejo.
Então se eu olho 10 plantinhas e 5 tem aranhas e 5 não tem aranhas, pronto 50% de ocupação, metade das plantas estão ocupadas por aranhas.
Mas e quanto a colonização? Ue, quando eu voltar para a segunda coleta, se tinha 5 plantas não ocupadas, e agora uma dessas 5 plantas foi colonizadas, são 20% de colonização, 1 colonização de 5 plantas.
E a extinção vai ser o mesmo esquema, se eu tenho 5 plantas que tinham aranhas, e vejo uma planta que tinha aranha e passou a não ter aranha, pronto 20% de chance de extinção, de 5 plantas que tinham aranha, tem uma agora que não tem.
Mas mesmo plantas que não mudam entre uma coleta e outra são informação sobre essa colonização/extinção.
Bem se eu tenho 5 plantas colonizadas, e um aparece na próxima coleta sem aranhas, pimba, 80% de sobrevivência, e o que 1-sobrevivencia é a extinção, ou seja 20% de extinção. A outra coisa é que se eu tinha 5 plantas disponíveis para serem colonizadas e somente uma o foi, eu tenho a chance de não-colonização, ou seja no final das contas eu retirei informação de todas as plantas no intervalo entre 2 coletas. Logos essas taxas valem para esse intervalo.

Então se a gente pode fazer um censo todo ano em várias plantinhas, e ir anotando para cada planta quem esta colonizado e quem não esta colonizado, quem esta sem aranhas. Temos toda a informação que precisamos. So manter a identidade de cada plantinha.

No final nossos dados vão se parecer com uma matriz assim:

> head(yy) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 1 0 0 1 1 1 1 0 0 [2,] 0 0 0 0 0 0 0 0 0 0 [3,] 1 1 1 1 1 1 1 1 1 1 [4,] 0 0 1 1 0 0 0 1 1 1 [5,] 0 1 1 1 1 1 1 1 1 0 [6,] 1 0 0 1 1 1 1 1 1 1

Linhas são as plantas, colunas as coletas, aqui só um pedacinho dos dados que simulamos. As 6 primeiras plantas.
Nesse caso simulamos os parâmetros que sempre viemos usando, 20% de chance de colonização e 10% de chance de extinção.
O que termina com 50% de ocupação.

Mas e agora o que fazer com os dados?
Uma opção no R é usar o pacote unmarked, ele tem várias funções para trabalhar com ocupação e metapopulações.
Na verdade transformar aquilo que falamos em números exige uma modelagem um pouco mais complexa que anovas da vida. Mas usando a distribuição binomial do cara ou coroa, isso é possível, afinal os dados dizem respeito a ver ou não as aranhas nas plantas. Alias os autores desse pacote são insanos, eu chutaria que metade dos livros de estatística bayesiana voltados para ecologia ou biologia são deles, mas aqui ainda não entramos na estatística bayesiana, ajustamos os parâmetros com maior verosimilhança, ou likelihood, que já vimos a algum tempo atrás.

Mas a partir dos dados, vamos estimar então a colonização e extinção de forma simples. Sem covariaveis que podem tornar tudo bem complexo de imaginar.

Primeiro transformamos nossa planilha de observações no formato que o pacote usa…

> simUMF <- unmarkedMultFrame(y = yy,numPrimary=T)

Nesse momento, nos também declararíamos covariaveis, como isolamento, tamanho das plantas e qualquer variável que pode afetar a colonização e extinção, já que esperamos que por exemplo o isolamento faça o acesso se tornar difícil, logo diminua a chance de colonização. Mas não trataremos disso aqui, imaginemos que todas as plantas são perfeitamente iguais.
Outra coisa que achei até maluco a primeira vez, foi como é possível lidar com a detecção imperfeita na coleta, por exemplo num arbusto, ele pode não estar com aranha porque elas não estão la ou ainda porque falhamos em ver elas. Estimar esse erro é possível aqui, o que traz estimativas de colonização e extinção muito melhores, mas ainda assim não simulamos isso.
A nossa simulação foi simplesmente colonização de 20%, extinção de 10% com detectabildiade das aranhas de 100%, vamos ver se realmente conseguimos recuperar esses parâmetros dos dados, usando a função colext() do pacote unmarked do jeito que falamos la em cima.

> summary(simUMF) unmarkedFrame Object 250 sites Maximum number of observations per site: 10 Mean number of observations per site: 10 Number of primary survey periods: 10 Number of secondary survey periods: 1 Sites with at least one detection: 234 Tabulation of y observations: 0 1 988 1512 0

Antes ainda podemos dar um summary pra ver o que criamos, e vemos com o que estamos lidando, são 250 plantas que coletamos durante 10 anos, so visitamos essas plantas 1 vez por coleta, e 234 plantas vimos aranhas pelo menos uma vez nelas.

> #estimando colonização e extinção > m1 <- colext(psiformula = ~1, # Ocupação do primeiro Ano + gammaformula = ~ 1, # Colonização + epsilonformula = ~ 1, # Extinção + pformula = ~ 1, # Deteção + data = simUMF)

Usamos o comando colext(), e nela poderíamos declarar tudo que achamos que pode influencias na colonização, extinção e detectabilidade, e na ocupação inicial, deixamos tudo 1, ou seja sem influencia nenhuma de covariaveis.
Mas assim obtemos o resultado.

> m1 Call: colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1, pformula = ~1, data = simUMF) Initial: Estimate SE z P(>|z|) 0.0164 0.127 0.13 0.897 Colonization: Estimate SE z P(>|z|) -1.43 0.0855 -16.8 5.56e-63 Extinction: Estimate SE z P(>|z|) -2.18 0.0936 -23.3 2.36e-120 Detection: Estimate SE z P(>|z|) 7.8 6.21 1.26 0.209 AIC: 2125.146

E ai estão as estimativas de Colonização e Extinção. Usando a função logística podemos obter os resultados em porcentagem, que é o que sempre falamos e que é mais fácil de pensar, então destransformando temos:

> plogis(m1@opt$par) [1] 0.5041035 0.1928461 0.1012372 0.9995916

Ohhh que lindo, ocupação inicial de 50%, praticamente 20% de colonização, levando em conta o erro para estabelecer o intervalo de confiança ali em cima, mas olhando também a extinção temos 10%, o que usamos na simulação, ohhhhhhhh muito bom, e como não simulamos falhas na detectabilidade, ai esta ela 100% a mão.
Realmente é muito legal ver dados saindo da forma bruta para esses parâmetros, que podem nos gerar tantos insights e conclusões de como essa espécie esta agora e qual o seu futuro. E é um exercício também para se aprender que nem toda analise de dados te da um valor p, aqui não tem nenhum.

Bem segue o código usado, essa simulação de dados esta no vignnete do pacote, apenas mudei os parâmetros e tirei as covariaveis.
Mas os vignnetes desse pacote são especialmente bem escritos diga-se de passagem, pelo menos bem mais simples de serem entendidos por mentes não inclinadas a matemática.

Bem e tem também material muito legal aqui nesse site da aula do Richard Chandler. Alias das apresentações dele que eu roubei aquele diagrama do inicio do post.

#######################################################
### Estimando Colonização e Extinção a partir de dados#
#######################################################
set.seed(666)
 
M <- 250 # Número de sítios
J <- 1 # Coletas dentro de um periodo
T <- 10 # Periodos de Coleta
 
psi <- rep(NA, T) # Probabilidade de Ocupação
muZ <- z <- array(dim = c(M, T)) # Ocorrencia observada e verdadeira
y <- array(NA, dim = c(M, J, T)) # Historico de Detecção
 
psi[1] <- 0.5 # Ocupação Inicial
p <- c(1,1,1,1,1,1,1,1,1,1) # Detectabilidade
phi <- runif(n=T-1, min=0.9, max=0.9) # Sobrevivencia (1-extinção)
gamma <- runif(n=T-1, min=0.2, max=0.2) # Colonização
 
z[,1] <- rbinom(M, 1, psi[1])
for(i in 1:M){
for(k in 2:T){
muZ[k] <- z[i, k-1]*phi[k-1] + (1-z[i, k-1])*gamma[k-1]
z[i,k] <- rbinom(1, 1, muZ[k])
}
}
for(i in 1:M){
for(k in 1:T){
prob <- z[i,k] * p[k]
for(j in 1:J){
y[i,j,k] <- rbinom(1, 1, prob)
}
}
}
for (k in 2:T){
psi[k] <- psi[k-1]*phi[k-1] + (1-psi[k-1])*gamma[k-1]
}
yy <- matrix(y, M, J*T)
 
head(yy)
#Transformando para a classe da dos usado pelo Unmarked
simUMF <- unmarkedMultFrame(y = yy,numPrimary=T)
summary(simUMF)
 
#estimando colonização e extinção
m1 <- colext(psiformula = ~1, # Ocupação do primeiro Ano
gammaformula = ~ 1, # Colonização
epsilonformula = ~ 1, # Extinção
pformula = ~ 1, # Deteção
data = simUMF)
 
#estimativas
m1
 
#transformando para probabildiade
plogis(m1@opt$par)

Explorando mais um pouco a teoria de metapopulações.

Após a primeira exploração ao modelo de Levins, podemos ainda tentar uma outra abordagem mais “visual” que é resolver o modelo dele para vários valores de colonização e extinção e ver o que acontece. Quando eu digo visual, é porque tudo pode ser deduzido diretamente da equação. Equações de primeiro grau são retas, de segundo grau são tortas e assim vai. Equações logísticas sempre tem essas forma além de outras características. Mas eu particularmente sou péssimo para entender isso, então eu tento resolver a equação para vários valores diferentes dos parâmetros dele e vejo o que acontece, como primeira exploração acho muito legal fazer isso.

Aqui começaremos sempre com a população com 50% de ocupação e usaremos os seguintes valores, 10%, 20%, 40% e 80% para ocupação e extinção, esses valores eu escolhi, eles não tem uma motivação especial. Bem isso da 16 possibilidades, com colonização 10% e extinção 10%, colonização 20% e extinção de 10% e assim até completar todas as 16 possibilidades, isso vai dar 16 gráficos que montaremos numa tabelinha de gráficos.

Recapitulando o que vimos no post anterior, vimos que quando a colonização é maior que a extinção não teremos uma extinção, já quando a colonização é igual a extinção a ocupação ira diminuir ao longo do tempo, e essa diminuição será mais intensa quando a extinção é maior que a colonização.

Então agora resolvendo tudo veremos:

Primeira vamos pensar como os gráficos foram organizados.
Da esquerda para direita, sempre a colonização esta aumentando, e de cima para baixo a extinção esta aumentando.
Logo, na diagonal são sempre valores de colonização e extinção iguais. E nessa diagonal vemos como a tendência ao longo do tempo é a extinção, mas conforme a extinção é maior, ela fica mais intensa, comparem os gráficos mais extremos onde a colonização (ci) e extinção (e) são 10% cada e 80% cada.

Outra coisa que vemos devido a essa organização dos gráficos é que fora dessa diagonal onde “ci” e “e” são sempre iguais, temos na parte superior casos onde a “ci” é maior que a “e”, na parte inferior casos onde “ci” é menor que “e”.

Na diagonal inferior vemos a confirmação da outra afirmação do post anterior que é sempre a tendência é extinção caso a “ci” seja menor que “e”. E vemos como essa tendência é mais forte, o declínio da ocupação ocorre mais rápida conforme a extinção é muito maior que a colonização.

E por último, vamos olhar a diagonal superior, onde sempre o valor de colonização é maior que a extinção.Primeira coisa legal é olhar como a ocupação se manteve em 50%, igual o valor inicial quando a colonização era o dobro da extinção, ou seja com colonização e extinção respectivamente com os valores de 20% e 10%, 40% e 20%, 80% e 40%, sempre a ocupação se manteve em 50%, ou seja para a gente ver como a ocupação depende de um “trade off” da Colonização versus extinção. Nos outros três casos vimos a ocupação aumentar rapidamente.

Bem essas são mais algumas características das metapopulações de Levins visto de um jeito, podemos dizer um pouco mais heurístico, de sair resolvendo a equação para tudo quanto é valor e ver o que acontece.

E no próximo post, chega de teoria e veremos como transformar observações e censos em nossas parâmetros vistos aqui, colonização e extinção.

###############################################################
#Explorando os possiveis resultados da metapopulação de Levins#
###############################################################
#Função de metapopulações de Levins
levins <- function(t, y, parms) {
p <- y[1]
with(as.list(parms), {
dp <- ci * p * (1 – p) – e * p
return(list(dp))
})
}
 
library(deSolve)
 
#Parametros usados
valores<-c(0.10,0.20,0.40,0.80)
#criando uma matrix com todas as combinações
prms<-expand.grid(ci=valores,e=valores)
 
#Ocupação Inicial (50%)
Initial.p <- 0.5
 
#Criando uma matriz para armazenar os resultados
out<-matrix(NA,nrow=100,ncol=16)
 
#Resolvendo para 100 eventos
for(i in 1:nrow(prms)) {
out[,i]<-ode(y = Initial.p, times = 1:100,func = levins,parms = prms[i,])[,2]
}
 
#gráfico
par(mfrow=c(4,4),mar=c(4,4,2,2))
for(i in 1:16) {
plot(out[,i]~c(1:100),type="l",ylim=c(0, 1),xlab="Tempo",ylab="Ocupação",lwd=3)
text(60,0.95,labels=c(paste("Ci=",prms[i,1]," e=",prms[i,2],sep="")),cex=1)
}