Biogeografia de ilhas

MacArthur e Wilson propuseram a teoria de biogeografia de ilhas onde o número de espécies de uma ilha oceânica seria uma função da taxa de imigração e extinção de espécies. O número de espécies em qualquer tempo é um equilíbrio dinâmico, resultando de ambos os lentos processos de extinção e continua chegada de novas espécies. Assim, é esperado que a composição de espécies pode mudar ao longo do tempo, gerando um “turnover” de espécies.

Vamos pensar da taxa de imigração, y, como uma função do número de espécies, x, que ja está em uma ilha. Essa relação tem uma inclinação negativa, porque conforme o número de espécies aumenta, a chance de um novo individuo chegar ser de uma nova espécie fica mais baixa, lembre-se de como é a distribuição das abundancias das espécies aqui. Assim, a imigração vai ter seu maior valor quando a ilha está vazia, quando o x é igual a zero, e vai cair para zero quando todas as espécies possíveis tiverem chegado la. Além disso, a inclinação deve desacelerar porque algumas espécies são muito mais prováveis que outras de imigrar. Isso significa que a taxa de imigração decai abruptamente assim que as espécies mais prováveis de imigrar chegam, e depois só aquelas que tem pouca chance de dispersar, ou baixas abundancias estão faltando chegar. Denovo pensando nas distribuições de abundancias de espécies, como Preston observou, algumas espécies são simplesmente mais abundantes que outras, produzindo mais indivíduos para dispersar, sendo “melhor dispersores”
Ou seja, esperamos algo desse jeito.

01

Agora vamos imaginar a taxa de extinção, y, como uma função do número de espécies, x. Essa relação vai ter uma relação positiva com o número de espécies, de tal forma que a probabilidade de extinção aumenta com o número de espécies, fácil pensar nisso, quando mais espécies, mais a gente vai adicionando a chance de extinção delas a taxa geral de extinção.

Essa relação também é predita de ter uma inclinação cada vez maior, mais ou menos pelos menos motivos do declínio cada vez maior da imigração. Algumas espécies são menos comuns que outras, e assim mais fácil de se tornarem extintas devido a estocasticidade demográfica e ambiental, tipo um drift genético, e segundo algumas espécies vão ter um menor fitness por muitas razões. Com o acumulo de espécies, maior sera o número de especies propicias a extinção (em relação a baixa abundancia e fitness) presentes, e que vão se extinguir. Lembrando que estamos falando de extinção localmente. Porque lembro que uma vez que apresentei um seminário sobre meta-populações,a palavra extinção soava agressiva para as pessoas, mas veja que extinção aqui não significa o fim da especie, só o fim local.
Então localmente, a taxa de extinção vai ser algo assim:

02

Assim, a taxa de mudança do número de especies em uma ilha vai ser \Delta R, que vai ser a diferença entre imigração I e extinção E, ou

\Delta R = I-E

Quando \Delta R=0, nos temos um equilíbrio. Se soubermos a forma quantitativa da imigração e extinção, nos podemos resolver para o equilíbrio. Esse equilíbrio vai ser um ponto no eixo x, onde as duas taxas se cruzam.

Na teoria do MacArthur e Wilson de biogeografia de ilhas, essas taxas podem ser determinadas pelo tamanho das ilhas, onde:

  • Ilhas maiores tem menor taxa de extinção por causa do maior tamanho médio das populações
  • Ilhas maiores tem maior taxa de colonização porque são alvos maiores para a dispersão

A distância entre as ilhas e as fontes de propágulos também influência essas taxas da seguinte forma:

  • Ilhas próximas do continente tem maiores taxas de colonização de novas espécies porque mais propágulos podem chegar, eles não morrem na jornada por exemplo
  • O efeito da área seria mais importante em ilhar longe do continente do que ilhas perto do continente, para reverter o continuo processo de extinção

Derivar uma curva de imigração ou extinção é um bocado complicado. Mas para ilustrar a teoria, a gente pode assumir que a taxa de imigração I, pode ser representada como uma função exponencial negativa e^{I_0-iR}, onde I_0 é a taxa de imigração em uma ilha vazia e i é o efeito negativo por espécies, que multiplica R que é quantas espécies tem la.

Então a gente pode criar um exemplo, com o número de espécies na ilha:

1
exemplo<-data.frame(especies=seq(0.1,50,0.2))

E ver a taxa de acordo com o número de espécies:

1
2
3
4
##Imigração
IO<-log(1)
b<-0.1
exemplo$imigracao<-exp(IO-b*exemplo$especies)

Agora para a extinção, E, está precisa ser zero se não houver espécies presentes. Imaginando que a extinção é uma função da densidade e que a densidade média declina conforme o número de espécies aumenta, então \bar{N}=\frac{1}{R}. então a gente vai ter algo mais ou menos assim

e^{dR}-1

Nos subtraimos esse 1 para ter certeiza que E=0 quando R=0.
O número de espécies, R, vai resultar de \Delta R = 0 = I-E, o ponto onde as linhas se cruzam.

I=e^{I_0-iR}
E=e^{dR}-1
\Delta R = 0 = I-E

03

O exato ponto onde essas linhas se cruzam podem ser encontrada de muitas formas, a gente tem visto isso em evolução aqui, mas vamos criar uma função no R para minimizar, nos vamos minimizar (I-E)^2, porque elevando ao quadrado a diferença nos temos uma quantidade com a conveniente propriedade que o mínimo ira ser aproximado tanto pelo lado positivo como o negativo.

1
2
3
deltaR<-function(R){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}

Nos usamos essa função de um parâmetro no optimize, e especificamos que queremos saber o ótimo em um intervalo de 1 a 50.

> optimize(f=deltaR,interval = c(1,50)) $minimum [1] 16.91321 $objective [1] 2.813237e-13

A saída nos diz que o mínimo é algo em torno de 16.9, e podemos colocar isso na figura

04

Agora suponha uma ilha um pouco mais longe, do mesmo tamanho, mas mais longe do continente, veja que a taxa de imigração diminuirá, porque pensa, as espécies vão morrer antes de alcançar a ilha e isso vai resultar um outro ponto mínimo, outro ponto de equilíbrio.

05

A beleza dess teoria é que ela foca nossa atenção em processos a nível de paisagem, geralmente fora dos limites espaciais e temporais das nossas capacidades de amostragem. Ele demostra que qualquer fator que ajude a determinar taxas de imigração ou extinção, incluindo a area de uma ilha ou a sua proximidade a uma fonte de propágulos, vai alterar o equilíbrio do número de espécies. Mas temos que lembrar que a identidade das espécies pode mudar ao longo do tempo, isso é, um turnover de espécies, mas um turnover devagar, porque as espécies que tem uma grande chance de imigrar e menor chance de extinguir vão se manter no ambiente muito provavelmente.

Bem é isso ai, incrível como as teorias mais legais são extremamente simples, mas com um efeito profundo na forma como pensamos sobre as coisas, 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:
M.H.H. Stevens 2011, A Primer of Ecology with R

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
library(ggplot2)
 
exemplo<-data.frame(especies=seq(0.1,50,0.2))
##Imigração
IO<-log(1)
b<-0.1
exemplo$imigracao<-exp(IO-b*exemplo$especies)
 
ggplot(data=exemplo,aes(x=especies,y=imigracao))+geom_line()+
    xlab("Número de espécies") +
    ylab("Taxa de imigração")
ggsave("01.jpg")
##curve(exp(IO-b*x),0,50,xlab="n. de espécies",ylab="Taxa Imigração")
 
##Extinção
d<-0.01
exemplo$extincao<-exp(d*exemplo$especies)-1
ggplot(data=exemplo,aes(x=especies,y=extincao))+geom_line()+
    xlab("Número de espécies") +
    ylab("Taxa de extinção")
ggsave("02.jpg")
##curve(exp(d*x)-1,0,50,xlab="n. de espécies",ylab="Taxa Extinção")
 
 
##Os dois juntos
exemplo2<-stack(exemplo[,c("imigracao","extincao")])
colnames(exemplo2)<-c("valores","tipo")
exemplo2$especies<-rep(exemplo$especies,2)
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração"))+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("03.jpg")
##
deltaR<-function(R){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}
estavel<-optimize(f=deltaR,interval = c(1,50))
segmento<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,2:3],1,min))))
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração"))+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento[2,],colour="black")+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("04.jpg")
 
exemplo<-data.frame(especies=seq(0.1,50,0.2))
exemplo$isp1<-exp(IO-0.1*exemplo$especies)
exemplo$isp2<-exp(IO-0.14*exemplo$especies)
exemplo$extincao<-exp(d*exemplo$especies)-1
exemplo2<-stack(exemplo[,c("isp1","isp2","extincao")])
colnames(exemplo2)<-c("valores","tipo")
exemplo2$especies<-rep(exemplo$especies,3)
 
##
deltaR<-function(R,b){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}
b<-0.1
estavel<-optimize(f=deltaR,interval = c(1,50),b=b)
segmento1<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,c(2,4)],1,min))))
b<-0.14
estavel<-optimize(f=deltaR,interval = c(1,50),b=b)
segmento2<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,c(3,4)],1,min))))
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração Sp1","Imigração Sp2"))+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento1,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento1[2,],colour="black")+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento2,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento2[2,],colour="black")+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("05.jpg")

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

Teoria neutra na evolução e na ecologia

A teoria neutra, na evolução, começou em 1960 com Motto Kimura, basicamente diz que a frequência de alguns alelos podem variar por acaso devido a esses alelos terem um sucesso reprodutivo, um fitness neutro, ambos tem o mesmo valor de fitness.

Bastante simples, mas com implicações fortes. Lembrando que alelos são as formas de expressão de um gene.
Se a gente pensar na natureza e na espécie modelo mais comum em genética, a mosquinha Drosophila melanogaster que tiver, uma mutação, por exemplo, asas curtas.

08
08

Uma mosca, na natureza, sem a capacidade de voar está ferrada, ela vai morrer de fome. Ela vai ter um fitness muito menor que a mosca com asas normais e a frequência desse gene vai tender a sumir, isso é conhecido como seleção natural.

Agora talvez, outra mutações não tenham impacto direto no fitness. Estou supondo aqui, a titulo de exemplo, ja que eu não entendo de Drosophilas, mas a cor dos olhos, se não houver seleção natural qualquer, as fêmeas não verem moscas de olhos vermelhos como mais bonitas que as de olhos laranjas, o que aconteceria?

08
08

Bem essa mutação poderia ser considerada neutra, ou seja ela num melhora nem piora o fitness, mas isso não implica em nenhum momento que se você acompanhar essa população, a frequência desse alelo sera constante. Isso porque temos uma população, na natureza sempre finita, a a frequência é sempre o número de indivíduos que exibe uma característica, o alelo em questão dividido pelo tamanho total da população.

Nesse caso que temos 2 alelos em um gene (situação fictícia, eu estou assumindo que essa mutação e neutra e não existe mais variantes, sei que tem, mas é só para o exemplo ser fácil de entender), os alelos são, olhos vermelhos e olhos laranjas.

Mas esse exemplo é só para visualizar a situação, a ideia é que se esse gene é neutro, existe um processo onde se você começa cada alelo presente em 50% da população, 50 indivíduos, vamos supor, a qualquer momento que alguém usar um mata moscas aleatoriamente e mate, por cagada, 9 moscas de olhos vermelhos e 1 mosca da olhos laranjas, sobrara 41 moscas de olhos vermelhos a 49 de olhos laranjas, e se cada uma tem sei la, 1 filho, na próxima geração, veremos 45,5% de moscas de olhos vermelhos e 54.4% de moscas de olhos laranjas.

O problema é que quando eu bato o mata moscas, eu tento mata moscas, eu não fico olhando qual mosca eu vou matar, e pior que isso, nenhum tipo de moscas é mais rápido, ou mais ágil para fugir nesse exemplo, o que gera uma situação onde essa mudança na frequência dos alelos é completamente imprevisível, pra que direção ela vai.

Essa é um processo conhecido como deriva genética, ou genetic Drift, um efeito da teoria neutra.

Existem várias formas de simular a deriva genética, lembrando que a herança pode ser variável mesmo na gente. Lembre-se que mitocôndrias a gente so ganha da nossa mãe enquanto para genes nucleares alelo vem do nosso pai e outro alelo vem da nossa mãe, e tem crossing-over. Se começarmos a pensar muito nisso a gente não vai pra frente. Mas vamos usar aqui só pra visualizar o modelo do Wright–Fisher, que considera gerações discretas, indivíduos diploides, e cada nova geração tem um tamanho fixo, basicamente cada nova geração a gente pega um número aleatória de indivíduos da frequência da geração anterior.

01

Veja que cada linha é uma simulação, no eixo x temos as gerações, e no eixo x estamos representando a frequência de um alelo, lembrando que as duas frequências vão dar 1, então se um alelo tem 100% o outro tem 0, e vice versa.

Veja que em alguns casos, laranja e verde por exemplo, esse alelo domina, e o outro some, em outras simulações, ele some e só fica o outro alelo na população.

Apesar de imprevisível, a gente pode dizer que a deriva genética acaba sempre mudando as frequências de alelos, ainda mais, quase sempre sumiu com algum dos alelos, ela tirou alguém da jogada, e assim diminuiu a diversidade genética, já que ter uma população com 2 alelos é mais diversa que uma população com 1 alelo, seja ele qual for.

Certo, mas ai vendo tudo isso, um cara chamado Stephen P. Hubbell pensou se isso não poderia acontecer em outra escala, na escala de espécies. Será que existe uma teoria neutra para riqueza de espécies? Então ele pensou na teoria neutra unificada da biodiversidade e biogeografia (the unified neutral theory of biodiversity and biogeography), mas vamos dizer teoria neutra que é menor, ninguém tem saco de escrever esse nome por muito tempo.
A sua importância é que ela possibilita mais que muitos modelos de comunidade, ela faz predições claras em muitos níveis de organização, tanto de processos evolucionários como ecológicos.

Comunidades e genes são computacionalmente e matematicamente relacionadas, geralmente usando modelos iguais aos usados na genética, logo as lições que aprendemos sobre a deriva genética se aplicam a modelos de comunidade. A dinâmica ecológica numa escala local é conhecida de como deriva ecológica, e populações mudam ao acaso, devido a processos estocásticos como o nosso mata-moscas.

Hubbell propôs a teoria neutra como um modelo nulo para a dinâmica de indivíduos, para explicar a abundância de árvores tropicais ao na costa rica, barro colorado. No pacote vegan, se você abrir o pacote e digitar data(BCI), você pode ver parte dos dados que o Hubbell trabalhava, as abundâncias na floresta da reserva do Barro Colorado na costa rica, alias, como curiosidade, se você digitar data(), sem nenhum argumento, o R lista todos os data-set de exemplos disponíveis, de todos os pacotes abertos.

Certo, mas agora o controverso é que o Hubbell propôs, diferente do que comumente associamos a modelos nulos (que são como as coisas acontecem ao acaso normalmente, para teste de hipoteses), que é realmente assim que as coisas acontecem na dinâmica de comunidades.

Mas isso é controverso? Porque?

Porque estamos constantemente avaliando fatores que influenciam na riqueza e diversidade de espécies, e vem um cara e fala que tudo isso num é assim, que quem define quais espécies estão num local é o acaso. Realmente uma afirmação que vai contra muitos trabalhos que estamos acostumados a ver.

Mas se a gente não for muito extremista (e começar a achar que tudo tem que ser ou 100% estocásticos ou 100% determinísticos, sem acordos), a palavra chave para a teoria neutra do Hubbell é metacomunidade, uma coleção de comunidades similares conectadas por dispersão.

A metacomunidade é povoada inteiramente por indivíduos que são funcionalmente idênticos, segundo a teoria neutra. A teoria neutra fica então uma teoria de dinâmica de indivíduos, modelando nascimentos, morte, migração e mutação desses indivíduos, sendo a riqueza de espécies uma medida derivada desses parametros. Assumimos que dentro de uma guilda, como a de árvores do final de sucessão em floresta tropical, as espécies são essencialmente neutras, todo mundo tem um fitness igual, equivalente. Isso significa que probabilidades de nascimento, morte, mutação e migração são próximas ou equivalentes para todos os indivíduos. Assim, mudanças nos tamanhos das populações ocorrem devido a “random walks”, eventos estocásticos. Como a gente viu nas frequências gênicas. A gente pode simular praticamente igual fizemos ao gene, mas vamos usar aqui o pacote untb, que já tem essa simulação pronta.

02

Mas calma, isso não implica necessariamente na ausência de competição ou outros efeitos indiretos mediados pela dependência de densidade. A competição é vista como difusa e igual entre indivíduos.
Efeitos dependentes da densidade aparecem através de formas especificas no número total de indivíduos da comunidade, ou como características dos indivíduos relacionado as probabilidades de natalidade, mortalidade e especiação.

Um paradoxo gerado na teoria neutra de Hubbell é que na ausência de migração ou mutação, a diversidade declina para zero, riqueza para um, uma monodominância. Uma caminha aleatória (random walk), devido a equivalência de fitness, vai eventualmente resultar na perda de todas as espécies, exceto uma. Entretanto a perda de diversidade em uma comunidade local é esperada de forma muito, muito devagar, e é contra-balanceada por imigração e especiação. Assim, espécies não aumentam deterministicamente quando raras, isso faz o conceito de coexistência diferente do conceito que a gente vê nos modelos de Lotka-Volterra. Coexistência no caso da teoria neutra é um efeito estocástico balanceado pelo aparecimento de novas especies localmente. Lembrando que estamos pensando dentro de uma guilda, um nicho, não estamos comparando gramas com árvores aqui, ou onças com aranhas.

Mas se todas as comunidades fazem um random walk para a monodominâncias, como a diversidade é mantida em qualquer região em particular? Dois fatores mantém espécies em qualquer comunidade local. Primeiro, imigração para a comunidade local de indivíduos da metacomunidade. Mesmo com cada comunidade sofrendo um random walk para a monodominância, cada comunidade pode vir a ser dominada por uma espécies diferente, formando o pool de espécies dessa metacomunidade, já que todas as comunidades apresentam fitness diferentes. Assim comunidades separadas são preditas de se tornarem dominadas por espécies diferentes, e essas diferenças entre espécies locais ajudam a manter a diversidade na paisagem ocupada pela metacomunidade.

Um segundo fator, numa escala de tempo muito maior, é a especiação que dentro da metacomunidade mantém a biodiversidade. Mutação e consequentemente, especiação prove uma forma de introdução de diversidade na metacomunidade. Random walks em direção a extinção em comunidades que são tao lentas que são contrabalanceadas por especiação.

E agora temos uma pequena ideia do que se trada a teoria neutra da ecologia, e como ela veio da teoria neutra da genética. As simulações foi mais pra fazer alguma figurinha que realmente explicar algo, mas post sem figurinha é chato, e ainda as figuras foram feitas no ggplot2, que faz figuras bem bonitas, mas que até hoje eu não aprendi a usar direito. o script está no repositório do github, além de aqui em baixo, e precisamos definitivamente de muitos posts ainda pra entender a teoria do Hubbell, até onde ela vai, e como estimar parâmetros em metacomunidades para confrontar dados com teoria nesse caso aqui, mas paciência é uma virtude, sempre.

Referencia:
M Henry H Stevens 2011 A primer of ecology with R. Springer

library(untb)
library(ggplot2)
library(reshape)
library(RColorBrewer)
set.seed(123)
 
###################################################
### Deriva Genética
###################################################
#Original - http://statisticalrecipes.blogspot.com.br/2012/02/simulating-genetic-drift.html
 
#Parametros da Simulação de Deriva Genetica
N = 20           # Número de individuos diploides
N.chrom = 2*N    # Número de cromossomos
N.gen = 100      # Número de gerações
N.sim = 10       # Número de simulações
p = 0.2
q = 1-p
 
# Simulation
saida <- matrix(0,nrow=N.gen,ncol=N.sim,dimnames = list(paste("Geraçao",1:N.gen),paste("S",1:N.sim)))
saida[1,] = rep(N.chrom*p,N.sim) # Inicializando o numero de alelos na primeira geração
for(j in 1:N.sim){
  for(i in 2:N.gen){
    saida[i,j]<-rbinom(1,N.chrom,prob=saida[i-1,j]/N.chrom)
    }
  }
#Transformando em frequencia
saida<-data.frame(saida/N.chrom)
 
# Modificando os dados para a estrutura que da para plotar
sim_data <- melt(saida)
ggplot(sim_data, aes(x = rep(c(1:100), N.sim), y = value, colour = variable)) +
    geom_line() +
    scale_colour_brewer(palette="Set3") +
    ggtitle("Simulação de deriva genética") +
    xlab("Geração") + ylab("Frequencia alélica") +
    ylim(0,1) +
    labs(colour = "Simulações")
#ggsave("01.jpg")
 
###################################################
### Deriva Ecologica
###################################################
saida <- untb(rep(1:10, 90), prob=0, gens=5000, D=9, keep=TRUE)
 
sim_data <- melt(data.frame(species.table(saida)))
 
ggplot(sim_data, aes(x = rep(c(1:5000), 10), y = value, colour = variable)) +
    geom_line() +
    scale_colour_brewer(palette="Set3") +
    ggtitle("Simulação de deriva populacional") +
    xlab("Geração") + ylab("Abundância da espécie") +
    labs(colour = "Simulações")
#ggsave("02.jpg")

As vezes, menos é melhor.

A gente ja viu vários posts com vários tipos de modelos, de poisson, binomial, mixtos e até expansões como overdispersion e zeros inflados.
Mas antes de mais nada, antes de sair tentando tudo e ver o que da mais certo. A gente tem que pensar, será que os parâmetros que estamos estimando fazem sentido. A gente consegue discutir aquele monte de número? Eles tem algum sentido biológico (para quem é da biologia)?

Vamos olhar um conjunto de dados, contagens de uma população ao longo do tempo, pensa que vemos o seguinte.

01

Não precisamos mentir. A maioria das pessoas (inclua eu nessa) pensaria em tacar uma regressão linear nesse conjunto de dados, um pensamento instantâneo, como se estivéssemos condicionados a fazer isso. Não que isso seja totalmente errado. Vamos ver como esse modelo ficaria.

02

Olha, não ficaria tão ruim, até que o modelo ajustaria. Mas e ai, o que estaríamos fazendo? Vendo como a população dessa espécie muda ao longo do tempo? E a inclinação é a intensidade como a população está mudando e o intercepto é como ela começou?

Eu acho que não é um interpretação tão ruim assim. Mas podemos melhorar mais e mais esse modelo, usando polinômios por exemplo.
Polinômios de segundo grau são uma parábola, de terceiro grau fica com uma curvinha, e cada vez o modelo vai ficando melhor, até que com 6 graus temos o seguinte.

03

Agora ta perfeito. Mas infelizmente, perfeito somente para esse grupo de dados, conforme a gente aumenta o número de parâmetros, a gente sempre melhora o ajuste, diminui o espaço para erro de medida, o problema é que o que os parâmetros fazem? O que significam? Nada mais, a gente so tem um ajuste bom para um modelo ininterpretável.
Olhe o que criamos, ele pode parecer bonito no grafico, mas não significa muito. E pense que fazemos modelos para pensar nos dados que coletamos, mas nos dados também que não coletamos, que outras pessoas vão coletar, e mais que isso, para entender como o mundo funciona, e será que temos muita capacidade de predição com esse modelo?
Se parar para pensar por um momento, a gente lembra que existem modelos de crescimento populacional que vão além de retas e linhas. Que cada letrinha, cada parametro tem um sentido especial que discutimos um bocado por aqui. Então se a gente passa tanto tempo lendo sobre eles, deveríamos gastar eles. até porque dai fica fácil discutir os dados.

04

Olha ae, talvez ele não seja tão bonito quanto ao ajuste, como o modelo polinimial, talvez até pior que o ajuste da regressão linear, mas é muito mais o que pensamos, e temos que pensar mais para frente, no tempo a frente.

05

Veja o que preveríamos para o futuro de acordo com cada modelo. O modelo de regressão linear preveria um crescimento constante, e o pior de tudo, o modelo polinomial preveria uma queda brusca na população. Agora você acha que tem algum evento obvio que faria a população diminuir. No fim das contas o modelo de crescimento populacional é o que faz a previsão mais razoável, que a população vai continuar crescendo. E é o que a gente acha. Os modelos de regressão são muito bons localmente, ali aonde temos dados, mas eles são péssimos fora da região onde estão os dados. Isso é algo que é bem comum, porque a regressão não assume muita coisa sobre o sistema biológico em si, o preço que se paga é isso.
Claro que se a gente não sabe nada, seria uma informação importante saber o que está acontecendo nesse período de tempo, e como está acontecendo, mas enquanto estamos construindo um conhecimento “bottom-up”, não sabemos nada, e estamos começando. Agora se a gente já tem modelos, a gente pode por eles a prova, algo mais “top-down”. A gente parte do que já sabe e vê se funciona na prática.

Agora quando falamos de chi-quadrado aqui, falamos também de como usar ele para seleção de modelos, mas aqui ferrou, porque para a seleção de modelos, precisamos de modelos aninhados. Por exemplo, modelos de Evolução de DNA temos vários modelos, mas um é uma versão mais simplificada de outro, um caso especial. A mesma coisa temos com os modelos de regressão linear, polinomial ou não. Agora o modelo de crescimento polinomial não aninha com os modelos de regressão, então acho que temos que escolher na argumentação, e não vai ter um valor de p dando sopa para justificar sem discussão essa escolha. Mas é a hora que pensar, apesar de no começo da curva de crescimento populacional (a vermelha) ser ruim no começo, lembramos que no começo a população é pequena, e coisas ao acaso acontecem mais comumente, so lembrar de genetic drift. Depois a população deve estabilizar mais e a curva fica melhor. Muito melhor porque claramente a população ta começando a crescer, e vai crescer bastante, e o modelo de crescimento populacional é o único que está dizendo isso.

Bem é isso ai, uma semana sem posts por conta de provas, mas agora voltando a atividades :), os scripts sempre no repositório recologia e é isso, até o próximo post.

tempo<-seq(0,12,by=2)
N<-c(230,1205,1150,2091,3010,6106,7323)
 
#figura 1
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,pch=19)
 
#modelo de regressão linear
modelo1<-lm(N~tempo)
modelo1
 
#figura 2
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,8000),pch=19)
curve(coef(modelo1)[1]+coef(modelo1)[2]*x,0,12,add=T,lwd=2,lty=1,col="green")
segments(x0=tempo, y0=N, x1 =tempo , y1 = predict(modelo1))
 
#modelo de regressão polinomial
modelo2<-lm(N~poly(tempo, degree=6, raw=TRUE))
modelo2
 
#figura 3
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,8000),pch=19)
curve(coef(modelo2)[1]+coef(modelo2)[2]*x+coef(modelo2)[3]*x^2+coef(modelo2)[4]*x^3+coef(modelo2)[5]*x^4+
      coef(modelo2)[6]*x^5+coef(modelo2)[7]*x^6,0,12,add=T,lwd=2,lty=1,col="blue")
segments(x0=tempo, y0=N, x1 =tempo , y1 = predict(modelo2))
 
#modelo não linear
modelo3<-nls(N~N0+exp(r*tempo),start=list(N0=1,r=1))
modelo3
 
#figura 4
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,8000),pch=19)
curve(coef(modelo3)[1]+exp(coef(modelo3)[2]*x),0,12,add=T,lwd=2,col="red")
segments(x0=tempo, y0=N, x1 =tempo , y1 = predict(modelo3))
 
#figura 5
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,10000),xlim=c(0,20))
curve(coef(modelo1)[1]+coef(modelo1)[2]*x,0,20,add=T,lwd=2,lty=1,col="green")
curve(coef(modelo2)[1]+coef(modelo2)[2]*x+coef(modelo2)[3]*x^2+coef(modelo2)[4]*x^3+coef(modelo2)[5]*x^4+
      coef(modelo2)[6]*x^5+coef(modelo2)[7]*x^6,0,20,add=T,lwd=2,lty=1,col="blue")
curve(coef(modelo3)[1]+exp(coef(modelo3)[2]*x),0,20,add=T,lwd=2,col="red")
points(N~tempo,pch=19,cex=1.2)
legend("topleft",col=c("green","blue","red"),lwd=2,bty="n",title ="Modelo",
       legend=c("Regressão Linear","Polinomial de grau 6","Crescimento populacional"))
 
logLik(modelo1)
logLik(modelo2)
logLik(modelo3)

Forrageamento ótimo: Como predadores ativos minimizam o tempo de busca?

Bem, vamos pensar um pouco sobre outra teoria bem legal em ecologia, a de como ser um bom predador.
Nem só de garras e presas vive um predador, a maioria dos predadores terrestres é limitada pela disponibilidade de alimentos a que está submetido. Isso significa que aqueles que conseguem mais que os outros da mesma espécie deixam mais descendentes (filhotinhos, prole, como quiser chamar) que aqueles indivíduos menos aptos a encontrar, capturar e matar suas presas.

Ao longo das gerações, o comportamento de forrageamento de uma população tende a ficar mais e mais eficiente. Varias adaptações podem contribuir para diminuir:

(1) o tempo gasto caçando
(2) a energia gasta caçando

Ou ambos, essa área de estudo é conhecida como Teoria do forrageamento ótimo (optimal foraging theory).

De forma simplificada, a gente pode definir duas categorias gerais de predadores: Os que fazem busca ativa (ficam perambulando procurando o que comer) e os que sentam e esperam (sit-and-wait).

Um predador do tipo busca ativa, se move através de seu habitat e captura as presas que ele quer quando precisa conseguir comida, como uma aranha que fica andando pela casa, pulando (Salticidae é conhecida como a aranha saltadora) atrás das presas.

Um predador que senta e espera, bem ele senta e espera, se posiciona em algum ponto de vantagem e espera por uma presa inocente passar ali por perto, como uma aranha parada numa flor (Uma aranha da família das Thomizidae por exemplo).

Aranhas são legais porque apresentam os dois tipos de estrategias entre suas várias famílias, mas vamos pensar primeiro num predador de busca ativa, cuja estratégia de forrageamento ideal é adquirir mais presas no menor espaço de tempo.

Considere uma aranha saltadora então, que se move através da serrapilheira a procura de insetos, larvas e outros pequenos invertebrados para comer. Sua atividade de forrageamento a torna mais vulnerável ​​a seus próprios predadores, como algumas aves por exemplo. Então quanto mais tempo uma aranha gasta forrageamento, maior a chance de morrer durante esse processo. Assim, uma aranha de maior sucesso (aquela que terá a barriga cheia para deixar mais filhotinhos) toma decisões que lhe dão mais presas na menor quantidade de tempo, já que cada segundo desperdiçado pode significar a morte.
Ela faz isso continuamente, decidir se deve parar e comer uma presa especial ou se vai seguir em frente em busca de algo melhor. A decisão é baseada em tempo:

Quanto tempo vai demorar para encontrar outra presa? Será uma presa diferente?
Quanto tempo vai demorar para capturar essa presa em comparação com o que ele já encontrou?

Estes dois componentes que formam o tempo total de forrageamento e são chamados de tempo de busca e manuseio.

Primeiro vamos pensar no tempo de busca, bem não é difícil visualizar que ele é inversamente proporcional à abundância das presas. Quanto mais presas, menos tempo se leva para encontrar uma presa.

A nossa amiga aranha, de certa forma, controla a abundância de suas presas, controlando o número de espécies de presas mantém no cardápio, um cardápio mais restrito pode ser mais raro enquanto um cardápio mais variado pode facilitar encontrar comida pela serrapilheira onde a dona aranha vive.

Se a dona aranha come apenas larvas, sua presa é bem menos abundantes e seu tempo de acaba sendo maior. Agora se no cardápio está incluso larvas e besouros por exemplo, a busca é mais rápida, dependendo da abundância de cada uma dessas presas na serapilheira.

Porém existe um tempo de manipulação da presa, para comer, e esse tempo depende do comportamento da presa. É mais rápido capturar uma larva que um besouro, podemos imaginar isso, porque um besouro está mais consciente da presença da aranha e pode se mover mais rápido, ou ter adaptações para sua defesa.

Mas a partir da discussão que começamos podemos fazer um modelo que represente o tempo médio que leva para encontrar uma presa, e esse será mais ou menos o inverso da abundância dessa presa, então:

{TeP} = {1\over{AP}}

Onde
TeP significa tempo para encontrar uma presa,
AP é abundância dessa presa.

Temos que usar umas siglas senão a formula fica muito grande e não cabe na tela. Mas está bem simples até aqui.

Certo, mas como gostamos mais de gráficos que de formulas, vamos visualizar isso num gráfico.

01

O que acontece, estamos simulando aqui o tempo de captura em relação a a abundância da presa, onde a abundância é expressa como ocorrência por minuto. Por exemplo, uma abundância média de larvas, 100 larvas por metro quadrado significa que a aranha leva 2 minutos mais ou menos para encontrar uma larva, isso da uma taxa de meia larva por minuto. Se tivermos 200 larvas por metro quadrado, bem abundânte, podemos ter uma taxa de quase 1 larva encontrada por minuto. O Contrario, sei la, 10 lavas por metro quadrado, pode dar uma taxa de 0.1 larvas por minuto, a aranha leva 10 minutos para achar uma larva. Esse valores eu estou inventando agora, veja que eu nem coloquei valores no eixo x do gráfico.
A ideia é que estamos modelando a vida da aranha, e estamos pensando sempre numa unidade para fazer comparação, e essa unidade é presas capturadas pelo tempo gasto. Como a gente sempre ta pensando numa unidade, podemos mais para frente fazer comparações, baseadas nessa unidade. Então essa unidade é como dinheiro, dinheiro sabemos que mais dinheiro é melhor, aqui, menos tempo é melhor. Temos que lembrar disso sempre, menos tempo é melhor.

Mas como dissemos, existe ainda outro componente, o tempo de manipulação da presa, e o tempo total de forrageamento torna-se então:

{TEC} = {{1\over{AP}} + TM}

Onde:
TEC é o tempo de encontrar e capturar a presa
AP é a abundância da presa
TM é o tempo que tem que ser gasto manipulando essa presa.

Que podemos olhar num gráfico assim:

01

Basicamente, somar aquele número ali, o tempo de manipulação aumenta o tempo geral para encontrar e capturar a presa, que no final das contas só pioras as coisas para a aranha, então um menor tempo de manipulação é melhor certo? So lembrar que nossa unidade é presa por minuto, e em qualquer abundância de presa ali, um tempo de manipulação faz gastarmos mais tempo total com aquela presa.

Mas estamos falando de um cardápio simples aqui, vamos supor dois itens de cardápio então, larvas e besouros.

Sabemos que, pelo que vimos até agora, temos que ter algo assim:

{TEC_l} = {{1\over{AP_l}} + TM_l}

Colocamos ali um elezinho subescrito, para lembrar que estamos falando de larva e

{TEC_b} = {{1\over{AP_b}} + TM_b}

um bezinho subescrito para lembrar que estamos falando do besouro.

Quando uma aranha decide comer apenas larvas, apenas besouros ou os dois?

Podemos expandir nosso modelo para considerar as duas presas juntas. Considere agora o tempo de forrageamento de uma aranha que decide deixar no cardápio larvas e besouros. A primeira coisa a considerar é que a abundância total vai ser a abundância de besouros mais a abundância de larva e consequentemente o tempo médio de busca vai ser menor.

{TEC} = {1\over{AP_l + AP_b}}

Certo, neste modelo temos o tempo que uma aranha leva para encontrar e capturar uma presa, mas qual a chance dessa presa ser uma larva ou um besouro?
Estamos perguntando mais especificamente: quando uma aranha encontra uma presa aceitável (larva ou besouro), qual é a probabilidade de que ela seja uma larva?

(A probabilidade, neste caso de duas presas, de que a presa seja um besouro vai ser apenas um menos a probabilidade de que ele é uma larva)

Estamos querendo saber a abundância de larvas em relação à abundância de todas as presas possíveis, larvas e besouros:

{p_l} = {AP_l\over{AP_l + AP_b}}

Achou que ia ser algo complicado né ^^. Não, a resposta é apenas a abundância de larvas divido pela abundância total. Para o besouro vai ser:

{p_b} = {AP_b\over{AP_l + AP_b}}

Ou , como dissemos la em cima:

{p_b} = 1 - p_l

Mas para que queremos essa informação?

Porque imaginamos que o tempo de manipulação seja diferente entre besouros e larvas, é mais fácil comer uma larva que um besouro, logo o tempo de manipulação para as larvas da dieta é:

 {TM_l} = { {AP_l\over{AP_l + AP_b}} \cdot{M_l}}

E nos podemos usar a mesma ideia para calcular o tempo para o besouro

 {TM_b} = { {AP_b\over{AP_l + AP_b}} \cdot{M_b}}

Agora é so juntar tudo isso, e temos o tempo médio nossa aranha leva para encontrar e capturar uma presa qualquer do seu cardápio.

{TEC} = {  {1\over{AP_l + AP_b}} + TM_l + TM_b }

E substituímos para cada caso, e ficamos com:

{TEC} = {  {1\over{AP_l + AP_b}} + { {AP_l\over{AP_l + AP_b}} \cdot{M_l}} + { {AP_b\over{AP_l + AP_b}} \cdot{M_b}} }

Pode parecer uma formula grande, mas continuamos ainda com a mesma coisa, podemos falar que o Tempo de encontrar e capturar uma presa depende da abundância os espécies do cardápio, aqui larvas e besouros, mais o tempo de manipular elas, mas ponderamos a manipulação pela abundância da presa em questão (larva ou besouro) em relação a todos as presa que a aranha pode encontrar.

Agora vamos fazer uma tentativa aqui, uma pequena simulação.
Vamos atribuir abundâncias aos vermes e besouros, sempre pensando a abundância de um grupo em relação ao outro, em vários cenários. Ai vamos procurar em qual cenário a aranha se daria melhor. Qual a melhor escolha de cardápio?

01

Bem, como o besouro tem um tempo de manipulação muito maior, a primeira coisa que vemos é que comer somente larvas é melhor que comer somente besouros. Mas algo talvez não tão obvio, é que em baixas abundâncias, manter larvas e besouros no cardápio é melhor que somente larvas, ja que a abundância de ambos é baixo, a aranha demora para encontrar algo para comer, então mesmo demorando para capturar o besouro, vale a pena parar e comer o besouro.

Mas veja que a linha vermelha, que representa somente larvas, se cruza com a linha preta, que inclui larvas e besouros no cardápio. Isso acontece porque nem sempre incluir o besouro será uma boa estratégia. Como assim? Se o tempo que você leva capturando um besouro for maior que o tempo para encontrar e manipular uma larva, a partir dai não vale mais a pena se dar o trabalho de comer besouros, porque no tempo de capturar um a gente consegue comer uma larva. E a partir dai, conforme aumentamos a densidade das larvas, passa sempre a valer a pena comer só larvas, veja que no final do eixo x, a linha vermelha está sempre mais perto de zero, que é o que queremos, mais presas por segundo.

Agora vamos manipular um pouco a abundância dos besouros, que está como a metade da abundância das larvas.

01

Agora fica mais fácil de ver o que estávamos comentando no gráfico anterior, veja que fomos aumentando a abundância do besouro para ver o que acontecia. Colocamos a abundância igual a abundância das larvas, depois o dobro, quatro vezes mais e ainda assim as linhas se cruzam no mesmo ponto. Que é o tempo de manipulação do besouro. Se a larva tem uma abundância muito alta, sempre vai acabar valendo a pena tirar o besouro do cardápio, agora veja que quanto maior a abundância do besouro, olhe a linha azul, passa a valer muito a pena incluir ele se a abundância de larva é media ou baixa. Quando a larva tem uma abundancia muito baixo, um besouro no prato garante a refeição de todos os dias.

Bem interessante, principalmente como não existe uma escolha simples, haverá sempre um balanço entre tempo que se passa manipulando uma presa e tempo para encontrar ela. E lembre que tempo de manipulação, a gente está falando de modo geral, mas em manipulação, inclua o tempo que você demora para digerir uma a presa, certos tecidos são bem mais difíceis de digerir que outros, e se seu estomago está cheio, não da para comer mais nada.

Ja a abundância pode mudar com a temperatura ambiental, besouros e larvas são animais ectotérmicos (de sangue frio), assim todo mundo vai ficar menos ativo no frio, mais ativo no calor, e podemos imaginar que a atividade dessas presas vai influenciar nas chances delas serem encontradas, então podemos incluir isso também como parte da “Abundância”.

E além de tudo isso, ainda temos o ganho energético proporcionado por cada tipo de presa, que não foi considerado aqui. Até aqui estamos dando pesos iguais para larvas e besouros, mas talvez um besouro seja mais gordinho e suculento, e isso compensa todo os problemas para capturar ele. Mas isso vai ficar para o próximo post.

Bem ficamos por aqui, o script está no repositório recologia, e até o próximo post.

Referência:
Bernstein, R. 2003 – Population Ecology – An Introduction to Computer Simulations. John Wiley & Sons, Ltd 159pp.

#Tempo para encontrar uma presa.
tep<-function(ap) {
    return (1/ap)
}
 
sequencia<-seq(0.2,1,0.01)
 
#
plot(sequencia,tep(sequencia),pch=19,type="b",frame=F,xaxt="n",
     ylab="Tempo para encontrar uma presa",xlab="Abundância da presa no ambiente")
axis(1,at=c(0.2,0.4,0.6,0.8,1),labels=F)
mtext("Muito raro", side = 1, line = 1, at = 0.2)
mtext("Raro", side = 1, line = 1, at = 0.4)
mtext("Mais ou menos", side = 1, line = 1, at = 0.6)
mtext("Abundante", side = 1, line = 1, at = 0.8)
mtext("Muito abundante", side = 1, line = 1, at = 1)
 
#
tepm<-function(ap,tm) {
    return ((1/ap)+tm)
}
 
#
curve(tep(x),0.1,1,frame=F,ylab="Tempo para encontrar uma presa",xlab="Abundância da presa no ambiente",xaxt="n")
axis(1,at=c(0.1,0.4,0.6,0.8,1),labels=F)
mtext("Muito raro", side = 1, line = 1, at = 0.1)
mtext("Raro", side = 1, line = 1, at = 0.4)
mtext("Mais ou menos", side = 1, line = 1, at = 0.6)
mtext("Abundante", side = 1, line = 1, at = 0.8)
mtext("Muito abundante", side = 1, line = 1, at = 1)
tm<-seq(1,5,1)
for(i in tm){
    curve(tepm(x,i),col=i,add=T,lwd=2,lty=2)
}
legend("topright",col=c(1,tm),lwd=c(1,rep(2,5)),lty=c(1,rep(2,5)),legend=paste("Tempo =",c(0,tm)),bty="n")
 
#
tiw<-function(aw,hw) {
    return(1/aw + hw)
    }
#Larva
hw<-1
aw<-seq(0.000,0.03,0.0005)
#Besouro
ab<-0.5*aw
hb<-60
 
plot(aw,tiw(aw,hw),pch=19,cex=0.5,type="b",col=2,frame=F,xlab="Abundância",ylab="Tempo total manipulação",
     ylim=c(0,200),xlim=c(0,0.03))
points(aw,tiw(ab,hb),pch=19,cex=0.5,type="b",col=3)
points(aw,tiwb(aw,ab,hw,hb),type="l",col=1,lwd=1,lty=3)
 
legend("topright",col=c(2,3,1),cex=1,pch=19,lty=c(1,1,3),bty="n",legend=c("Larva","Besouro","Ambos"))
 
tiwb<- function(aw,ab,hw,hb) {
    return(1/(aw+ab)+hw*aw/(aw+ab)+hb*ab/(aw+ab))
}
#
plot(aw,tiwb(aw,ab,hw,hb),type="l",col=1,lwd=2,lty=2,frame=F,xlab="Abundância",ylab="Tempo total manipulação",
     ylim=c(0,200),xlim=c(0,0.03),main="Tempo de manipulação:\n Larvas = 1\n Besouros=60")
#
ab<-aw;
points(aw,tiwb(aw,ab,hw,hb),type="l",col=2,lwd=2,lty=2)
#
ab<-2*aw;
points(aw,tiwb(aw,ab,hw,hb),type="l",col=3,lwd=2,lty=2)
#
ab<-4*aw;
points(aw,tiwb(aw,ab,hw,hb),type="l",col=4,lwd=2,lty=2)
#
points(aw,tiw(aw,hw),type="l",col=1,lwd=1,lty=3)
legend("topright",col=c(1,2,3,4,1),lwd=c(2,2,2,2,1),lty=c(2,2,2,2,3),bty="n",title="Relações de abundância",
       legend=c("larvas=0.5*besouros","larvas=besouros","larvas=2*besouros","larvas=4*besouros","Somente larvas"))